Wandering Star - Challenge

Question

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

Link to the problem
Problem:
Suppose that the preliminary processing of the preliminary image is already completed, and you have data in the form of the coordinates of the stars in two pictures. These images are about 100x100 millimeters, and the coordinates are also indicated in millimeters relative to their center. Take a look at the schematic diagram below: You see that in both pictures the stars are shown in approximately a circular area (think of it as the aperture of our telescope), and you can find out that they represent the same part of the sky - thoulgh are slightly rotated and slightly shifted.

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

Your task is to find such a "stray star", since it can be very likely cometary or asteroidal.

Please note that some stars that are close to the edge may not appear in one of the pictures (due to a shift), but the “wandering star” is not so far from the center and therefore is presented in both images.

The input contains two sections corresponding to two images. Each sequence begins with one integer - the number of stars listed. Then follow the coordinates (X and Y) of the stars.

The answer should give two indexes (based on 0) of the wandering star in the first and second sections, respectively.

An example is the same as above. Star number 29 from the first section with coordinates (-18.2, 11.1) coincides with star No. 3 from the second section with coordinates (-19.7, 6.9). Sample input:
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
...
...

Problem i am facing
The problem is that the vectors do not correspond to each other and do not even have the same size. So, for example, the first vector in the first section does not correspond to the first in the second section, so I can not calculate the rotation matrix based on this. I also tried to calculate it based on the centroids of each section, but some points on the edge may be missing, so they will have different centroids (I tried to include only those vectors whose length is <40, the size still does not match).

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

What I did was implement functions to find the rotation matrix between two given vectors. The formula I use:
transform_vector = R * original_vector + t
where R is the rotation matrix, and since the vector also moves along the axis a little, I also add t. Now all I need is the two vectors on which my calculations are based.

Edit: I should probably mention that in fact they give me two arrays of vectors, one for each image, in fact they do not give me images. I need to find a star that moved based on these vectors.

Thanks!

Answer 1

[edit2] complete reedit

Find some time / mood to make it more reliable.

• let xy0[],xy1[] - list of input stars
• let max_r - closest search area treshld
• let max_err - maximum allowable cluster matching error

so here is the Algorithm:

• sort xy0 [] by x asc
  • it speeds up and simplifies the search
• Find star clusters in xy0[]
  • through all the stars
  • and cross-referencing them with stars nearby
  • due to the sorting of nearby stars, it will also be next to the current stellar index
  • so just find the closing area before and after this star in the array
  • until x-distance crosses max_r
  • add a cluster to cl0[] list of the cluster, if found
  • (cluster equal to 2 and closer stars)
  • before adding a new cluster
  • check if there is a cluster next to it
  • and merge if it is too close to another cluster
• full list of found clusters
  • average speed
  • the distances between all the stars inside
  • sort them by distance asc
• do 1., 2., 3. also for xy1[],cl1[]
• find matches between clusters
  • check that the list of distances inside is the same
  • (remember the smallest amount of abs error)
  • if the error is greater than max_err rejects this cluster match
  • this is a strong match that has been tested on many clusters (large max_r) without missing matches in this dataset.
• take 2 clusters from found clusters in cl0[] , which found a match
• take also mapped clusters
• calculate the conversion between xy0[],xy1[] from these 4 points
  • I use the middle coordinate of the cluster and fits perfectly

This is the visual result:

• 
• left side xy0[] set
• middle xy1[] set
• and to the right, where the blue big points are xy0[]
• and green smaller dots are converted xy1[]
• digits are a cluster matching error (-1 means no match found)

[notes]

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

• it just dynamically redistributes the linear array
• List<int> x; matches int x[];
• where x[i] access to the item is available
• x.num - the number of elements in the array
• x.add(5); matches x[x.num]=5; x.num++; x[x.num]=5; x.num++;

From now on, you can check the correspondence between xy0 and transformed xy1

• so that the flags match the stars to avoid duplication of their use.
• use something like max_err for this
• what stars are left
• find the two closest to each other
• It must be a wandering star ... have fun
• (you can sort the converted xy1 again)
• do not forget to use the original stellar indices ix0[],ix1[] to display the result

[edit3] also other work

 //--------------------------------------------------------------------------- // 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; dd=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; cx=xy0[i0][0]; for (i1=0;i1<dist.num;i1++) c.x+=xy0[dist[i1].ix][0]; cx/=dist.num+1; cy=xy0[i0][1]; for (i1=0;i1<dist.num;i1++) c.y+=xy0[dist[i1].ix][1]; cy/=dist.num+1; for (e=1,i1=0;i1<cl0.num;i1++) { pc=&cl0[i1]; x=cx-pc->x; x*=x; y=cy-pc->y; y*=y; a=x+y; if (a<max_rr) // merge if too close to another cluster { pc->x=0.5*(pc->x+cx); pc->y=0.5*(pc->y+cy); 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; dd=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; cx=xy1[i0][0]; for (i1=0;i1<dist.num;i1++) c.x+=xy1[dist[i1].ix][0]; cx/=dist.num+1; cy=xy1[i0][1]; for (i1=0;i1<dist.num;i1++) c.y+=xy1[dist[i1].ix][1]; cy/=dist.num+1; for (e=1,i1=0;i1<cl1.num;i1++) { pc=&cl1[i1]; x=cx-pc->x; x*=x; y=cy-pc->y; y*=y; a=x+y; if (a<max_rr) // merge if too close to another cluster { pc->x=0.5*(pc->x+cx); pc->y=0.5*(pc->y+cy); 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 is the final result:

• 
• as you can see that the result matches the answer from the original link

Answer 2

First, for simplicity, let it first pretend that there is no wandering star, and that each star located in one image is in another - that is, that image B was created from image A using only (possibly zero) shift and ( possibly zero) rotation.

Captions

I think that the best way to approach this is to “forget about” the position of each star in each image and instead write for each star only a signature: some information that will not be changed by a change (translation) or rotation. Then you can look at each signature in image A and try to match it with the signatures in image B. If we make the signatures in a reasonable way, then they should distinguish the stars well, so that no two stars in the same image get the same or very similar signatures .

For a given star, what does not change when shifting or rotating is the distance to any other star in the same image. Thus, you can use the full range of distances for all other n-1 stars as a star signature (stored as an unordered set, since we don’t know what stars there are). But this is probably enough (and faster, and ultimately more reliable) just to use a subset of them, in which case the distances to the nearest neighbors are probably the most informative (unless you have repeating patterns of stars appearing). How much to use? I will start by calculating the distance of each star to its nearest neighbor in the same image; if on image A all of these distances are “different” (that is, they have a difference below some arbitrary threshold that you choose), as well as on image B, then you’ll finish - that is, the signatures that you have calculated already distinguish quite well stars - otherwise go back and for each star, paste the distance to your second nearest neighbor in your signature, repeating this until the “uniqueness” of the signatures is achieved.

Signature Comparison

This raises the question of how to decide whether two signatures are “different” when the signature contains more than one distance. Now you are comparing two sets of numbers instead of two numbers. One way is to sort the distances in each signature, and then compare them, calculating some measure of the distance between points, such as the sum of the squared differences (for example, if each signature contains 3 distances, this will be (u [1] -v [1]) ^ 2 + (u [2] -v [2]) ^ 2 + (u [3] -v [3]) ^ 2 for two stars u and v, and u [i] is the distance from the star u to its i to the nearest neighbor in the same image, as well as for v.) In your example dataset, I would suggest that probably 3-4 distances per star would be enough.

Reliable and efficient signature comparison

However, this will not be particularly stable when wandering and missed stars are considered: suppose that some star u in image A has a wandering star v as its closest neighbor, and in image B, v moves away so that now it’s the 5th nearest neighbor. Then, if we compare the u-signature on image A with its signature on image B, we get an incorrectly large distance, because we will “blindly” compare uA [1] with uB [1], etc. Instead of comparing uA [2] with uB [1], uA [3] with uB [2], etc., I would like (since these distances will be equal). Thus, the best way to compare the two signatures is to take an ordered list of distances for each of them and let the numbers “slide until they match” the numbers in the other.

Fortunately, this can be done in linear time by merging lists . Given two sorted lists of numbers X and Y, we carry out the usual strategy of merging the pick list of the smallest element from any list, deleting it from there and writing it down, but we also track two possible partial solutions: Bm, the estimate of the best partial solution, in which the last exit number is matched with an earlier number, and Bf is the estimate of the best partial solution in which it is free (so far it is not matched with anything). Choose the cost of the fine p to assign to numbers that do not match any number in another list, and the counting function of the estimate (a, b), which assigns the value 0 when the values ​​a = b and higher have more different values ​​of a and b ( for example, the estimate (a, b) = (a - b) ^ 2). Then the optimal counter sigDiff (X, Y) can be calculated using:

• Set Bm = Bf = 0.
• If X and Y are both empty, then stop. The end result: min (Bf + p, Bm).
• Call the element at the beginning of X x and the element in front of Y y Let z be smaller than x and y, and set the list identifier for this list (X or Y) from which z came.
• Remove the first element (z) from the list named whichList.
• If ifList == prevWhichList (i.e. if the previous smallest number, prev, was also from the same list as z), or if this is the first iteration, set newBm = min (Bf + p, Bm) + p.
• Otherwise, it is possible to match z with the previous number, so set newBm = Bf + score (z, prev).
• Set Bf = min (Bf + p, Bm).
• Set Bm = newBm.
• Set prev = z and prevWhichList = whichList.
• Go to 2.

For example, suppose we have two distance lists (that is, two signatures) X = [10, 20, 30, 40] and Y = [11, 18, 41, 50], the penalty is 20, and the estimate ( ) is the sum of the squared differences, as indicated 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 

Unlike the naive square sum method:

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

Matching captions between two images

We built all this equipment so that we can try to match the stars between the two images by matching their signatures. Probably the “right” way to do this would be to solve the two-sided problem with maximum weight correspondence on the graph in which the vertices on one side are the signatures of stars in image A, the vertices on the other side are the signatures of stars in image B, and on each vertex X in each side has an edge of each vertex Y on the other side having a weight of -sigDiff (X, Y) (negated because we want to minimize the total difference).

However, solving this may take some time, and it is also possible that this algorithm will find some incorrect “approximate” matches in an attempt to get the total minimum cost. We would prefer to focus only on those pairs of stars that we definitely respond to each other. This can be done much faster using a heuristic based on the idea that if for star X in image A, the closest star in image B (based on sigDiff ()) is Y, it turns out that the nearest star Y in image A is also X ( i.e., if the best match of the best match of X is X), then we can be sure that X and Y really match. These confidently matching pairs can be found quickly and used to evaluate the affine transformation from image A to image B using least squares .

Ignoring Missing Stars

First, we calculate the convex hull of the images of B stars. Then, using a transformation determined from confident star pairs, we will transform each star in image A to its supposed corresponding location on image B and take the convex hull of this converted image. Then we intersect these two convex hulls: each star in any image falling inside this intersection is a star that we expect to find in both images. Finally, we throw away all the stars from any image located outside this intersecting convex hull, and start over!

After starting everything all over again (maybe it’s not necessary to repeat everything really), we should find that on image A there are exactly 1 star in image A and 1 star in image B, which have a bad resemblance to all other stars on another image (as measured by sigDiff (), as usual). These "two" stars, of course, are the only "wandering" stars :)