查找连续的未掩盖值

Question

我有一个很大的三维数组，这个数组的维度分别是时间、经度和纬度。大部分数据都是被遮罩的。我需要找出那些在特定时间步长内，遮罩状态为假（也就是没有被遮住）的数据，这个特定的时间步长我称之为 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 再转换成一个遮罩。

from scipy.ndimage import measurements
structure = np.zeros((3, 3, 3))
structure[:, 1, 1] = 1
labels, nr_labels = measurements.label(1 - mask, structure=structure)
_, counts = np.unique(labels, return_counts=True)
labels_selected = [i_count for i_count, count in enumerate(counts)
                   if count >= threshold]

相关文章：

• 暂无相关问题

Answer 1

为了完整起见，这里也提供了我在编辑中提到的解决方案。虽然在性能上比Divakar的两个方案差得多（大约慢10倍，跟numpy_binary_closing比），但它可以处理三维数组。此外，它还提供了输出聚类位置的功能（虽然这不是问题的一部分，但可能是有趣的信息）。

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 模块，特别是 scipy.ndimage.morphology.binary_closing，在这里我们需要提供一个合适的全是 1 的一维核，长度为 threshold。另外，Scipy 的 binary closing 函数只会给我们闭合后的掩膜。所以，要得到想要的输出，我们需要将其与输入掩膜进行 或运算。因此，具体的实现看起来会像这样 -

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

那我们来看看 NumPy 版本的闭合操作吧!

闭合操作其实就是 图像膨胀 和 图像腐蚀，所以我们可以用可靠的卷积操作来模拟这个过程，而在 NumPy 中我们可以使用 np.convolve。和 scipy 的二元闭合操作一样，我们在这里也需要用到相同的核，并且在膨胀和腐蚀时都使用它。具体的实现会是 -

def numpy_binary_closing(mask,threshold):

    # Define kernel
    K = np.ones(threshold)

    # Perform dilation and threshold at 1
    dil = np.convolve(mask,K,mode='same')>=1

    # Perform erosion on the dilated mask array and threshold at given threshold
    dil_erd = np.convolve(dil,K,mode='same')>= threshold
    return dil_erd

示例运行 -

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

案例 #2 : 更稀疏且阈值更大

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，而且差距很大!

---

边界条件

看起来如果 1 离边界 足够近，边界也需要进行闭合。要解决这些情况，你可以在输入布尔数组的开头和结尾各加一个 1，然后使用之前的代码，最后再去掉第一个和最后一个元素。因此，使用 scipy 的 binary_closing 方法的完整实现会是 -

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)