| 1361 | | \section{Our approach} |
| 1362 | | \label{ourapproach} |
| 1363 | | |
| 1364 | | \subsection{Preprocessing} |
| 1365 | | |
| 1366 | | Astronomers have developed computational methods for detecting and |
| 1367 | | localizing stars, galaxies, and other objects in astronomical images. |
| 1368 | | In astronomy this process is known as \emph{source extraction}. |
| 1369 | | Telescopes are designed to produce a predictable point-spread function |
| 1370 | | which covers several pixels, and this allows stars, which are |
| 1371 | | effectively point sources, to be localized with subpixel accuracy. |
| 1372 | | Sources that have image profiles larger than the point spread function |
| 1373 | | must be extended objects, and this allows stars to be differentiated |
| 1374 | | from galaxies and other objects. Cosmic rays (single, energetic |
| 1375 | | particles) deposit energy in single pixels; since no distant object |
| 1376 | | can produce such an image profile, this source of noise can easily be |
| 1377 | | filtered from the image. |
| 1378 | | |
| 1379 | | |
| 1380 | | The result of source extraction is a list of object positions (with |
| 1381 | | subpixel resolution), along with the brightness of each object. If |
| 1382 | | the telescope has been properly calibrated, this brightness is |
| 1383 | | proportional to the electromagnetic \emph{flux} from the object, but |
| 1384 | | in general we simply treat the brightness as a means of ordering the |
| 1385 | | objects. It is also possible to classify each source as a star or |
| 1386 | | galaxy, and even classify the galaxy type, but we currently do not use |
| 1387 | | this information. |
| 1388 | | |
| 1389 | | |
| 1390 | | For the purposes of this paper, we will assume that source extraction |
| 1391 | | has been performed as a preprocessing step, so we are given a list of |
| 1392 | | source positions sorted by brightness. |
| 1393 | | |
| 1394 | | |
| 1395 | | |
| 1396 | | \subsection{Quad features} |
| 1397 | | |
| 1398 | | \begin{figure} |
| 1399 | | \begin{center} |
| 1400 | | \includegraphics[width=0.5\textwidth]{quad-fig} |
| 1401 | | \end{center} |
| 1402 | | \caption{The local coordinate frame of a quad. The most distant pair of sources, $A$ and $B$, |
| 1403 | | define the origin and $(1,1)$ of the coordinate frame, respectively. |
| 1404 | | The remaining two sources $C$ and $D$ are required to lie inside the circle with diameter $AB$. |
| 1405 | | To resolve ambiguities, the points must be labelled such that $C_x \le \smallhalf$ and |
| 1406 | | $C_x + D_x \le 1$. The geometric feature descriptor is $(C_x, C_y, D_x, D_y)$.} |
| 1407 | | \label{quad} |
| 1408 | | \end{figure} |
| 1409 | | |
| 1410 | | We have taken a geometric feature-matching approach to this problem. |
| 1411 | | Since individual sources can be well localized but have no distinctive |
| 1412 | | properties, we build distinctive features by describing the relative |
| 1413 | | geometry of small sets of objects. |
| 1414 | | |
| 1415 | | |
| 1416 | | To solve the blind astrometry problem, we need to index a large number |
| 1417 | | of features: we want to solve images that cover $10^{-6}$ or less of |
| 1418 | | the area of the sky, and we want to have tens of features per image in |
| 1419 | | order to increase our robustness to distractors, dropouts, and |
| 1420 | | occlusion. |
| 1421 | | |
| 1422 | | |
| 1423 | | We are matching two-dimensional objects, and we assume that |
| 1424 | | astronomical images can be well described by similarity |
| 1425 | | transformations: scaling, rotation and translation. Therefore, we |
| 1426 | | require two sources to define a local coordinate system. |
| 1427 | | |
| 1428 | | |
| 1429 | | The standard geometric hashing recipe would have us use triangle |
| 1430 | | features. However, the positional noise in astronomical image is |
| 1431 | | sufficiently large that triangles are not distinctive: with the number |
| 1432 | | of features we require, \emph{any} triangle will match an excessively |
| 1433 | | large number of triangles in our index due to feature aliasing. In |
| 1434 | | order to building a more distinctive feature, we consider sets of |
| 1435 | | \emph{four} sources-- \emph{quads}. Two of the sources are the basis |
| 1436 | | set that defines the local coordinate system, and our geometric |
| 1437 | | feature descriptor is composed of the positions of the remaining two |
| 1438 | | sources in this coordinate system. |
| 1439 | | |
| 1440 | | Using a quad feature can be seen as a way of incorporating a voting |
| 1441 | | scheme directly into the feature-matching phase: a quad is essentially |
| 1442 | | the conjunction of two triangles. We are effectively squaring the |
| 1443 | | distinctiveness of our geometric features by demanding that \emph{two} |
| 1444 | | triangles be matched simultaneously. |
| 1445 | | |
| 1446 | | One disadvantage of using quads rather than triangles is that we |
| 1447 | | become more sensitive to distractors and dropouts: if the proportion |
| 1448 | | of distractors and dropouts is $d$, the proportion of overlapping |
| 1449 | | sources is $(1-d)$, and we expect about $(1-d)^3$ triangles but |
| 1450 | | $(1-d)^4$ quads to overlap. |
| 1451 | | |
| 1452 | | |
| 1453 | | In standard geometric feature-matching, one enumerates all basis sets, |
| 1454 | | then lists the positions of all the remaining interest points within |
| 1455 | | that basis. Since we have a two-dimensional basis and we list the |
| 1456 | | positions of pairs of sources, this would yield $\mathcal{O}(n^4)$ |
| 1457 | | features. Astronomical images typically contain hundreds to thousands |
| 1458 | | of sources, and our reference catalog contains $10^9$ sources, so this |
| 1459 | | would be prohibitively expensive. We therefore impose two constraints |
| 1460 | | on the properties of the quads we create. See Figure \ref{quad}. |
| 1461 | | |
| 1462 | | Given a set of four sources, we first find the two sources most |
| 1463 | | distant from each other and label them $A$ and $B$. If the distance |
| 1464 | | from $A$ to $B$ is outside a specified range, we reject the quad. |
| 1465 | | Next, we demand that the two remaining sources lie inside the circle |
| 1466 | | that has $AB$ as a diameter. We place $A$ at the origin of the local |
| 1467 | | coordinate system, and place $B$ at position $(x, y) = (1, 1)$. We |
| 1468 | | then compute the $(x,y)$ coordinates of the two remaining sources. |
| 1469 | | The source with the smaller $x$ coordinate is labelled $C$, and the |
| 1470 | | last source is labelled $D$, so we have $C_x \le D_x$. The geometric |
| 1471 | | feature descriptor is $(C_x, C_y, D_x, D_y)$. |
| 1472 | | |
| 1473 | | Observe that this feature has one redundancy: we labelled $A$ and $B$ |
| 1474 | | arbitrarily, and if we swap their labels then the coordinate system is |
| 1475 | | rotated $180$ degrees about $(\half, \half)$. The feature descriptor |
| 1476 | | becomes $(1 - C_x, 1 - C_y, 1 - D_x, 1 - D_y)$. We can deal with this |
| 1477 | | redundancy by demanding that $A$ and $B$ be labelled such that $C_x + |
| 1478 | | D_x \le 1$. |
| 1479 | | |
| 1480 | | |
| 1481 | | Recall the idea mentioned earlier in this paper of solving images of unknown |
| 1482 | | scale by building a cascade of indices, each covering a range of scales. |
| 1483 | | We use this approach with a scale factor of $2$ or $\sqrt{2}$. |
| 1484 | | This allows us to assume that the image scale is known to |
| 1485 | | within a small factor. When deciding what range of quad sizes ($AB$ distances) we wish to |
| 1486 | | allow, we are faced with a tradeoff: if we use quads with large |
| 1487 | | diameter, then the probability that all four stars will be within the bounds of |
| 1488 | | the image decreases. Meanwhile, the number of quads of diameter $d$ that |
| 1489 | | can be created is $\mathcal{O}(d^5)$, assuming a uniform distribution |
| 1490 | | of sources and ignoring edge effects.\footnote{% |
| 1491 | | For a given star $A$, the number of |
| 1492 | | stars that can become star $B$ goes as $d$, while the number of stars |
| 1493 | | $C$ and $D$ goes as $d^2$, so the number of $BCD$ combinations goes as $d^5$.} |
| 1494 | | On the other hand, if we build small quads, then the relative positional error |
| 1495 | | will be larger, which means that when we search for matching quads in the index |
| 1496 | | we must use a larger search radius, which means that many more false matches |
| 1497 | | will be found. We have not yet found an analytic solution to optimize this |
| 1498 | | tradeoff; we currently set the desired quad scale to about |
| 1499 | | a third of the size of the images we wish to solve. |
| 1500 | | |
| 1501 | | |
| 1502 | | Our quad features have good noise properties: a small change of |
| 1503 | | position of any one of sources $ABCD$ leads to a small change in the |
| 1504 | | feature descriptor. The noise is nearly isotropic: movement of a |
| 1505 | | source in any direction yields an almost equal change in the feature |
| 1506 | | descriptor. The feature has some edge effects: $A$ and $B$ are |
| 1507 | | defined to be the most distant pair, but if $C$ and $D$ are nearly |
| 1508 | | diametrically opposite each other then a small amount of positional |
| 1509 | | noise can push them further apart so that they become the most distant |
| 1510 | | pair. As a result, $C$ and $D$ will be labelled $A$ and $B$, and vice |
| 1511 | | versa, and a completely different geometric feature will be produced. |
| 1512 | | This is a minor effect since it only occurs when both $C$ and $D$ have |
| 1513 | | extreme values; we ignore it. Another edge effect is that $C$ and $D$ |
| 1514 | | are required to lie inside the circle with diameter $AB$; at test time |
| 1515 | | we simply expand the allowed diameter of the quad enough to compensate |
| 1516 | | for a given level of positional noise. Finally, if the distance $AB$ |
| 1517 | | is near the limit of the allowed range of scales, then the quad can |
| 1518 | | become invalid due to noise. This effect is equally easy to correct |
| 1519 | | by expanding the range slightly at test time. |
| 1520 | | |
| 1521 | | |
| 1522 | | Notice that our constraints on the quads we accept has the effect of |
| 1523 | | restricting the regions of feature space in which valid quads live. |
| 1524 | | Since $C$ and $D$ are required to fall inside a circle with $AB$ on |
| 1525 | | the diameter, we have: |
| 1526 | | \begin{eqnarray*} |
| 1527 | | \left(C_x - \smallhalf\right)^2 + \left(C_y - \smallhalf\right)^2 & \le & \smallhalf \\ |
| 1528 | | \left(D_x - \smallhalf\right)^2 + \left(D_y - \smallhalf\right)^2 & \le & \smallhalf |
| 1529 | | \end{eqnarray*} |
| 1530 | | We also have the constraints $C_x + D_x \le 1$ and $C_x \le D_x$ to |
| 1531 | | eliminate the redundancy of the labelling of the source pairs $AB$ and |
| 1532 | | $CD$. As a result, our feature space is rather strangely shaped, so a |
| 1533 | | grid-based hashing scheme would result in poor utilization of the hash |
| 1534 | | table storage. Also note that since our noise is nearly isotropic and |
| 1535 | | our quads are restricted to a known range of scales, we can use the |
| 1536 | | Euclidean distance metric with a fixed search radius $r$ to do quad |
| 1537 | | feature matching. |
| 1538 | | |
| 1539 | | |
| 1540 | | |
| 1541 | | |
| 1542 | | \subsection{Indexing} |
| 1543 | | |
| 1544 | | |
| 1545 | | We want our index to satisfy the following constraints. It should |
| 1546 | | have uniform source and quad density over the whole sky, because we |
| 1547 | | want to be equally likely to solve images over the whole sky. We want |
| 1548 | | to build quads from bright stars, because we expect these to be more |
| 1549 | | likely to be visible in test images, yet because of distractors, |
| 1550 | | dropouts, and occlusions, we don't want to depend too heavily on any |
| 1551 | | one star; we want to trade off using bright stars and using a variety |
| 1552 | | of stars. |
| 1553 | | |
| 1554 | | |
| 1555 | | We achieve this by first selecting a spatially uniform subset of the |
| 1556 | | brightest sources in the reference catalog. We place a grid over the |
| 1557 | | sky whose cells are sized so that a typical image is a few cells wide, |
| 1558 | | then we take the $n$ brightest sources in each cell. |
| 1559 | | |
| 1560 | | |
| 1561 | | Next, we place another similarly-sized grid over the sky and in each |
| 1562 | | cell we try to build a quad whose center is within the cell. We try |
| 1563 | | building quads in order of source brightness, but we also keep track |
| 1564 | | of how many times each source has been used. In our first pass |
| 1565 | | through the cells we only allow each star to be used once, and in each |
| 1566 | | subsequent pass we increase the number of times a source can be used. |
| 1567 | | |
| 1568 | | |
| 1569 | | Once we have built all the quads, we build a \kdtree in feature space. |
| 1570 | | This allows us to rapidly find all the quads in our index whose |
| 1571 | | geometric feature descriptor is close to a query feature. |
| 1572 | | |
| 1573 | | We also build a \kdtree over the spatially uniform subset of sources; |
| 1574 | | this is used for our verification scheme. |
| 1575 | | |
| 1576 | | |
| 1577 | | \subsection{Solving} |
| 1578 | | |
| 1579 | | At test time, we begin enumerating all valid quads in the image, |
| 1580 | | starting with those that can be built from the brightest sources. For |
| 1581 | | each quad we perform a range search in the \kdtree in order to find |
| 1582 | | matching quads. Typically each image quad results in a few matches. |
| 1583 | | We run a verification process on each match. |
| 1584 | | |
| 1585 | | The verification process first computes a least-squares transformation |
| 1586 | | between image coordinates and celestial coordinates. We then search |
| 1587 | | in the \kdtree of sources for all sources that should be visible in |
| 1588 | | the image. We apply the transformation to bring them into image |
| 1589 | | coordinates, and for each one find the nearest source in the image. |
| 1590 | | We use a probabilistic model of the positional noise to compute the |
| 1591 | | odds ratio between the match being a true positive and a false |
| 1592 | | positive. If this ratio exceeds a threshold we declare the image |
| 1593 | | solved. By using an odds ratio threshold, we allow the user to trade |
| 1594 | | off the probability of false negatives and false positives. |
| 1595 | | |
| 1596 | | |
| 1597 | | |
| 1598 | | \comment{ |
| 1599 | | We build an \emph{index} which maps these features |
| 1600 | | back to their locations on the sky.\footnote{% |
| 1601 | | In the database community (and some early computer science literature) |
| 1602 | | this would be called an \emph{inverted} index, but we feel this |
| 1603 | | is contrary to common usage: the index of a book maps from concepts |
| 1604 | | (features) to pages (locations); the inverse mapping is called a |
| 1605 | | \emph{table of contents}.} |
| 1606 | | |
| 1607 | | At test time, we extract features from the image, then search our |
| 1608 | | index for matching features. Each matched feature is a hypothesized |
| 1609 | | astrometric solution. We use a separate probabilistic verification |
| 1610 | | stage to test these hypotheses. |
| 1611 | | |
| 1612 | | |
| 1613 | | Since our input images have unknown scale, rotation, and |
| 1614 | | translation, we want our features to be invariant under these |
| 1615 | | transformations. One option is to consider triangles of stars and |
| 1616 | | use the interior angles (or relative edge lengths) as a |
| 1617 | | two-dimensional feature descriptor, as in \cite{pal2006}. Since |
| 1618 | | both the image and reference catalog are noisy, we must allow each |
| 1619 | | feature in the image to match a region in feature space. When using |
| 1620 | | geometric features, we are faced with a tradeoff: if we consider |
| 1621 | | sets of stars that span a large portion of the image, then the |
| 1622 | | relative positional error is smaller (and hence each feature matches |
| 1623 | | a smaller region of feature space - features are more distinctive), |
| 1624 | | but there are \emph{many} features to test. If we choose smaller |
| 1625 | | features, there are fewer features to test, but the relative error |
| 1626 | | is larger so each feature matches a larger region in feature space; |
| 1627 | | features can be confused more easily. |
| 1628 | | |
| 1629 | | |
| 1630 | | Although triangles are appealing features from a theoretical |
| 1631 | | perspective, in practice they are problematic due to the magnitude |
| 1632 | | of positional noise in typical images we wish to solve. Any |
| 1633 | | triangle we find in an image will generate a large number of false |
| 1634 | | matches to triangles in our index because triangles aren't |
| 1635 | | distinctive enough. We therefore use sets of four stars -- |
| 1636 | | \emph{quads} -- as our features. The two most widely separated |
| 1637 | | stars (A and B) define the corners of a coordinate system, and we |
| 1638 | | list the positions of the remaining two stars (C and D) in that |
| 1639 | | coordinate system. We demand that stars C and D lie inside the |
| 1640 | | circle centered at the midpoint of stars A and B and diameter AB. |
| 1641 | | We impose an arbitrary order to decide which star is labelled C and |
| 1642 | | which is labelled D, based on their positions in the local |
| 1643 | | coordinate system. This yields a four-dimensional feature |
| 1644 | | descriptor. If we exchange the labels of A and B, then the local |
| 1645 | | coordinate system is flipped and a different feature descriptor is |
| 1646 | | produced. We store one of these descriptors in our index, and at |
| 1647 | | test time we search for matches to both descriptors. |
| 1648 | | |
| 1649 | | |
| 1650 | | This choice of feature and feature descriptor (\emph{code}) has |
| 1651 | | several desirable properties. If stars are distributed uniformly in |
| 1652 | | space, then codes are distributed uniformly in \emph{code space}. |
| 1653 | | The error propagation properties are good: small positional errors |
| 1654 | | in any of the four stars yield small errors in the code value, and |
| 1655 | | isotropic errors in the positions of the stars produce nearly |
| 1656 | | isotropic errors in the code value. This allows us to find matching |
| 1657 | | codes based on Euclidean distance in code space. |
| 1658 | | |
| 1659 | | Our index is composed of a large number of code values and mappings |
| 1660 | | back to locations on the sky. When given a new image to solve, we |
| 1661 | | extract many quads, compute their codes, and for each one search for |
| 1662 | | matching codes in the index. Due to the isotropic error properties |
| 1663 | | of codes, searching for matching codes is simply \emph{range search} |
| 1664 | | in four-dimensional code space. We use \kdtrees to accelerate this |
| 1665 | | process. |
| 1666 | | |
| 1667 | | |
| 1668 | | Each time we match a feature in an image to a feature in our index, |
| 1669 | | we generate a hypothesis about the location (pointing, scale, and |
| 1670 | | rotation) of the image on the sky. We are typically operating in |
| 1671 | | the regime where \emph{any} feature will produce several matches, so |
| 1672 | | we use a probabilistic model to determine whether a given match is a |
| 1673 | | false positive or a true solution. We do this by computing a |
| 1674 | | transformation between image coordinates and the sky based on the |
| 1675 | | matching quad. We then search for stars in our reference catalog |
| 1676 | | that should be visible in the image if the hypothesis is correct. |
| 1677 | | By combining our knowledge of the positional errors in the reference |
| 1678 | | catalog and the image with an estimate of the proportion of stars |
| 1679 | | that should be visible in both, we can assign a probability to the |
| 1680 | | hypothesis. If it is above some threshold, we accept the match. |
| 1681 | | |
| 1682 | | |
| 1683 | | \emph{ |
| 1684 | | \begin{itemize} |
| 1685 | | \item Brightness ordering to speed search only. |
| 1686 | | \item More about the constraints and goals of our system to differentiate it |
| 1687 | | from the spacecraft attitude estimation problem. |
| 1688 | | \item A quad is like the conjunction of two triangles: we're doing implicit agreement of two triangle features! |
| 1689 | | \item Applications of our system? |
| 1690 | | \end{itemize} |
| 1691 | | } |
| 1692 | | |
| 1693 | | |
| 1694 | | To summarize, the key elements of our approach are: |
| 1695 | | \begin{itemize} |
| 1696 | | \item \textbf{geometric features} composed of sets of point features; |
| 1697 | | \item \textbf{continuous feature descriptors} that describe the relative positions of the point features; |
| 1698 | | \item \textbf{fast range-search} to match features descriptors; |
| 1699 | | \item \textbf{probabilistic validation} to eliminate false positives. |
| 1700 | | \end{itemize} |
| 1701 | | |
| 1702 | | |
| 1703 | | |
| 1704 | | %rummage through the toolbox of computer vision and object |
| 1705 | | %recognition techniques that are related to our approach or would seem |
| 1706 | | %to be applicable to our problem. |
| 1707 | | } |
| 1708 | | |
| 1709 | | |
| 1710 | | |