更高效的计算交点的方法？

Question

我有一个包含300000个列表（纤维束轨迹）的列表，其中每个轨迹是由(x,y,z)元组/坐标组成的列表：

tracks=
[[(1,2,3),(3,2,4),...]
 [(4,2,1),(5,7,3),...]
 ...
]

我还有一个面具组，每个面具都被定义为(x,y,z)元组/坐标的列表：

mask_coords_list=
[[(1,2,3),(8,13,4),...]
 [(6,2,2),(5,7,3),...]
 ...
]

我正在尝试为所有可能的掩模对找到以下内容：

1. 相交于每个掩模-掩模对的轨道数量（以创建连通性矩阵）
2. 与每个掩模相交的轨道子集，以便为子集中的每个轨道的（x，y，z）坐标加1（以创建“密度”图像）

我目前正在这样做第一部分：

def mask_connectivity_matrix(tracks,masks,masks_coords_list):
    connect_mat=zeros((len(masks),len(masks)))
    for track in tracks:
        cur=[]
        for count,mask_coords in enumerate(masks_coords_list):
            if any(set(track) & set(mask_coords)):
                cur.append(count)
            for x,y in list(itertools.combinations(cur,2)):
                connect_mat[x,y] += 1

并且第二部分如下：

def mask_tracks(tracks,masks,masks_coords_list):
    vox_tracks_img=zeros((xdim,ydim,zdim,len(masks)))
    for track in tracks:
        for count,mask in enumerate(masks_coords_list):
            if any(set(track) & set(mask)):
                for x,y,z in track:
                    vox_tracks_img[x,y,z,count] += 1

使用集合来查找交集显着加快了此过程，但是当我有70个或更多掩码的列表时，这两部分仍需要超过一小时。是否有比为每个轨道迭代更有效的方法？

- jbrown

---

所有的答案似乎只是些微小的改进，但我认为你需要更多。 Translated text: - McPherrinM

如果您能够在Pastebin上发布一个样本数据集和正确答案，那么您可能会得到更多的帮助。 - Jason Orendorff

我看到的是，交叉点仅被定义为两个坐标元组相同，而不是坐标之间相交的线段，这是正确的吗？ - Svante

没错，轨迹和掩码中的相同元组构成一个交集。 - jbrown

Answer 1

将体素坐标线性化，并将它们放入两个scipy.sparse.sparse.csc矩阵中。

假设有v个体素，m个掩模和t条轨迹。 令M为掩模csc矩阵，大小为（m x v），其中（i，j）处的1表示掩模i与体素j重叠。 令T为轨迹csc矩阵，大小为（t x v），其中（k，j）处的1表示轨迹k与体素j重叠。

Overlap = (M * T.transpose() > 0)  # track T overlaps mask M  
Connected = (Overlap * Overlap.tranpose() > 0) # Connected masks
Density[mask_idx] = numpy.take(T, nonzero(Overlap[mask_idx, :])[0], axis=0).sum(axis=0)

我可能在最后一个问题上错了，而且我不确定css_matrices能否被非零和take操作。你可能需要在循环中提取每一列并将其转换为完整矩阵。

---

我运行了一些实验，试图模拟我认为合理的数据量。下面的代码在一台两年前的MacBook上大约需要2分钟。如果使用csr_matrices，则需要大约4分钟。这可能取决于每个轨道的长度。

from numpy import *
from scipy.sparse import csc_matrix

nvox = 1000000
ntracks = 300000
nmask = 100

# create about 100 entries per track
tcoords = random.uniform(0, ntracks, ntracks * 100).astype(int)
vcoords = random.uniform(0, nvox, ntracks * 100).astype(int)
d = ones(ntracks * 100)
T = csc_matrix((d,  vstack((tcoords, vcoords))), shape=(ntracks, nvox), dtype=bool)

# create around 10000 entries per mask
mcoords = random.uniform(0, nmask, nmask * 10000).astype(int)
vcoords = random.uniform(0, nvox, nmask * 10000).astype(int)
d = ones(nmask * 10000)
M = csc_matrix((d, vstack((mcoords, vcoords))), shape=(nmask, nvox), dtype=bool)

Overlap = (M * T.transpose()).astype(bool) # mask M overlaps track T
Connected = (Overlap * Overlap.transpose()).astype(bool) # mask M1 and M2 are connected
Density = Overlap * T.astype(float) # number of tracks overlapping mask M summed across voxels

- user227667

---

如果矩阵的数据类型设置为布尔型，我认为“> 0”位就不再必要了。 - user227667

实际上，这并不正确。至少对于稀疏矩阵来说，乘法会将它们提升为一个字节。我希望这并不意味着还存在环绕问题。 - user227667

谢谢，这让我在平均曲长约为10和平均掩码大小约为500的情况下加速到不到一分钟。 - jbrown

很高兴能帮忙。你有验证过慢方法以确保环绕不是问题吗？我猜没有，因为numpy非常稳定，但我喜欢对这些事情小心谨慎（而且我想知道答案以备将来参考）。 - user227667

实际上，这确实导致了环绕。在发现连接矩阵中出现负数并需要弄清它们来自何处之前，我不熟悉这个概念。我不确定csc_matrix乘法的最大大小上限是多少，但似乎确实存在限制-我的T矩阵为140000x2000000，而我的M矩阵为150 x 2000000。 - jbrown

Answer 2

好的，我认为我终于找到了一些可以减少复杂性的方法。相比你现在拥有的代码，这个代码应该会运行得非常快。

似乎首先你需要知道哪些轨道与哪些掩模重叠，也就是关联矩阵。

import numpy
from collections import defaultdict

def by_point(sets):
    d = defaultdict(list)
    for i, s in enumerate(sets):
        for pt in s:
            d[pt].append(i)
    return d

def calc(xdim, ydim, zdim, mask_coords_list, tracks):
    masks_by_point = by_point(mask_coords_list)
    tracks_by_point = by_point(tracks)

    a = numpy.zeros((len(mask_coords_list), len(tracks)), dtype=int)
    for pt, maskids in masks_by_point.iteritems():
        for trackid in tracks_by_point.get(pt, ()):
            a[maskids, trackid] = 1
    m = numpy.matrix(a)

你需要的邻接矩阵是 m * m.T。

你已经编写的代码只计算了上三角。你可以使用triu函数来仅获取该部分。

    am = m * m.T  # calculate adjacency matrix
    am = numpy.triu(am, 1)  # keep only upper triangle
    am = am.A  # convert matrix back to array

体素计算也可以使用关联矩阵。

    vox_tracks_img = numpy.zeros((xdim, ydim, zdim, len(mask_coords_list)), dtype=int)
    for trackid, track in enumerate(tracks):
        for x, y, z in track:
            vox_tracks_img[x, y, z, :] += a[:,trackid]
    return am, vox_tracks_img

对于我来说，处理数百个掩模和轨迹的数据集只需要不到一秒钟的时间。

如果您有许多出现在掩模中但不在任何轨迹上的点，则在进入循环之前从 masks_by_point 中删除这些点的条目可能是值得的。

- Jason Orendorff

Answer 3

你可以尝试将这两个函数结合起来，同时生成两个结果。此外，在循环之前无需创建组合列表，因为它已经是一个生成器，这可能会节省你一些时间。

def mask_connectivity_matrix_and_tracks(tracks,masks,masks_coords_list):
    connect_mat=zeros((len(masks),len(masks)))
    vox_tracks_img=zeros((xdim,ydim,zdim,len(masks)))
    for track in tracks:
        cur=[]
        for count,mask_coords in enumerate(masks_coords_list):
            if any(set(track) & set(mask_coords)):
                cur.append(count)
                for x,y,z in track:
                    vox_tracks_img[x,y,z,count] += 1
            for x,y in itertools.combinations(cur,2):
                connect_mat[x,y] += 1

此外，这可能永远不会像“在我们死之前完成”那样“快”，因此最好的方法是最终使用Cython将其编译为Python的C模块。

- Tor Valamo

Answer 4

如果你将每个掩码点集(1,2,3),(1,2,4),(1,3,1)存储为如下字典：{1: [{2: set([3, 4])}, {3: set([1])}]}，也许可以更快地检查匹配...但也可能不行。

- Brian

Answer 5

通过删除冗余操作，可以进行轻微的优化（相同的大O符号，略小的乘数）：

1. 不要在每个轨道和掩码上调用set这么多次：每个轨道和每个掩码只调用一次，以设置辅助“并行”集合列表，然后对其进行操作
2. if any(someset):在语义上与if someset:相同，但速度稍慢

不会有太大的差异，但可能会稍微有所帮助。

- Alex Martelli

Answer 6

我知道这听起来很无聊，但是：

使用Python的长整型可以将小整数集建模为位向量。假设您用小整数ID替换每个元组，然后将每个轨道和每组掩码坐标转换为这些小ID的集合。您可以将这些集合表示为长整数，从而使交集操作变得更快（但不是渐近地更快）。

- Jason Orendorff