我有一个300000列表(光纤轨道)的列表,其中每个轨道是(x,y,z)元组/坐标的列表：
tracks=
[[(1,2,3),(3,4),...]
[(4,1),(5,7,...]
...
]
我还有一组蒙版,其中每个蒙版被定义为(x,z)元组/坐标的列表：
mask_coords_list=
[[(1,(8,13,...]
[(6,2),...]
...
]
我试图找到所有可能的面具对：
>与每个掩码 – 掩码对相交的轨道数(以创建连接矩阵)
>与每个蒙版相交的轨道子集,以便为子集中的每个轨道的每个(x,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_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,z in track:
vox_tracks_img[x,z,count] += 1
使用集合查找交叉点已大大加快了这一过程,但当我有70个或更多掩码的列表时,这两个部分仍然需要一个多小时.有没有比为每个轨道迭代更有效的方法？
最佳答案
线性化体素坐标,并将它们放入两个scipy.sparse.sparse.csc矩阵中.
原文链接：https://www.f2er.com/python/439001.html设v是体素的数量,m是掩模的数量,t是轨道的数量.
令M为掩模csc矩阵,大小(m×v),其中a at(i,j)表示掩模i与体素j重叠.
令T为轨道csc矩阵,大小(t×v),其中a(k,j)表示轨道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可以通过非零和&采取.您可能需要在循环中拉出每列并将其转换为完整矩阵.
我做了一些试验,试图模拟我认为合理数量的数据.在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,nmask * 10000).astype(int)
d = ones(nmask * 10000)
M = csc_matrix((d,vstack((mcoords,shape=(nmask,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