python - 将 argsort 与掩码相结合以在移动窗口中获取最接近的值
Posted
技术标签:
【中文标题】python - 将 argsort 与掩码相结合以在移动窗口中获取最接近的值【英文标题】:python - combining argsort with masking to get nearest values in moving window 【发布时间】:2014-12-11 23:59:53 【问题描述】:我有一些代码用于根据 2D 圆形窗口中的相邻值计算图像中的缺失值。它还使用来自相同位置的一个或多个时间相邻图像的值(即在第 3 维中移动的相同 2D 窗口)。
对于每个缺失的位置,我需要计算的值不一定基于整个窗口中可用的所有值,而仅基于空间上最近的 n 个具有值的单元格(在两个图像/Z 轴中)位置),其中 n 是小于 2D 窗口中单元格总数的某个值。
此时,计算窗口中的所有内容要快得多,因为我的排序方法是使用数据获取最近的 n 个单元格是函数中最慢的部分,因为它每次都必须重复,即使距离很远在窗口坐标方面不改变。我不确定这是必要的,并且觉得我必须能够获得一次排序的距离,然后在仅选择可用单元格的过程中屏蔽这些距离。
这是我用于在间隙单元位置的窗口中选择要使用的数据的代码:
# radius will in reality be ~100
radius = 2
y,x = np.ogrid[-radius:radius+1, -radius:radius+1]
dist = np.sqrt(x**2 + y**2)
circle_template = dist > radius
# this will in reality be a very large 3 dimensional array
# representing daily images with some gaps, indicated by 0s
dataStack = np.zeros((2,5,5))
dataStack[1] = (np.random.random(25) * 100).reshape(dist.shape)
dataStack[0] = (np.random.random(25) * 100).reshape(dist.shape)
testdata = dataStack[1]
alternatedata = dataStack[0]
random_gap_locations = (np.random.random(25) * 30).reshape(dist.shape) > testdata
testdata[random_gap_locations] = 0
testdata[radius, radius] = 0
# in reality we will go through every gap (zero) location in the data
# for each image and for each gap use slicing to get a window of
# size (radius*2+1, radius*2+1) around it from each image, with the
# gap being at the centre i.e.
# testgaplocation = [radius, radius]
# and the variables testdata, alternatedata below will refer to these
# slices
locations_to_exclude = np.logical_or(circle_template, np.logical_or
(testdata==0, alternatedata==0))
# the places that are inside the circular mask and where both images
# have data
locations_to_include = ~locations_to_exclude
number_available = np.count_nonzero(locations_to_include)
# we only want to do the interpolation calculations from the nearest n
# locations that have data available, n will be ~100 in reality
number_required = 3
available_distances = dist[locations_to_include]
available_data = testdata[locations_to_include]
available_alternates = alternatedata[locations_to_include]
if number_available > number_required:
# In this case we need to find the closest number_required of elements, based
# on distances recorded in dist, from available_data and available_alternates
# Having to repeat this argsort for each gap cell calculation is slow and feels
# like it should be avoidable
sortedDistanceIndices = available_distances.argsort(kind = 'mergesort',axis=None)
requiredIndices = sortedDistanceIndices[0:number_required]
selected_data = np.take(available_data, requiredIndices)
selected_alternates = np.take(available_alternates , requiredIndices)
else:
# we just use available_data and available_alternates as they are...
# now do stuff with the selected data to calculate a value for the gap cell
这可行,但函数总时间的一半以上用于掩码空间距离数据的 argsort。 (~900uS 总共 1.4mS - 这个函数会运行数百亿次,所以这是一个重要的区别!)
我确信我必须能够在函数外部执行此 argsort 一次,当空间距离窗口最初设置时,然后将这些排序索引包含在掩码中,以获得第一个 howManyToCalculate 索引而无需重新排序。答案可能涉及将我们从中提取的各种位放入记录数组中 - 但我不知道如何,如果是这样。谁能看到我怎样才能使这部分过程更有效率?
【问题讨论】:
您的代码真的很难阅读...您可能想阅读 PEP8 并关注它:它有助于与其他 Python 程序员共享代码。 我同意 Jaime 的观点,即这很难阅读,尤其是代码,但描述也留下了一些解释空间。所以我不会冒险提供答案,但如果我愿意,我会检查一些工具(如果我至少模糊地正确理解了你的问题)。 sklearn.feature_extraction.image.extract_patches 为您提供补丁的视图,您可以将其屏蔽。它将创建一个副本,因此请注意内存问题。 您可能还对看似完全不相关的函数感兴趣,该函数使用膨胀来估算缺失值。它不会为您提供确切的结果,但可能是一个很好的代理:nilearn.masking._extrapolate_out_img 公平评论,我已经澄清了问题并用一个更简单的位替换了代码示例,该位不显示移动窗口的设置/边界的处理 记录数组或结构化数组可以方便地在单个操作中对多个数据集进行(布尔)索引操作。对于大型索引操作,这种方式也可能更快。 【参考方案1】:所以你想在循环之外进行排序:
sorted_dist_idcs = dist.argsort(kind='mergesort', axis=None)
然后使用原始代码中的一些变量,这就是我能想到的,尽管它仍然感觉像是一个主要的往返..
loc_to_incl_sorted = locations_to_include.take(sorted_dist_idcs)
sorted_dist_idcs_to_incl = sorted_dist_idcs[loc_to_incl_sorted]
required_idcs = sorted_dist_idcs_to_incl[:number_required]
selected_data = testdata.take(required_idcs)
selected_alternates = alternatedata.take(required_idcs)
注意required_idcs
指的是testdata
中的位置,而不是原始代码中的available_data
。而这个 sn-p 我使用take
是为了方便地索引扁平数组。
【讨论】:
谢谢!这回答了我在编辑后的版本中所说的问题,所以我接受了答案。但是,当间隙距离数据边缘 【参考方案2】:@moarningsun - 感谢您的评论和回答。这些让我走上了正确的轨道,但是当间隙距离数据边缘
不幸的是,当我澄清原始问题时,我编辑了我的那部分代码,但结果证明它是相关的。
我现在意识到,如果您在 argsort
的输出上再次使用 argsort
,那么您将获得排名;即对整个数组进行排序时每个项目的位置。我们可以安全地屏蔽这些,然后取其中最小的number_required
(并在结构化数组上执行此操作以同时获取相应的数据)。
这意味着循环中的另一种排序,但实际上我们可以使用分区而不是完整排序,因为我们只需要最小的num_required
项。如果num_required
大大少于数据项的数量,那么这比argsort
快得多。
例如,num_required
= 80 和 num_available
= 15000 完整的 argsort
需要 ~900µs 而argpartition
后跟索引和切片来获得前 80 个需要 ~110µs。我们仍然需要在开始时执行argsort
来获得排名(而不是仅仅基于距离进行分区),以便获得合并排序的稳定性,从而在距离不唯一时获得“正确的”。
如下所示,我的代码现在可以在大约 610uS 的实际数据上运行,包括此处未显示的实际计算。我现在对此感到满意,但似乎还有其他几个明显次要的因素会影响难以理解的运行时。
例如,将circle_template
与dist
、ranks
以及此处未显示的另一个字段一起放入结构化数组中,会使整个函数的运行时间加倍(即使我们不访问circle_template
in循环!)。更糟糕的是,与使用 np.argpartition
相比,在带有 order=['ranks']
的结构化数组上使用 np.partition
会使整个函数运行时间增加几乎 两个数量级,如下所示!
# radius will in reality be ~100
radius = 2
y,x = np.ogrid[-radius:radius+1, -radius:radius+1]
dist = np.sqrt(x**2 + y**2)
circle_template = dist > radius
ranks = dist.argsort(axis=None,kind='mergesort').argsort().reshape(dist.shape)
diam = radius * 2 + 1
# putting circle_template in this array too doubles overall function runtime!
fullWindowArray = np.zeros((diam,diam),dtype=[('ranks',ranks.dtype.str),
('thisdata',dayDataStack.dtype.str),
('alternatedata',dayDataStack.dtype.str),
('dist',spatialDist.dtype.str)])
fullWindowArray['ranks'] = ranks
fullWindowArray['dist'] = dist
# this will in reality be a very large 3 dimensional array
# representing daily images with some gaps, indicated by 0s
dataStack = np.zeros((2,5,5))
dataStack[1] = (np.random.random(25) * 100).reshape(dist.shape)
dataStack[0] = (np.random.random(25) * 100).reshape(dist.shape)
testdata = dataStack[1]
alternatedata = dataStack[0]
random_gap_locations = (np.random.random(25) * 30).reshape(dist.shape) > testdata
testdata[random_gap_locations] = 0
testdata[radius, radius] = 0
# in reality we will loop here to go through every gap (zero) location in the data
# for each image
gapz, gapy, gapx = 1, radius, radius
desLeft, desRight = gapx - radius, gapx + radius+1
desTop, desBottom = gapy - radius, gapy + radius+1
extentB, extentR = dataStack.shape[1:]
# handle the case where the gap is < search radius from the edge of
# the data. If this is the case, we can't use the full
# diam * diam window
dataL = max(0, desLeft)
maskL = 0 if desLeft >= 0 else abs(dataL - desLeft)
dataT = max(0, desTop)
maskT = 0 if desTop >= 0 else abs(dataT - desTop)
dataR = min(desRight, extentR)
maskR = diam if desRight <= extentR else diam - (desRight - extentR)
dataB = min(desBottom,extentB)
maskB = diam if desBottom <= extentB else diam - (desBottom - extentB)
# get the slice that we will be working within
# ranks, dist and circle are already populated
boundedWindowArray = fullWindowArray[maskT:maskB,maskL:maskR]
boundedWindowArray['alternatedata'] = alternatedata[dataT:dataB, dataL:dataR]
boundedWindowArray['thisdata'] = testdata[dataT:dataB, dataL:dataR]
locations_to_exclude = np.logical_or(boundedWindowArray['circle_template'],
np.logical_or
(boundedWindowArray['thisdata']==0,
boundedWindowArray['alternatedata']==0))
# the places that are inside the circular mask and where both images
# have data
locations_to_include = ~locations_to_exclude
number_available = np.count_nonzero(locations_to_include)
# we only want to do the interpolation calculations from the nearest n
# locations that have data available, n will be ~100 in reality
number_required = 3
data_to_use = boundedWindowArray[locations_to_include]
if number_available > number_required:
# argpartition seems to be v fast when number_required is
# substantially < data_to_use.size
# But partition on the structured array itself with order=['ranks']
# is almost 2 orders of magnitude slower!
reqIndices = np.argpartition(data_to_use['ranks'],number_required)[:number_required]
data_to_use = np.take(data_to_use,reqIndices)
else:
# we just use available_data and available_alternates as they are...
pass
# now do stuff with the selected data to calculate a value for the gap cell
【讨论】:
以上是关于python - 将 argsort 与掩码相结合以在移动窗口中获取最接近的值的主要内容,如果未能解决你的问题,请参考以下文章
c_cpp 以位为单位显示十六进制;与掩码一起使用按位AND。