更有效的交叉路口计算方法？

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(以创建"密度"图像)

我现在正在做第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

和第2部分一样:

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个或更多掩码的列表时,这两个部分仍然需要一个多小时.有没有比为每个轨道迭代更有效的方法？

Answer 1

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

令 v 为体素数，m 为掩模数，t 为轨道数。
设 M 为掩模 csc 矩阵，大小为 (mxv)，其中 (i,j) 处的 1 表示掩模 i 与体素 j 重叠。
令 T 为轨迹 csc 矩阵，大小为 (txv)，其中 (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 是否可以通过 nonzero & take 进行操作。您可能需要在循环中取出每一列并将其转换为完整矩阵。

我进行了一些实验，试图模拟我认为合理的数据量。下面的代码在使用了 2 年的 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