Wandering star - codeabbey task

Question

I'm trying to solve this problem and I'm not sure what to do next.

Link to the problem
Problem statement:
Suppose that some preliminary image preprocessing was already done and you have data in form of coordinates of stars on two pictures. These pictures are about 100x100 millimeters, and coordinates are also given in millimeters relative to their center. Look at the schematic reprsentation below: You see that in both pictures stars are shown in roughly circular area (think of it as of apperture of our telescope), and you can find out that they represent the same piece of the sky - thoulgh slightly rotated and slightly shifted.

You also can see that one of the stars (marked with red arrows) have changed its position relative to others.

Your task is to find out such a "wandering star" for it could be comet or asteroid with high probability.

Note that some stars which are close to the edge could be absent from one of pictures (due to shift) - but the "wandering star" is not that far from the center and therefore is presented on both of images.

Input data contains two sections corresponding to two images. Each sequence starts with a single integer - amount of stars listed. Then the coordinates (X and Y) of stars follow.

Answer should give two indexes (0-based) of the wandering star in the first and second section respectively.

The example is the same as pictures above. The star #29 from the first section with coordinates (-18.2, 11.1) is the same as the star #3 from the second section with coordinates (-19.7, 6.9). Example input data:
94 # section 1 contains 94 stars
-47.5 -10.4
19.1 25.9
18.9 -10.4
-2.1 -47.6
...
...
92 # section 2 contains 92 stars
-14.8 10.9
18.8 -0.1
-11.3 5.7
-19.7 6.9
-11.5 -16.7
-45.4 -15.3
6.0 -46.9
-24.1 -26.3
30.2 27.4
...
...

The problem I'm facing
The problem is that the vectors do no match, and don't even have the same size. So for example the first vector in the first section doesn't match the first in the second section, so I can't calculate the rotation matrix based on that. I also tried calculating it based on the centroids of each section, but some points on the edge might be absent, so they will have different centroids(I tried including only the vectors who's length is < 40, the size still doesn't match).

So my question is what should I base my calculations on? How do I find some matching vectors so I can calculate the rotation matrix on them? I need a push in the right direction.

What I did is implement the functions to find the rotation matrix between two given vectors. The formula I'm using:
transformed_vector = R * original_vector + t
where R is the rotation matrix and since the vector also moves along the axis a bit I'm also adding t
Now all I need is two vectors to base my calculations on.

Edit: I should probably mention that I'm actually given two arrays of vectors, one for each image, I'm not actually given the images. I need to find the star that moved based on those vectors.

Thanks!

Answer 1

First, for simplicity let's initially pretend that there is no wandering star, and that every star that is in one image is in the other -- that is, that image B has been produced from image A by applying just a (possibly zero) shift and a (possibly zero) rotation.

Signatures

I think the best way to approach this is to "forget about" the position of each star in each image, and instead record for each star only a signature: some information that will not be changed by a shift (translation) or rotation. Then you can look at every signature in image A, and try to match it with a signature in image B. If we make the signatures in a sensible way, then they should do a good job of distinguishing the stars, so that no two stars in the same image get the same or highly similar signatures.

For a given star, something that is not changed by shifting or rotation is its distance to any other star in the same image. So you could use the complete set of distances to all the other n-1 stars as a star's signature (stored as an unordered set, since we don't know "which stars are which"). But it's probably sufficient (and faster, and ultimately more robust) to just use a subset of these, and in that case the distances to its nearest neighbours are probably the most informative (unless you have repeating patterns of stars appearing). How many to use? I would start by calculating the distance of each star to its nearest neighbour in the same image; if in image A, all of these these distances are "different enough" (i.e., have difference below some arbitrary threshold that you choose), and likewise in image B, then you're done -- that is, the signatures that you have calculated already distinguish the stars well enough -- otherwise go back and for each star insert the distance to its second-nearest neighbour into its signature, repeating this until "uniqueness" of signatures is achieved.

Comparing signatures

This brings up the question of how to decide whether two signatures are "different enough" when a signature contains more than one distance. Now you are comparing two sets of numbers instead of just two numbers. One way would be to sort the distances in each signature, and then compare them by calculating some measure of distance between points like the sum of squared differences (e.g., if each signature contains 3 distances, this would be (u[1]-v[1])^2 + (u[2]-v[2])^2 + (u[3]-v[3])^2 for two stars u and v, with u[i] being the distance from star u to its ith-nearest neighbour in the same image, and likewise for v.) On your example dataset, I would guess that probably 3-4 distances per star would be enough.

Robust and efficient signature comparison

However this will turn out to be not particularly robust once wandering and missing stars are considered: Suppose that some star u in image A has the wandering star v as its nearest neighbour, and in image B, v has moved away so that it is now the 5th-nearest neighbour to u. Then if we compare u's signature in image A with its signature in image B, we will get an incorrectly large distance, because we will be "blindly" comparing uA[1] with uB[1] and so on, instead of comparing uA[2] with uB[1], uA[3] with uB[2] and so on as we would like to (since these distances will be equal). So a better way to compare two signatures would be by taking the ordered list of distances of each one and allowing the numbers to "slide along until they fit" the numbers in the other one.

Happily this can be done in linear time with a list merge. Given two sorted lists of numbers X and Y, we perform the usual list merge strategy of choosing the smallest element from either list, removing it from there and writing it out, but we also keep track of two candidate partial solutions: Bm, the score of the best partial solution in which the most recently output number is matched with an earlier number, and Bf, the score of the best partial solution in which it is free (not matched with anything yet). Choose a penalty cost p to assign to numbers that are not matched with any number in the other list, and a scoring function score(a, b) that assigns a value of 0 when a = b and higher values the more different a and b are (e.g., score(a, b) = (a - b)^2). Then the optimal score sigDiff(X, Y) can be computed using:

1. Set Bm = Bf = 0.
2. If X and Y are both empty, then halt. The final score is min(Bf + p, Bm).
3. Call the element at the front of X x, and the element at the front of Y y. Let z be the smaller of x and y, and set whichList to the ID of the list (X or Y) that z came from.
4. Remove the first element (z) from the list named by whichList.
5. If whichList == prevWhichList (i.e. if the previous smallest number, prev, was also from the same list as z), or if this is the first iteration, then set newBm = min(Bf + p, Bm) + p.
6. Otherwise, it's possible to match z with the previous number written out, so set newBm = Bf + score(z, prev).
7. Set Bf = min(Bf + p, Bm).
8. Set Bm = newBm.
9. Set prev = z and prevWhichList = whichList.
10. Goto 2.

For example, suppose we have the two distance lists (i.e., two signatures) X = [10, 20, 30, 40] and Y = [11, 18, 41, 50], the penalty is 20, and score() is the sum of squared differences as above. The above algorithm will "align" the two sequences as follows:

X        10          20   30   40
matches    \        /            \
Y           11    18              41   50
cost         1     4      20       1   20      Total cost: 46

By contrast, the "naive" sum-of-squares method would give:

X        10   20   30   40
matches   |    |    |    |
Y        11   18   41   50
cost      1    4  121  100                     Total cost: 226

Matching signatures between the two images

We have built all this machinery so that we can try to match stars between the two images by matching their signatures. Arguably the "right" way to do this would be to solve the bipartite maximum weighted matching problem on the graph in which the vertices in one side are the signatures of the stars in image A, the vertices in the other side are the signatures of the stars in image B, and there is an edge from every vertex X in one side to every vertex Y in the other side having weight -sigDiff(X, Y) (negated because we want to minimise the total difference).

However, solving this might take a while, and it's also possible that this algorithm will find some incorrect "approximate" matches in an attempt to get an overall minimum cost. We would prefer to focus on just those pairs of stars which we are certain correspond to one another. That can be done much faster with a heuristic based on the idea that if, for a star X in image A whose closest star in image B (based on sigDiff()) is Y, it turns out that Y's closest star in image A is also X (i.e., if the best match of X's best match is X), then we can be confident that X and Y really match each other. These confident matching pairs can be found quickly, and used to estimate an affine transform from image A to image B using least squares.

Ignoring missing stars

We first compute the convex hull of image B's stars. Then, using the transform determined from the confident star pairs, we transform every star in image A to its estimated corresponding location in image B, and take the convex hull of this transformed image. We then intersect these two convex hulls: every star in either image that falls inside this intersection is a star that we expect to find in both images. Finally we throw away all stars from either image that are outside this intersected convex hull, and start again!

After running everything again (maybe it's not necessary to rerun everything, actually), we should find that there is exactly 1 star in image A and 1 star in image B that have poor similarity to all other stars in the other image (as measured by sigDiff(), as usual). These "two" stars are of course the single "wandering" star :)

Answer 2

[edit2] complete reedit

Have found some time/mood for this to make it more robust

• let xy0[],xy1[] be the input star lists
• let max_r be the nearby search area treshld
• let max_err be the max acceptable cluster match error

so here is the Algorithm:

1. sort xy0[] by x asc
   • this makes the search faster and easier
2. find star clusters in xy0[]
   • loop through all stars
   • and cross reference them with stars nearby
   • because of sorting the nearby stars will be also near current star index
   • so just search close area before and after this star in array
   • until x-distance cross the max_r
   • add cluster to cl0[] cluster list if found
   • (cluster is 2 and more close stars)
   • before adding new cluster
   • check if no cluster is near it
   • and merge if too close to another cluster instead
3. full recompute found clusters
   • avg coordinate
   • distances between all stars inside
   • sort them by distance asc
4. do the 1.,2.,3. also for xy1[],cl1[]
5. find matches between clusters
   • so check if the distances list inside are the same
   • (remember the lowest sum of abs error)
   • if error is bigger then max_err reject this cluster match
   • this is strong matching have tested on many clusters (big max_r) without miss match on this dataset
6. take 2 clusters from found clusters in cl0[] that have match found
7. take also the matched clusters
8. compute transformation between xy0[],xy1[] from this 4 points
   • I use avg coordinate of cluster and it matches perfectly

This is the visual result:

• 
• left side is xy0[] set
• middle is xy1[] set
• and on the right both where blue bigger dots are xy0[]
• and green smaller dots are transformed xy1[]
• the numbers are the error of cluster match (-1 means no match found)

[notes]

I use my own List<T> template ...

• it is just dynamically reallocating linear array
• List<int> x; is the same as int x[];
• where x[i] is item access
• x.num is number of items in the array
• x.add(5); is the same as x[x.num]=5; x.num++;

From this point you can check for matches between xy0 and transformed xy1

• so flag matched stars to avoid duplicate use of them
• use some treshold for this like max_err
• from what stars are left
• find two that are closest to each other
• this should be the wandering star ... have fun
• (you can sort the transformed xy1 again)
• do not forget to use the original star indexes ix0[],ix1[] for result output

[edit3] also the rest works

//---------------------------------------------------------------------------
// answer: 29 3
// input data:
const int n0=94; double xy0[n0][2]=
    {
    -47.5,-10.4,19.1,25.9,18.9,-10.4,-2.1,-47.6,41.8,-12.1,-15.7,12.1,-11.0,-0.6,
    -15.6,-7.6,14.9,43.5,16.6,0.1,3.6,-33.5,-14.2,20.8,17.8,-29.8,-2.2,-12.8,
    44.6,19.7,17.9,-41.3,24.6,37.0,43.9,14.5,23.8,19.6,-4.2,-40.5,32.0,17.2,
    22.6,-26.9,9.9,-33.4,-13.6,6.6,48.5,-3.5,-9.9,-39.9,-28.2,20.7,7.1,15.5,
    -36.2,-29.9,-18.2,11.1,-1.2,-13.7,9.3,9.3,39.2,15.8,-5.2,-16.2,-34.9,5.0,
    -13.4,-31.8,24.7,-29.1,1.4,24.0,-24.4,18.0,11.9,-29.1,36.3,18.6,30.3,38.4,
    4.8,-20.5,-46.8,12.1,-44.2,-6.0,-1.4,-39.7,-1.0,-13.7,13.3,23.6,37.4,-7.0,
    -22.3,37.8,17.6,-3.3,35.0,-9.1,-44.5,13.1,-5.1,19.7,-12.1,1.7,-30.9,-1.9,
    -19.4,-15.0,10.8,31.9,19.7,3.1,29.9,-16.6,31.7,-26.8,38.1,30.2,3.5,25.1,
    -14.8,19.6,2.1,29.0,-9.6,-32.9,24.8,4.9,-2.2,-24.7,-4.3,-37.4,-3.0,37.4,
    -34.0,-21.2,-18.4,34.6,9.3,-45.2,-21.1,-10.3,-19.8,29.1,31.3,37.7,27.2,19.3,
    -1.6,-45.6,35.3,-23.5,-39.9,-19.8,-3.8,40.6,-15.7,12.5,-0.8,-16.3,-5.1,13.1,
    -13.8,-25.7,43.8,5.6,9.2,38.6,42.2,0.2,-10.0,-48.6,14.1,-6.5,34.6,-26.8,
    11.1,-6.7,-6.1,25.1,-38.3,8.1,
    };
const int n1=92; double xy1[n1][2]=
    {
    -14.8,10.9,18.8,-0.1,-11.3,5.7,-19.7,6.9,-11.5,-16.7,-45.4,-15.3,6.0,-46.9,
    -24.1,-26.3,30.2,27.4,21.4,-27.2,12.1,-36.1,23.8,-38.7,41.5,5.3,-8.7,25.5,
    36.6,-5.9,43.7,-14.6,-9.7,-8.6,34.7,-19.3,-15.5,19.3,21.4,3.9,34.0,29.8,
    6.5,19.5,28.2,-21.7,13.4,-41.8,-25.9,-6.9,37.5,27.8,18.1,44.7,-43.0,-19.9,
    -15.7,18.0,2.4,-31.6,9.6,-37.6,15.4,-28.8,43.6,-11.2,4.6,-10.2,-8.8,38.2,
    8.7,-34.6,-4.7,14.1,-1.7,31.3,0.6,27.9,26.3,13.7,-1.2,26.3,32.1,-17.7,
    15.5,32.6,-14.4,-12.6,22.3,-22.5,7.0,48.5,-6.4,20.5,-42.9,4.2,-23.0,31.6,
    -24.6,14.0,-30.2,-26.5,-29.0,15.7,6.0,36.3,44.3,13.5,-27.6,33.7,13.4,-43.9,
    10.5,28.9,47.0,1.4,10.2,14.0,13.3,-15.9,-3.4,-25.6,-14.7,10.5,21.6,27.6,
    21.8,10.6,-37.8,-14.2,7.6,-21.8,-8.6,1.3,6.8,-13.3,40.9,-15.3,-10.3,41.1,
    6.0,-10.8,-1.5,-31.4,-35.6,1.0,2.5,-14.3,24.4,-2.6,-24.1,-35.3,-29.9,-34.7,
    15.9,-1.0,19.5,7.0,44.5,19.1,39.7,2.7,2.7,42.4,-23.0,25.9,25.0,28.2,31.2,-32.8,
    3.9,-38.4,-44.8,2.7,-39.9,-19.3,-7.0,-0.6,5.8,-10.9,-44.5,19.9,-31.5,-1.2,
    };
//---------------------------------------------------------------------------
struct _dist                        // distance structure
    {
    int ix;                         // star index
    double d;                       // distance to it
    _dist(){}; _dist(_dist& a){ *this=a; }; ~_dist(){}; _dist* operator = (const _dist *a) { *this=*a; return this; }; /*_dist* operator = (const _dist &a) { ...copy... return this; };*/
    };
struct _cluster                     // star cluster structure
    {
    double x,y;                     // avg coordinate
    int iy;                         // ix of cluster match in the other set or -1
    double err;                     // error of cluster match
    List<int> ix;                   // star ix
    List<double> d;                 // distances of stars ix[] against each other
    _cluster(){}; _cluster(_cluster& a){ *this=a; }; ~_cluster(){}; _cluster* operator = (const _cluster *a) { *this=*a; return this; }; /*_cluster* operator = (const _cluster &a) { ...copy... return this; };*/
    };
const double max_r=5.0;             // find cluster max radius
const double max_err=0.2;           // match cluster max distance error treshold
const double max_rr=max_r*max_r;
const double max_errr=max_err*max_err;
int wi0,wi1;                        // result wandering star ix ...
int ix0[n0],ix1[n1];                // original star indexes
List<_cluster> cl0,cl1;             // found clusters

double txy1[n1][2];                 // transformed xy1[]
//---------------------------------------------------------------------------
double atanxy(double x,double y)
    {
    const double pi=M_PI;
    const double pi2=2.0*M_PI;
    int sx,sy;
    double a;
    const double _zero=1.0e-30;
    sx=0; if (x<-_zero) sx=-1; if (x>+_zero) sx=+1;
    sy=0; if (y<-_zero) sy=-1; if (y>+_zero) sy=+1;
    if ((sy==0)&&(sx==0)) return 0;
    if ((sx==0)&&(sy> 0)) return 0.5*pi;
    if ((sx==0)&&(sy< 0)) return 1.5*pi;
    if ((sy==0)&&(sx> 0)) return 0;
    if ((sy==0)&&(sx< 0)) return pi;
    a=y/x; if (a<0) a=-a;
    a=atan(a);
    if ((x>0)&&(y>0)) a=a;
    if ((x<0)&&(y>0)) a=pi-a;
    if ((x<0)&&(y<0)) a=pi+a;
    if ((x>0)&&(y<0)) a=pi2-a;
    return a;
    }
//---------------------------------------------------------------------------
void compute()
    {
    int i0,i1,e,f;
    double a,x,y;
    // original indexes (to keep track)
    for (e=0;e<n0;e++) ix0[e]=e;
    for (e=0;e<n1;e++) ix1[e]=e;
    // sort xy0[] by x asc
    for (e=1;e;) for (e=0,i0=0,i1=1;i1<n0;i0++,i1++)
     if (xy0[i0][0]>xy0[i1][0])
        {
        e=ix0[i0]   ; ix0[i0]   =ix0[i1]   ; ix0[i1]   =e; e=1;
        a=xy0[i0][0]; xy0[i0][0]=xy0[i1][0]; xy0[i1][0]=a;
        a=xy0[i0][1]; xy0[i0][1]=xy0[i1][1]; xy0[i1][1]=a;
        }
    // sort xy1[] by x asc
    for (e=1;e;) for (e=0,i0=0,i1=1;i1<n1;i0++,i1++)
     if (xy1[i0][0]>xy1[i1][0])
        {
        e=ix1[i0]   ; ix1[i0]   =ix1[i1]   ; ix1[i1]   =e; e=1;
        a=xy1[i0][0]; xy1[i0][0]=xy1[i1][0]; xy1[i1][0]=a;
        a=xy1[i0][1]; xy1[i0][1]=xy1[i1][1]; xy1[i1][1]=a;
        }
    _dist d;
    _cluster c,*pc,*pd;
    List<_dist> dist;
    // find star clusters in xy0[]
    for (cl0.num=0,i0=0;i0<n0;i0++)
        {
        for (dist.num=0,i1=i0+1;(i1<n0)&&(fabs(xy0[i0][0]-xy0[i1][0])<=max_r);i1++) // stars nearby
            {
            x=xy0[i0][0]-xy0[i1][0]; x*=x;
            y=xy0[i0][1]-xy0[i1][1]; y*=y; a=x+y;
            if (a<=max_rr) { d.ix=i1; d.d=a; dist.add(d); }
            }
        if (dist.num>=2)                                                            // add/compute cluster if found
            {
            c.ix.num=0; c.err=-1.0;
            c.ix.add(i0);   for (i1=0;i1<dist.num;i1++) c.ix.add(dist[i1].ix); c.iy=-1;
            c.x=xy0[i0][0]; for (i1=0;i1<dist.num;i1++) c.x+=xy0[dist[i1].ix][0]; c.x/=dist.num+1;
            c.y=xy0[i0][1]; for (i1=0;i1<dist.num;i1++) c.y+=xy0[dist[i1].ix][1]; c.y/=dist.num+1;
            for (e=1,i1=0;i1<cl0.num;i1++)
                {
                pc=&cl0[i1];
                x=c.x-pc->x; x*=x;
                y=c.y-pc->y; y*=y; a=x+y;
                if (a<max_rr)   // merge if too close to another cluster
                    {
                    pc->x=0.5*(pc->x+c.x);
                    pc->y=0.5*(pc->y+c.y);
                    for (e=0;e<c.ix.num;e++)
                        {
                        for (f=0;f<pc->ix.num;f++)
                         if (pc->ix[f]==c.ix[e]) { f=-1; break; }
                        if (f>=0) pc->ix.add(c.ix[e]);
                        }
                    e=0; break;
                    }
                }
            if (e) cl0.add(c);
            }
        }
    // full recompute clusters
    for (f=0,pc=&cl0[f];f<cl0.num;f++,pc++)
        {
        // avg coordinate
        pc->x=0.0;  for (i1=0;i1<pc->ix.num;i1++) pc->x+=xy0[pc->ix[i1]][0]; pc->x/=pc->ix.num;
        pc->y=0.0;  for (i1=0;i1<pc->ix.num;i1++) pc->y+=xy0[pc->ix[i1]][1]; pc->y/=pc->ix.num;
        // distances
        for (pc->d.num=0,i0=   0;i0<pc->ix.num;i0++)
        for (            i1=i0+1;i1<pc->ix.num;i1++)
            {
            x=xy0[pc->ix[i1]][0]-xy0[pc->ix[i0]][0]; x*=x;
            y=xy0[pc->ix[i1]][1]-xy0[pc->ix[i0]][1]; y*=y;
            pc->d.add(sqrt(x+y));
            }
        // sort by distance asc
        for (e=1;e;) for (e=0,i0=0,i1=1;i1<pc->d.num;i0++,i1++)
         if (pc->d[i0]>pc->d[i1])
            {
            a=pc->d[i0]; pc->d[i0]=pc->d[i1]; pc->d[i1]=a; e=1;
            }
        }

    // find star clusters in xy1[]
    for (cl1.num=0,i0=0;i0<n1;i0++)
        {
        for (dist.num=0,i1=i0+1;(i1<n1)&&(fabs(xy1[i0][0]-xy1[i1][0])<=max_r);i1++) // stars nearby
            {
            x=xy1[i0][0]-xy1[i1][0]; x*=x;
            y=xy1[i0][1]-xy1[i1][1]; y*=y; a=x+y;
            if (a<=max_rr) { d.ix=i1; d.d=a; dist.add(d); }
            }
        if (dist.num>=2)                                                            // add/compute cluster if found
            {
            c.ix.num=0; c.err=-1.0;
            c.ix.add(i0);   for (i1=0;i1<dist.num;i1++) c.ix.add(dist[i1].ix); c.iy=-1;
            c.x=xy1[i0][0]; for (i1=0;i1<dist.num;i1++) c.x+=xy1[dist[i1].ix][0]; c.x/=dist.num+1;
            c.y=xy1[i0][1]; for (i1=0;i1<dist.num;i1++) c.y+=xy1[dist[i1].ix][1]; c.y/=dist.num+1;
            for (e=1,i1=0;i1<cl1.num;i1++)
                {
                pc=&cl1[i1];
                x=c.x-pc->x; x*=x;
                y=c.y-pc->y; y*=y; a=x+y;
                if (a<max_rr)   // merge if too close to another cluster
                    {
                    pc->x=0.5*(pc->x+c.x);
                    pc->y=0.5*(pc->y+c.y);
                    for (e=0;e<c.ix.num;e++)
                        {
                        for (f=0;f<pc->ix.num;f++)
                         if (pc->ix[f]==c.ix[e]) { f=-1; break; }
                        if (f>=0) pc->ix.add(c.ix[e]);
                        }
                    e=0; break;
                    }
                }
            if (e) cl1.add(c);
            }
        }
    // full recompute clusters
    for (f=0,pc=&cl1[f];f<cl1.num;f++,pc++)
        {
        // avg coordinate
        pc->x=0.0;  for (i1=0;i1<pc->ix.num;i1++) pc->x+=xy1[pc->ix[i1]][0]; pc->x/=pc->ix.num;
        pc->y=0.0;  for (i1=0;i1<pc->ix.num;i1++) pc->y+=xy1[pc->ix[i1]][1]; pc->y/=pc->ix.num;
        // distances
        for (pc->d.num=0,i0=   0;i0<pc->ix.num;i0++)
        for (            i1=i0+1;i1<pc->ix.num;i1++)
            {
            x=xy1[pc->ix[i1]][0]-xy1[pc->ix[i0]][0]; x*=x;
            y=xy1[pc->ix[i1]][1]-xy1[pc->ix[i0]][1]; y*=y;
            pc->d.add(sqrt(x+y));
            }
        // sort by distance asc
        for (e=1;e;) for (e=0,i0=0,i1=1;i1<pc->d.num;i0++,i1++)
         if (pc->d[i0]>pc->d[i1])
            {
            a=pc->d[i0]; pc->d[i0]=pc->d[i1]; pc->d[i1]=a; e=1;
            }
        }
    // find matches
    for (i0=0,pc=&cl0[i0];i0<cl0.num;i0++,pc++) if  (pc->iy<0){ e=-1; x=0.0;
    for (i1=0,pd=&cl1[i1];i1<cl1.num;i1++,pd++) if (pc->d.num==pd->d.num)
            {
            for (y=0.0,f=0;f<pc->d.num;f++) y+=fabs(pc->d[f]-pd->d[f]);
            if ((e<0)||(x>y)) { e=i1; x=y; }
            }
        x/=pc->d.num;
        if ((e>=0)&&(x<max_err))
            {
            if (cl1[e].iy>=0) cl0[cl1[e].iy].iy=-1;
            pc->iy =e; cl1[e].iy=i0;
            pc->err=x; cl1[e].err=x;
            }
        }
    // compute transform
    double tx0,tx1,ty0,ty1,tc,ts;
    tx0=0.0; tx1=0.0; ty0=0.0; ty1=0.0; tc=1.0; ts=0.0; i0=-1; i1=-1;
    for (e=0,f=0,pc=&cl0[e];e<cl0.num;e++,pc++) if (pc->iy>=0)  // find 2 clusters with match
        {
        if (f==0)   i0=e;
        if (f==1) { i1=e; break; }
        f++;
        }
    if (i1>=0)
        {
        pc=&cl0[i0];        // translation and offset from xy0 set
        pd=&cl0[i1];
        tx1=pc->x;
        ty1=pc->y;
        a =atanxy(pd->x-pc->x,pd->y-pc->y);
        pc=&cl1[pc->iy];    // translation and offset from xy1 set
        pd=&cl1[pd->iy];
        tx0=pc->x;
        ty0=pc->y;
        a-=atanxy(pd->x-pc->x,pd->y-pc->y);
        tc=cos(a);
        ts=sin(a);
        }
    // transform xy1 -> txy1 (in xy0 coordinate system)
    for (i1=0;i1<n1;i1++)
        {
        x=xy1[i1][0]-tx0;
        y=xy1[i1][1]-ty0;
        txy1[i1][0]=x*tc-y*ts+tx1;
        txy1[i1][1]=x*ts+y*tc+ty1;
        }
    // sort txy1[] by x asc (after transfrm)
    for (e=1;e;) for (e=0,i0=0,i1=1;i1<n1;i0++,i1++)
     if (txy1[i0][0]>txy1[i1][0])
        {
        e= ix1[i0]   ;  ix1[i0]   = ix1[i1]   ;  ix1[i1]   =e; e=1;
        a=txy1[i0][0]; txy1[i0][0]=txy1[i1][0]; txy1[i1][0]=a;
        a=txy1[i0][1]; txy1[i0][1]=txy1[i1][1]; txy1[i1][1]=a;
        }
    // find match between xy0,txy1 (this can be speeded up by exploiting sorted order)
    int ix01[n0],ix10[n1];
    for (i0=0;i0<n0;i0++) ix01[i0]=-1;
    for (i1=0;i1<n1;i1++) ix10[i1]=-1;
    for (i0=0;i0<n0;i0++){ a=-1.0;
    for (i1=0;i1<n1;i1++)
        {
        x=xy0[i0][0]-txy1[i1][0]; x*=x;
        y=xy0[i0][1]-txy1[i1][1]; y*=y; x+=y;
        if (x<max_errr)
         if ((a<0.0)||(a>x)) { a=x; ix01[i0]=i1; ix10[i1]=i0; }
        }}
    // find the closest stars from unmatched stars
    a=-1.0; wi0=-1; wi1=-1;
    for (i0=0;i0<n0;i0++) if (ix01[i0]<0)
    for (i1=0;i1<n1;i1++) if (ix10[i1]<0)
        {
        x=xy0[i0][0]-txy1[i1][0]; x*=x;
        y=xy0[i0][1]-txy1[i1][1]; y*=y; x+=y;
        if ((wi0<0)||(a>x)) { a=x; wi0=i0; wi1=i1; }
        }
    }
//---------------------------------------------------------------------------
void draw()
    {
    bmp->Canvas->Font->Charset=OEM_CHARSET;
    bmp->Canvas->Font->Name="System";
    bmp->Canvas->Font->Pitch=fpFixed;
    bmp->Canvas->Font->Color=0x00FFFF00;
    bmp->Canvas->Brush->Color=0x00000000;
    bmp->Canvas->FillRect(TRect(0,0,xs,ys));
    _cluster *pc;
    int i,x0,y0,x1,y1,x2,y2,xx,yy,r,_r=4;
    double x,y,m;
    x0=xs/6; x1=3*x0; x2=5*x0;
    y0=ys/2; y1=  y0; y2=  y0;
    x=x0/60.0; y=y0/60.0; if (x<y) m=x; else m=y;
    // clusters match
    bmp->Canvas->Pen  ->Color=clAqua;
    bmp->Canvas->Brush->Color=0x00303030;
    for (i=0,pc=&cl0[i];i<cl0.num;i++,pc++)
     if (pc->iy>=0)
        {
        x=pc->x*m; xx=x0+x;
        y=pc->y*m; yy=y0-y;
        bmp->Canvas->MoveTo(xx,yy);
        x=cl1[pc->iy].x*m; xx=x1+x;
        y=cl1[pc->iy].y*m; yy=y1-y;
        bmp->Canvas->LineTo(xx,yy);
        }
    // clusters area
    for (i=0,pc=&cl0[i];i<cl0.num;i++,pc++)
        {
        x=pc->x*m; xx=x0+x;
        y=pc->y*m; yy=y0-y;
        r=pc->d[pc->d.num-1]*m*0.5+_r;
        bmp->Canvas->Ellipse(xx-r,yy-r,xx+r,yy+r);
        bmp->Canvas->TextOutA(xx+r,yy+r,AnsiString().sprintf("%.3lf",pc->err));
        }
    for (i=0,pc=&cl1[i];i<cl1.num;i++,pc++)
        {
        x=pc->x*m; xx=x1+x;
        y=pc->y*m; yy=y1-y;
        r=pc->d[pc->d.num-1]*m*0.5+_r;
        bmp->Canvas->Ellipse(xx-r,yy-r,xx+r,yy+r);
        bmp->Canvas->TextOutA(xx+r,yy+r,AnsiString().sprintf("%.3lf",pc->err));
        }
    // stars
    r=_r;
    bmp->Canvas->Pen  ->Color=clAqua;
    bmp->Canvas->Brush->Color=clBlue;
    for (i=0;i<n0;i++)
        {
        x=xy0[i][0]*m; xx=x0+x;
        y=xy0[i][1]*m; yy=y0-y;
        bmp->Canvas->Ellipse(xx-r,yy-r,xx+r,yy+r);
        }
    for (i=0;i<n1;i++)
        {
        x=xy1[i][0]*m; xx=x1+x;
        y=xy1[i][1]*m; yy=y1-y;
        bmp->Canvas->Ellipse(xx-r,yy-r,xx+r,yy+r);
        }
    // merged sets
    r=_r;
    bmp->Canvas->Pen  ->Color=clBlue;
    bmp->Canvas->Brush->Color=clBlue;
    for (i=0;i<n0;i++)
        {
        x=xy0[i][0]*m; xx=x2+x;
        y=xy0[i][1]*m; yy=y2-y;
        bmp->Canvas->Ellipse(xx-r,yy-r,xx+r,yy+r);
        }
    r=_r-2;
    bmp->Canvas->Pen  ->Color=clGreen;
    bmp->Canvas->Brush->Color=clGreen;
    for (i=0;i<n1;i++)
        {
        x=txy1[i][0]*m; xx=x2+x;
        y=txy1[i][1]*m; yy=y2-y;
        bmp->Canvas->Ellipse(xx-r,yy-r,xx+r,yy+r);
        }
    // wandering star
    r=_r+5;
    bmp->Canvas->Pen  ->Color=0x00FF8000;
    bmp->Canvas->Font ->Color=0x00FF8000;
    bmp->Canvas->Brush->Style=bsClear;
    x=xy0[wi0][0]*m; xx=x2+x;
    y=xy0[wi0][1]*m; yy=y2-y;
    bmp->Canvas->Ellipse(xx-r,yy-r,xx+r,yy+r);
    bmp->Canvas->TextOutA(xx+r,yy+r,ix0[wi0]);

    bmp->Canvas->Pen  ->Color=0x0040FF40;
    bmp->Canvas->Font ->Color=0x0040FF40;
    x=txy1[wi1][0]*m; xx=x2+x;
    y=txy1[wi1][1]*m; yy=y2-y;
    bmp->Canvas->Ellipse(xx-r,yy-r,xx+r,yy+r);
    bmp->Canvas->TextOutA(xx+r,yy+r,ix1[wi1]);
    bmp->Canvas->Brush->Style=bsSolid;

    Form1->Canvas->Draw(0,0,bmp);
    }
//---------------------------------------------------------------------------

And here the final result:

• 
• as you can see the result matches the answer from source link