Search code examples
pythonalgorithmvectorpuzzle

Wandering star - codeabbey task


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: enter image description here 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!


Solution

  • [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:

    • output example
    • 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:

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