python - 获取每个坐标的对角线长度的更有效方法

Question

我有一个代表匹配的 x 和 y 值（坐标）数组，对于这些 x,y 中的每一个，我想知道它所属的对角线的长度。例如，让我们取这些坐标

资料说明

coords = np.asarray([[0,0], [0,7], [1,1], [1,6], [2,2], [2,5], [3,3],[3,4], [4,4]])
# [[0 0]
#  [0 7]
#  [1 1]
#  [1 6]
#  [2 2]
#  [2 5]
#  [3 3]
#  [3 4]
#  [4 4]]

我们可以将其转换为矩阵，但在我使用大量表格的情况下，这太低效了（例如，scipytodia()会抛出低效警告；见下文）。无论如何，让我们制作矩阵以使问题更清楚：

[[1 0 0 0 0 0 0 1]
 [0 1 0 0 0 0 1 0]
 [0 0 1 0 0 1 0 0]
 [0 0 0 1 1 0 0 0]
 [0 0 0 0 1 0 0 0]]

目标
查看上表，我们看到两条对角线（或一条对角线和一条对角线）。对于对角线的每个位置，我想知道它所在的对角线的长度，所以像这样的表格：

# x, y, diag length
[[0 0 5]
 [1 1 5]
 [2 2 5]
 [3 3 5]
 [4 4 5]
 [3 4 4]
 [2 5 4]
 [1 6 4]
 [0 7 4]]

低效的解决方案我认为我可以在稀疏的 scipy 矩阵
中表示这些数据，而这给出了将稀疏矩阵转换为对角坐标矩阵的所需结果对于 100 个对角线来说已经是低效的，更不用说我拥有的数千个对角线了。

from scipy.sparse import dia_matrix, coo_matrix
coords = np.asarray([[0,0], [0,7], [1,1], [1,6], [2,2], [2,5], [3,3],[3,4], [4,4]])

# Create the scipy coord matrix
x = coords[:,0]
y = coords[:,1]
tot_elem = coords.shape[0]*2
data = np.repeat(1, len(x))
co_mat = coo_matrix( (data, (x, y)), shape=(max(x)+1, max(y)+1))

# Get the diagonal matrix
dia_mat = dia_matrix(co_mat).tocoo()
diag_coords = np.column_stack((dia_mat.row, dia_mat.col))

# Get the consecutive values to put them to lengths
difs = np.diff(diag_coords[:, 1])
cuts = [0] + list(np.where(difs != 1)[0] + 1) + [diag_coords.shape[0]]
sizes = np.diff(cuts)
sizes = np.repeat(sizes, sizes)

# Combine with the original coords
dia_sizes = np.column_stack((dia_mat.row, dia_mat.col, sizes))
print(dia_sizes)

*刚刚意识到坐标可以是对角线和对角线的一部分，在这种情况下，我可以报告两者或只报告最长对角线的长度 - 我的解决方案没有处理:(

编辑： 更有效的解决方案

在这里查看todia() 代码，我注意到他们使用了一个聪明的技巧来查看点是否在对角线上，即x-y对于同一对角线上的点应该是相同的。然而，这对于反对角线是不正确的。所以我假设相反，x + y确实给了我们在同一个对角线上的观点。使用这个我想出了已经比使用 scipy 快得多的代码。

import numpy as np

coords = np.asarray([[0,0], [0,7], [1,1], [1,6], [2,2], [2,5], [3,3],[3,4], [4,4]])
x = coords[:,0]
y = coords[:,1]

# Get the diagonal (inspired by scripy todia code)
ks1 = y - x

# Unlike scipy, I think we can do the same by summing to get the anti-diagonal
ks2 = y + x

# Sort these to get the groups in the same diagonal
idx = np.argsort(ks1)
anti_idx = np.argsort(ks2)

def get_dia_len(arr,ori):
    sizes = np.diff([0] + list(np.where(np.diff(arr)!= ori)[0] + 1) + [arr.shape[0]])
    size_arr = np.repeat(sizes, sizes)
    return size_arr

# Get the diagonal lengths, i.e. cut at changing values and get the gaps between them
norm_sizes = get_dia_len(x[idx],1)
anti_sizes = get_dia_len(y[anti_idx],-1)

# Gather this in a table
norm = np.column_stack([x[idx], y[idx], norm_sizes])
anti = np.column_stack([x[anti_idx], y[anti_idx], anti_sizes])
dia_coord = np.concatenate((norm, anti))

# We only have a diagonal when we have >1 value
dia_coord = dia_coord[dia_coord[:, -1] > 1]
print(dia_coord)

一段时间以来我一直在低头，很想知道是否有人有聪明的方法来解决这个问题:)

标签： pythonnumpymatrixscipydiagonal

Answer 1

一种方法是遍历坐标并通过每个点构建45 度线（假设这就是“对角线”的含义），然后从coords列表中删除位于这条线上的任何点 -

此函数计算固定点的45 度线上的点，并仅返回coords列表中的点

coords = [[0,0], [0,7], [1,1], [1,6], [2,2], [2,5], [3,3],[3,4], [4,4]]
coords = [tuple(_) for _ in coords]

def get_y(x, fixed_point, allowed_slopes=(1, -1), coords=coords.copy()):
    coords = [tuple(_) for _ in coords]
    x_fixed, y_fixed = fixed_point
    possible_y = [y_fixed + slope*(x - x_fixed) for slope in allowed_slopes]
    possible_coords = [(x, y) for y in possible_y]
    available_coords = list(set(possible_coords) & set(coords))
    return available_coords
print(get_y(1, (0,0)))
#[(1, 1)]
print(get_y(6, (0,0)))
#[] because (6, 6) is not on coords

然后我们可以循环遍历coords，同时删除同一行上的所有点。使用list.pop确保我们不必为同一组点多次不必要地计算对角线

idx = 0
grouped_points = list()
while coords:
    group = list()
    fixed_point = coords.pop()
    print(f'fixed_point is now {fixed_point}')
    group.append(fixed_point)
    print(f'group is now {group}')
    available_x = set([x for (x, y) in coords])
    print(f'available_x is now {available_x}')
    for x in available_x:
        pt, *_ = get_y(x, fixed_point)
        print(f'pt is now {pt}')
        if pt and pt in coords:
            group.append(pt)
            coords.remove(pt)
        print(f'coords is now {coords}')
        print(f'group is now {group}')
    print(idx, group, sep='\t')
    grouped_points.append(group)
    idx += 1

然后将长度附加到输出以获得所需的结果

grouped_points = [(*pt, len(group)) for group in grouped_points for pt in group]
print(*grouped_points, sep='\n')
#(4, 4, 5)
#(0, 0, 5)
#(1, 1, 5)
#(2, 2, 5)
#(3, 3, 5)
#(3, 4, 4)
#(0, 7, 4)
#(1, 6, 4)
#(2, 5, 4)

计时这个使用timeit表明这个解决方案比这组快 10 倍coords