流星-codeabbey任务
Wandering star - codeabbey task
我正在尝试解决这个问题，但不确定下一步该怎么做。
Link to the problem
问题陈述：
假设一些初步的图像预处理已经完成，并且您有两张图片上星星坐标形式的数据。这些图片大约为 100x100 毫米，坐标也以相对于它们中心的毫米为单位给出。请看下面的示意图：
你会看到在两张图片中星星都显示在大致圆形的区域中（可以把它想象成我们 telescope 的光圈），你会发现它们代表同一片天空 - 尽管略微旋转和略微移动.
您还可以看到其中一颗星星（标有红色箭头）相对于其他星星改变了位置。
你的任务是找出这样一个 "wandering star" 因为它很有可能是彗星或小行星。
请注意，其中一张图片中可能没有一些靠近边缘的星星（由于偏移）- 但 "wandering star" 离中心不远，因此在两张图片上都出现了.
输入数据 包含对应于两个图像的两个部分。
每个序列都以一个整数开始——列出的星星数量。
然后是星星的坐标（X 和 Y）。
回答应该分别给出第一段和第二段流星的两个索引（从0开始）
示例同上图。第一部分中坐标为 (-18.2, 11.1) 的星号 #29 与第二部分中坐标为 (-19.7, 6.9) 的星号 #3 相同。
示例输入数据：
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
...
...
我面临的问题
问题是向量不匹配，甚至大小都不相同。因此，例如，第一部分中的第一个向量与第二部分中的第一个向量不匹配，因此我无法基于此计算旋转矩阵。我也尝试根据每个部分的质心计算它，但是边缘上的一些点可能不存在，所以它们会有不同的质心（我尝试只包括长度 < 40 的向量，大小仍然不匹配） .
所以我的问题是我应该根据什么进行计算？如何找到一些匹配向量，以便计算它们的旋转矩阵？我需要朝着正确的方向前进。
我所做的是实现找到两个给定向量之间的旋转矩阵的函数。我使用的公式：
transformed_vector = R * original_vector + t
其中 R 是旋转矩阵，并且由于矢量也沿轴移动了一点，所以我还添加了 t
现在我只需要两个向量作为我计算的基础。
编辑：我可能应该提一下，实际上我得到了两个矢量数组，每个图像一个，我实际上并没有得到图像。我需要根据这些矢量找到移动的星星。
谢谢！
[edit2] 完成重新编辑
已经找到一些 time/mood 来让它更健壮
- 设
xy0[],xy1[]为输入星表
- 设
max_r为附近搜索区域treshld
- 让
max_err 成为可接受的最大簇匹配误差
所以这里是算法：
- 按 x 升序对 xy0[] 排序
- 这使得搜索更快更容易
- 在
xy0[]中找到星团
- 遍历所有星星
- 并将它们与附近的星星交叉引用
- 因为排序附近的星星也会靠近当前的星星索引
- 所以只搜索阵列中这颗星前后的附近区域
- 直到 x 距离越过
max_r
- 如果找到，将集群添加到
cl0[] 集群列表
- （星团是 2 颗或更近的恒星）
- 添加新集群之前
- 检查附近是否没有集群
- 如果离另一个集群太近则合并
- 完全重新计算发现集群
- 平均坐标
- 里面所有星星之间的距离
- 按距离升序排序
- 执行 1.,2.,3。也为
xy1[],cl1[]
- 查找集群之间的匹配项
- 所以检查里面的距离列表是否相同
- (记住abs误差的最小总和)
- 如果误差较大，则 max_err 拒绝该簇匹配
- 这是强匹配，已经在许多集群（大 max_r）上测试过，在此数据集上没有遗漏匹配
- 从
cl0[] 中找到匹配项的簇中取出 2 个簇
- 也取匹配的簇
- 从这 4 个点计算
xy0[],xy1[] 之间的转换
- 我使用集群的平均坐标，它完美匹配
这是视觉结果：
- 左边是
xy0[]集
- 中间是
xy1[]集
- 右边两个蓝色大点是
xy0[]
- 并且绿色小点被变换
xy1[]
- 数字是簇匹配的错误（-1表示没有找到匹配）
[注释]
我使用自己的List<T>模板...
- 只是动态重新分配线性数组
List<int> x; 等同于 int x[];- 其中
x[i] 是项目访问权限
x.num 是数组中的项目数x.add(5); 等同于 x[x.num]=5; x.num++;
从这一点开始，您可以检查 xy0 and transformed xy1
之间的匹配项
- 所以标记匹配的星星以避免重复使用它们
- 为此使用一些阈值，例如
max_err
- 剩下的星星
- 找到彼此最接近的两个
- 这应该是流星吧……玩得开心
- (可以对变换后的xy1重新排序)
- 不要忘记使用原始星号索引
ix0[],ix1[]作为结果输出
[edit3] 其余的也有效
//---------------------------------------------------------------------------
// 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);
}
//---------------------------------------------------------------------------
这里是最终结果：
- 如您所见，结果与来源 link
的答案相匹配
首先，为简单起见，我们先假设没有流星，一张图像中的每颗星星都在另一张图像中——也就是说，图像 B 是从图像 A 中仅应用一个（可能为零）移位和（可能为零）旋转。
签名
我认为解决这个问题的最好方法是 "forget about" 每张图片中每颗星星的位置，而不是只为每颗星星记录一个 签名 ：一些不会因移位（平移）或旋转而改变的信息。然后你可以看看图A中的每一个签名，并尝试将其与图B中的签名进行匹配。如果我们以合理的方式制作签名，那么它们应该做好区分星星的工作，这样就不会出现两颗星在同一张图片中得到相同或高度相似的签名。
对于给定的恒星，不会因移动或旋转而改变的是它与同一图像中任何其他恒星的 距离。因此，您可以使用到所有其他 n-1 颗恒星的完整距离集作为恒星的签名（存储为无序集，因为我们不知道 "which stars are which"）。但仅使用其中的一个子集可能就足够了（而且速度更快，最终更稳健），在这种情况下，与其最近的邻居的距离可能是最有信息的（除非你有重复出现的星星模式）。使用多少？我将从计算同一图像中每颗恒星与其最近邻居的距离开始；如果在图像 A 中，所有这些距离都是 "different enough"（即，差异低于您选择的某个任意阈值），同样在图像 B 中，那么您就完成了——也就是说，您的签名已经计算出足够好地区分星星 - 否则返回并为每颗星星将到其第二近邻的距离插入其签名，重复此操作直到获得 "uniqueness" 个签名。
比较签名
这就提出了一个问题，当一个签名包含一个以上的距离时，如何判断两个签名是否"different enough"。现在您正在比较两组数字，而不仅仅是两个数字。一种方法是对每个签名中的距离进行排序，然后通过计算点之间的一些距离度量（例如平方差之和）来比较它们（例如，如果每个签名包含 3 个距离，这将是 (u[1]-v [1])^2 + (u[2]-v[2])^2 + (u[3]-v[3])^2 对于两颗星 u 和 v，其中 u[i] 是到的距离在同一图像中将 u 标记为它的第 i 个最近的邻居，对于 v 也是如此。）在您的示例数据集上，我猜想每颗星可能有 3-4 个距离就足够了。
稳健高效的签名比较
然而，一旦考虑到流星和缺失星，这将变得不是特别稳健：假设图像 A 中的某颗星 u 将流星 v 作为其最近的邻居，而在图像 B 中，v 已经移开所以它现在是你的第 5 个最近的邻居。那么如果我们比较你在图像 A 中的签名和它在图像 B 中的签名，我们会得到一个不正确的大距离，因为我们将 "blindly" 比较 uA[1] 和 uB[1] 等等，而不是比较uA[2] 与 uB[1]、uA[3] 与 uB[2] 等等（因为这些距离相等）。因此，比较两个签名的更好方法是获取每个签名的有序距离列表，并允许数字 "slide along until they fit" 另一个签名中的数字。
令人高兴的是，这可以通过 列表合并 在线性时间内完成。给定两个已排序的数字 X 和 Y 列表，我们执行通常的列表合并策略，从任一列表中选择最小的元素，将其从那里移除并写出，但我们还跟踪两个候选部分解决方案：Bm，分数最近的输出数字与较早的数字匹配的最佳部分解决方案的分数，以及 Bf，最佳部分解决方案的分数，其中它是免费的（尚未与任何匹配）。选择一个惩罚成本 p 分配给与另一个列表中的任何数字都不匹配的数字，以及一个评分函数 score(a, b) 当 a = b 时分配一个值 0 并且值越高，a 和 b 越不同是（例如，score(a, b) = (a - b)^2）。然后可以使用以下公式计算最佳分数 sigDiff(X, Y)：
- 设置 Bm = Bf = 0。
- 如果X和Y都为空，则停止。最终得分为min(Bf + p, Bm).
- 把X前面的元素叫做x，把Y前面的元素叫做y。令 z 为 x 和 y 中的较小者，并将 whichList 设置为 z 来自的列表（X 或 Y）的 ID。
- 从 whichList 命名的列表中删除第一个元素 (z)。
- 如果 whichList == prevWhichList（即如果前一个最小数字 prev 也来自与 z 相同的列表），或者如果这是第一次迭代，则设置 newBm = min(Bf + p, Bm) + p.
- 否则可能会把z和之前写出来的数字匹配起来，所以设newBm = Bf + score(z, prev)。
- 设置 Bf = min(Bf + p, Bm)。
- 设置 Bm = newBm。
- 设置 prev = z 和 prevWhichList = whichList。
- 转到 2。
例如，假设我们有两个距离列表（即两个签名）X = [10, 20, 30, 40] 和 Y = [11, 18, 41, 50]，惩罚是 20，和 score() 是上面的平方差之和。上述算法将"align"两个序列如下：
X 10 20 30 40
matches \ / \
Y 11 18 41 50
cost 1 4 20 1 20 Total cost: 46
相比之下，"naive" 平方和方法将给出：
X 10 20 30 40
matches | | | |
Y 11 18 41 50
cost 1 4 121 100 Total cost: 226
两个图像之间的匹配签名
我们已经建立了所有这些机器，以便我们可以通过匹配它们的签名来尝试匹配两张图像之间的星星。可以说 "right" 方法可以解决图上的 bipartite maximum weighted matching problem 问题，其中一侧的顶点是图像 A 中星星的签名，另一侧的顶点是签名图像 B 中的星星，并且从一侧的每个顶点 X 到另一侧的每个顶点 Y 都有一条边，其权重为 -sigDiff(X, Y)（取反，因为我们想最小化总差）。
但是，解决这个问题可能需要一段时间，而且此算法也有可能会找到一些不正确的 "approximate" 匹配项，以尝试获得最低的总体成本。我们更愿意只关注那些我们确定彼此对应的恒星对。这可以通过基于以下想法的启发式方法更快地完成，即如果对于图像 A 中的恒星 X，其图像 B 中最近的恒星（基于 sigDiff()）是 Y，则结果是图像 A 中 Y 最接近的恒星是也 X（即，如果 X 的最佳匹配中的最佳匹配是 X），那么我们可以确信 X 和 Y 确实相互匹配。可以快速找到这些可靠的匹配对，并使用 least squares.
来估计从图像 A 到图像 B 的仿射变换
忽略缺失的星星
我们首先计算图像 B 的星星的凸包。然后，使用从置信星对确定的变换，我们将图像 A 中的每颗星变换到其在图像 B 中的估计对应位置，并采用该变换图像的凸包。然后我们将这两个凸包相交：落在该交点内的任一图像中的每颗星都是我们期望在两个图像中找到的星。最后，我们将两幅图像中位于相交凸包外的所有星星都扔掉，然后重新开始！
再次运行 everything之后（实际上可能不需要重新运行everything），我们应该会发现图像A中恰好有1颗星和1颗星在图像 B 中，与另一幅图像中的所有其他恒星的相似性很差（像往常一样由 sigDiff() 测量）。这些 "two" 颗星当然是单颗 "wandering" 颗星 :)
我正在尝试解决这个问题，但不确定下一步该怎么做。
Link to the problem
问题陈述：
假设一些初步的图像预处理已经完成，并且您有两张图片上星星坐标形式的数据。这些图片大约为 100x100 毫米，坐标也以相对于它们中心的毫米为单位给出。请看下面的示意图：
您还可以看到其中一颗星星（标有红色箭头）相对于其他星星改变了位置。
你的任务是找出这样一个 "wandering star" 因为它很有可能是彗星或小行星。
请注意，其中一张图片中可能没有一些靠近边缘的星星（由于偏移）- 但 "wandering star" 离中心不远，因此在两张图片上都出现了.
输入数据 包含对应于两个图像的两个部分。 每个序列都以一个整数开始——列出的星星数量。 然后是星星的坐标（X 和 Y）。
回答应该分别给出第一段和第二段流星的两个索引（从0开始）
示例同上图。第一部分中坐标为 (-18.2, 11.1) 的星号 #29 与第二部分中坐标为 (-19.7, 6.9) 的星号 #3 相同。
示例输入数据：
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
...
...
我面临的问题
问题是向量不匹配，甚至大小都不相同。因此，例如，第一部分中的第一个向量与第二部分中的第一个向量不匹配，因此我无法基于此计算旋转矩阵。我也尝试根据每个部分的质心计算它，但是边缘上的一些点可能不存在，所以它们会有不同的质心（我尝试只包括长度 < 40 的向量，大小仍然不匹配） .
所以我的问题是我应该根据什么进行计算？如何找到一些匹配向量，以便计算它们的旋转矩阵？我需要朝着正确的方向前进。
我所做的是实现找到两个给定向量之间的旋转矩阵的函数。我使用的公式：
transformed_vector = R * original_vector + t
其中 R 是旋转矩阵，并且由于矢量也沿轴移动了一点，所以我还添加了 t
现在我只需要两个向量作为我计算的基础。
编辑：我可能应该提一下，实际上我得到了两个矢量数组，每个图像一个，我实际上并没有得到图像。我需要根据这些矢量找到移动的星星。
谢谢！
[edit2] 完成重新编辑
已经找到一些 time/mood 来让它更健壮
- 设
xy0[],xy1[]为输入星表 - 设
max_r为附近搜索区域treshld - 让
max_err成为可接受的最大簇匹配误差
所以这里是算法：
- 按 x 升序对 xy0[] 排序
- 这使得搜索更快更容易
- 在
xy0[]中找到星团- 遍历所有星星
- 并将它们与附近的星星交叉引用
- 因为排序附近的星星也会靠近当前的星星索引
- 所以只搜索阵列中这颗星前后的附近区域
- 直到 x 距离越过
max_r - 如果找到，将集群添加到
cl0[]集群列表 - （星团是 2 颗或更近的恒星）
- 添加新集群之前
- 检查附近是否没有集群
- 如果离另一个集群太近则合并
- 完全重新计算发现集群
- 平均坐标
- 里面所有星星之间的距离
- 按距离升序排序
- 执行 1.,2.,3。也为
xy1[],cl1[] - 查找集群之间的匹配项
- 所以检查里面的距离列表是否相同
- (记住abs误差的最小总和)
- 如果误差较大，则 max_err 拒绝该簇匹配
- 这是强匹配，已经在许多集群（大 max_r）上测试过，在此数据集上没有遗漏匹配
- 从
cl0[]中找到匹配项的簇中取出 2 个簇 - 也取匹配的簇
- 从这 4 个点计算
xy0[],xy1[]之间的转换- 我使用集群的平均坐标，它完美匹配
这是视觉结果：
- 左边是
xy0[]集 - 中间是
xy1[]集 - 右边两个蓝色大点是
xy0[] - 并且绿色小点被变换
xy1[] - 数字是簇匹配的错误（-1表示没有找到匹配）
[注释]
我使用自己的List<T>模板...
- 只是动态重新分配线性数组
List<int> x;等同于int x[];- 其中
x[i]是项目访问权限 x.num是数组中的项目数x.add(5);等同于x[x.num]=5; x.num++;
从这一点开始，您可以检查 xy0 and transformed xy1
- 所以标记匹配的星星以避免重复使用它们
- 为此使用一些阈值，例如
max_err - 剩下的星星
- 找到彼此最接近的两个
- 这应该是流星吧……玩得开心
- (可以对变换后的xy1重新排序)
- 不要忘记使用原始星号索引
ix0[],ix1[]作为结果输出
[edit3] 其余的也有效
//---------------------------------------------------------------------------
// 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);
}
//---------------------------------------------------------------------------
这里是最终结果：
- 如您所见，结果与来源 link 的答案相匹配
首先，为简单起见，我们先假设没有流星，一张图像中的每颗星星都在另一张图像中——也就是说，图像 B 是从图像 A 中仅应用一个（可能为零）移位和（可能为零）旋转。
签名
我认为解决这个问题的最好方法是 "forget about" 每张图片中每颗星星的位置，而不是只为每颗星星记录一个 签名 ：一些不会因移位（平移）或旋转而改变的信息。然后你可以看看图A中的每一个签名，并尝试将其与图B中的签名进行匹配。如果我们以合理的方式制作签名，那么它们应该做好区分星星的工作，这样就不会出现两颗星在同一张图片中得到相同或高度相似的签名。
对于给定的恒星，不会因移动或旋转而改变的是它与同一图像中任何其他恒星的 距离。因此，您可以使用到所有其他 n-1 颗恒星的完整距离集作为恒星的签名（存储为无序集，因为我们不知道 "which stars are which"）。但仅使用其中的一个子集可能就足够了（而且速度更快，最终更稳健），在这种情况下，与其最近的邻居的距离可能是最有信息的（除非你有重复出现的星星模式）。使用多少？我将从计算同一图像中每颗恒星与其最近邻居的距离开始；如果在图像 A 中，所有这些距离都是 "different enough"（即，差异低于您选择的某个任意阈值），同样在图像 B 中，那么您就完成了——也就是说，您的签名已经计算出足够好地区分星星 - 否则返回并为每颗星星将到其第二近邻的距离插入其签名，重复此操作直到获得 "uniqueness" 个签名。
比较签名
这就提出了一个问题，当一个签名包含一个以上的距离时，如何判断两个签名是否"different enough"。现在您正在比较两组数字，而不仅仅是两个数字。一种方法是对每个签名中的距离进行排序，然后通过计算点之间的一些距离度量（例如平方差之和）来比较它们（例如，如果每个签名包含 3 个距离，这将是 (u[1]-v [1])^2 + (u[2]-v[2])^2 + (u[3]-v[3])^2 对于两颗星 u 和 v，其中 u[i] 是到的距离在同一图像中将 u 标记为它的第 i 个最近的邻居，对于 v 也是如此。）在您的示例数据集上，我猜想每颗星可能有 3-4 个距离就足够了。
稳健高效的签名比较
然而，一旦考虑到流星和缺失星，这将变得不是特别稳健：假设图像 A 中的某颗星 u 将流星 v 作为其最近的邻居，而在图像 B 中，v 已经移开所以它现在是你的第 5 个最近的邻居。那么如果我们比较你在图像 A 中的签名和它在图像 B 中的签名，我们会得到一个不正确的大距离，因为我们将 "blindly" 比较 uA[1] 和 uB[1] 等等，而不是比较uA[2] 与 uB[1]、uA[3] 与 uB[2] 等等（因为这些距离相等）。因此，比较两个签名的更好方法是获取每个签名的有序距离列表，并允许数字 "slide along until they fit" 另一个签名中的数字。
令人高兴的是，这可以通过 列表合并 在线性时间内完成。给定两个已排序的数字 X 和 Y 列表，我们执行通常的列表合并策略，从任一列表中选择最小的元素，将其从那里移除并写出，但我们还跟踪两个候选部分解决方案：Bm，分数最近的输出数字与较早的数字匹配的最佳部分解决方案的分数，以及 Bf，最佳部分解决方案的分数，其中它是免费的（尚未与任何匹配）。选择一个惩罚成本 p 分配给与另一个列表中的任何数字都不匹配的数字，以及一个评分函数 score(a, b) 当 a = b 时分配一个值 0 并且值越高，a 和 b 越不同是（例如，score(a, b) = (a - b)^2）。然后可以使用以下公式计算最佳分数 sigDiff(X, Y)：
- 设置 Bm = Bf = 0。
- 如果X和Y都为空，则停止。最终得分为min(Bf + p, Bm).
- 把X前面的元素叫做x，把Y前面的元素叫做y。令 z 为 x 和 y 中的较小者，并将 whichList 设置为 z 来自的列表（X 或 Y）的 ID。
- 从 whichList 命名的列表中删除第一个元素 (z)。
- 如果 whichList == prevWhichList（即如果前一个最小数字 prev 也来自与 z 相同的列表），或者如果这是第一次迭代，则设置 newBm = min(Bf + p, Bm) + p.
- 否则可能会把z和之前写出来的数字匹配起来，所以设newBm = Bf + score(z, prev)。
- 设置 Bf = min(Bf + p, Bm)。
- 设置 Bm = newBm。
- 设置 prev = z 和 prevWhichList = whichList。
- 转到 2。
例如，假设我们有两个距离列表（即两个签名）X = [10, 20, 30, 40] 和 Y = [11, 18, 41, 50]，惩罚是 20，和 score() 是上面的平方差之和。上述算法将"align"两个序列如下：
X 10 20 30 40
matches \ / \
Y 11 18 41 50
cost 1 4 20 1 20 Total cost: 46
相比之下，"naive" 平方和方法将给出：
X 10 20 30 40
matches | | | |
Y 11 18 41 50
cost 1 4 121 100 Total cost: 226
两个图像之间的匹配签名
我们已经建立了所有这些机器，以便我们可以通过匹配它们的签名来尝试匹配两张图像之间的星星。可以说 "right" 方法可以解决图上的 bipartite maximum weighted matching problem 问题，其中一侧的顶点是图像 A 中星星的签名，另一侧的顶点是签名图像 B 中的星星，并且从一侧的每个顶点 X 到另一侧的每个顶点 Y 都有一条边，其权重为 -sigDiff(X, Y)（取反，因为我们想最小化总差）。
但是，解决这个问题可能需要一段时间，而且此算法也有可能会找到一些不正确的 "approximate" 匹配项，以尝试获得最低的总体成本。我们更愿意只关注那些我们确定彼此对应的恒星对。这可以通过基于以下想法的启发式方法更快地完成，即如果对于图像 A 中的恒星 X，其图像 B 中最近的恒星（基于 sigDiff()）是 Y，则结果是图像 A 中 Y 最接近的恒星是也 X（即，如果 X 的最佳匹配中的最佳匹配是 X），那么我们可以确信 X 和 Y 确实相互匹配。可以快速找到这些可靠的匹配对，并使用 least squares.
来估计从图像 A 到图像 B 的仿射变换忽略缺失的星星
我们首先计算图像 B 的星星的凸包。然后，使用从置信星对确定的变换，我们将图像 A 中的每颗星变换到其在图像 B 中的估计对应位置，并采用该变换图像的凸包。然后我们将这两个凸包相交：落在该交点内的任一图像中的每颗星都是我们期望在两个图像中找到的星。最后，我们将两幅图像中位于相交凸包外的所有星星都扔掉，然后重新开始！
再次运行 everything之后（实际上可能不需要重新运行everything），我们应该会发现图像A中恰好有1颗星和1颗星在图像 B 中，与另一幅图像中的所有其他恒星的相似性很差（像往常一样由 sigDiff() 测量）。这些 "two" 颗星当然是单颗 "wandering" 颗星 :)