匹配x,y指向另一个缩放,旋转,平移和缺少元素的集合
Gab*_*iel 10 python performance geometry astronomy affinetransform
(为什么我这样做？见下面的解释)
考虑两个点集,A并B如下图所示
它可能看起来不像,但是set A在集合中是"隐藏的" B.它不容易被看到,因为相对于点B的缩放,旋转和平移.更糟糕的是,存在的一些点缺失,并且包含许多不存在的点.(x, y)AABBA
我需要找到必须应用于B集合的适当缩放,旋转和平移,以便将其与集合匹配A.在上面显示的情况中,正确的值是:
scale = 0.14, rot_angle = 0.0, x_transl = 35.0, y_transl = 2.0
产生(足够好)匹配
(红色,仅显示匹配的B点;这些点位于1000<x<2000, y~2000右侧第一个图中的扇区中).鉴于很多自由度(DoF:缩放+旋转+ 2D平移)我知道错过匹配的可能性,但是点的坐标不是随机的(尽管它们可能看起来像它)所以这个概率发生的事情很小.
我编写的代码(见下文)使用强力循环遍历从每个预定义范围取得的所有可能的DoF值.代码的核心是基于最小化每个点A到任何点的距离B
代码有效(它实际上生成了上面提到的解决方案),但由于解决方案的数量(即每个DoF的可接受值的组合)与更大的范围进行扩展,因此它可能会变得非常快速(也会使所有的我系统中的RAM)
如何提高代码的性能？我愿意接受任何解决方案,包括numpy和/或scipy.或许类似于Basing-Hopping来搜索最佳匹配(或相对接近的匹配)而不是我目前使用的强力方法？
import numpy as np
from scipy.spatial import distance
import math
def scalePoints(B_center, delta_x, delta_y, scale):
"""
Scales xy points.
http://codereview.stackexchange.com/q/159183/35351
"""
x_scale = B_center[0] - scale * delta_x
y_scale = B_center[1] - scale * delta_y
return x_scale, y_scale
def rotatePoints(center, x, y, angle):
"""
Rotates points in 'xy' around 'center'. Angle is in degrees.
Rotation is counter-clockwise
http://stackoverflow.com/a/20024348/1391441
"""
angle = math.radians(angle)
xy_rot = x - center[0], y - center[1]
xy_rot = (xy_rot[0] * math.cos(angle) - xy_rot[1] * math.sin(angle),
xy_rot[0] * math.sin(angle) + xy_rot[1] * math.cos(angle))
xy_rot = xy_rot[0] + center[0], xy_rot[1] + center[1]
return xy_rot
def distancePoints(set_A, x_transl, y_transl):
"""
Find the sum of the minimum distance of points in set_A to points in set_B.
"""
d = distance.cdist(set_A, zip(*[x_transl, y_transl]), 'euclidean')
# Sum of all minimal distances.
d_sum = np.sum(np.min(d, axis=1))
return d_sum
def match_frames(
set_A, B_center, delta_x, delta_y, tol, sc_min, sc_max, sc_step,
ang_min, ang_max, ang_step, xmin, xmax, xstep, ymin, ymax, ystep):
"""
Process all possible solutions in the defined ranges.
"""
N = len(set_A)
# Ranges
sc_range = np.arange(sc_min, sc_max, sc_step)
ang_range = np.arange(ang_min, ang_max, ang_step)
x_range = np.arange(xmin, xmax, xstep)
y_range = np.arange(ymin, ymax, ystep)
print("Total solutions: {:.2e}".format(
np.prod([len(_) for _ in [sc_range, ang_range, x_range, y_range]])))
d_sum, params_all = [], []
for scale in sc_range:
# Scaled points.
x_scale, y_scale = scalePoints(B_center, delta_x, delta_y, scale)
for ang in ang_range:
# Rotated points.
xy_rot = rotatePoints(B_center, x_scale, y_scale, ang)
# x translation
for x_tr in x_range:
x_transl = xy_rot[0] + x_tr
# y translation
for y_tr in y_range:
y_transl = xy_rot[1] + y_tr
# Find minimum distance sum.
d_sum.append(distancePoints(set_A, x_transl, y_transl))
# Store solutions.
params_all.append([scale, ang, x_tr, y_tr])
# Condition to break out if given tolerance for match
# is achieved.
if d_sum[-1] <= tol * N:
print("Match found:", scale, ang, x_tr, y_tr)
i_min = d_sum.index(min(d_sum))
return i_min, params_all
# Print best solution found so far.
i_min = d_sum.index(min(d_sum))
print("d_sum_min = {:.2f}".format(d_sum[i_min]))
return i_min, params_all
# Data.
set_A = [[2015.81, 1981.665], [1967.31, 1960.865], [1962.91, 1951.365],
[1964.91, 1994.565], [1894.41, 1957.065]]
set_B = [
[2689.28, 3507.04, 2895.67, 1051.3, 1929.49, 1035.97, 752.44, 130.62,
620.06, 2769.06, 1580.77, 281.76, 224.54, 3848.3],
[2061.19, 3700.27, 2131.2, 1837.3, 2017.52, 80.96, 3524.61, 3821.22,
3711.53, 1812.12, 1868.33, 3865.77, 3273.77, 2100.71]]
# This is necessary to apply the scaling.
x, y = np.asarray(set_B)
# Center of B points, defined as the center of the minimal rectangle that
# contains all points.
B_center = [(min(x) + max(x)) * .5, (min(y) + max(y)) * .5]
# Difference between the center coordinates and the xy points.
delta_x, delta_y = B_center[0] - x, B_center[1] - y
# Tolerance in pixels for match.
tol = 1.
# Ranges for each DoF.
sc_min, sc_max, sc_step = .01, .2, .01
ang_min, ang_max, ang_step = -30., 30., 1.
xmin, xmax, xstep = -150., 150., 1.
ymin, ymax, ystep = -150., 150., 1.
# Find proper scaling + rotation + translation for set_B.
i_min, params_all = match_frames(
set_A, B_center, delta_x, delta_y, tol, sc_min, sc_max, sc_step,
ang_min, ang_max, ang_step, xmin, xmax, xstep, ymin, ymax, ystep)
# Best match found
print(params_all[i_min])
我为什么要这样做？
当天文学家观测到一个恒星场时,它还需要观察所谓的"恒星标准场".这需要能够将观测恒星的"仪器量值"(亮度的对数度量)转换为通用的通用比例,因为这些量值取决于所使用的光学系统(望远镜+ CCD阵列).在此处显示的示例中,标准字段可以在左下方看到,而观察到的字段在右侧.
请注意,集合中的点A(上面使用的)是标准字段中标记的星号,集合B是在观察区域中检测到的那些星星(上面用红色标记)
即使在今天,识别观察区域中与标准区域中标记的恒星相对应的那些恒星的过程也是通过眼睛进行的.这是由于任务的复杂性.
在上面观察到的图像中,有相当多的缩放,但没有旋转和很少的平移.这构成了一个相当有利的方案; 它可能会更糟糕.我正在尝试开发一种简单的算法,以避免必须逐个手动识别观察视野中的恒星作为标准场中的恒星.
由litepresence提出的解决方案
这是我按照litepresence的回答制作的脚本.
import itertools
import numpy as np
import matplotlib.pyplot as plt
def getTriangles(set_X, X_combs):
"""
Inefficient way of obtaining the lengths of each triangle's side.
Normalized so that the minimum length is 1.
"""
triang = []
for p0, p1, p2 in X_combs:
d1 = np.sqrt((set_X[p0][0] - set_X[p1][0]) ** 2 +
(set_X[p0][1] - set_X[p1][1]) ** 2)
d2 = np.sqrt((set_X[p0][0] - set_X[p2][0]) ** 2 +
(set_X[p0][1] - set_X[p2][1]) ** 2)
d3 = np.sqrt((set_X[p1][0] - set_X[p2][0]) ** 2 +
(set_X[p1][1] - set_X[p2][1]) ** 2)
d_min = min(d1, d2, d3)
d_unsort = [d1 / d_min, d2 / d_min, d3 / d_min]
triang.append(sorted(d_unsort))
return triang
def sumTriangles(A_triang, B_triang):
"""
For each normalized triangle in A, compare with each normalized triangle
in B. find the differences between their sides, sum their absolute values,
and select the two triangles with the smallest sum of absolute differences.
"""
tr_sum, tr_idx = [], []
for i, A_tr in enumerate(A_triang):
for j, B_tr in enumerate(B_triang):
# Absolute value of lengths differences.
tr_diff = abs(np.array(A_tr) - np.array(B_tr))
# Sum the differences
tr_sum.append(sum(tr_diff))
tr_idx.append([i, j])
# Index of the triangles in A and B with the smallest sum of absolute
# length differences.
tr_idx_min = tr_idx[tr_sum.index(min(tr_sum))]
A_idx, B_idx = tr_idx_min[0], tr_idx_min[1]
print("Smallest difference: {}".format(min(tr_sum)))
return A_idx, B_idx
# Data.
set_A = [[2015.81, 1981.665], [1967.31, 1960.865], [1962.91, 1951.365],
[1964.91, 1994.565], [1894.41, 1957.065]]
set_B = [
[2689.28, 3507.04, 2895.67, 1051.3, 1929.49, 1035.97, 752.44, 130.62,
620.06, 2769.06, 1580.77, 281.76, 224.54, 3848.3],
[2061.19, 3700.27, 2131.2, 1837.3, 2017.52, 80.96, 3524.61, 3821.22,
3711.53, 1812.12, 1868.33, 3865.77, 3273.77, 2100.71]]
set_B = zip(*set_B)
# All possible triangles.
A_combs = list(itertools.combinations(range(len(set_A)), 3))
B_combs = list(itertools.combinations(range(len(set_B)), 3))
# Obtain normalized triangles.
A_triang, B_triang = getTriangles(set_A, A_combs), getTriangles(set_B, B_combs)
# Index of the A and B triangles with the smallest difference.
A_idx, B_idx = sumTriangles(A_triang, B_triang)
# Indexes of points in A and B of the best match triangles.
A_idx_pts, B_idx_pts = A_combs[A_idx], B_combs[B_idx]
print 'triangle A %s matches triangle B %s' % (A_idx_pts, B_idx_pts)
# Matched points in A and B.
print "A:", [set_A[_] for _ in A_idx_pts]
print "B:", [set_B[_] for _ in B_idx_pts]
# Plot
A_pts = zip(*[set_A[_] for _ in A_idx_pts])
B_pts = zip(*[set_B[_] for _ in B_idx_pts])
plt.scatter(*A_pts, s=10, c='k')
plt.scatter(*B_pts, s=10, c='r')
plt.show()
该方法几乎是即时的,并产生正确的匹配:
Smallest difference: 0.0314154749597
triangle A (0, 1, 4) matches triangle B (3, 4, 10)
A: [[2015.81, 1981.665], [1967.31, 1960.865], [1894.41, 1957.065]]
B: [(1051.3, 1837.3), (1929.49, 2017.52), (1580.77, 1868.33)]
1) 我会通过标记所有点并从每组中找到 3 个点的所有可能组合来解决这个问题。
# normalize B data to same format as A
set_Bx, set_By = (set_B)
set_B = []
for i in range(len(set_Bx)):
set_B.append([set_Bx[i],set_By[i]])
'''
set_B = [[2689.28, 2061.19], [3507.04, 3700.27], [2895.67, 2131.2],
[1051.3, 1837.3], [1929.49, 2017.52], [1035.97, 80.96], [752.44,
3524.61], [130.62, 3821.22], [620.06, 3711.53], [2769.06, 1812.12],
[1580.77, 1868.33], [281.76, 3865.77], [224.54, 3273.77], [3848.3,
2100.71]]
'''
list(itertools.combinations(range(len(set_A)), 3))
list(itertools.combinations(range(len(set_B)), 3))
2）对于每3个点组，计算各自三角形的边；对 A 组和 B 组重复该过程。
dist = sqrt( (x2 - x1)**2 + (y2 - y1)**2 )
3）然后减少每个三角形的边比，使每个三角形的最小边等于1；其他双方的比例均有所减少。
在两个相似三角形中：
两个三角形的周长与边长之比相同。相应的边、中线和高度都将采用相同的比例。
http://www.mathopenref.com/similartrianglesparts.html
4) 最后，对于 A 组中的每个三角形，与 B 组中的每个三角形进行逐元素减法比较。然后将所得元素相加，并从 A 和 B 中找出总和最小的三角形。
list(numpy.array(list1)-numpy.array(list2))
5）给定匹配的三角形；就 CPU/RAM 而言，找到适当的缩放、平移和旋转应该相对简单。
ETA1：草稿脚本
ETA2：修补了评论中讨论的错误：使用 sum(abs()) 而不是 abs(sum())。现在可以用了，速度也很快！
'''
known correct solution
A = [[1894.41, 1957.065],[1967.31, 1960.865],[2015.81, 1981.665]]
B = [[1051.3, 1837.3],[1580.77, 1868.33],[1929.49, 2017.52]]
'''
import numpy as np
import itertools
import math
import operator
set_A = [[2015.81, 1981.665], [1967.31, 1960.865], [1962.91, 1951.365],
[1964.91, 1994.565], [1894.41, 1957.065]]
set_B = [[2689.28, 3507.04, 2895.67, 1051.3, 1929.49, 1035.97, 752.44, 130.62,
620.06, 2769.06, 1580.77, 281.76, 224.54, 3848.3],
[2061.19, 3700.27, 2131.2, 1837.3, 2017.52, 80.96, 3524.61, 3821.22,
3711.53, 1812.12, 1868.33, 3865.77, 3273.77, 2100.71]]
# normalize set B data to set A format
set_Bx, set_By = (set_B)
set_B = []
for i in range(len(set_Bx)):
set_B.append([set_Bx[i],set_By[i]])
'''
set_B = [[2689.28, 2061.19], [3507.04, 3700.27], [2895.67, 2131.2],
[1051.3, 1837.3], [1929.49, 2017.52], [1035.97, 80.96], [752.44, 3524.61],
[130.62, 3821.22], [620.06, 3711.53], [2769.06, 1812.12], [1580.77, 1868.33],
[281.76, 3865.77], [224.54, 3273.77], [3848.3, 2100.71]]
'''
print set_A
print set_B
print len(set_A)
print len(set_B)
set_A_tri = list(itertools.combinations(range(len(set_A)), 3))
set_B_tri = list(itertools.combinations(range(len(set_B)), 3))
print set_A_tri
print set_B_tri
print len(set_A_tri)
print len(set_B_tri)
'''
set_A = [[2015.81, 1981.665], [1967.31, 1960.865], [1962.91, 1951.365],
[1964.91, 1994.565], [1894.41, 1957.065]]
set_A_tri = [(0, 1, 2), (0, 1, 3), (0, 1, 4), (0, 2, 3), (0, 2, 4), (0, 3, 4),
(1, 2, 3), (1, 2, 4), (1, 3, 4), (2, 3, 4)]
'''
def distance(x1,y1,x2,y2):
return math.sqrt((x2 - x1)**2 + (y2 - y1)**2 )
def tri_sides(set_x, set_x_tri):
triangles = []
for i in range(len(set_x_tri)):
point1 = set_x_tri[i][0]
point2 = set_x_tri[i][1]
point3 = set_x_tri[i][2]
point1x, point1y = set_x[point1][0], set_x[point1][1]
point2x, point2y = set_x[point2][0], set_x[point2][1]
point3x, point3y = set_x[point3][0], set_x[point3][1]
len1 = distance(point1x,point1y,point2x,point2y)
len2 = distance(point1x,point1y,point3x,point3y)
len3 = distance(point2x,point2y,point3x,point3y)
min_side = min(len1,len2,len3)
len1/=min_side
len2/=min_side
len3/=min_side
t=[len1,len2,len3]
t.sort()
triangles.append(t)
return triangles
A_triangles = tri_sides(set_A, set_A_tri)
B_triangles = tri_sides(set_B, set_B_tri)
print A_triangles
'''
[[1.0, 5.0405616860744304, 5.822935502560814],
[1.0, 1.5542012854321234, 1.5619803879976761],
[1.0, 1.3832883678507584, 2.347214708755337],
[1.0, 1.2141910838179129, 1.4096730529373076],
[1.0, 1.1275138587537166, 2.0318412465223665],
[1.0, 1.5207417600732074, 2.3589630093994876],
[1.0, 3.2270326342163584, 4.13069930678442],
[1.0, 6.565440477766354, 6.972550347780966],
[1.0, 2.1606693015281997, 2.3635387983160885],
[1.0, 1.589425903498476, 1.846471085870448]]
'''
print B_triangles
def list_subtract(list1,list2):
return np.absolute(np.array(list1)-np.array(list2))
sums = []
threshold = 1
for i in range(len(A_triangles)):
for j in range(len(B_triangles)):
k = sum(list_subtract(A_triangles[i], B_triangles[j]))
if k < threshold:
sums.append([i,j,k])
# sort by smallest sum
sums = sorted(sums, key=operator.itemgetter(2))
print sums
print 'winner %s' % sums[0]
print sums[0][0]
print sums[0][1]
match_A = set_A_tri[sums[0][0]]
match_B = set_B_tri[sums[0][1]]
print 'triangle A %s matches triangle B %s' % (match_A, match_B)
match_A_pts = []
match_B_pts = []
for i in range(3):
match_A_pts.append(set_A[match_A[i]])
match_B_pts.append(set_B[match_B[i]])
print 'triangle A has points %s' % match_A_pts
print 'triangle B has points %s' % match_B_pts
'''
winner [2, 204, 0.031415474959738399]
2
204
triangle A (0, 1, 4) matches triangle B (3, 4, 10)
triangle A has points [[2015.81, 1981.665], [1967.31, 1960.865], [1894.41, 1957.065]]
triangle B has points [[1051.3, 1837.3], [1929.49, 2017.52], [1580.77, 1868.33]]
'''
| 归档时间： |
|
| 查看次数： |
1174 次 |
| 最近记录： |