How (the heck) does Astrometry.net work?

Question

What magic is this?

This answer to my question Astronomy detective question: what part of the sky are these photos of? What are a few of the stars? nails it, and the output (linked there) contains the following images along with a whole lot more information.

The images are analyzed by astrometry.net

How does it do what it does? If I understand correctly you can upload random star images with a wide variety of fields of view (from telescopes to cell phones) and random orientations and no further information, and it figures it out.

Is there a way to explain how the algorithm works with more detail than "a computer matches the dots to known stars"?

Related question about these cell phone images: Milky Way with a cell phone, how is this possible?

Answer 1

I've been trying to figure out the technical details of astrometry.net for quite some time. As others already pointed out, the main input to the whole process is a list of stars. I will not go into details on how astrometry.net does it, just note that you can either use its internal simplexy algorithm or use SExtractor. In the end you need a list of coordinates for stars (plus optional flux/intensity/brightness and background).

The .xyls files are used to store these input values (small example):

 x-coord(px)   y-coord(px)        flux  background
 1008.911987,   557.925659,  10.556271,  32.320175,
 1449.509277,   643.280212,   6.580036,  27.963276,
  185.978119,  1253.869751,   5.525373,  27.713015,

Background is approximated before/during star extraction (SExtractor has a function for that).

Now the algorithm goes into creating quads. The cited paper always talks about 4 stars that form a quad. But by my own experiments they seem to be just triangles. Maybe they should have better named them "Asterism" (the code has a lot of very bad naming, as we soon will discover, took me hours to figure out).

So we have a triangle and want to know if a similar triangle exists anywhere in the selected astrometry.net index. So this is where the real magic happens IMO. The papers talks about geometrically hashed lookup. IMO the word "hash" is a bit misleading here as it has not much in common with real hashes. Real hashes should normally give very different results for small changes in the input (e.g. checksum or hash-table bucket distributions). I would rather say that they perform a geometric transformation where the result must be a one-dimensional value (e.g. a double value).

Abstractly speaking we want to map numerous properties to a single value. The value should vary very little if the overall properties also vary little, so from looking at the delta of both translation results we can deduct how similar the two objects were. One such property could e.g. be the angle between two lines.

From here it should be "obvious" how one can use this approach to drive a search. In this simple one-dimensional case we could e.g. use a binary search. For astrometry.net this problem gets technically more complicated, as we need to search for two "hashes". This is done with a KD-Tree. The one used by astrometry.net is optimized for the pre-built index files, so they can be accessed very fast and without much (memory) overhead. Basically a KD-Tree can optimize the question "give me the closest point(s) to X/Y".

In the picture below I tried to visualize how we can search for similar triangles by reducing the question to two numbers. We basically search for the normalized blue vectors (or ones that are pretty close). It should be obvious that this eliminates any rotation and scaling in the question asked, so any similarly shaped triangle will match, regardless of orientation or size. IMO it simply boils down to the fact that the shape of a triangle is defined by two parameters, e.g. not sure why they didn't use two angles.

In reality we have to make the lookup multiple times to also search for flipped and/or inverted variations. After a similar triangle is found (this will happen a lot), the process goes into the verification step. The resulting triangle will give the program the hint on how all other stars must be translated to match the given triangle. With that it can try to match all other stars with its know star catalog.

Query images may contain some extra stars that are not in your index catalogue, and some catalogue stars may be missing from the image.

These can be seen as distractors and conflicts? in the debug output of the program:

verify: logodds -1.38629, 0 matches, 0 conflicts, 1 distractors after 0 field objects.
verify: logodds 333.123, 99 matches, 0 conflicts, 115 distractors after 213 field objects.
108 matches, 424 distractors, 2 conflicts (all sources)

I believe these numbers are not "really" accurate as the checker will probably bail out once it thinks it is impossibly a match (e.g. after certain distractors are found without any matches). Anyway, this is more or less to whole process in detail. I left out some of the magic, e.g. there are a lot of probability checks for speed. Also the way potential quads are chosen from all stars is quite elaborate in reality (and should probably match closely to how this was done when the index files were created).

Side note on indexes: As they contain x/y positions for known stars and stars wander a little over time (a few faster, most very slow), the index can get outdated and start to not match (no idea if this would be in a decade or a millenium). Regular star catalogs give x/y position at a given epoch time plus x/y-speeds to calculate the actual position in any given point of time. IMO with modern CPUs this can be done in seconds for millions of stars (so time independent indexes could be doable, although probably not for the initial triangle matching, but surely for the verify phase). Also with the new ESO Gaia data releases there shouldn't be any gaps anymore (as noted by astrometry.net as small holes in the USNO-B1.0 catalog). But it seems the official available indexes haven't been updated yet to use the new gaia catalog.

Disclaimer: This knowledge was acquired mostly by reading and testing with the astrometry.net source code. So any conclusions I made could be wrong. But I would say it all makes sense if put together. Below I'll give a few more details into the actual implementation inside astrometry.net.

Edit: After reading https://iopscience.iop.org/article/10.1088/0004-6256/139/5/1782 I came to the conclusion that real quads with four stars are probably used in indexes with smaller view angles (zoomed in). I used a picture taken with a regular 55mm DSLR lens for my tests. This would make sense and it basically is exactly what an n-dimensional tree is made for (the question for the closest neighbors now involves 4 parameters).

Each index_t contains two KD-Trees, namely codekd and starkd. The tree codekd contains the information of all "quads", while starkd contains regular star coordinates for later verification.

The "hashing" bits can interestingly be found in "solver.c" in the function "check_inbox".

static void check_inbox(pquad* pq, int start, solver_t* solver) {
    int i;
    double Ax, Ay;
    field_getxy(solver, pq->fieldA, &Ax, &Ay);
    // check which C, D points are inside the circle.
    for (i = start; i < pq->ninbox; i++) {
        double r;
        double Cx, Cy, xxtmp;
        double tol = solver->codetol;
        if (!pq->inbox[i])
            continue;
        field_getxy(solver, i, &Cx, &Cy);
        Cx -= Ax;
        Cy -= Ay;
        xxtmp = Cx;
        Cx = Cx * pq->costheta + Cy * pq->sintheta;
        Cy = -xxtmp * pq->sintheta + Cy * pq->costheta;

        // make sure it's in the circle centered at (0.5, 0.5)
        // with radius 1/sqrt(2) (plus codetol for fudge):
        // (x-1/2)^2 + (y-1/2)^2   <=   (r + codetol)^2
        // x^2-x+1/4 + y^2-y+1/4   <=   (1/sqrt(2) + codetol)^2
        // x^2-x + y^2-y + 1/2     <=   1/2 + sqrt(2)*codetol + codetol^2
        // x^2-x + y^2-y           <=   sqrt(2)*codetol + codetol^2
        r = (Cx * Cx - Cx) + (Cy * Cy - Cy);
        if (r > (tol * (M_SQRT2 + tol))) {
            pq->inbox[i] = FALSE;
            continue;
        }
        setx(pq->xy, i, Cx);
        sety(pq->xy, i, Cy);
    }
}

The pquads are the potential "asterisms" we are currently creating, and this function's job is to set the "query" parameters via "set[xy]" on the bottom. Those are the actual values later looked up in the codekd tree.

A pquad is basically a line between two points (also called backbone-stars) called fieldA and fieldB (which are actually indexes to get an x/y position). Additionally it must have at least one additioanl xy point (exactly one in our triangular case). The pquad also contains sintheta and costheta (set in check_scale).

    double dx, dy;
    dx = field_getx(s, pq->fieldB) - field_getx(s, pq->fieldA);
    dy = field_gety(s, pq->fieldB) - field_gety(s, pq->fieldA);
    pq->scale = dx*dx + dy*dy;
    pq->costheta = (dy + dx) / pq->scale;
    pq->sintheta = (dy - dx) / pq->scale;

As we see both code parts translate and scale the vector AC in relation to the line AB.

I hope this info is useful to somebody, even if it got a bit long!

Answer 2

Finding an astrometric solution from an image with Astrometry.net is usually called plate solving. As mentioned in the comments, it is based on pattern matching, using a large set of databases that are pre-computed for various field of view and plate (or pixel) scale. The ArXiv paper Astrometry.net: Blind astrometric calibration of arbitrary astronomical images provides some details of how it works. From their abstract ...

After robust source detection is performed in the input image, asterisms (sets of four or five stars) are geometrically hashed and compared to pre-indexed hashes to generate hypotheses about the astrometric calibration. A hypothesis is only accepted as true if it passes a Bayesian decision theory test against a null hypothesis.

The blind solution described in the paper is obtained directly from a photo without any additional input, but even a rough guess for the field of view can provide considerable speedup by narrowing the number of databases to be searched. Additional speedup is possible if the scale (arcsec/pixel), FOV and approximate coordinates in the sky are provided.

Astrometry.net code and databases can also be downloaded and run locally on unix/linux computers, and on Windows using Cygwin. There are several other plate solvers available, often included with, or designed to interface with, software that controls a telescope, mount, camera, etc. for astrophotography.

Edit: Some additional info related to the comments.

In finding a solution, they must deal with imperfections of real images, including defects, optical and atmospheric distortion, etc. along with small positional changes of stars. Roughly speaking, they obtain a solution by using a best fit to several asterisms, each containing a small number of stars.

Noise in the image and distortion caused by the atmosphere and telescope optics lead to noise in the measured positions of stars in the image. In general this noise causes the stars in a quad to move slightly with respect to each other, which yields small changes in the hash code (i.e., position in code space) of the quad. Therefore, we must always match the image hash code with a neighborhood of hash codes in the index.

Also, they are able to ignore defects and moving objects like planets, asteroids, etc. in finding matches.

However, we also know that some fraction of the stars in the query image will have no counterpart in the index, due to occlusions or artifacts in the images, errors in star detection or localization, differences in the spectral bandpass, or because the query image "star" is actually a planet, satellite, comet, or some other non-star, non-galaxy object. True stars can be lost, and false stars can be added. Our foreground model is therefore a mixture of a uniform probability that a star will be found anywhere in the image--a query star that has no counterpart in the index--plus a blob of probability around each star in the index, where the size of the blob is determined by the combined positional variances of the index and query stars.

Their goal was not to achieve high astrometric precision.

There are several automated calibration systems that refine the astrometric calibration of an image to produce a high-precision alignment to a reference catalog given a good first guess (for example, Valdes et al. 1995; Mink 2006; Bertin 2005). These systems are reliable and robust, but they require a reasonable first guess about the image pointing, orientation, and scale. Our system can be used to create that good first guess.