查找连续的无掩码值

Question

我有一个大的三维（时间，经度，纬度）输入数组。大多数条目都被屏蔽了。我需要找到那些掩码为False的条目，其时间超过了特定的连续时间步数（这里我称之为threshold）。结果应该是一个与输入掩码形状相同的掩码。

下面是一些伪代码，希望能更清楚地说明我的意思：

new_mask = find_consecutive(mask, threshold=3)
mask[:, i_lon, i_lat]
# [1, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0]
new_mask[:, i_lon, i_lat]
# [1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1]

编辑：

我不确定到目前为止我的方法是否有意义。它在性能方面做得很好，给了我一个标签数组和我想要的标签的知识。我只是想不出一个有效的方法把labels再转换成一个掩码。

^{pr2}$

Answer 1

为了完整起见，这里还有我在编辑中概述的方法的解决方案。在性能方面，它确实比两个Divakars解决方案差得多（与numpy_binary_closing相比大约是10倍），但允许处理3D数组。此外，它还提供了写出集群位置的可能性（这不是问题的一部分，但它可以是有趣的信息）

import numpy as np
from scipy.ndimage import measurements

def select_consecutive(mask, threshold):

    structure = np.zeros((3, 3, 3))
    structure[:, 1, 1] = 1
    labels, _ = measurements.label(1 - mask, structure=structure)

    # find positions of all unmasked values
    # object_slices = measurements.find_objects(labels)

    _, counts = np.unique(labels, return_counts=True)
    labels_selected = [i_count for i_count, count in enumerate(counts)
                       if count >= threshold and i_count != 0]
    ind = np.in1d(labels.flatten(), labels_selected).reshape(mask.shape)

    mask_new = np.ones_like(mask)
    mask_new[ind] = 0

    return mask_new

Answer 2

这是^{}的一个经典例子。要解决这个问题，您可以从scipy模块获得帮助，特别是-^{}，在我们提供了一个适当的1D内核，该内核的长度为ONES，长度为threshold。{Scipy>也只给出了Scipy}函数。因此，为了得到所需的输出，我们需要使用输入掩码OR。因此，实现将如下所示-

from scipy.ndimage import binary_closing
out = mask | binary_closing(mask, structure=np.ones(threshold))

一个新版本的二进制关闭怎么样！

现在，一个结束操作基本上是^{} and ^{}，所以我们可以用可信卷积操作来模拟这种行为，在NumPy中我们有^{}。与scipy的二进制关闭操作类似，我们在这里也需要相同的内核，我们将使用它来进行膨胀和腐蚀。实施将是-

^{pr2}$

样本运行-

In [133]: mask
Out[133]: 
array([ True, False, False, False, False,  True,  True, False, False,
        True, False], dtype=bool)

In [134]: threshold = 3

In [135]: binary_closing(mask, structure=np.ones(threshold))
Out[135]: 
array([False, False, False, False, False,  True,  True,  True,  True,
        True, False], dtype=bool)

In [136]: numpy_binary_closing(mask,threshold)
Out[136]: 
array([False, False, False, False, False,  True,  True,  True,  True,
        True, False], dtype=bool)

In [137]: mask | binary_closing(mask, structure=np.ones(threshold))
Out[137]: 
array([ True, False, False, False, False,  True,  True,  True,  True,
        True, False], dtype=bool)

In [138]: mask| numpy_binary_closing(mask,threshold)
Out[138]: 
array([ True, False, False, False, False,  True,  True,  True,  True,
        True, False], dtype=bool)

运行时测试（Scipy vs Numpy！）

情形1：一致稀疏

In [163]: mask = np.random.rand(10000) > 0.5

In [164]: threshold = 3

In [165]: %timeit binary_closing(mask, structure=np.ones(threshold))
1000 loops, best of 3: 582 µs per loop

In [166]: %timeit numpy_binary_closing(mask,threshold)
10000 loops, best of 3: 178 µs per loop

In [167]: out1 = binary_closing(mask, structure=np.ones(threshold))

In [168]: out2 = numpy_binary_closing(mask,threshold)

In [169]: np.allclose(out1,out2) # Verify outputs
Out[169]: True

案例二：稀疏度越高阈值越大

In [176]: mask = np.random.rand(10000) > 0.8

In [177]: threshold = 11

In [178]: %timeit binary_closing(mask, structure=np.ones(threshold))
1000 loops, best of 3: 823 µs per loop

In [179]: %timeit numpy_binary_closing(mask,threshold)
1000 loops, best of 3: 331 µs per loop

In [180]: out1 = binary_closing(mask, structure=np.ones(threshold))

In [181]: out2 = numpy_binary_closing(mask,threshold)

In [182]: np.allclose(out1,out2) # Verify outputs
Out[182]: True

胜利者是Numpy，而且差距很大！在

---

似乎边界也需要关闭，如果1s离花盆足够近。要解决这些情况，可以在输入布尔数组的开头和结尾各填充一个1，使用发布的代码，然后在末尾取消选择第一个和最后一个元素。因此，使用scipy的二进制关闭方法的完整实现将是-

mask_ext = np.pad(mask,1,'constant',constant_values=(1))
out = mask_ext | binary_closing(mask_ext, structure=np.ones(threshold))
out = out[1:-1]

样本运行-

In [369]: mask
Out[369]: 
array([False, False,  True, False, False, False, False,  True,  True,
       False, False,  True, False], dtype=bool)

In [370]: threshold = 3

In [371]: mask_ext = np.pad(mask,1,'constant',constant_values=(1))
     ...: out = mask_ext | binary_closing(mask_ext, structure=np.ones(threshold))
     ...: out = out[1:-1]
     ...: 

In [372]: out
Out[372]: 
array([ True,  True,  True, False, False, False, False,  True,  True,
        True,  True,  True,  True], dtype=bool)