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 $\mathcal{L}_{ij}$ or its intersection points
$\\{\textbf{\emph{c}}^{(1)}_{i},\textbf{\emph{c}}^{(2)}_{i},\textbf{\emph{c}}^{(1)}_{j},\textbf{\emph{c}}^{(2)}_{j}\\}$
with the two 3D conics. This is because the two conics do not actually
intersect one another. That no projective invariant exists from a pair of non-
coplanar conics agrees with the findings of Section 5.4.3. The non-existence
of projective invariants for three or more arbitrarily placed non-coplanar
conics follows by similar arguments.

In the special case where the conics do intersect one another, their
intersection is preserved when imaged from an arbitrary pose by a projective
camera. One of the most flexible ways to ensure two 3D conics intersect over
the complex numbers is to constrain them to lie on a nondegenerate quadric
surface (e.g., sphere, ellipsoid). This is now shown.

####  6.1.2 Pair of Non-Coplanar Conics on a Nondegenerate Quadric Surface

Consider craters modeled as conics lying on a nondegenerate quadric surface
$\mathcal{S}$. For the case of lunar craters, we observe that the Moon is very
nearly a sphere. Therefore, to a good approximation, define the surface of
Moon by the quadric locus parameterized by $\textbf{\emph{Q}}_{M}$,

$\mathcal{S}=\left\\{\bar{\bm{\xi}}\in\mathbb{P}^{3}\;|\;\bar{\bm{\xi}}^{T}\textbf{\emph{Q}}_{M}\bar{\bm{\xi}}=0\right\\}$
(65)

For a spherical moon of radius $R_{M}$,

$\textbf{\emph{Q}}_{M}=\left[\begin{array}[]{c c
c}R^{-2}_{M}\textbf{\emph{I}}_{3\times 3}&\textbf{0}_{3\times 1}\\\
\textbf{0}_{1\times 3}&-1\end{array}\right]$ (66)

Now, consider crater $i$ lying on plane $\mathcal{P}_{i}$. For the 3D rim of a
crater $i$ to also be on the surface of the Moon, it must lie on the
intersection of this plane with the $\mathcal{S}$. Therefore, the 3D conic
locus for crater $i$ must be the intersection

$\mathcal{C}_{i}=\mathcal{P}_{i}\cap\mathcal{S}$ (67)

which produces a circular crater if $\mathcal{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. 39, the crater may also be
viewed as the conic section produced by intersecting plane $\mathcal{P}_{i}$
with the Moon-centered quadric cone $\mathcal{X}_{i}$. However, since the
crater must lie on both $\mathcal{X}_{i}$ and $\mathcal{S}$, it is clear that
we can also write the 3D crater conic as,

$\mathcal{C}_{i}=\mathcal{X}_{i}\cap\mathcal{S}$ (68)

Therefore, we see that

$\mathcal{C}_{i}=\mathcal{X}_{i}\cap\mathcal{S}=\mathcal{P}_{i}\cap\mathcal{S}$
(69)

Suppose we have two craters: $\mathcal{C}_{i}$ and $\mathcal{C}_{j}$. Because
$\mathcal{C}_{i}$ and $\mathcal{C}_{j}$ are formed by the intersections of
planes $\mathcal{P}_{i}$ and $\mathcal{P}_{j}$ with the quadric surface
$\mathcal{S}$, we observe that,

$\mathcal{C}_{i}\cap\mathcal{C}_{j}=\left(\mathcal{P}_{i}\cap\mathcal{S}\right)\cap\left(\mathcal{P}_{j}\cap\mathcal{S}\right)=\left(\mathcal{P}_{i}\cap\mathcal{P}_{j}\right)\cap\mathcal{S}=\mathcal{L}_{ij}\cap\mathcal{S}$
(70)

where $\mathcal{L}_{ij}$ is the 3D line formed by the intersection of the two
crater planes (see Eq. 61).

If the two craters physically intersect, the line $\mathcal{L}_{ij}$ pierces
$\mathcal{S}$ and the intersection points are real. If the craters don’t
physically intersect, the intersection occurs over the complex numbers. Let
$\bar{\textbf{\emph{s}}}^{(1)}_{ij}$ and $\bar{\textbf{\emph{s}}}^{(2)}_{ij}$
be the two intersection points

$\\{\bar{\textbf{\emph{s}}}^{(1)}_{ij},\bar{\textbf{\emph{s}}}^{(2)}_{ij}\\}=\mathcal{L}_{ij}\cap\mathcal{S}$
(71)

Therefore, define $\mathcal{P}_{ij}$ as the plane passing through the center
of the Moon and the line $\mathcal{L}_{ij}$ (or, equivalently passing through
the center of the Moon, $\bar{\textbf{\emph{s}}}^{(1)}_{ij}$, and
$\bar{\textbf{\emph{s}}}^{(2)}_{ij}$). Clearly, the numerical stability of
computing $\mathcal{P}_{ij}$ in this manner is poor if the craters are close
to one another and $\mathcal{P}_{i}$ and $\mathcal{P}_{j}$ are nearly
coplanar. Better numeric stability may be achieved by taking a different
approach.

Figure 14: Two non-coplanar conics $\mathcal{C}_{i}$ and $\mathcal{C}_{j}$
lying on a nondegenerate quadric surface must intersect in two points (which
happen to be complex valued in this example). The two intersection points lie
on the line $\mathcal{L}_{ij}$ formed by the intersection of the two crater
planes. The intersection points also lie in the plane $\mathcal{P}_{ij}$
formed by the join of $\mathcal{L}_{ij}$ and the body center.

The conic intersection from Eq. 70 may also be written in terms of the quadric
cones instead of the planes,

$\mathcal{C}_{i}\cap\mathcal{C}_{j}=\left(\mathcal{X}_{i}\cap\mathcal{S}\right)\cap\left(\mathcal{X}_{j}\cap\mathcal{S}\right)=\left(\mathcal{X}_{i}\cap\mathcal{X}_{j}\right)\cap\mathcal{S}$
(72)

Without constraints, the intersection of two cones
$\mathcal{X}_{i}\cap\mathcal{X}_{j}$ is a set of four lines passing through
the center of the Moon that pierce $\mathcal{S}$ at eight different points.
Four of the points will be on the wrong side and come from conics that don’t
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,

$\left(\mathcal{X}_{i}\cap\mathcal{X}_{j}\right)\cap\mathcal{S}=\mathcal{C}_{i}\cap\mathcal{C}_{j}=\mathcal{L}_{ij}\cap\mathcal{S}=\\{\bar{\textbf{\emph{s}}}^{(1)}_{ij},\bar{\textbf{\emph{s}}}^{(2)}_{ij}\\}$
(73)

The advantage of using the quadric cone intersections to describe
$\bar{\textbf{\emph{s}}}^{(1)}_{ij}$ and $\bar{\textbf{\emph{s}}}^{(2)}_{ij}$
is that this remains well defined and more numerically stable as the craters
become close and their planes become nearly parallel. This is especially
important when viewing craters at a local level, where craters are nearly
coplanar.

Since the points $\bar{\textbf{\emph{s}}}^{(1)}_{ij}$ and
$\bar{\textbf{\emph{s}}}^{(2)}_{ij}$ lie on $\mathcal{C}_{i}$, the projection
of these points lie on the projection of $\mathcal{C}_{i}$ (which we called
$\mathcal{A}_{i}$). The same holds for $\mathcal{C}_{j}$, and its projection
$\mathcal{A}_{j}$. Therefore,the projection of
$\bar{\textbf{\emph{s}}}^{(1)}_{ij}$ and $\bar{\textbf{\emph{s}}}^{(2)}_{ij}$
is recoverable as two of the intersection points of the image conics. That is,

$\bar{\textbf{\emph{u}}}^{(k)}_{ij}=\textbf{\emph{P}}^{M}_{C}\bar{\textbf{\emph{s}}}^{(k)}_{ij}\subset\mathcal{A}_{i}\cap\mathcal{A}_{j},\quad
k=1,2$ (74)

where $\textbf{\emph{P}}^{M}_{C}$ is the projection matrix from Eq. 47 and the
pixel coordinate $\bar{\textbf{\emph{u}}}^{(k)}_{ij}$ is the projection of
$\bar{\textbf{\emph{s}}}^{(k)}_{ij}$ into the image.

The challenge, therefore, is to determine which of the four intersection
points of $\mathcal{A}_{i}$ and $\mathcal{A}_{j}$ are the projection of
$\bar{\textbf{\emph{s}}}^{(1)}_{ij}$ and $\bar{\textbf{\emph{s}}}^{(2)}_{ij}$.
There are six possible combinations of these four points, corresponding to six
possibilities for the projection of line $\mathcal{L}_{ij}$ into the line
$\ell_{ij}$. We now discuss how to uniquely compute $\ell_{ij}$.

Consider the two conics $\mathcal{A}_{i}$ and $\mathcal{A}_{j}$ described by
the conic locus matrices $\textbf{\emph{A}}_{i}$ and $\textbf{\emph{A}}_{j}$.
We recall that

$\bar{\textbf{\emph{u}}}^{T}\textbf{\emph{A}}_{i}\bar{\textbf{\emph{u}}}=0\text{
and
}\bar{\textbf{\emph{u}}}^{T}\textbf{\emph{A}}_{j}\bar{\textbf{\emph{u}}}=0$
(75)

It follows that we can form a pencil of conics parameterized by the scalar
ratio $\lambda:\mu$ that also pass through the same four intersection points
as $\textbf{\emph{A}}_{i}$ and $\textbf{\emph{A}}_{i}$,

$\bar{\textbf{\emph{u}}}^{T}\left(\lambda\textbf{\emph{A}}_{i}+\mu\textbf{\emph{A}}_{j}\right)\bar{\textbf{\emph{u}}}=0$
(76)

If we force the matrix $\lambda\textbf{\emph{A}}_{i}+\mu\textbf{\emph{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 $\ell_{ij}$ that we seek.

Letting $\mu=1$ for easy computation, we force the conic to be degenerate by
setting the determinant to zero,

$\text{det}\left(\lambda\textbf{\emph{A}}_{i}+\textbf{\emph{A}}_{j}\right)=0$
(77)

Since both $\textbf{\emph{A}}_{i}$ and $\textbf{\emph{A}}_{j}$ are full rank
and well-conditioned, we may rewrite this as

$\text{det}\left[\textbf{\emph{A}}_{j}\left(-\textbf{\emph{A}}_{i}\right)^{-1}-\lambda\textbf{\emph{I}}_{3\times
3}\right]=0$ (78)

which is a simple $3\times 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 [112]. The specific line we
seek is $\ell_{ij}$, which is the projection of $\mathcal{L}_{ij}$. As the
intersection of two non-coplanar planes, $\mathcal{L}_{ij}$ is real valued by
construction. Thus, its projection $\ell_{ij}$ is also real valued. Therefore
we always seek the eigenvalue leading to a degenerate conic of real valued
lines, discarding the four complex valued line possibilities. It will soon
become apparent which eigenvalue to choose.

Therefore, suppose we choose one of the eigenvalues to construct the
degenerate conic formed by

$\textbf{\emph{B}}_{ij}=\lambda\textbf{\emph{A}}_{i}+\textbf{\emph{A}}_{j}$
(79)

where we observe that $\textbf{\emph{B}}_{ij}$ has rank 2. The remaining steps
are now to find the lines described by $\textbf{\emph{B}}_{ij}$, determine if
they are real valued, and then determine which one is $\ell_{ij}$.

A degenerate conic made of two lines, g and h, is defined by the symmetric
matrix

$\textbf{\emph{B}}_{ij}=\textbf{\emph{g}}\textbf{\emph{h}}^{T}+\textbf{\emph{h}}\textbf{\emph{g}}^{T}$
(80)

where our task is now to find g and h given $\textbf{\emph{B}}_{ij}$. There
are a variety of ways to accomplish this task. Our specific approach is
similar to the framework outlined in [112]. Therefore, proceeding in this way,
we first recall that the intersection of the two lines in homogeneous
coordinates is computed as [48]

$\bar{\textbf{\emph{z}}}=\textbf{\emph{g}}\times\textbf{\emph{h}}$ (81)

Further, observing that

$\left[\bar{\textbf{\emph{z}}}\times\right]=\textbf{\emph{g}}\textbf{\emph{h}}^{T}-\textbf{\emph{h}}\textbf{\emph{g}}^{T}$
(82)

we find that that

$\textbf{\emph{B}}_{ij}+\left[\bar{\textbf{\emph{z}}}\times\right]=\left(\textbf{\emph{g}}\textbf{\emph{h}}^{T}+\textbf{\emph{h}}\textbf{\emph{g}}^{T}\right)+\left(\textbf{\emph{g}}\textbf{\emph{h}}^{T}-\textbf{\emph{h}}\textbf{\emph{g}}^{T}\right)=2\textbf{\emph{g}}\textbf{\emph{h}}^{T}$
(83)

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

$\textbf{\emph{B}}^{\ast}_{ij}=\left(\textbf{\emph{g}}\textbf{\emph{h}}^{T}+\textbf{\emph{h}}\textbf{\emph{g}}^{T}\right)^{\ast}=-\left(\textbf{\emph{g}}\times\textbf{\emph{h}}\right)\left(\textbf{\emph{g}}\times\textbf{\emph{h}}\right)^{T}=-\bar{\textbf{\emph{z}}}\bar{\textbf{\emph{z}}}^{T}$
(84)

Or, in terms of the original ellipses,

$\textbf{\emph{B}}^{\ast}_{ij}=\left(\lambda\textbf{\emph{A}}_{i}+\textbf{\emph{A}}_{j}\right)^{\ast}=-\bar{\textbf{\emph{z}}}\bar{\textbf{\emph{z}}}^{T}$
(85)

The diagonals of $-\textbf{\emph{B}}^{\ast}_{ij}$ are the squares of the
elements of $\bar{\textbf{\emph{z}}}$. Thus, defining the columns and elements
of $\textbf{\emph{B}}^{\ast}_{ij}$ as

$\textbf{\emph{B}}^{\ast}_{ij}=\left[\begin{array}[]{c c
c}\textbf{\emph{b}}^{\ast}_{1}&\textbf{\emph{b}}^{\ast}_{2}&\textbf{\emph{b}}^{\ast}_{3}\end{array}\right]=\left[\begin{array}[]{c
c c}b^{\ast}_{11}&b^{\ast}_{12}&b^{\ast}_{13}\\\
b^{\ast}_{21}&b^{\ast}_{22}&b^{\ast}_{23}\\\
b^{\ast}_{31}&b^{\ast}_{32}&b^{\ast}_{33}\end{array}\right]$ (86)

any column of $\textbf{\emph{B}}^{\ast}_{ij}$ may be scaled to find
$\bar{\textbf{\emph{z}}}$ according to

$\bar{\textbf{\emph{z}}}=-\textbf{\emph{b}}^{\ast}_{k}/\sqrt{-b^{\ast}_{kk}}$
(87)

where any $k\in\\{1,2,3\\}$ with $b_{kk}\neq 0$ will work. Best numerical
performance is achieved by selecting the value for $k$ that yields
$\max_{k}\|b^{\ast}_{kk}\|$. If the ellipses do not intersect, only one
eigenvalue from Eq. 78 will produce a real valued estimate of z. The
eigenvalue to choose is the one that makes the diagonal of
$\textbf{\emph{B}}^{\ast}_{ij}$ negative.

Therefore, with $\bar{\textbf{\emph{z}}}$ known, compute D

$\textbf{\emph{D}}=\textbf{\emph{B}}_{ij}+\left[\bar{\textbf{\emph{z}}}\times\right]=2\textbf{\emph{g}}\textbf{\emph{h}}^{T}$
(88)

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 $\ell_{ij}$. The
other does not. Each of the lines g and h divide ${\mathbb{P}}^{2}$ into two
regions. Since the line we seek comes from the projection of
$\mathcal{L}_{ij}=\mathcal{P}_{i}\cap\mathcal{P}_{j}$, the two conics must not
lie in the same region. That is, we seek a line that passes between the two
conics. Only one line will satisfy this condition and this is the line we
choose for $\ell_{ij}$.

If we were to compute the intersection of $\ell_{ij}$ with either
$\mathcal{A}_{i}$ or $\mathcal{A}_{j}$ we would obtain two points that are the
projection $\bar{\textbf{\emph{s}}}^{(1)}_{ij}$ and
$\bar{\textbf{\emph{s}}}^{(2)}_{ij}$. Fortunately such a computation is not
necessary since we work directly with $\ell_{ij}$ in subsequent discussions.

Briefly, we observe that $\mathcal{C}_{i}$ and the two points
$\bar{\textbf{\emph{s}}}^{(1)}_{ij}$ and $\bar{\textbf{\emph{s}}}^{(2)}_{ij}$
all lie in the plane $\mathcal{P}_{i}$. While a projective invariant exists
for two points and a conic (all coplanar) [39], the points must not lie on the
conic. Since $\bar{\textbf{\emph{s}}}^{(1)}_{ij}$ and
$\bar{\textbf{\emph{s}}}^{(2)}_{ij}$ lie on $\mathcal{C}_{i}$ by construction,
we seek an alternative way of constructing an invariant.

Consider instead the 2D line $\ell_{ij}$ formed by the projection of 3D line
$\mathcal{L}_{ij}$ (or, equivalently, by the join of image points
$\bar{\textbf{\emph{u}}}^{(1)}_{ij}$ and
$\bar{\textbf{\emph{u}}}^{(2)}_{ij}$). While there is not an invariant for a
single line and a conic, there is an invariant for two lines and a conic
[118]. This motivates the study of invariants for a triad of craters.

####  6.1.3 Triad of Non-Coplanar Conics on a Nondegenerate Quadric Surface

Suppose a projective camera observes three craters: $\mathcal{C}_{i}$,
$\mathcal{C}_{j}$, and $\mathcal{C}_{k}$. The intersection of the
corresponding planes ($\mathcal{P}_{i}$, $\mathcal{P}_{j}$, and
$\mathcal{P}_{k}$) produces the three 3D lines $\mathcal{L}_{ij}$,
$\mathcal{L}_{ik}$, and $\mathcal{L}_{jk}$. Three planes intersect in a point,
which is also the location where the 3D line formed by two planes intersects
the third plane. This point, defined as $a_{ijk}$, is the apex of the pyramid
formed by the three planes,

$a_{ijk}=\mathcal{P}_{i}\cap\mathcal{P}_{j}\cap\mathcal{P}_{k}=\mathcal{L}_{ij}\cap\mathcal{P}_{k}=\mathcal{L}_{ik}\cap\mathcal{P}_{j}=\mathcal{L}_{jk}\cap\mathcal{P}_{i}$
(89)

Now, let the three 3D craters $\mathcal{C}_{i}$, $\mathcal{C}_{j}$, and
$\mathcal{C}_{k}$ project to 2D conics $\mathcal{A}_{i}$, $\mathcal{A}_{j}$,
and $\mathcal{A}_{k}$ in the image. For any given pair of 3D conics (e.g.,
$\mathcal{C}_{i}$ and $\mathcal{C}_{j}$), the 3D line through their
intersection points (e.g., $\mathcal{L}_{ij}$) is coplanar with both of these
3D conics and it is possible to recover its projection (e.g., $\ell_{ij}$)
from just the projected image conics (e.g., $\mathcal{A}_{i}$ and
$\mathcal{A}_{j}$). Thus, for a triad of observed craters in an image, we may
compute the three lines $\ell_{ij}$, $\ell_{ik}$, and $\ell_{jk}$. Since
$\mathcal{C}_{i}$, $\mathcal{L}_{ij}$, $\mathcal{L}_{ik}$ are all coplanar, an
invariant exists in the image using $\mathcal{A}_{i}$, $\ell_{ij}$, and
$\ell_{ik}$. The same is true for craters $\mathcal{C}_{j}$ and
$\mathcal{C}_{k}$.

It is easy to see 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 $\mathbb{P}^{2}$, we can also form a cross-ratio of four
lines passing through a point. Therefore, as shown in Fig. 15, consider a
conic $\mathcal{A}_{i}$. Let the two lines $\ell_{ij}$ and $\ell_{ik}$ be
described by the $3\times 1$ vectors $\bm{\ell}_{ij}$ and $\bm{\ell}_{ik}$,
and let these two lines intersect at the point
$\bar{\textbf{\emph{p}}}=\bm{\ell}_{ij}\times\bm{\ell}_{ik}$. Using the pole-
polar relation for a conic, we may find the two lines from
$\bar{\textbf{\emph{p}}}$ that are tangent to the conic $\mathcal{A}_{i}$,
which we call $\textbf{\emph{w}}_{1}$ and $\textbf{\emph{w}}_{2}$. Since these
four lines go through the point $\bar{\textbf{\emph{p}}}$, we may form a cross
ratio that is a projective invariant

$\text{Cr}\left(\bm{\ell}_{ij},\textbf{\emph{w}}_{1},\textbf{\emph{w}}_{2},\bm{\ell}_{ik}\right)$
(90)

Figure 15: Visualization of contributing components to invariant for a conic
$\mathcal{A}_{i}$ and two lines $\bm{\ell}_{ij}$ and $\bm{\ell}_{jk}$. The
lines $\textbf{\emph{w}}_{1}$ and $\textbf{\emph{w}}_{2}$ pass through p and
are tangent to $\mathcal{A}_{i}$. The line m is the polar of point
$\bar{\textbf{\emph{p}}}$ with respect to $\mathcal{A}_{i}$.

We can now form a classical Cayley-Klein metric as

$d(\bm{\ell}_{ij},\bm{\ell}_{ik})=\frac{1}{2}\ln\left[\text{Cr}\left(\bm{\ell}_{ij},\textbf{\emph{w}}_{1},\textbf{\emph{w}}_{2},\bm{\ell}_{ik}\right)\right]$
(91)

which, from hyperbolic geometry, is known to be equivalent to

$\alpha_{i}=\cosh\left[d(\bm{\ell}_{ij},\bm{\ell}_{ik})\right]=\frac{\bm{\ell}^{T}_{ij}\textbf{\emph{A}}^{\ast}_{i}\bm{\ell}_{ik}}{\sqrt{\left(\bm{\ell}^{T}_{ij}\textbf{\emph{A}}^{\ast}_{i}\bm{\ell}_{ij}\right)\left(\bm{\ell}^{T}_{ik}\textbf{\emph{A}}^{\ast}_{i}\bm{\ell}_{ik}\right)}}$
(92)

where the adjugate matrix $\textbf{\emph{A}}^{\ast}_{i}$ describes the conic
envelope of $\mathcal{A}_{i}$. We may also write, therefore,

$\alpha^{2}_{i}=\cosh^{2}\left[d(\bm{\ell}_{ij},\bm{\ell}_{ik})\right]=\frac{\left(\bm{\ell}^{T}_{ij}\textbf{\emph{A}}^{\ast}_{i}\bm{\ell}_{ik}\right)^{2}}{\left(\bm{\ell}^{T}_{ij}\textbf{\emph{A}}^{\ast}_{i}\bm{\ell}_{ij}\right)\left(\bm{\ell}^{T}_{ik}\textbf{\emph{A}}^{\ast}_{i}\bm{\ell}_{ik}\right)}$
(93)

Since the cross ratio from Eq. 90 is a projective invariant, it follows
immediately that the Cayley-Klein metric from Eq. 91 and its hyperbolic cosine
in Eq. 92 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 $\alpha^{2}_{i}$ (the typical invariant found in
the literature for a coplanar set of a conic and two lines [118]) loses
descriptiveness for small distances $d(\bm{\ell}_{ij},\bm{\ell}_{ik})$. The
spacecraft navigator will immediately see that this is analogous to how small
inter-star angles $\theta$ lose descriptiveness when cataloged as $\cos\theta$
instead of by $\theta$ (i.e., we need to tighten the matching tolerance as
$\theta\rightarrow 0$ when using a $\cos\theta$ index). Thus, we prefer to
index on $d$ instead of $\alpha^{2}$.

If we choose to only consider craters that do not physically intersect (i.e.
do not overlap), then we are guaranteed that the cross ratio is positive. This
implies that $d\geq 0$ and is real, and that the denominator on the right-hand
side of Eq. 92 is also real. Since $d\geq 0$, we know that

$d=\text{arcosh}\,\alpha_{i}=\ln\left[\alpha_{i}+\sqrt{\alpha_{i}^{2}-1}\right]$
(94)

Therefore, we may relate the cross-ratio, the Cayley-Klein metric
$(\bm{\ell}_{ij},\bm{\ell}_{ik})$, and the invariant $\alpha_{i}$ according to

$\ln\left[\alpha_{i}+\sqrt{\alpha_{i}^{2}-1}\right]=d(\bm{\ell}_{ij},\bm{\ell}_{ik})=\frac{1}{2}\ln\left[\text{Cr}\left(\bm{\ell}_{ij},\textbf{\emph{w}}_{1},\textbf{\emph{w}}_{2},\bm{\ell}_{ik}\right)\right]$
(95)

where we choose to construct the invariant to index with the Cayley-Klein
metric,

$J_{i}=d(\bm{\ell}_{ij},\bm{\ell}_{ik})=\text{acosh}\left\\{\frac{\bm{\ell}^{T}_{ij}\textbf{\emph{A}}^{\ast}_{i}\bm{\ell}_{ik}}{\sqrt{\left(\bm{\ell}^{T}_{ij}\textbf{\emph{A}}^{\ast}_{i}\bm{\ell}_{ij}\right)\left(\bm{\ell}^{T}_{ik}\textbf{\emph{A}}^{\ast}_{i}\bm{\ell}_{ik}\right)}}\right\\}$
(96)

An identical procedure produces corresponding invariants for craters $j$ and
$k$

$J_{j}=d(\bm{\ell}_{ij},\bm{\ell}_{jk})=\text{acosh}\left\\{\frac{\bm{\ell}^{T}_{ij}\textbf{\emph{A}}^{\ast}_{j}\bm{\ell}_{jk}}{\sqrt{\left(\bm{\ell}^{T}_{ij}\textbf{\emph{A}}^{\ast}_{j}\bm{\ell}_{ij}\right)\left(\bm{\ell}^{T}_{jk}\textbf{\emph{A}}^{\ast}_{j}\bm{\ell}_{jk}\right)}}\right\\}$
(97)
$J_{k}=d(\bm{\ell}_{ik},\bm{\ell}_{jk})=\text{acosh}\left\\{\frac{\bm{\ell}^{T}_{ik}\textbf{\emph{A}}^{\ast}_{k}\bm{\ell}_{jk}}{\sqrt{\left(\bm{\ell}^{T}_{ik}\textbf{\emph{A}}^{\ast}_{k}\bm{\ell}_{ik}\right)\left(\bm{\ell}^{T}_{jk}\textbf{\emph{A}}^{\ast}_{k}\bm{\ell}_{jk}\right)}}\right\\}$
(98)

We know from our review of the literature that many previous attempts at
crater identification have accidentally employed feature descriptors with
elements that are not independent. Thus, before proceeding further, it is
necessary to ensure that the invariants $J_{i}$, $J_{j}$, and $J_{k}$ are
independent. Since we know from Section 5.5 that a $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 $J_{i}$, $J_{j}$, and $J_{k}$
are independent then we have found all the algebraically independent
invariants that exist. This is now shown.

We observe that the invariants $\alpha_{i}^{2},\alpha_{j}^{2},\alpha_{k}^{2}$
(and hence also $\alpha_{i},\alpha_{j},\alpha_{k}$ and $J_{i}$, $J_{j}$,
$J_{k}$) are algebraically independent. We calculate these invariants for a
particular model, since it suffices to show that the invariants are
independent when restricted to special subclasses of models. We assume that
the lunar surface is the sphere $x^{2}+y^{2}+z^{2}=1$, and that the three
craters are the circles defined by intersecting the sphere with the planes
given by $x=t_{1}$, $y=t_{2}$ and $z=t_{3}$ respectively. Within the plane
$z=t_{3}$ we have the crater given by $x^{2}+y^{2}=1-t_{3}^{2}$ and the lines
$x=t_{1}$ and $y=t_{2}$. So we get

$\textbf{\emph{A}}_{3}=\begin{pmatrix}1&0&0\\\ 0&1&0\\\
0&0&t_{3}^{2}-1\end{pmatrix}\mbox{ and
}\textbf{\emph{A}}^{*}_{3}=\begin{pmatrix}t_{3}^{2}-1&0&0\\\
0&t_{3}^{2}-1&0\\\ 0&0&1\end{pmatrix}$ (99) $\ell_{13}=\begin{pmatrix}1\\\
0\\\ -t_{1}\end{pmatrix}\mbox{ and }\ell_{23}=\begin{pmatrix}0\\\ 1\\\
-t_{2}\end{pmatrix}$ (100)

Thus, we may compute $\alpha_{3}^{2}$

$\alpha_{3}^{2}=\frac{(\ell_{13}\textbf{\emph{A}}_{3}^{*}\ell_{23})^{2}}{(\ell^{T}_{12}\textbf{\emph{A}}_{3}^{*}\ell_{12})(\ell^{T}_{13}\textbf{\emph{A}}_{3}^{*}\ell_{13})}=\frac{t_{1}^{2}t_{2}^{2}}{(t_{1}^{2}+t_{3}^{2}-1)(t_{2}^{2}+t_{3}^{2}-1)}$
(101)

and, by symmetry, we also find

$\alpha_{2}^{2}=\frac{t_{1}^{2}t_{3}^{2}}{(t_{1}^{2}+t_{2}^{2}-1)(t_{2}^{2}+t_{3}^{2}-1)}$
(102)
$\alpha_{1}^{2}=\frac{t_{2}^{2}t_{3}^{2}}{(t_{1}^{2}+t_{2}^{2}-1)(t_{1}^{2}+t_{3}^{2}-1)}$
(103)

It is now easy to verify that the Jacobi matrix

$J(\alpha_{1}^{2},\alpha_{2}^{2},\alpha_{3}^{2})=\Big{(}\frac{\partial\alpha_{i}^{2}}{\partial
t_{j}}\Big{)}_{1\leq i,j\leq 3}$ (104)

is invertible for some choice of $t_{1},t_{2},t_{3}$ (e.g.,
$t_{1}=t_{2}=t_{3}=1/2$). This implies that
$\alpha_{1}^{2},\alpha_{2}^{2},\alpha_{3}^{2}$ (and, therefore,
$J_{i},J_{j},J_{k}$) are algebraically independent.

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 $\mathcal{A}_{j}$ and $\mathcal{A}_{k}$, thus necessitating
the consideration of non-coplanar invariants. Two substantially different
views of the same crater pattern are shown. The image on the left is nearly
nadir pointing and the image on the right is pointed over 30 deg off nadir. As
can be seen, the non-coplanar invariants $J_{1}$, $J_{2}$, and $J_{3}$ are
identical in both images since they are formal projective invariants. This
makes evident the power of such invariants for crater pattern recognition.

Figure 16: Example of invariants $J_{i}$, $J_{j}$, and $J_{k}$ (blue) for the
same triad of craters ($\mathcal{A}_{i}$, $\mathcal{A}_{j}$, and
$\mathcal{A}_{k}$) seen from two different vantage points. The invariants
remain unchanged under perspective projection from any viewing geometry. The
triangle interior angles (yellow) are not a projective invariant. Synthetic
images of the lunar surface were produced using PANGU Planet Surface
Simulation Software developed by the Space Technology Centre at the University
of Dundee, Scotland [106].

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.
[46]. 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.

###  6.2 Geometry of Invariants for Coplanar Conics

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

####  6.2.1 Pairs of Coplanar Conics

A pair of coplanar conics has two projective invariants [39]. 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 $\mathbb{P}^{2}$, we can also form a projective invariant by the
cross-ratio of four lines passing through a point. Chasles’ theorem [123]
states that the cross-ratio is a constant for the pencil of four lines from
four points on a conic to a fifth point also on the conic (Fig. 17).
Therefore, we may directly form two invariants: the first by finding the cross
ratio of the lines from the four intersection points to any other point on
first conic, and the second by finding the cross ratio of the lines from the
same four intersection points to any other point on the second conic.

Figure 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
$\\{\bar{\textbf{\emph{a}}}_{k}\\}_{k=1}^{4}=\mathcal{A}_{i}\cap\mathcal{A}_{j}$
be the four intersection points of the conics in the image plane and let
$\bar{\textbf{\emph{p}}}$ be an arbitrary point on conic $\mathcal{A}_{i}$.
Since the points are written in homogeneous coordinates, we may form the line
joining $\bar{\textbf{\emph{p}}}$ and $\bar{\textbf{\emph{a}}}_{k}$ as
$\textbf{\emph{w}}_{k}=\bar{\textbf{\emph{p}}}\times\bar{\textbf{\emph{a}}}_{k}$.
Thus, the cross ratio
$\text{Cr}\left(\textbf{\emph{w}}_{1},\textbf{\emph{w}}_{2},\textbf{\emph{w}}_{3},\textbf{\emph{w}}_{4}\right)$
for ellipse $\mathcal{A}_{i}$ 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 $\mathcal{A}_{j}$. This
cross-ratio may also be rewritten in terms of the eigenvalues of the matrix
$\textbf{\emph{A}}_{i}^{-1}\textbf{\emph{A}}_{j}$ [83].

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 $\textbf{\emph{A}}_{i}$
and $\textbf{\emph{A}}_{j}$ describing the image conics. To see this, consider
the pencil of conics $\lambda\textbf{\emph{A}}_{i}+\mu\textbf{\emph{A}}_{j}$
parameterized by the scalars $\lambda$ and $\mu$. Recall now that 3D crater
$\mathcal{C}_{i}$ projects to image ellipse $\mathcal{A}_{i}$ according to the
homography from Eq. 59. Since the craters are assumed coplanar their
appearance in an image is related by a common homography,

$\textbf{\emph{H}}_{C}^{T}\left(\lambda\textbf{\emph{A}}_{i}+\mu\textbf{\emph{A}}_{j}\right)\textbf{\emph{H}}_{C}=s_{i}\lambda\textbf{\emph{C}}_{i}+s_{j}\mu\textbf{\emph{C}}_{j}$
(105)

Taking the determinant yields

$|\textbf{\emph{H}}_{C}|^{2}\,\left|\lambda\textbf{\emph{A}}_{i}+\mu\textbf{\emph{A}}_{j}\right|=|s_{i}\lambda\textbf{\emph{A}}^{\prime}_{i}+s_{j}\mu\textbf{\emph{A}}^{\prime}_{j}|$
(106)

which may be expanded as

$\displaystyle|\textbf{\emph{H}}_{C}|^{2}\,|\left|\lambda\textbf{\emph{A}}_{i}+\mu\textbf{\emph{A}}_{j}\right|$
$\displaystyle=|\textbf{\emph{H}}_{C}|^{2}\,\left(\Theta_{1}\lambda^{3}+\Theta_{2}\lambda^{2}\mu+\Theta_{3}\lambda\mu^{2}+\Theta_{4}\mu^{3}\right)$
$\displaystyle=\Theta^{\prime}_{1}\lambda^{3}+\Theta^{\prime}_{2}\lambda^{2}\mu+\Theta^{\prime}_{3}\lambda\mu^{2}+\Theta^{\prime}_{4}\mu^{3}$
(107)
$\displaystyle=|s_{i}\lambda\textbf{\emph{A}}^{\prime}_{i}+s_{j}\mu\textbf{\emph{A}}^{\prime}_{j}|$

Thus we find that

$\displaystyle|\textbf{\emph{H}}_{C}|^{2}\Theta_{1}=\Theta^{\prime}_{1}\quad|\textbf{\emph{H}}_{C}|^{2}\Theta_{4}=\Theta^{\prime}_{4}$
(108)
$\displaystyle|\textbf{\emph{H}}_{C}|^{2}\Theta_{2}=\Theta^{\prime}_{2}\quad|\textbf{\emph{H}}_{C}|^{2}\Theta_{3}=\Theta^{\prime}_{3}$
(109)

The coefficients ${\Theta_{1},\Theta_{2},\Theta_{3},\Theta_{4}}$ are only
_semi-invariants_ since their numerical value changes with
$|\textbf{\emph{H}}_{C}|$ and the arbitrary choices of $s_{i}$ and $s_{j}$.
These may easily be turned into _rational invariants_ by considering
appropriate ratios. To see this, first express the coefficients in terms of
$\textbf{\emph{A}}_{i}$ and $\textbf{\emph{A}}_{j}$,

$\Theta_{1}=\left|\textbf{\emph{A}}_{i}\right|\quad\quad\Theta_{4}=\left|\textbf{\emph{A}}_{j}\right|$
(110)
$\Theta_{2}=\text{Tr}\left[\textbf{\emph{A}}^{\ast}_{i}\textbf{\emph{A}}_{j}\right]=\left|\textbf{\emph{A}}_{i}\right|\,\text{Tr}\left[\textbf{\emph{A}}^{-1}_{i}\textbf{\emph{A}}_{j}\right]=\Theta_{1}\,\text{Tr}\left[\textbf{\emph{A}}^{-1}_{i}\textbf{\emph{A}}_{j}\right]$
(111)
$\Theta_{3}=\text{Tr}\left[\textbf{\emph{A}}^{\ast}_{j}\textbf{\emph{A}}_{i}\right]=\left|\textbf{\emph{A}}_{j}\right|\,\text{Tr}\left[\textbf{\emph{A}}^{-1}_{j}\textbf{\emph{A}}_{i}\right]=\Theta_{4}\,\text{Tr}\left[\textbf{\emph{A}}^{-1}_{j}\textbf{\emph{A}}_{i}\right]$
(112)

Likewise, we find that

$\Theta^{\prime}_{1}=s^{3}_{i}\left|\textbf{\emph{C}}_{i}\right|\quad\quad\Theta^{\prime}_{4}=s^{3}_{j}\left|\textbf{\emph{C}}_{j}\right|$
(113)
$\Theta^{\prime}_{2}=\frac{s_{j}}{s_{i}}\Theta^{\prime}_{1}\,\text{Tr}\left[\textbf{\emph{C}}_{i}^{-1}\textbf{\emph{C}}_{j}\right]$
(114)
$\Theta^{\prime}_{3}=\frac{s_{i}}{s_{j}}\Theta^{\prime}_{4}\,\text{Tr}\left[\textbf{\emph{C}}_{j}^{-1}\textbf{\emph{C}}_{i}\right]$
(115)

As shown in [123], the simplest pair of independent rational invariants are

$\frac{\Theta_{1}\Theta_{3}}{\Theta_{2}^{2}}=\frac{\Theta^{\prime}_{1}\Theta^{\prime}_{3}}{{\Theta^{\prime}_{2}}^{2}}\quad\text{and}\quad\frac{\Theta_{2}\Theta_{4}}{\Theta_{3}^{2}}=\frac{\Theta^{\prime}_{2}\Theta^{\prime}_{4}}{{\Theta^{\prime}_{3}}^{2}}$
(116)

which we quickly verify produce a pair of non-trivial scalars that are not
dependent on $\textbf{\emph{H}}_{C}$, $s_{i}$, or $s_{j}$ (hence, they are the
two fully generic invariants for a pair of conics):

$\frac{|\textbf{\emph{A}}_{j}|}{|\textbf{\emph{A}}_{i}|}\frac{\text{Tr}\left[\textbf{\emph{A}}_{j}^{-1}\textbf{\emph{A}}_{i}\right]}{\left(\text{Tr}\left[\textbf{\emph{A}}_{i}^{-1}\textbf{\emph{A}}_{j}\right]\right)^{2}}=\frac{\Theta_{1}\Theta_{3}}{\Theta_{2}^{2}}=\frac{\Theta^{\prime}_{1}\Theta^{\prime}_{3}}{{\Theta^{\prime}_{2}}^{2}}=\frac{|\textbf{\emph{C}}_{j}|}{|\textbf{\emph{C}}_{i}|}\frac{\text{Tr}\left[\textbf{\emph{C}}_{j}^{-1}\textbf{\emph{C}}_{i}\right]}{\left(\text{Tr}\left[\textbf{\emph{C}}_{i}^{-1}\textbf{\emph{C}}_{j}\right]\right)^{2}}$
(117)
$\frac{|\textbf{\emph{A}}_{i}|}{|\textbf{\emph{A}}_{j}|}\frac{\text{Tr}\left[\textbf{\emph{A}}_{i}^{-1}\textbf{\emph{A}}_{j}\right]}{\left(\text{Tr}\left[\textbf{\emph{A}}_{j}^{-1}\textbf{\emph{A}}_{i}\right]\right)^{2}}=\frac{\Theta_{2}\Theta_{4}}{\Theta_{3}^{2}}=\frac{\Theta^{\prime}_{2}\Theta^{\prime}_{4}}{{\Theta^{\prime}_{3}}^{2}}=\frac{|\textbf{\emph{C}}_{i}|}{|\textbf{\emph{C}}_{j}|}\frac{\text{Tr}\left[\textbf{\emph{C}}_{i}^{-1}\textbf{\emph{C}}_{j}\right]}{\left(\text{Tr}\left[\textbf{\emph{C}}_{j}^{-1}\textbf{\emph{C}}_{i}\right]\right)^{2}}$
(118)

Since the left-hand side of these expressions are independent of
$\textbf{\emph{H}}_{C}$, the same value is computed for any arbitrary
projection of the conics $\mathcal{C}_{i}$ and $\mathcal{C}_{j}$. These
results are in agreement with [110], which also follows [123].

Following the convention of [39], we now observe that the equations simplify
considerably by choosing the scale such that
$|\textbf{\emph{A}}_{i}|=|\textbf{\emph{A}}_{j}|=1$, making the coefficients
$\Theta_{1}$ and $\Theta_{4}$ trivial. Using this simplification, substitute
Eq. 110 and 113 into Eq. 108 to find

$\displaystyle|\textbf{\emph{H}}|^{2}=s_{i}^{3}=s_{j}^{3}$ (119)

which, in turn, leads to

$\Theta_{2}=\Theta^{\prime}_{2}\quad\rightarrow\quad\text{Tr}\left[\textbf{\emph{A}}^{-1}_{i}\textbf{\emph{A}}_{j}\right]=\text{Tr}\left[\textbf{\emph{C}}_{i}^{-1}\textbf{\emph{C}}_{j}\right]$
(120)
$\Theta_{3}=\Theta^{\prime}_{3}\quad\rightarrow\quad\text{Tr}\left[\textbf{\emph{A}}^{-1}_{j}\textbf{\emph{A}}_{i}\right]=\text{Tr}\left[\textbf{\emph{C}}_{j}^{-1}\textbf{\emph{C}}_{i}\right]$
(121)

Therefore, when conics are pre-scaled to
$|\textbf{\emph{A}}_{i}|=|\textbf{\emph{A}}_{j}|=1$, the two unique invariants
are simply

$I_{ij}=\Theta_{2}=\text{Tr}\left[\textbf{\emph{A}}^{-1}_{i}\textbf{\emph{A}}_{j}\right]$
(122)
$I_{ji}=\Theta_{3}=\text{Tr}\left[\textbf{\emph{A}}^{-1}_{j}\textbf{\emph{A}}_{i}\right]$
(123)

We note that [39] follows a completely different approach to finding $I_{ij}$
and $I_{ji}$, and even more geometrically-inspired derivations may be found in
[92]. 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 [19, 17]. We find, however, that better
performance is often achieved by considering a triad of craters.

####  6.2.2 Triads of Coplanar Conics

It is established in Section 5.2 that there exist $5d-8$ projective invariants
for a $d$-tuple of coplanar conics ($d\geq 2$). Thus, for a triad of coplanar
conics, there exists $15-8=7$ projective invariants. This fact has been known
for some time [111, 50].

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 [123]. As observed in [111], this evaluates to

$\displaystyle|\lambda\textbf{\emph{A}}_{i}+\mu\textbf{\emph{A}}_{j}+\sigma\textbf{\emph{A}}_{k}|=$
$\displaystyle\Theta_{1}\lambda^{3}+\Theta_{2}\lambda^{2}\mu+\Theta_{3}\lambda\mu^{2}+\Theta_{4}\mu^{3}+\Theta_{5}\lambda^{2}\sigma$
(124)
$\displaystyle+\Theta_{6}\lambda\sigma^{2}+\Theta_{7}\sigma^{3}+\Theta_{8}\mu^{2}\sigma+\Theta_{9}\mu\sigma^{2}+\Theta_{10}\lambda\mu\sigma$

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

$\Theta_{1}=\left|\textbf{\emph{A}}_{i}\right|,\quad\Theta_{4}=\left|\textbf{\emph{A}}_{j}\right|,\quad\Theta_{7}=\left|\textbf{\emph{A}}_{k}\right|$
(125)

and that

$\Theta_{2}=\Theta_{1}\,\text{Tr}\left[\textbf{\emph{A}}^{-1}_{i}\textbf{\emph{A}}_{j}\right],\quad\Theta_{3}=\Theta_{4}\,\text{Tr}\left[\textbf{\emph{A}}^{-1}_{j}\textbf{\emph{A}}_{i}\right]$
(126)
$\Theta_{5}=\Theta_{1}\,\text{Tr}\left[\textbf{\emph{A}}^{-1}_{i}\textbf{\emph{A}}_{k}\right],\quad\Theta_{6}=\Theta_{7}\,\text{Tr}\left[\textbf{\emph{A}}^{-1}_{k}\textbf{\emph{A}}_{i}\right]$
(127)
$\Theta_{8}=\Theta_{4}\,\text{Tr}\left[\textbf{\emph{A}}^{-1}_{j}\textbf{\emph{A}}_{k}\right],\quad\Theta_{9}=\Theta_{7}\,\text{Tr}\left[\textbf{\emph{A}}^{-1}_{k}\textbf{\emph{A}}_{j}\right]$
(128)

Note that $\Theta_{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

$\Theta_{10}=\frac{1}{2}\text{Tr}\left\\{\left[\left(\textbf{\emph{A}}_{j}+\textbf{\emph{A}}_{k}\right)^{\ast}-\left(\textbf{\emph{A}}_{j}-\textbf{\emph{A}}_{k}\right)^{\ast}\right]\textbf{\emph{A}}_{i}\right\\}$
(129)

where the reader is briefly reminded that $\textbf{\emph{A}}^{\ast}$ is the
adjugate of A. Computation of this result (which, interestingly, does not
appear in [111], [50], 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. 124 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 $\textbf{\emph{A}}_{i}$,
$\textbf{\emph{A}}_{j}$, and $\textbf{\emph{A}}_{k}$ such that
$|\textbf{\emph{A}}_{i}|=|\textbf{\emph{A}}_{j}|=|\textbf{\emph{A}}_{k}|=1$.
Thus the coefficients $\Theta_{1}$, $\Theta_{4}$, and $\Theta_{7}$ become
trivial, and the remaining seven coefficients become the unique invariants.
That is, with
$|\textbf{\emph{A}}_{i}|=|\textbf{\emph{A}}_{j}|=|\textbf{\emph{A}}_{k}|=1$,
we find that

$I_{ij}=\text{Tr}\left[\textbf{\emph{A}}^{-1}_{i}\textbf{\emph{A}}_{j}\right]\quad
I_{ji}=\text{Tr}\left[\textbf{\emph{A}}^{-1}_{j}\textbf{\emph{A}}_{i}\right]$
(130)
$I_{ik}=\text{Tr}\left[\textbf{\emph{A}}^{-1}_{i}\textbf{\emph{A}}_{k}\right]\quad
I_{ki}=\text{Tr}\left[\textbf{\emph{A}}^{-1}_{k}\textbf{\emph{A}}_{i}\right]$
(131)
$I_{jk}=\text{Tr}\left[\textbf{\emph{A}}^{-1}_{j}\textbf{\emph{A}}_{k}\right]\quad
I_{kj}=\text{Tr}\left[\textbf{\emph{A}}^{-1}_{k}\textbf{\emph{A}}_{j}\right]$
(132)
$I_{ijk}=\text{Tr}\left\\{\left[\left(\textbf{\emph{A}}_{j}+\textbf{\emph{A}}_{k}\right)^{\ast}-\left(\textbf{\emph{A}}_{j}-\textbf{\emph{A}}_{k}\right)^{\ast}\right]\textbf{\emph{A}}_{i}\right\\}$
(133)

Therefore, if three coplanar craters $\mathcal{C}_{i}$, $\mathcal{C}_{j}$, and
$\mathcal{C}_{k}$ project into three image conics $\mathcal{A}_{i}$,
$\mathcal{A}_{j}$, and $\mathcal{A}_{k}$ (described by the matrices
$\textbf{\emph{A}}_{i}$, $\textbf{\emph{A}}_{j}$, and
$\textbf{\emph{A}}_{k}$), then the seven scalar values in Eqs. 130 to 133 will
remain exactly the same under perspective projection regardless of the camera
position and attitude.

##  7 Creating a Crater Pattern Descriptor from Invariants

Having identified projective invariants that may be computed from the observed
crater rim geometry in an image (three for a triad of non-coplanar conics on a
quadric surface and seven for a triad of coplanar conics), we wish to use this
information to construct a pattern descriptor. This descriptor must have a
structure that can be consistently reproduced to allow for matching against a
static index, which requires some care since the invariants from Section 6
depend on the ordering of the observed craters. We note that the sensitivity
of invariants to ordering is not inherently a bad thing (after all, we must
specifically match each individual crater to the catalog)—though it does
result in a few different methods for building descriptors.

The index is built using descriptors for crater triads. We must choose,
therefore, whether we want to impose ordering of the three craters within the
triad before or after matching to the index. If we choose the former, then the
invariants developed before may be used directly. If we choose the latter, we
must either sort the projective invariants from Section 6 into a prescribed
order or transform them into _projective and permutation_ ($p^{2}$)
invariants. These conventions each have their advantages, which we now
discuss.

###  7.1 Crater Pattern Descriptors with Projective Invariants

Suppose we observe a triad of craters that we want to match to a precomputed
index in an order-dependent fashion. Since the crater pattern must be observed
from above (i.e., the camera must be looking downwards from above the lunar
surface because it cannot be inside the Moon), there are three possible
orderings of the craters. If we choose to arrange these in a clockwise fashion
within the image (see Fig. 18) then the three possibilities are $(1,2,3)$,
$(3,1,2)$, and $(2,3,1)$. For a match to be successful, the observed
descriptor must match the index descriptor. Therefore, if the descriptor
entries are order dependent, then a match only occurs when the observation
ordering matches index ordering.

There are at least two options for matching to an index with order-dependent
descriptors. First, if each crater triad has a single index entry, the index
must be searched three times (once with each possible observation ordering).
Alternatively, the index could include three entries for each crater triad
(one for each ordering) and the larger index is searched only once. The first
method (one index per triad) performs best in the presence of measurement
noise.

Figure 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.

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 descriptor:

$\textbf{\emph{J}}^{T}_{ijk}=\left[\begin{array}[]{c c c c c c
c}J_{i}&J_{j}&J_{k}\end{array}\right]$ (134)

where the individual elements are from Eqs. 96 to 98.

Likewise, for a triad of conics on a plane, we build a seven-element
descriptor:

$\textbf{\emph{I}}^{T}_{ijk}=\left[\begin{array}[]{c c c c c c
c}I_{ij}&I_{jk}&I_{ki}&I_{ji}&I_{kj}&I_{ik}&I_{ijk}\end{array}\right]$ (135)

where the individual elements are from Eqs. 130 to 133.

###  7.2 Crater Pattern Descriptors with Sorted Projective Invariants

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
$\textbf{\emph{J}}^{\prime
T}_{ijk}=[J_{i}^{\prime},J_{j}^{\prime},J_{k}^{\prime}]$, where
$J_{i}^{\prime}\leq J_{j}^{\prime}\leq J_{k}^{\prime}$ is a rearrangement of
$(J_{i},J_{j},J_{k})$, i.e., we sort the numbers $J_{i},J_{j},J_{k}$ from
small to large and let the result be the descriptor. This arrangement was
employed by Hanak [47, 46] 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., [82]
presort the entries in their index of crater triangles (which consists of two
triangle interior angles) such that the smallest interior angle is first.
Other examples of this idea abound within the space navigation community,
ranging from crater pattern identification to star pattern identification.

Returning to the problem at hand, we know from our labeling convention that
the crater pattern indices may only undergo a cyclic permutation. Thus, we may
cyclically permute the entries until smallest valued invariant is in the first
position—that is, where $(J_{i}^{\prime},J_{j}^{\prime},J_{k}^{\prime})$ is a
cyclic permutation of $(J_{i},J_{j},J_{k})$ with
$J^{\prime}_{i}=\min(J_{i},J_{j},J_{k})$.

Thus, for a triad of conics on a quadric surface, we build a three-element
descriptor:

$\textbf{\emph{J}}^{T}_{ijk}=\left[\begin{array}[]{c c c c c c
c}J^{\prime}_{i}&J^{\prime}_{j}&J^{\prime}_{k}\end{array}\right]$ (136)

where we label observed crater $i$ such that
$J^{\prime}_{i}=\min(J_{i},J_{j},J_{k})$ and the rest of the primed invariants
follow by cyclic permutation. The individual elements are from Eqs. 96 to 98.

Likewise, for a triad of conics on a plane, we build a seven-element
descriptor:

$\textbf{\emph{I}}^{T}_{ijk}=\left[\begin{array}[]{c c c c c c
c}I^{\prime}_{ij}&I^{\prime}_{jk}&I^{\prime}_{ki}&I^{\prime}_{ji}&I^{\prime}_{kj}&I^{\prime}_{ik}&I^{\prime}_{ijk}\end{array}\right]$
(137)

where we label observed crater $i$ such that
$I^{\prime}_{ij}=\min(I_{ij},I_{jk},I_{ki})$ and the rest of the primed
invariants follow by cyclic permutation of the crater labels $i,j,k$ (and, to
be clear, not by cyclic permutation of the elements of the descriptor vector).
The individual elements are from Eqs. 130 to 133.

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.

###  7.3 Crater Pattern Descriptors with Projective and Permutation ($p^{2}$)
Invariants

Our objective is to construct a descriptor that is _invariant_ to the order in
which the craters are observed. It should be no surprise, therefore, that
invariant theory is (once again) a useful mathematical framework for achieving
this objective. There is some precedence for this approach using the so-called
projective and permutation ($p^{2}$) invariants. The $p^{2}$ invariants are
discussed for a set of five points in [72] and for a pair of conics in [83].
An attempt was made to apply this to the problem of crater identification in
[105], 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.

####  7.3.1 Non-Coplanar Crater Patterns (Three-Element Descriptor)

We first develop the $p^{2}$ invariants for a triad of conics lying on a
quadric surface. We know there exists exactly three algebraically independent
invariants for this case, and we are able to compute three such invariants
$J_{i},J_{j},J_{k}$ directly from an image using Eqs. 96 to 98. Our goal is to
compute rational functions of $J_{i},J_{j},J_{k}$ that are invariant to a
cyclic permutation of the crater labels.

Consider, therefore, the triple of numbers $\\{x,y,z\\}$. The standard
generators of the ring of polynomial invariants for the symmetric group,
$S_{3}$, of all coordinate permutations are the elementary symmetric functions
(see [99, Chapter 4])

$\displaystyle e_{1}$ $\displaystyle=x+y+z$ (138) $\displaystyle e_{2}$
$\displaystyle=xy+yz+zx$ (139) $\displaystyle e_{3}$ $\displaystyle=xyz$ (140)

which are clearly invariant to permutations of $\\{x,y,z\\}$. To obtain the
generators of the invariant ring for the cyclic group
${\mathbb{Z}}/3{\mathbb{Z}}$, which we may also view as the alternating group
$A_{3}$, we require one more invariant:

$\Delta=(x-y)(y-z)(z-x)$ (141)

(see for example [99, Example 4.24]). We observe that $\Delta$ is invariant
under cyclic permutations, but that odd permutations change its sign. While
these polynomials could certainly be used to construct $p^{2}$ invariants and
build a pattern descriptor, we note that the polynomials are of different
degree which often leads to undesirable numerical properties for indexing.
Thus, we instead look for three rational functions, where the difference in
the degrees of the numerator and denominator are the same for all three
functions.

To do this, we begin by defining a polynomial map

$F(x,y,z)=(F_{1}(x,y,z),F_{2}(x,y,z),F_{3}(x,y,z))$ (142)

where $F_{1},F_{2},F_{3}$ are continuous rational functions with the property
that $F(x,y,z)=F(x^{\prime},y^{\prime},z^{\prime})$ if and only if
$(x^{\prime},y^{\prime},z^{\prime})$ is a cyclic permutation of $(x,y,z)$. Let
the cyclic group ${\mathbb{Z}}/3{\mathbb{Z}}$ act on ${\mathbb{R}}^{3}$ by
permuting the coordinates. Then we desire $F_{1},F_{2},F_{3}$ to be invariant
functions on ${\mathbb{R}}^{3}$.

Whenever a finite group acts by a linear transformation, there exists a
particular coordinate change where the action of the group becomes diagonal.
That is, in some coordinate system, the matrix describing the action of each
group element is a diagonal matrix. This follows from two facts from the
representation theory of finite groups (see [40, Chapter 1]). First, every
representation of a finite group is a direct sum of irreducible
representations. Second, every irreducible representation of an finite abelian
group is 1-dimensional. For the specific problem under consideration here (the
action of ${\mathbb{Z}}/3{\mathbb{Z}}$ on ${\mathbb{R}}^{3}$), we can make
this very explicit: namely, the change of coordinates making the action
diagonal is a discrete Fourier transform (DFT) (which involves complex
numbers). To make these ideas explicit, let $a,b,c$ be

$\displaystyle a$ $\displaystyle=x+y+z$ (143) $\displaystyle b$
$\displaystyle=x+\zeta y+\zeta^{2}z$ (144) $\displaystyle c$
$\displaystyle=x+\zeta^{2}y+\zeta z$ (145)

where $\zeta=e^{2\pi i/3}=-\frac{-1+\sqrt{3}i}{2}$. Then the generator
$\sigma$ of ${\mathbb{Z}}/3{\mathbb{Z}}$ acts by the diagonal matrix

$\left[\begin{array}[]{ccc}1&0&0\\\ 0&\zeta&0\\\
0&0&\zeta^{2}\end{array}\right]$

We define $F_{1}$, $F_{2}$, and $F_{3}$ as rational functions of $a,b,c$,

$\displaystyle F_{1}(x,y,z)$ $\displaystyle=a=x+y+z$ (146) $\displaystyle
F_{2}(x,y,z)$ $\displaystyle=\frac{b^{2}}{c}+\frac{c^{2}}{b}$ (147)
$\displaystyle=\frac{2(x^{3}+y^{3}+z^{3})+12xyz-3(x^{2}y+y^{2}x+y^{2}z+z^{2}y+z^{2}x+x^{2}z)}{x^{2}+y^{2}+z^{2}-(xy+yz+zx)}$
$\displaystyle F_{3}(x,y,z)$
$\displaystyle=\frac{1}{i}(\frac{b^{2}}{c}-\frac{c^{2}}{b})$ (148)
$\displaystyle=\frac{-3\sqrt{3}(x-y)(y-z)(z-x)}{x^{2}+y^{2}+z^{2}-(xy+yz+zx)}$

The denominator of $F_{2}$ and $F_{3}$ is equal to
$\frac{1}{2}[(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

$F_{2}^{2}+F_{3}^{2}=2[(x-y)^{2}+(y-z)^{2}+(z-x)^{2}]=2bc$

This implies that if $(x,y,z)$ is a sequence that converges to a point
$(t,t,t)$, then $F_{2}(x,y,z)$ and $F_{3}(x,y,z)$ converge to $0$.
Consequently, we can extend $F_{2}$ and $F_{3}$ to continuous functions on all
of ${\mathbb{R}}^{3}$ by defining them to be $0$ whenever $x=y=z$.

######  Lemma 3.

$F(x,y,z)=F(x^{\prime},y^{\prime},z^{\prime})$ if and only if
$(x^{\prime},y^{\prime},z^{\prime})$ is a cyclic permutation of $(x,y,z)$ for
all $(x,y,z),(x^{\prime},y^{\prime},z^{\prime})\in{\mathbb{R}}^{3}$.

###### Proof.

Suppose now that $F(x,y,z)=F(x^{\prime},y^{\prime},z^{\prime})$. Let $a,b,c$
be as given in Eqs. 146–148 and define similarly
$a^{\prime}=x^{\prime}+y^{\prime}+z^{\prime}$, etc. To begin, observe that if
$F_{2}(x,y,z)=F_{2}(x^{\prime},y^{\prime},z^{\prime})=F_{3}(x,y,z)=F_{3}(x^{\prime},y^{\prime},z^{\prime})=0$
then we must have $x=y=z$ and $x^{\prime}=y^{\prime}=z^{\prime}$ and from
$F_{1}(x,y,z)=3x=3x^{\prime}=F_{1}(x^{\prime},y^{\prime},z^{\prime})$ it
follows that $(x,y,z)=(x^{\prime},y^{\prime},z^{\prime})$. Now suppose that
$F_{2}(x,y,z)^{2}+F_{3}(x,y,z)^{2}$ (and, therefore,
$F_{2}(x^{\prime},y^{\prime},z^{\prime})^{2}+F_{3}(x^{\prime},y^{\prime},z^{\prime})^{2}$)
is nonzero. Then $b,b^{\prime},c,c^{\prime}$ are all nonzero. Since
$F_{2}(x,y,z)=F_{2}(x^{\prime},y^{\prime},z^{\prime})$ and
$F_{3}(x,y,z)=F_{3}(x^{\prime},y^{\prime},z^{\prime})$ we get

$\frac{b^{2}}{c}=\frac{F_{2}(x,y,z)+iF_{3}(x,y,z)}{2}=\frac{F_{2}(x^{\prime},y^{\prime},z^{\prime})+iF_{3}(x^{\prime},y^{\prime},z^{\prime})}{2}=\frac{(b^{\prime})^{2}}{c^{\prime}}$

and

$2bc=F_{2}(x,y,z)^{2}+F_{3}(x,y,z)^{2}=F_{2}(x^{\prime},y^{\prime},z^{\prime})^{2}+F_{3}(x^{\prime},y^{\prime},z^{\prime})^{2}=2b^{\prime}c^{\prime}.$

This implies that

$b^{3}=\frac{b^{2}}{c}(bc)=\frac{(b^{\prime})^{2}}{c^{\prime}}b^{\prime}c^{\prime}=(b^{\prime})^{3}.$

Therefore $b^{\prime}=\zeta^{s}b$ for some $s\in\\{0,1,2\\}$. Taking complex
conjugation, we get that $c^{\prime}=\zeta^{-s}c$. We conclude that $(a,b,c)$
and $(a^{\prime},b^{\prime},c^{\prime})$ lie in the same orbit
${\mathbb{Z}}/3{\mathbb{Z}}$, which means that
$(x^{\prime},y^{\prime},z^{\prime})$ is a permutation of $(x,y,z)$. ∎

One can also verify that the field of rational invariants
${\mathbb{R}}(x,y,z)^{{\mathbb{Z}}/3{\mathbb{Z}}}$ is generated by
$F_{1},F_{2},F_{3}$.

####  7.3.2 Coplanar Crater Patterns (Seven-Element Descriptor)

The development of $p^{2}$ invariants for a triad of conics lying on a plane
follows similar arguments as for conics on a quadric, though it requires more
care due to additional relationships between the seven elements in the
descriptor. Therefore, recall from Eqs. 130–133 that the seven invariants
consists of three invariant pairs,
$(I_{ij},I_{ji}),(I_{jk},I_{kj}),(I_{ki},I_{ik})$, and a 7th invariant for the
entire triad $I_{ijk}$. Because they’re constructed from a common set of
crater observations, the two invariants within any given pair always undergo a
common permutation. We also observe that $I_{ijk}$ is already symmetric and
invariant to the ellipse label order.

Simply by applying the results of the previous section, we can find six
${\mathbb{Z}}/3{\mathbb{Z}}$-invariant functions from
$F_{\ell}(I_{ij},I_{jk},I_{ki})$ and $F_{\ell}(I_{ji},I_{kj},I_{ik})$. The
descriptiveness of these six invariants, however, is not optimal because
members of a common pair are not forced to undergo the same permutation.
Therefore, we introduce two new invariants.

Again making use of a DFT, define $a_{\ell},b_{\ell},c_{\ell}$ as

$\displaystyle a_{\ell}$ $\displaystyle=x_{\ell}+y_{\ell}+z_{\ell}$ (149)
$\displaystyle b_{\ell}$ $\displaystyle=x_{\ell}+\zeta
y_{\ell}+\zeta^{2}z_{\ell}$ (150) $\displaystyle c_{\ell}$
$\displaystyle=x_{\ell}+\zeta^{2}y_{\ell}+\zeta z_{\ell}$ (151)

for $\ell=1,2$. Note that $c_{\ell}=\overline{b}_{\ell}$ where
$x_{\ell},y_{\ell},z_{\ell}$ are real. A generator of
${\mathbb{Z}}/3{\mathbb{Z}}$ acts on $(x_{1},y_{1},z_{1})$ and
$(x_{2},y_{2},z_{2})$ by the same cyclic permutation of coordinates, and it
acts on $b_{1},b_{2}$ with a scalar $\zeta$ and on $c_{1},c_{2}$ with a scalar
$\zeta^{2}$. We see, therefore, that $b_{1}c_{2}$ is invariant. We define
$G_{1}(x_{1},y_{1},z_{1},x_{2},y_{2},z_{2})$ and
$G_{2}(x_{1},y_{1},z_{1},x_{2},y_{2},z_{2})$ be the real and imaginary part of
$b_{1}c_{2}$, which we compute as

$b_{1}c_{2}=(x_{1}+\zeta y_{1}+\zeta^{2}z_{1})(x_{2}+\zeta^{2}y_{2}+\zeta
z_{2})=\\\
=(x_{1}x_{2}+y_{1}y_{2}+z_{1}z_{2})+\zeta(x_{1}z_{2}+y_{1}x_{2}+z_{1}y_{2})+\zeta^{2}(x_{1}y_{2}+y_{1}z_{2}+z_{1}x_{2})$

We may isolate the real part to obtain $G_{1}$

$G_{1}(x_{1},y_{1},z_{1},x_{2},y_{2},z_{2})=\\\
=(x_{1}x_{2}+y_{1}y_{2}+z_{1}z_{2})-{\textstyle\frac{1}{2}}(x_{1}z_{2}+y_{1}x_{2}+z_{1}y_{2}+x_{1}y_{2}+y_{1}z_{2}+z_{1}x_{2})\\\
={\textstyle\frac{3}{2}}(x_{1}x_{2}+y_{1}y_{2}+z_{1}z_{2})-{\textstyle\frac{1}{2}}(x_{1}+y_{1}+z_{1})(x_{2}+y_{2}+z_{2})$

and the imaginary part to obtain $G_{2}$

$G_{2}(x_{1},y_{1},z_{1},x_{2},y_{2},z_{2})=\frac{\sqrt{3}}{2}\big{[}(x_{1}z_{2}+y_{1}x_{2}+z_{1}y_{2})-(x_{1}y_{2}+y_{1}z_{2}+z_{1}x_{2})\big{]}\\\
=-\frac{\sqrt{3}}{2}\left|\begin{array}[]{ccc}1&1&1\\\ x_{1}&y_{1}&z_{1}\\\
x_{2}&y_{2}&z_{2}\end{array}\right|$

######  Lemma 4.

Suppose that
$x_{\ell},y_{\ell},z_{\ell},x_{\ell}^{\prime},y_{\ell}^{\prime},z_{\ell}^{\prime}\in{\mathbb{R}}$
for $\ell=1,2$ such that

$\displaystyle F(x_{1},y_{1},z_{1})$ $\displaystyle=$ $\displaystyle
F(x_{1}^{\prime},y_{1}^{\prime},z_{1}^{\prime})$ $\displaystyle
F(x_{2},y_{2},z_{2})$ $\displaystyle=$ $\displaystyle
F(x_{2}^{\prime},y_{2}^{\prime},z_{2}^{\prime})$ $\displaystyle
G(x_{1},y_{1},z_{1},x_{2},y_{2},z_{2})$ $\displaystyle=$ $\displaystyle
G(x_{1}^{\prime},y_{1}^{\prime},z_{1}^{\prime},x_{2}^{\prime},y_{2}^{\prime},z_{2}^{\prime})$

Then

$(x_{1},x_{2}),(y_{1},y_{2}),(z_{1},z_{2})$

is a cyclic permutation of

$(x_{1}^{\prime},x_{2}^{\prime}),(y_{1}^{\prime},y_{2}^{\prime}),(z_{1}^{\prime},z_{2}^{\prime})$

###### Proof.

By Lemma 3 we know that $(x_{1}^{\prime},y_{1}^{\prime},z_{1}^{\prime})$ is a
cyclic permutation of $(x_{1},y_{1},z_{1})$ and
$(x_{2}^{\prime},y_{2}^{\prime},z_{2}^{\prime})$ is a cyclic permutation of
$(x_{2},y_{2},z_{2})$, i.e., $b_{1}^{\prime}=\zeta^{s}b_{1}$ and
$b_{2}^{\prime}=\zeta^{t}b_{2}$ for some $s,t\in\\{0,1,2\\}$. It is not yet
clear that this is the same permutation, i.e., whether $s=t$. From
$G(x_{1},y_{1},z_{1},x_{2},y_{2},z_{2})=G(x_{1}^{\prime},y_{1}^{\prime},z_{1}^{\prime},x_{2}^{\prime},y_{2}^{\prime},z_{2}^{\prime})$
follows that
$b_{1}\overline{b}_{2}=b_{1}c_{2}=b_{1}^{\prime}c_{2}^{\prime}=b_{1}^{\prime}\overline{b}_{2}^{\prime}$.
On the other hand,
$b_{1}\overline{b}_{2}=b_{1}^{\prime}\overline{b}_{2}^{\prime}=(\zeta^{s}b_{1})(\zeta^{-t}\overline{b}_{2})=\zeta^{s-t}b_{1}\overline{b}_{2}$.
So either $b_{1}=0$, $b_{2}=0$, or $s=t$. If $t=s$, then we are done. If
$b_{1}=0$ then
$x_{1}=y_{1}=z_{1}=x_{1}^{\prime}=y_{1}^{\prime}=z_{1}^{\prime}$ and we are
also done because $(x_{2}^{\prime},y_{2}^{\prime},z_{2}^{\prime})$ is a cyclic
permutation of $(x_{2},y_{2},z_{2})$. The case $b_{2}=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

$[F(I_{ij},I_{jk},I_{ki}),F(I_{ji},I_{kj},I_{ik}),I_{ijk}]$

As was noted above, the descriptiveness of such a scheme is not optimal,
because this descriptor remains the same when the values of
$I_{ij},I_{jk},I_{ki}$ cyclically rotate while the values of
$I_{ji},I_{kj},I_{ik}$ stay the same. This means that 3 different
configurations might have the same descriptor.

A second option is to form a nine-element descriptor

$[F(I_{ij},I_{jk},I_{ki}),F(I_{ji},I_{kj},I_{ik}),\widetilde{G}(I_{ij},I_{jk},I_{ki},I_{ji},I_{kj},I_{ik}),I_{ijk}]$

This descriptor is more descriptive but uses more space. Here, we have
replaced the invariant $G$ with $\widetilde{G}$, which we define as

$\widetilde{G}(x_{1},y_{1},z_{1},x_{2},y_{2},z_{2})=\\\
\frac{G(x_{1},y_{1},z_{1},x_{2},y_{2},z_{2})}{\sqrt[4]{\Big{(}(x_{1}-y_{1})^{2}+(y_{1}-z_{1})^{2}+(z_{1}-x_{1})^{2}\Big{)}\Big{(}(x_{2}-y_{2})^{2}+(y_{2}-z_{2})^{2}+(z_{2}-x_{2})^{2}\Big{)}}}$

This is an essential step, since $G$ scales quadratically with
$(x_{1},y_{1},z_{1},x_{2},y_{2},z_{2})$ while $F$ scales linearly. Using $G$
instead of $\widetilde{G}$ results in a poorly scaled descriptor that
complicates (and sometimes precludes) nearest neighbor searches with efficient
data structures. Conversely, the function $\widetilde{G}$ scales linearly. We
note that $\widetilde{G}$ is not a rational function because of the $4$-th
root, but it does extend to a continuous function on ${\mathbb{R}}^{6}$ by
defining the function to be 0 whenever $x_{1}=y_{1}=z_{1}$ or
$x_{2}=y_{2}=z_{2}$.

A third option (and the one we suggest) is to use the seven-element descriptor

$[F(I_{ij},I_{jk},I_{ki}),F_{1}(I_{ji},I_{kj},I_{ik}),\widetilde{G}(I_{ij},I_{jk},I_{ki},I_{ji},I_{kj},I_{ik}),I_{ijk}]$

For generic configurations, this descriptor is as good as the nine-element
descriptor. The subfield of
$L={\mathbb{C}}(x_{1},y_{1},z_{1},x_{2},y_{2},z_{2})$ generated by
$F_{\ell}(x_{1},y_{1},z_{1})$, $\ell=1,2,3$, $F_{1}(x_{2},y_{2},z_{2})$ and
$G_{\ell}(x_{1},y_{1},z_{1},x_{2},y_{2},z_{2})$, for $\ell=1,2$ is equal to
$K={\mathbb{C}}(a_{1},a_{2},b_{1}^{3},b_{1}c_{1},b_{1}c_{2},c_{1}b_{2})$. It
is easy to see that $L=K(b_{1})$ and $b_{1}$ has degree three over $K$. So the
degree of the field extension $L/K$ is at most three, $K$ is contained in the
fixed field $L^{{\mathbb{Z}}/3{\mathbb{Z}}}$ and the degree of the extension
$L/L^{{\mathbb{Z}}/3{\mathbb{Z}}}$ is three. This implies that
$K=L^{{\mathbb{Z}}/3{\mathbb{Z}}}$. In other words, every
${\mathbb{Z}}/3{\mathbb{Z}}$-invariant rational function in
$I_{ij},I_{jk},I_{ki},I_{ji},I_{kj},I_{ik}$ is a rational function in the
first 6 elements of the descriptor.

However, in the degenerate case $I_{ij}=I_{jk}=I_{ki}=t$ for some fixed $t$ we
have $G(I_{ij},I_{jk},I_{ki},I_{ji},I_{kj},I_{ik})=0$. In that case, the
descriptor cannot distinguish between different configurations of the same
pattern (e.g., a pattern and it’s mirror). This is to be expected, since the
degenerate case corresponds to equally sized craters around an equilateral
triangle. Here, while there may be enough information to uniquely identify the
triad (i.e., match the triad descriptor to the database), there is not
sufficient information to disambiguate the specific crater labels with the
projected crater rims alone. No special action is required, since this
ambiguity is naturally handled in the pattern verification process (Section
9).

###  7.4 Remarks on Choosing a Descriptor Convention

In the sections above we introduce descriptors based on the projective
invariants, sorted projective invariants, and on the $p^{2}$ invariants. Which
one is best is often application dependent.

Descriptors built directly on the projective invariants (Section 7.1) have
better matching performance when presented with noisy data (see left-hand plot
of Fig. 19), though they require the index be searched three times per
observed pattern (one for each possible permutation).

Conversely, the descriptors built with the sorted invariants (Section 7.2) or
$p^{2}$ invariants (Section 7.3) only require one index search per observed
pattern, but are more sensitive to measurement noise. As a consequence the
sorted invariant and $p^{2}$ invariant descriptors are faster for low noise
situations, but slower for high noise situations (see right-hand plot of Fig.
19). This may seem counter-intuitive. The explanation, however, is
straightforward. As measurement noise increases, the sorted invariant and
$p^{2}$ invariant descriptors must attempt more triads (on average, as
compared to the ordered projective invariant descriptor) before finding a
correct nearest neighbor match. Thus, with large amounts of measurement noise,
the ordered projective invariant descriptor 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.

Figure 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.

##  8 Building a Global Crater Index

Global lunar crater identification requires multi-scale indexing and careful
catalog curation. When close to the Moon, a spacecraft sees small craters that
are nearly coplanar. Conversely, when far from the Moon, a spacecraft sees
larger craters distributed over the lunar sphere (not coplanar). This
immediately suggests at least two crater indexes be built—one for small
patterns of nearly coplanar craters and another for larger patterns of non-
coplanar craters. In practice, we find that the most efficient real-time
performance is achieved by creating a single index for the planar case and 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 [114] contains 1.3M craters over about 1–2 km
in diameter. It is immediately obvious that consideration of all the
$\binom{1.3\text{M}}{3}=3.66\times 10^{17}$ combinations of possible crater
triads is not reasonable. This further motivates construction of a hierarchy
of indexes with different scales to better manage the combinations of craters
that could plausibly be observed at the same time.

We propose the Hierarchical Equal Area isoLatitude Pixelization
(HEALPix)333Software implementations of HEALPix in C++, FORTRAN, and Python
are freely available online at https://healpix.sourceforge.io. framework [42]
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 [68]. We believe our work to
be the first application of HEALPix to the management of lunar surface
features.

All indexes are built using the same fundamental approach, but at different
scales. The procedure is as follows:

First, the lunar surface is tiled into $N_{pix}=12(2^{2k})$ surface pixels of
equal area. Next, a list is constructed for each pixel containing the catalog
entries for craters within a specified size range (e.g., minimum/maximum
diameter) and whose catalog fit was constructed using at least 90% of the
rim’s circumference.

Second, given a list of usable craters in each HEALPix surface pixel, we loop
through all the pixels to create crater triads. At each pixel, we consider
craters from the $3\times 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. 20 shows an example $3\times 3$
HEALPix grid (region of support for surface pixel 6318) overlayed on an
example Metric Camera image from Apollo 17.

Figure 20: 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 images: NASA/JSC/ASU. See [116, 70].

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 [89, 90] has seen extensive use for star pattern matching
in space applications [91, 129], and the $n$-d $k$-vector has been proposed
for space applications as well [71]. The specific choice of data structure is
not of primary concern in the present analysis so long as a reasonable
selection is made.

Figure 21: 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. Figure 22: 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.

###  8.1 Indexing Local Crater Patterns with Coplanar Invariants

To construct a global index of local-scale (small) crater patterns, we choose
a HEALPix resolution of $k=5$ for craters of diameter 4–30 km. This yields
$N_{pix}=12(2^{10})=12,288$ surface pixels, with each surface pixel having a
surface area of approximately $3,086\text{ 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
$\binom{20,737}{3}=1.49\times 10^{12}$ combinations, the HEALPix grouping
strategy only produces the 4.8M crater triads that are nearly coplanar.

The local crater patterns are assumed coplanar and craters are allowed to
follow an arbitrary elliptical shape. Thus, we construct a seven-element
descriptor using an appropriate method from Section 7.

###  8.2 Indexing Global Crater Patterns with Non-Coplanar Invariants

For crater patterns that are regional or global in extent, we find it helpful
to construct at least two indexes—though specific mission needs may require
more or less. The extension to a single global index or many global indexes is
trivial.

The Moon is very nearly a sphere at the global level. Therefore, because the
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\leq 1.1$. Since the global patterns are assumed to lie on a
nondegenerate quadric surface, we construct a three-element descriptor using
an appropriate method from Section 7.

The first global crater index has a HEALPix resolution of $k=3$ for craters of
diameter 25–125 km. This yields $N_{pix}=12(2^{6})=768$ surface pixels, with
each surface pixel having a surface area of approximately $49,390\text{
km}^{2}$. The index is good for matching crater patterns at the regional
level. There are 904 nearly circular craters ($a/b\leq 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
$\binom{904}{3}=1.23\times 10^{8}$ combinations, the HEALPix grouping strategy
only produces the $140,929$ crater triads. The example crater triad in Fig. 16
(see Section 6.1) was recognized using this index.

The second global crater index has a HEALPix resolution of $k=1$ for craters
having a diameter over $100$ km. This yields $N_{pix}=12(2^{2})=48$ surface
pixels, with each surface pixel having a surface area of approximately
$790,300\text{ 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\leq 1.1$) with catalog parameters
supported by $>90\%$ of the crater rim circumference and having a diameter
over $100$ km. Rather than producing $\binom{31}{3}=4,495$ combinations, the
HEALPix grouping strategy only produces the $707$ crater triads.

###  8.3 Remarks on the Utility of Crater Index Hierarchies

The primary purpose of having more than one index is to apply coplanar
invariants at the local scale (allowing for arbitrary elliptical crater shape)
and non-coplanar invariants at the regional/global scale (where curvature of
the Moon makes craters lie in substantially different planes).

In practice, the orbital regime is often known ahead of time (e.g., LLO,
cislunar), allowing us to only search crater patterns in the index of
appropriate scale. If, however, we are truly “lost-in-space” and don’t 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.

##  9 Matching Observed Crater Patterns to a Pre-Built Index

When supplied an image of the lunar surface, the objective is to recognize a
pattern of observed craters using the image conics corresponding to the
projection of the crater rims.

This is accomplished using a straightforward matching and and verification
procedure. First, a crater detection algorithm (CDA) produces an ellipse fit
to a crater rims in an image. Second, we consider triads of these image
conics, cycling through combinations using the Enhanced Pattern Shifting (EPS)
method [5]. For each crater triad, we compute the appropriate descriptor (see
Section 7) and query the corresponding pre-built index to find potential
matches (e.g., best $N\geq 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. 48. 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.

###  9.1 Pose from Non-Coplanar Conics in Correspondence

The usual approach for computing pose from matched craters only uses the
coordinates of the crater center. Here, it is essential to remember that the
center of an image ellipse does not generally produce a line-of-sight vector
to the center of a 3D ellipse (or circle). Some past crater identification
pipelines consider this effect (e.g., [105]), 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, 73]. 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 [130, 64, 136].
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 [81], which is generally cited as the solution to this problem.

The approach of [81], 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 [77, 78]). 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 [129, 113]. Regardless, we presume the attitude
transformation matrix $\textbf{\emph{T}}^{M}_{C}$ is known, such that the pose
estimation problem only requires a solution for the selenographic camera
location $\textbf{\emph{r}}_{M}$.

Therefore, let us begin with consideration of the projection of a single conic
as described by Eq. 59. Now, substitute for $\textbf{\emph{H}}_{C_{i}}$ from
Eq. 51 and expand to find

$\textbf{\emph{H}}_{M_{i}}^{T}\textbf{\emph{B}}_{i}\textbf{\emph{H}}_{M_{i}}-\textbf{\emph{H}}_{M_{i}}^{T}\textbf{\emph{B}}_{i}\textbf{\emph{r}}_{M}\textbf{\emph{k}}^{T}-\textbf{\emph{k}}\textbf{\emph{r}}_{M}^{T}\textbf{\emph{B}}_{i}\textbf{\emph{H}}_{M_{i}}+\textbf{\emph{r}}_{M}^{T}\textbf{\emph{B}}_{i}\textbf{\emph{r}}_{M}\textbf{\emph{k}}\textbf{\emph{k}}^{T}=s_{i}\textbf{\emph{C}}_{i}$
(152)

where

$\textbf{\emph{B}}_{i}=\textbf{\emph{T}}^{C}_{M}\textbf{\emph{K}}^{T}\textbf{\emph{A}}_{i}\textbf{\emph{K}}\textbf{\emph{T}}^{M}_{C}$
(153)

Recall here that $\textbf{\emph{A}}_{i}$ describes the measured image conic,
$\textbf{\emph{C}}_{i}$ describes the 3D crater conic in the $i$th crater’s
ENU frame, $\textbf{\emph{T}}^{M}_{C}$ is presumed known, and
$\textbf{\emph{k}}^{T}=[0\;0\;1]$. Thus, the only unknown in Eq. 152 is
$\textbf{\emph{r}}_{M}$.

Now, recalling Eq. 35 for $\textbf{\emph{H}}_{M_{i}}$, rearrange the left-hand
side into matrix form

$\left[\begin{array}[]{c c
c}\textbf{\emph{S}}^{T}\textbf{\emph{T}}_{E_{i}}^{M}\textbf{\emph{B}}_{i}\textbf{\emph{T}}^{E_{i}}_{M}\textbf{\emph{S}}&\textbf{\emph{S}}^{T}\textbf{\emph{T}}_{E_{i}}^{M}\textbf{\emph{B}}_{i}\left(\textbf{\emph{p}}_{M_{i}}-\textbf{\emph{r}}_{M}\right)\\\
\left(\textbf{\emph{p}}_{M_{i}}-\textbf{\emph{r}}_{M}\right)^{T}\textbf{\emph{B}}_{i}^{T}\textbf{\emph{T}}^{E_{i}}_{M}\textbf{\emph{S}}&\left(\textbf{\emph{p}}_{M_{i}}-\textbf{\emph{r}}_{M}\right)^{T}\textbf{\emph{B}}_{i}\left(\textbf{\emph{p}}_{M_{i}}-\textbf{\emph{r}}_{M}\right)\end{array}\right]=s_{i}\textbf{\emph{C}}_{i}$
(154)

where S is from Eq. 36. It is immediately evident that the upper-left $2\times
2$ submatrix is independent of $\textbf{\emph{r}}_{M}$, the lower-right
element is quadratic in $\textbf{\emph{r}}_{M}$, and the remaining terms are
linear in $\textbf{\emph{r}}_{M}$. Therefore, we first find the unknown scalar
$s_{i}$ using only the upper-left $2\times 2$ submatrix. Using the Frobenius
norm, compute $\hat{s}_{i}$ as

$\hat{s}_{i}=\arg\min_{s_{i}}\left\|\textbf{\emph{S}}^{T}\textbf{\emph{T}}_{E_{i}}^{M}\textbf{\emph{B}}_{i}\textbf{\emph{T}}^{E_{i}}_{M}\textbf{\emph{S}}-s_{i}\textbf{\emph{S}}^{T}\textbf{\emph{C}}_{i}\textbf{\emph{S}}\right\|^{2}_{F}$
(155)

which has the least squares solution

$\hat{s}_{i}=\frac{\text{vec}(\textbf{\emph{S}}^{T}\textbf{\emph{C}}_{i}\textbf{\emph{S}})^{T}\text{vec}(\textbf{\emph{S}}^{T}\textbf{\emph{T}}_{E_{i}}^{M}\textbf{\emph{B}}_{i}\textbf{\emph{T}}^{E_{i}}_{M}\textbf{\emph{S}})}{\text{vec}(\textbf{\emph{S}}^{T}\textbf{\emph{C}}_{i}\textbf{\emph{S}})^{T}\text{vec}(\textbf{\emph{S}}^{T}\textbf{\emph{C}}_{i}\textbf{\emph{S}})}$
(156)

With the scale $s_{i}$ known, we can now take the upper-right $2\times 1$
submatrix to form the linear system for all of the observed conics

$\left[\begin{array}[]{c c
c}\textbf{\emph{S}}^{T}\textbf{\emph{T}}_{E_{1}}^{M}\textbf{\emph{B}}_{1}\\\
\vdots\\\
\textbf{\emph{S}}^{T}\textbf{\emph{T}}_{E_{n}}^{M}\textbf{\emph{B}}_{n}\end{array}\right]\textbf{\emph{r}}_{M}=\left[\begin{array}[]{c
c
c}\textbf{\emph{S}}^{T}\textbf{\emph{T}}_{E_{1}}^{M}\textbf{\emph{B}}_{i}\textbf{\emph{p}}_{M_{1}}-\hat{s}_{i}\textbf{\emph{S}}^{T}\textbf{\emph{C}}_{1}\textbf{\emph{k}}\\\
\vdots\\\
\textbf{\emph{S}}^{T}\textbf{\emph{T}}_{E_{n}}^{M}\textbf{\emph{B}}_{n}\textbf{\emph{p}}_{M_{n}}-\hat{s}_{n}\textbf{\emph{S}}^{T}\textbf{\emph{C}}_{n}\textbf{\emph{k}}\end{array}\right]$
(157)

This may be solved for the selenographic camera location
$\textbf{\emph{r}}_{M}$ in the least squares sense.

###  9.2 Comparing Two Image Conics

In order to verify a crater match hypothesis, we require a measure of the
distance between two conics. For crater pattern verification, this distance of
interest is usually between the image ellipse we expect from the projection of
a crater’s rim ($\mathcal{A}_{i}$) and the image ellipse we measure of the
same crater’s rim ($\tilde{\mathcal{A}}_{i}$). For the moment, however, we
will briefly discuss how to compute the distance between two arbitrary
ellipses: $\mathcal{A}_{i}$ and $\mathcal{A}_{j}$.

Given two ellipses in an image, $\mathcal{A}_{i}$ and $\mathcal{A}_{j}$, we
seek a scalar distance metric $d(\mathcal{A}_{i},\mathcal{A}_{j})$ that
satisfies the three usual axioms for a distance metric [30]

Minimality: $d(\mathcal{A}_{i},\mathcal{A}_{j})=0\text{ iff
}\mathcal{A}_{i}=\mathcal{A}_{j}$. That is, the distance between an ellipse
and itself is zero.

Symmetry:
$d(\mathcal{A}_{i},\mathcal{A}_{j})=d(\mathcal{A}_{j},\mathcal{A}_{i})$. That
is, the distance from $\mathcal{A}_{i}$ to $\mathcal{A}_{j}$ is the same as
the distance from $\mathcal{A}_{j}$ to $\mathcal{A}_{i}$.

Triangle Inequality: $d(\mathcal{A}_{i},\mathcal{A}_{j})\leq
d(\mathcal{A}_{i},\mathcal{A}_{k})+d(\mathcal{A}_{k},\mathcal{A}_{j})$ as well
as a fourth axiom unique to this application

Similarity Invariance:
$d(\mathcal{A}_{i},\mathcal{A}_{j})=d(S[\mathcal{A}_{i}],S[\mathcal{A}_{j}])$,
where $S[\cdot]$ is a similarity transformation. That is, the distance between
$\mathcal{A}_{i}$ to $\mathcal{A}_{j}$ should not change if the two ellipses
undergo a common translation, rotation, or scaling in the image (i.e., undergo
a common similarity transformation). See Fig. 23. After reviewing the
literature, we found numerous approaches for comparing ellipses, but all fail
to meet one (or more) of the above four axioms. Thus, after a brief review of
existing techniques, a novel method is proposed.

Figure 23: Examples of an ellipse pair undergoing a simiarity transformation.
The left-most ellipse pair may be translated ($S_{1}$); translated and rotated
($S_{2}$); or translated, rotated, and scaled ($S_{3}$). Under the similarity
invariance axiom for an ellipse pair distance metric, all four of these
ellipse pairs should produce the same value for distance since they all have
the same relative geometry. Observe that all four of these ellipse pairs have
the same Jaccard distance and Gaussian angle. The area of ellipse intersection
(i.e., ellipse overlap) is shaded in light blue.

There exist various _ad hoc_ methods based on the explicit ellipse parameters
($a$, $b$, $x_{c}$, $y_{c}$, and $\psi$) [108, 28] or by measuring the
distance only at one specific point (e.g., via the Hausdorff distance [131,
127]), 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., $\text{det}(\textbf{\emph{A}}_{i})=1$)
and then compute the Frobenius norm of their difference [56]. There are
various normalizations one could choose for the matrix A, as summarized in
[63]. Nevertheless, such a comparison metric is given by

$d_{F}(\textbf{\emph{A}}_{i},\textbf{\emph{A}}_{j})=\|\textbf{\emph{A}}_{i}-\textbf{\emph{A}}_{j}\|^{2}_{F}$
(158)

though the specific geometric meaning of this is not clear. This metric also
fails to meet the similarity invariance axiom.

We also considered computing the distance as the line integral about one
ellipse of the Euclidean distance to the other ellipse. While this admits an
elegant and efficient solution with many interesting properties, it fails to
meet both the symmetry and similarity invariance axioms. To meet the symmetry
requirement, it would be possible to compute this metric in both directions
and sum the result, though this is expensive. Regardless, a distance computed
in this fashion will still fail to meet the similarity invariance axiom. The
method was thus abandoned for the present application.

The deficiencies of the above methods motivate the need for another metric
that satisfies all four distance axioms, has a geometrically consistent
meaning, and possesses well-understood statistics. There are at least two such
distance metrics: one using the Jaccard distance and another by interpreting
the ellipse as a multivariate Gaussian.

####  9.2.1 Ellipse Pair Distance Metric with the Jaccard Index

The Jaccard index is the ratio of an intersection to a union. In this case, we
take it to be the ratio of the areas of the ellipse intersection to the
ellipse union,

$J(\mathcal{A}_{i},\mathcal{A}_{j})=\frac{\text{Area}(\mathcal{A}_{i}\cap\mathcal{A}_{j})}{\text{Area}(\mathcal{A}_{i}\cup\mathcal{A}_{j})}$
(159)

The Jaccard distance, which is well-known to satisfy all three of the
classical distance metric axioms [75], is computed as

$d_{J}(\mathcal{A}_{i},\mathcal{A}_{j})=1-J(\mathcal{A}_{i},\mathcal{A}_{j})=1-\frac{\text{Area}(\mathcal{A}_{i}\cap\mathcal{A}_{j})}{\text{Area}(\mathcal{A}_{i}\cup\mathcal{A}_{j})}$
(160)

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 [57], 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 $\mathcal{A}_{i}$ by finding those pixels that satisfy the
inequality
$\bar{\textbf{\emph{u}}}^{T}\textbf{\emph{A}}_{i}\bar{\textbf{\emph{u}}}<0$,
which is a simple and parallelizable computation. Once we count the pixels in
this way, computing the Jaccard distance is trivial.

####  9.2.2 Ellipse Pair Distance Metric with the Gaussian Angle

It is possible to interpret the 2D ellipse $\mathcal{A}_{i}$ as a bivariate
Gaussian probability density function. If that interpretation is used for both
$\mathcal{A}_{i}$ and $\mathcal{A}_{j}$, the distance between two ellipses may
be computed using any number of measures for comparing Gaussian PDFs. The
authors of [139] 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
$\mathcal{A}_{i}$ as a bivariate Gaussian. Therefore, begin by subdividing the
conic locus matrix $\textbf{\emph{A}}_{i}$ describing the image ellipse
$\mathcal{A}_{i}$ as

$\textbf{\emph{A}}_{i}\propto\left[\begin{array}[]{c c
c}\textbf{\emph{Y}}_{i}&-\textbf{\emph{Y}}_{i}\textbf{\emph{y}}_{i}\\\
-\textbf{\emph{y}}_{i}^{T}\textbf{\emph{Y}}_{i}^{T}&(\textbf{\emph{y}}_{i}^{T}\textbf{\emph{Y}}_{i}\textbf{\emph{y}}_{i}-1)\end{array}\right]$
(161)

and it follows that that
$\bar{\textbf{\emph{u}}}^{T}_{i}\textbf{\emph{A}}_{i}\bar{\textbf{\emph{u}}}_{i}=0$
and that

$(\textbf{\emph{u}}_{i}-\textbf{\emph{y}}_{i})^{T}\textbf{\emph{Y}}_{i}(\textbf{\emph{u}}_{i}-\textbf{\emph{y}}_{i})=1$
(162)

With this particular scaling of $\textbf{\emph{A}}_{i}$, we immediately
recognize the $2\times 2$ submatrix $\textbf{\emph{Y}}_{i}$ as the shape of
the ellipse

$\textbf{\emph{Y}}_{i}=\left[\begin{array}[]{c c
c}\cos\psi_{i}&-\sin\psi_{i}\\\
\sin\psi_{i}&\cos\psi_{i}\end{array}\right]\left[\begin{array}[]{c c
c}1/a^{2}_{i}&0\\\ 0&1/b_{i}^{2}\end{array}\right]\left[\begin{array}[]{c c
c}\cos\psi_{i}&\sin\psi_{i}\\\ -\sin\psi_{i}&\cos\psi_{i}\end{array}\right]$
(163)

and $\textbf{\emph{y}}^{T}_{i}=[u_{c_{i}},v_{c_{i}}]$ as the image coordinates
of the ellipse center. If we interpret the ellipse as the $1\sigma$
isofrequency contour of a bivariate Gaussian distribution (i.e., the so-called
$1\sigma$ error ellipse), the ellipse $\mathcal{A}_{i}$ corresponds to the
distribution $N(\textbf{\emph{y}}_{i},\textbf{\emph{Y}}^{-1}_{i})$ with
probability density

$p_{\mathcal{A}_{i}}(\textbf{\emph{u}})=\frac{|\textbf{\emph{Y}}_{i}|}{2\pi}\exp\left[-\frac{1}{2}(\textbf{\emph{u}}-\textbf{\emph{y}}_{i})^{T}\textbf{\emph{Y}}_{i}(\textbf{\emph{u}}-\textbf{\emph{y}}_{i})\right]$
(164)

Therefore, given two ellipses $\mathcal{A}_{i}$ and $\mathcal{A}_{j}$, we may
compute the distance between these two ellipses as the angle between
$p_{\mathcal{A}_{i}}(\textbf{\emph{u}})$ and
$p_{\mathcal{A}_{j}}(\textbf{\emph{u}})$. Such an angle is computed as the
inner product

$\cos\theta=\frac{\int
p_{\mathcal{A}_{i}}(\textbf{\emph{u}})\,p_{\mathcal{A}_{j}}(\textbf{\emph{u}})}{\sqrt{\int
p^{2}_{\mathcal{A}_{i}}(\textbf{\emph{u}})\int
p^{2}_{\mathcal{A}_{j}}(\textbf{\emph{u}})}}$ (165)

and we compute a distance metric as

$d_{GA}(\mathcal{A}_{i},\mathcal{A}_{j})=\theta=\arccos\left[\frac{\int
p_{\mathcal{A}_{i}}(\textbf{\emph{u}})\,p_{\mathcal{A}_{j}}(\textbf{\emph{u}})}{\sqrt{\int
p^{2}_{\mathcal{A}_{i}}(\textbf{\emph{u}})\int
p^{2}_{\mathcal{A}_{j}}(\textbf{\emph{u}})}}\right]$ (166)

From inspection, this clearly satisfies the minimality and symmetry
conditions. Satisfaction of the triangle inequality follows from the
Cauchy–Schwarz inequality.

We compute

$\int
p^{2}_{\mathcal{A}_{i}}(\textbf{\emph{u}})=\frac{|\textbf{\emph{Y}}_{i}|^{2}}{4\pi^{2}}\int\exp\left[-\frac{1}{2}(\textbf{\emph{u}}-\textbf{\emph{y}}_{i})^{T}\textbf{\emph{Y}}_{i}(\textbf{\emph{u}}-\textbf{\emph{y}}_{i})\right]=\frac{|\textbf{\emph{Y}}_{i}|^{2}}{4\pi^{2}}\cdot\frac{2\pi}{|2\textbf{\emph{Y}}_{i}|}=\frac{|\textbf{\emph{Y}}_{i}|}{8\pi}$
(167)

To find the numerator, we first write

$\int
p_{\mathcal{A}_{i}}(\textbf{\emph{u}})\,p_{\mathcal{A}_{j}}(\textbf{\emph{u}})=\frac{|\textbf{\emph{Y}}_{i}||\textbf{\emph{Y}}_{j}|}{4\pi^{2}}\int\exp\left[-\frac{1}{2}\left((\textbf{\emph{u}}-\textbf{\emph{y}}_{i})^{T}\textbf{\emph{Y}}_{i}(\textbf{\emph{u}}-\textbf{\emph{y}}_{i})+(\textbf{\emph{u}}-\textbf{\emph{y}}_{j})^{T}\textbf{\emph{Y}}_{j}(\textbf{\emph{u}}-\textbf{\emph{y}}_{j})\right)\right]$
(168)

Now, considering the term inside the exponential,

$(\textbf{\emph{u}}-\textbf{\emph{y}}_{i})^{T}\textbf{\emph{Y}}_{i}(\textbf{\emph{u}}-\textbf{\emph{y}}_{i})+(\textbf{\emph{u}}-\textbf{\emph{y}}_{j})^{T}\textbf{\emph{Y}}_{j}(\textbf{\emph{u}}-\textbf{\emph{y}}_{j})=(\textbf{\emph{u}}-\textbf{\emph{z}})^{T}(\textbf{\emph{Y}}_{i}+\textbf{\emph{Y}}_{j})(\textbf{\emph{u}}-\textbf{\emph{z}})+c$
(169)

where

$\textbf{\emph{z}}=(\textbf{\emph{Y}}_{i}+\textbf{\emph{Y}}_{j})^{-1}(\textbf{\emph{Y}}_{i}\textbf{\emph{y}}_{i}+\textbf{\emph{Y}}_{j}\textbf{\emph{y}}_{j})$
(170)
$c=(\textbf{\emph{y}}_{i}-\textbf{\emph{y}}_{j})^{T}\textbf{\emph{Y}}_{i}(\textbf{\emph{Y}}_{i}+\textbf{\emph{Y}}_{j})^{-1}\textbf{\emph{Y}}_{j}(\textbf{\emph{y}}_{i}-\textbf{\emph{y}}_{j})$
(171)

Thus, noting that $c$ is a constant, we rewrite the integral as

$\displaystyle\int\exp\left[-\frac{1}{2}\left((\textbf{\emph{u}}-\textbf{\emph{y}}_{i})^{T}\textbf{\emph{Y}}_{i}(\textbf{\emph{u}}-\textbf{\emph{y}}_{i})+(\textbf{\emph{u}}-\textbf{\emph{y}}_{j})^{T}\textbf{\emph{Y}}_{j}(\textbf{\emph{u}}-\textbf{\emph{y}}_{j})\right)\right]$
(172)
$\displaystyle=\int\exp\left[-\frac{1}{2}(\textbf{\emph{u}}-\textbf{\emph{z}})^{T}(\textbf{\emph{Y}}_{i}+\textbf{\emph{Y}}_{j})(\textbf{\emph{u}}-\textbf{\emph{z}})+c\right]$
(173)
$\displaystyle=\exp\left[-\frac{c}{2}\right]\int\exp\left[-\frac{1}{2}(\textbf{\emph{u}}-\textbf{\emph{z}})^{T}(\textbf{\emph{Y}}_{i}+\textbf{\emph{Y}}_{j})(\textbf{\emph{u}}-\textbf{\emph{z}})\right]$
(174)
$\displaystyle=\exp\left[-\frac{c}{2}\right]\frac{2\pi}{|\textbf{\emph{Y}}_{i}+\textbf{\emph{Y}}_{j}|}$
(175)

Consequently,

$\displaystyle\int
p_{\mathcal{A}_{i}}(\textbf{\emph{u}})\,p_{\mathcal{A}_{j}}(\textbf{\emph{u}})$
$\displaystyle=\frac{|\textbf{\emph{Y}}_{i}||\textbf{\emph{Y}}_{j}|}{2\pi|\textbf{\emph{Y}}_{i}+\textbf{\emph{Y}}_{j}|}\exp\left[-\frac{c}{2}\right]$
(176)
$\displaystyle=\frac{|\textbf{\emph{Y}}_{i}||\textbf{\emph{Y}}_{j}|}{2\pi|\textbf{\emph{Y}}_{i}+\textbf{\emph{Y}}_{j}|}\exp\left[-\frac{1}{2}(\textbf{\emph{y}}_{i}-\textbf{\emph{y}}_{j})^{T}\textbf{\emph{Y}}_{i}(\textbf{\emph{Y}}_{i}+\textbf{\emph{Y}}_{j})^{-1}\textbf{\emph{Y}}_{j}(\textbf{\emph{y}}_{i}-\textbf{\emph{y}}_{j})\right]$

We may therefore compute the distance metric $d_{GA}$ by substitution of Eq.
167 and 176 into Eq. 166,

$d_{GA}(\mathcal{A}_{i},\mathcal{A}_{j})=\arccos\left\\{\frac{4\sqrt{|\textbf{\emph{Y}}_{i}||\textbf{\emph{Y}}_{j}|}}{|\textbf{\emph{Y}}_{i}+\textbf{\emph{Y}}_{j}|}\exp\left[-\frac{1}{2}(\textbf{\emph{y}}_{i}-\textbf{\emph{y}}_{j})^{T}\textbf{\emph{Y}}_{i}(\textbf{\emph{Y}}_{i}+\textbf{\emph{Y}}_{j})^{-1}\textbf{\emph{Y}}_{j}(\textbf{\emph{y}}_{i}-\textbf{\emph{y}}_{j})\right]\right\\}$
(177)

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 $\mathcal{A}_{i}$ and
$\mathcal{A}_{j}$.

####  9.2.3 A Fast Method for Assessing Ellipse Correspondence

We propose to use the Gaussian angle distance metric as the criteria for
evaluating a potential crater correspondence. In the case where
$\mathcal{A}_{j}$ is a perturbed version of $\mathcal{A}_{i}$, we find that,
to an excellent approximation,

$\frac{d^{2}_{GA}}{\sigma^{2}}\sim\chi^{2}_{4}$ (178)

where

$\sigma\approx\frac{0.85}{\sqrt{a_{i}b_{i}}}\sigma_{img}$ (179)

The validity of this approximation is illustrated in Fig. 24 for both a nearly
circular crater and a very elliptical crater.

Figure 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
histograms at the bottom show the statistic
$d^{2}/\sigma^{2}\sim\chi^{2}_{4}$. The red overlay on the histogram is the
analytic $\chi^{2}_{4}$ probability density function (PDF). (Credit for raw
scans of Apollo flight film images: NASA/JSC/ASU. See [116, 70].)

Since the distance metric is observed to follow a $\chi^{2}_{4}$ distribution,
we may construct a statistically-informed criterion for accepting a crater
match hypothesis. This is accomplished by evaluating the quantile function for
the $\chi^{2}_{4}$ distribution. For example, the 99th percentile for the
$\chi^{2}_{4}$ distribution occurs at 13.277, and would lead to an acceptance
criteria of

$\text{Hypothesis}(\mathcal{A}_{i}=\tilde{\mathcal{A}}_{i})\leftarrow\left\\{\begin{array}[]{l
l c}\text{Accept}&d^{2}_{{GA}_{i}}/\sigma^{2}\leq 13.276\\\
\text{Reject}&\text{otherwise}\end{array}\right.$ (180)

The appropriate percentile for accepting or rejecting craters is application
dependent and left as a design variable left to the analyst.

##  10 Numerical Results

###  10.1 Sensitivity to Crater Rim Fit and Viewing Geometry

We performed a Monte Carlo analysis to evaluate the sensitivity of the
proposed crater matching methodology to differing levels of measurement noise
and under different viewing geometries.

Therefore, as an example, consider a spacecraft placed randomly (with a
uniform distribution) around the lunar sphere, thus creating a situation where
we randomly image with equal probability any part of the lunar globe. We
assume a camera with a full field-of-view (FOV) of about $73.7\times 73.7$ deg
(same as the Apollo metric camera [36, 138]) and with a $2,200\times 2,200$
pixel focal plane array (creating a 4.8 MPixel image).

We now consider two parametric scenarios: (1) a nadir pointing camera with
varying errors in crater rim localization and (2) a fixed level of error in
crater rim localization and increasing off-nadir pointing angle. For the first
scenario (nadir pointing images), let the error in the ellipse parameters $a$,
$b$, $u_{c}$, $v_{c}$ all vary by $\sigma_{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.

####  10.1.1 Local Crater Matching Results

Suppose the randomly placed spacecraft has an altitude of 150 km above the
lunar surface. At this altitude, the camera views primarily local crater
patterns. In this case, results for the two Monte Carlo scenarios are shown in
Table 2 and Table 3.

Table 2: Monte Carlo results (100 trials) for local crater matching using seven-element descriptor. A nadir-pointing camera is placed randomly around the lunar globe at an altitude of 150 km, with ellipse fit error varying from 0 to 3 pixel. Ellipse Fit | Correct | Incorrect | No Match | Less Than | RSS Camera
---|---|---|---|---|---
Error | Matches | Matches | Found | Three Craters | Position Error
0.0 pixel | 100 | 0 | 0 | 0 | 2.5$\times 10^{-6}$ m
0.5 pixel | 96 | 0 | 1 | 3 | 116 m
1.0 pixel | 96 | 0 | 4 | 0 | 285 m
1.5 pixel | 94 | 0 | 4 | 2 | 428 m
2.0 pixel | 91 | 0 | 7 | 2 | 620 m
2.5 pixel | 93 | 0 | 7 | 0 | 696 m
3.0 pixel | 83 | 0 | 15 | 2 | 923 m
Table 3: Monte Carlo results (100 trials) for local crater matching using seven-element descriptor. A camera is placed randomly around the lunar globe at an altitude of 150 km, with an off-nadir pointing varying from 0 to 30 deg. Off-Nadir | Correct | Incorrect | No Match | Less Than | RSS Camera
---|---|---|---|---|---
Angle | Matches | Matches | Found | Three Craters | Position Error
0 deg | 98 | 0 | 1 | 1 | 140 m
10 deg | 97 | 0 | 3 | 0 | 147 m
20 deg | 99 | 0 | 1 | 0 | 134 m
30 deg | 96 | 0 | 4 | 0 | 178 m

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. 180), which would suggest that
$1-0.99^{3}\approx 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.

####  10.1.2 Global Crater Matching Results

Suppose the randomly placed spacecraft has an altitude of 600 km above the
lunar surface. At this altitude, the camera views primarily regional/global
crater patterns. In this case, results for the two Monte Carlo scenarios are
shown in Table 4 and Table 5.

Table 4: Monte Carlo results (100 trials) for global crater matching using three-element descriptor. A nadir-pointing camera is placed randomly around the lunar globe at an altitude of 600 km, with ellipse fit error varying from 0 to 3 pixel. Ellipse Fit | Correct | Incorrect | No Match | Less Than | RSS Camera
---|---|---|---|---|---
Error | Matches | Matches | Found | Three Craters | Position Error
0.0 pixel | 100 | 0 | 0 | 0 | 1.1$\times 10^{-6}$ m
0.5 pixel | 100 | 0 | 0 | 0 | 485 m
1.0 pixel | 100 | 0 | 0 | 0 | 1,294 m
1.5 pixel | 99 | 0 | 1 | 0 | 1,790 m
2.0 pixel | 99 | 0 | 1 | 0 | 2,605 m
2.5 pixel | 96 | 0 | 4 | 0 | 3,549 m
3.0 pixel | 95 | 0 | 5 | 0 | 4,500 m
Table 5: Monte Carlo results (100 trials) for global crater matching using three-element descriptor. A camera is placed randomly around the lunar globe at an altitude of 600 km, with an off-nadir pointing varying from 0 to 30 deg. Off-Nadir | Correct | Incorrect | No Match | Less Than | RSS Camera
---|---|---|---|---|---
Angle | Matches | Matches | Found | Three Craters | Position Error
0 deg | 100 | 0 | 0 | 0 | 549 m
10 deg | 99 | 0 | 1 | 0 | 539 m
20 deg | 99 | 0 | 1 | 0 | 652 m
30 deg | 97 | 0 | 2 | 1 | 866 m

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.

###  10.2 Local Matching on Clementine Images

We demonstrate the seven-element local matching descriptor on real images
collected by the Clementine spacecraft’s UV/Visible (UV/Vis) camera. During
it’s 71 days in lunar orbit in 1994 [100], 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\times 5.6$ deg FOV camera
with a $288\times 384$ pixel CCD focal plane array [67, 55]. 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 obtained from the image) were
perturbed by a Gaussian distribution with a standard deviation of
$\sigma_{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 [114] using the local crater index (HEALPix with $k=5$) described in
Section 8.

Figure 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 [96]. 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.

###  10.3 Global Matching on Synthetic Images on Reference Lunar Trajectory

We demonstrate the three-element global matching descriptor on synthetic
images along a reference lunar trajectory.

The synthetic images were generated using PANGU, and have a resolution of
$2,048\times 2,592$ pixels with a 30 deg horizontal FOV. Three example images
are shown in Fig. 26. For each image, the parameters of reference crater
ellipse ($a$, $b$, $u_{c}$, and $v_{c}$) were then perturbed by with
uncorrelated Gaussian noise of $\sigma_{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 ($\sigma_{img}=1.0$ pixel on $a$,
$b$, $u_{c}$, and $v_{c}$), we see the histograms widen as the spacecraft gets
farther away from the Moon (top to bottom in figure). This occurs because the
size of the pattern shrinks relative to the fixed magnitude of ellipse
localization error—hence the fixed error causes a greater perturbation of the
crater rim geometry.

Figure 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-$\sigma$ 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 [106].

##  11 Conclusion

The identification of craters in a digital image of the lunar surface is a
critical capability, both for terrain relative navigation (TRN) and for the
registration of scientific images. In this work, we provide the first
comprehensive and mathematically rigorous approach for image-based crater
identification using concepts from invariant theory. We first show that crater
rims tend to be elliptical in shape, suggesting that they are well modeled as
a conic (e.g., circle, ellipse) lying on the lunar surface. We then show that
there are no projective invariants for conics lying on the surface of an
arbitrarily shaped celestial body, but there are invariants for regularly
shaped bodies like the Moon. Specifically, for a $d$-tuple of conics, we show
that there are $3d-6$ algebraically independent projective invariants for
conics lying on a nondegenerate quadric surface and there are $5d-8$
algebraically independent projective invariants for conics lying on a plane.

With the number of algebraically independent projective invariants known, we
then develop a practical means of computing these invariants from just the
apparent crater rims in a digital image (i.e., from the perspective projection
of the 3D conics). These invariants may be used to form a pattern descriptor
that is always the same, regardless of camera viewpoint. Thus, these
descriptors may be precomputed for known crater patterns and stored in a
searchable index. We find that descriptors for a triad of craters provide good
matching results.

When given an image with observed crater rims, a match hypothesis for the
crater rim pattern may be found by a simple nearest neighbor search of the
index. With the resulting image-to-map crater correspondence hypothesis, we
show how to estimate the location of the camera in the least-squares sense
using the three crater rims directly (instead of just the crater center
coordinates). Finally, to verify the crater match and pose hypothesis, we
reproject the expected crater rim contours into the images and compare these
with the measured crater rims. This is accomplished with a new distance metric
that compares the entire ellipse fit (instead of just the crater center
coordinates). By verifying with a distance metric considering the full crater
rim contours, we are able to verify a pattern using fewer craters than methods
using only center coordinates.

These techniques are demonstrated in a variety of numerical experiments—both
on synthetic and real images. For a camera placed randomly around the Moon
with no _a priori_ position knowledge, successful matches generally exceed 90%
of the test cases. The remaining cases return “no match,” either because there
were fewer than three craters in a particular image or because measurement
noise corrupted the crater rim fits too much for unambiguous recognition. In
no case was an incorrect match returned, a sign that the new distance metric
can indeed verify match hypotheses with only three craters (instead of the
usual $\geq 5$ craters).

While this work represents a significant advance in methods for crater-based
TRN, a great deal of work remains to be done. What follows are a few thoughts
regarding the next steps:

While the repurposing of scientific crater catalogs for TRN is convenient,
these catalogs were built with different objectives and are generally
suboptimal for TRN applications. Thus, we suggest that the lunar exploration
community would benefit greatly from a purpose-built TRN crater catalog.

We provide a complete set of independent projective invariants for crater
patterns on the surface of the Moon in Section 6 (three invariants for a triad
of conics on a quadric surface, seven invariants for a triad of conics on a
plane). While there are no additional algebraically independent invariants, it
is possible to construct other algebraically _dependent_ invariants as
functions of the ones provided in this work (e.g., the $p^{2}$ invariants in
Section 7 are algebraic functions of the projective invariants from Section
6). There may exist other forms of these invariants that make lunar crater
patterns more distinct or that have superior stability in the presence of
measurement noise. Searching for such alternative formulations of the
invariants from Section 6 is an interesting problem.

This work develops a global index of crater descriptors at three different
scales using HEALPix. Most real missions likely don’t need a global index—but,
instead, just an index of craters along a reference trajectory. Although
HEALPix may still be used to manage varying scales, there is a great deal of
forward work in customizing the proof-of-concept outlined here to a mission-
specific application. The best design for an actual in-flight index is not
self evident and is the topic of ongoing work.

Not addressed in this work are crater detection algorithms (CDAs). Though many
different CDAs exist, it is not clear which specific algorithms are best
paired with the crater identification framework outlined here. Our crater
identification algorithm is critically dependent on the localization of the
best-fit ellipse to the crater rim, which is often not the metric used to
design or evaluate CDAs. Thus, finding (or developing) an appropriate CDA is
an obvious topic of follow-on work.

In conclusion, we find that crater-based TRN holds great promise for lost-in-
space navigation near the Moon. We introduce a robust method for crater
identification in this manuscript, though we note that additional work is
required to achieve a flight-capable system. Crater-based TRN is a rich field
for future study.

##  12 Acknowledgements

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
manuscript: Stuart Robbins, Joseph Mundy, Richard Hartley, Jeffrey Banks,
James McCabe, Kristen Schell, Jacob Hinkle, and Christopher D’Souza. 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.).

## References

  * [1] Acton, C.: Ancillary data services of nasa’s navigation and ancillary information facility.  Planetary and Space Sciences 44(1), 65–70 (1996).  DOI 10.1016/0032-0633(95)00107-7
  * [2] Acton, C., Bachman, N., Semenov, B., Wright, E.: A look towards the future in the handling of space science mission geometry.  Planetary and Space Sciences 150, 9–12 (2018).  DOI 10.1016/j.pss.2017.02.013
  * [3] Adams, D., Criss, T., Shankar, U.: Passive optical terrain relative navigation using APLNav.  In: IEEE Aerospace Conference (2008).  DOI 10.1109/AERO.2008.4526303
  * [4] Ansar, A., Daniilidis, K.: Linear pose estimation from points or lines.  IEEE Transactions on Pattern Analysis and Machine Intelligence 25(5), 578–589 (2003).  DOI 10.1109/TPAMI.2003.1195992
  * [5] Arnas, D., Fialho, M.A.A., Mortari, D.: Fast and robust kernel generators for star trackers.  Acta Astronautica 134, 291–302 (2017).  DOI 10.1016/j.actaastro.2017.02.016
  * [6] Arnas, D., Leake, C., Mortari, D.: The n-dimensional k-vector and its application to orthogonal range searching.  Applied Mathematics and Computation 372, 125010 (2020).  DOI 10.1016/j.amc.2019.125010
  * [7] Bay, H., Ess, A., Tuytelaars, T., Van Gool, L.: Speeded-up robust features (surf).  Computer Vision and Image Understanding 110(3), 346–359 (2008).  DOI 10.1016/j.cviu.2007.09.014
  * [8] Bellantoni, J.: Unidentified landmark navigation.  AIAA Journal 5(8), 1478–1483 (1967).  DOI 10.2514/3.4222
  * [9] Benedix, G., Norman, C., Bland, P., Towner, M., Paxman, J., Tan, T.: Automated detection of martian craters using a convolutional neural network.  In: 49th Lunar and Planetary Science Conference 2018, 2083 (2018)
  * [10] Bennett, C.L., et. al: The microwave anisotropy probe mission.  The Astrophysical Journal 583, 1–23 (2003).  DOI 10.1086/345346
  * [11] Bentley, J.L.: Multidimensional binary search trees used for associative searching.  Communications of the ACM 18(9), 509–517 (1975).  DOI 10.1145/361002.361007
  * [12] Borkowski, K.M.: Accurate algorithms to transform geocentric to geodetic coordinates.  Bulletin géodésique 63, 50–56 (1989).  DOI 10.1007/BF02520228
  * [13] Bottke, W., Love, S., Tytell, D., Glotch, T.: Interpreting the elliptical crater populations on mars, venus, and the moon.  Icarus 145, 108–121 (2000).  DOI 10.1006/icar.1999.6323
  * [14] Bowring, B.R.: The accuracy of geodetic latitude and height equations.  Survey Review 28(218), 202–206 (1985).  DOI 10.1179/sre.1985.28.218.202
  * [15] Burns, J.B., Weiss, R.S., Riseman, E.M.: Geometric Invariance in Computer Vision, chap. The non-existence of general-case view-invariants, pp. 120–131.  MIT Press, Cambridge, MA (1992)
  * [16] Chen, H.F., Belhumeur, P.N., Jacobs, D.W.: In search of illumination invariants.  In: IEEE Conference on Computer Vision and Pattern Recognition (2000).  DOI 10.1109/CVPR.2000.855827
  * [17] Cheng, Y., Ansar, A.: Landmark based position estimation for pinpoint landing on mars.  In: IEEE International Conference on Robotics and Automation (2005).  DOI 10.1109/ROBOT.2005.1570338
  * [18] Cheng, Y., Goguen, J., Johnson, A., Leger, C., Matthies, L., San Martin, M., Wilson, R.: The mars exploration rovers decent image motion estimation system.  IEEE Intelligent Systems 19(3), 13–21 (2004).  DOI 10.1109/MIS.2004.18
  * [19] Cheng, Y., Miller, J.: Autonomous landmark based spacecraft navigation.  In: AAS/AIAA Astrodynamics Specialist Conference, AAS 03-223 (2003)
  * [20] Christian, J.: Accurate planetary limb localization for image-based spacecraft navigation.  Journal of Spacecraft and Rockets 54(3), 708–730 (2017).  DOI 10.2514/1.A33692
  * [21] Christian, J.: Perspective projection of ellipses and ellipsoids with applications to spacecraft navigation.  In: AAS/AIAA Astrodynamics Specialist Conference, AAS 17-800 (2017)
  * [22] Christian, J.: Autonomous initial orbit determination with optical observations of unknown planetary landmarks.  Journal of Spacecraft and Rockets 56(1), 211–220 (2019).  DOI 10.2514/1.A34259
  * [23] Christian, J., Benhacine, L., Hikes, J., D’Souza, C.: Geometric calibration of the orion optical navigation camera using star field images.  The Journal of the Astronautical Sciences 63(4), 335–353 (2016).  DOI 10.1007/s40295-016-0091-3
  * [24] Christian, J.A.: Optical navigation for a spacecraft in a planeatry system.  Ph.D. thesis, The University of Texas at Austin (2010)
  * [25] Christian, J.A., Hong, L., McKee, P., Christensen, R., Crain, T.P.: Image-based lunar terrain relative navigation without a map: Measurements.  Journal of Spacecraft and Rockets (2020)
  * [26] Clemens, D.T., Jacobs, D.W.: Space and time bounds on indexing 3d models from 2d images.  IEEE Transactions on Pattern Analysis and Machine Intelligence 13(10), 1007–1017 (1991).  DOI 10.1109/34.99235
  * [27] Collins, G.S., Elbeshausen, D., Davison, T.M., Robbins, S.J., Hynek, B.M.: The size-frequency distribution of elliptical impact craters.  Earth and Planetary Science Letters 310, 1–8 (2011).  DOI 10.1016/j.epsl.2011.07.023
  * [28] Cuevas, E., González, M., Zaldívar, D., Pérez-Cisneros, M.: Multi-ellipses detection on images inspired by collective animal behavior.  Neural Computing and Applications volume 24, 1019–1033 (20014).  DOI 10.1007/s00521-012-1332-4
  * [29] Cui, P., Gao, X., Zhu, S., Shao, W.: Visual navigation using edge curve matching for pinpoint planetary landing.  Acta Astronautica 146, 171–180 (2018).  DOI 10.1016/j.actaastro.2018.02.033
  * [30] Cullinane, M.J.: Metric axioms and distance.  The Mathematical Gazette 95, 414–419 (2011).  DOI 10.1017/S0025557200003508
  * [31] DeLatte, D., Crites, S., Guttenberg, N., Yairi, T.: Automated crater detection algorithms from a machine learning perspective in the convolutional neural network era.  Advances in Space Research 64(8), 1615–1628 (2019).  DOI 10.1016/j.asr.2019.07.017
  * [32] Derksen, H., Kemper, G.: Computational invariant theory, _Encyclopaedia of Mathematical Sciences_ , vol. 130, enlarged edn.  Springer, Heidelberg (2015).  DOI 10.1007/978-3-662-48422-7.  With two appendices by Vladimir L. Popov, and an addendum by Norbert A’Campo and Popov, Invariant Theory and Algebraic Transformation Groups, VIII
  * [33] Dieudonné, J.A., Carrell, J.B.: Invariant theory, old and new.  Advances in Math. 4, 1–80 (1970) (1970).  DOI 10.1016/0001-8708(70)90015-0
  * [34] Dieudonné, J.A., Carrell, J.B.: Invariant theory, old and new.  Academic Press, New York-London (1971)
  * [35] Downes, L., Steiner, T., How, J.: Deep learning crater detection for lunar terrain relative navigation.  In: AIAA SciTech Forum (2020).  DOI 10.2514/6.2020-1838
  * [36] Edmundson, K., Alexandrov, O., Archinal, B., Becker, K., Kirk, R., Moratto, Z., Nefian, A., Richie, J., Robinson, M.: Photogrammetric processing of apollo 15 metric camera oblique images.  In: The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLI-B4 (2016).  DOI 10.5194/isprs-archives-XLI-B4-375-2016
  * [37] Elbeshausen, D., Wünnemann, K., Collins, G.S.: The transition from circular to elliptical impact craters.  Journal of Geophysical Research: Planets 118, 2295–2309 (2013).  DOI 10.1002/2013JE004477
  * [38] Emami, E., Bebis, G., Nefian, A., Fong, T.: Automatic crater detection using convex grouping and convolutional neural networks.  In: 11th International Symposium on Visual Computing (2015).  DOI 10.1007/978-3-319-27863-6˙20
  * [39] Forsyth, D., Mundy, J., Zisserman, A., Coelho, C., Heller, A., Rothwell, C.: Invariant descriptors for 3-D object recognition and pose.  IEEE Transactions on Pattern Analysis and Machine Intelligence 13(10), 971–991 (1991).  DOI 10.1109/34.99233
  * [40] Fulton, W., Harris, J.: Representation theory, _Graduate Texts in Mathematics_ , vol. 129.  Springer-Verlag, New York (1991).  DOI 10.1007/978-1-4612-0979-9.  A first course, Readings in Mathematics
  * [41] Gaskell, R.W., Barnouin-Jha, O.S., Scheers, D.J., et al.: Characterizing and navigating small bodies with imaging data.  Meteoritics & Planetary Science 43(6), 1049–1061 (2008).  DOI 10.1111/j.1945-5100.2008.tb00692.x
  * [42] Górski, K., Hivon, E., Banday, A., Wandelt, B., Hansen, F., Reinecke, M., Bartelmann, M.: Healpix: A framework for high-resolution discretization and fast analysis of data distributed on the sphere.  The Astrophysical Journal 622(2), 759–771 (2005).  DOI 10.1086/427976
  * [43] Group, C.A.T.W.: Standard techniques for presentation and analysis of crater size-frequency data.  Icarus 37, 467–474 (1979).  DOI 10.1016/0019-1035(79)90009-5
  * [44] Gupta, R., Hartley, R.I.: Linear pushbroom cameras.  IEEE Transactions on Pattern Analysis and Machine Intelligence 19(9), 963–975 (1997).  DOI 10.1109/34.615446
  * [45] Haboush, W.J.: Reductive groups are geometrically reductive.  Ann. of Math. (2) 102(1), 67–83 (1975).  DOI 10.2307/1970974
  * [46] Hanak, C., Crain, T., Bishop, R.: Crater identification algorithm for the lost in low lunar orbit scenario.  In: AAS Guidance & Control Conference, AAS 10-052 (2008)
  * [47] Hanak, F.: Lost in low lunar orbit crater pattern detection and identification.  Ph.D. thesis, The University of Texas at Austin (2009)
  * [48] Hartley, R., Zisserman, A.: Multiple View Geometry, 2nd Ed.  Cambridge University Press, Cambridge, UK (2003)