游星 - codeabbey任务任务、游星、codeabbey

2023-09-11 05:36:28 作者:烟花耀眼灼伤我的回忆丶

我试图解决这个问题,我不知道下一步该怎么做。

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

链接到这个问题 问题陈述: 假设一些preliminary图像preprocessing中已经完成了,你有数据分上的两个图像坐标的形式。这些照片是约100×100毫米,坐标相对于它们的中心毫米给出。看看下面的示意图reprsentation: 你看,在这两个图片明星都显示在大致呈圆形的区域(认为它作为apperture我们的望远镜的),你可以发现,他们重新present同一片天空 - thoulgh稍微旋转和轻微移动

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.

请注意,某些恒星这是靠近边缘可能会缺席照片(因转移)之一 - 但流浪的星星是不是远离中心,因此是psented对两幅图像的$ P $ 。

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.

输入数据包含两个部分,分别对应于两个图像。 每个序列开始于一个整数 - 量的上市恒星。 然后分的坐标(X和Y)跟随。

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.

答案应该给两个索引(从0开始)徘徊的明星分别是第一和第二部分。

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

的例子是一样的画面的上方。星形#29从坐标(-18.2,11.1)的第一个部分是一样的星#3从第二部分与坐标(-19.7,6.9)。 示例输入数据: 94#第1条包含94星 -47.5 -10.4 19.1 25.9 18.9 -10.4 -2.1 -47.6 ... ... 92#第2节包含92星 -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 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 ... ...

我现在面临的问题 问题是,向量做的不匹配,甚至不具有相同的尺寸。因此,例如,在首节的第一个向量不匹配,首先在第二部分,所以我不能计算了基于该旋转矩阵。我也试过根据各部分的重心计算,但在边缘的一些点可能不存在,所以他们将有不同的质心(我试过,包括只有谁的长度向量为< 40,尺寸还是不匹配)。

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.

我所做的就是实现功能找到两个给定的向量之间的旋转矩阵。我使用的公式: transformed_vector = R * original_vector + T 其中R是旋转矩阵,并且由于载体还沿轴移动一点,我也加入吨 现在,我需要的是两个向量立足我的计算上。

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.

编辑:我也许应该提及的是,我居然给两个数组载体,每一个形象,我没有真正考虑到的图像。我需要找到感动基于这些载体的明星。

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.

谢谢!

推荐答案

完整的重新编辑

已经发现了一些时间/心情为这使其更加坚固

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

xy0 [],XY1 [] 是输入明星名单 让 max_r 是附近搜索区域treshld 让 max_err 是最大可接受的集群匹配误差 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:

排序xy0 []用x递增 这使搜索更快,更容易 sort xy0[] by x asc this makes the search faster and easier 在通过所有的恒星环 和交叉引用它们与恒星附近 ,因为排序附近恒星也将是接近当前星级指数 所以只搜索紧密区域之前和之后这个星阵 直到的x距离越过 max_r 如果发现添加集群 CL0 [] 集群列表 (集群2多收星) 添加新的群集之前 检查,如果没有集群是接近它 如果太靠近另一个集群,而不是合并 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 平均坐标 在里面的所有恒星之间的距离 将距离ASC进行排序 在如此检查,如果距离列表里面都是一样的 (记住绝​​对误差最低的总和) 如果误差更大然后max_err拒绝这一组的比赛 这是强大的配套都对这个数据集测试的许多簇(大max_r)没有错过比赛 在我使用的AVG集群的协调和它完美匹配

这是视觉效果:

在左侧 xy0 [] 设置 在中间的是 XY1 [] 设置 和右边两个,其中蓝色的大点是 xy0 [] 和绿色的小点被转化 XY1 [] 的数字集群匹配误差(-1表示没有找到匹配) 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)

[注意事项]

我用我自己的名单,其中,T&GT; 模板...

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

它只是重新分配动态线性阵列 名单,其中,INT&GT; X; 是一样的 INT X []; 其中, X [I] 是项目获得 x.num 是项目数数组中 x.add(5); 是一样的 X [x.num] = 5; x.num ++; 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++;

从这个角度,你可以检查的匹配 xy0并转化XY1

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

在这样的标志匹配的明星,以避免重复使用它们 使用一些treshold这就像 max_err 从什么星星都留下 找到两个是最接近对方 这应该是流浪的星星......玩得开心 (您可以重新排序的转化XY1) 请不要忘记使用原来的星级指数 IX0 [],IX1 [] 的结果输出 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

还剩下的工作原理

//---------------------------------------------------------------------------
// 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:

正如你所看到的结果,答案从源链接匹配