【问题标题】:Removing completely isolated cells from Python array?从 Python 数组中删除完全隔离的单元格?
【发布时间】:2015-04-01 04:24:32
【问题描述】:

我试图通过删除所有完全隔离的单个单元格来减少二进制 python 数组中的噪声,即如果它们完全被其他“0”包围,则将“1”值单元格设置为 0。我已经能够通过使用循环删除大小等于 1 的 blob 来获得一个可行的解决方案,但这对于大型数组来说似乎是一个非常低效的解决方案:

import numpy as np
import scipy.ndimage as ndimage
import matplotlib.pyplot as plt    

# Generate sample data
square = np.zeros((32, 32))
square[10:-10, 10:-10] = 1
np.random.seed(12)
x, y = (32*np.random.random((2, 20))).astype(np.int)
square[x, y] = 1

# Plot original data with many isolated single cells
plt.imshow(square, cmap=plt.cm.gray, interpolation='nearest')

# Assign unique labels
id_regions, number_of_ids = ndimage.label(square, structure=np.ones((3,3)))

# Set blobs of size 1 to 0
for i in xrange(number_of_ids + 1):
    if id_regions[id_regions==i].size == 1:
        square[id_regions==i] = 0

# Plot desired output, with all isolated single cells removed
plt.imshow(square, cmap=plt.cm.gray, interpolation='nearest')

在这种情况下,侵蚀和扩张我的数组将不起作用,因为它也会删除宽度为 1 的特征。我觉得解决方案位于 scipy.ndimage 包中的某个位置,但是所以到目前为止,我还无法破解它。任何帮助将不胜感激!

【问题讨论】:

    标签: python numpy scipy python-2.6 ndimage


    【解决方案1】:

    非常感谢 Jaime 和 Kazemakase 的回复。手动邻居检查方法确实删除了所有孤立的补丁,但也删除了一个角落(即样本阵列中正方形的右上角)连接到其他补丁的补丁。面积求和表运行良好,在小样本阵列上非常快,但在较大阵列上速度较慢。

    我最终采用了一种使用 ndimage 的方法,该方法似乎对非常大和稀疏的数组有效(5000 x 5000 数组为 0.91 秒,而面积表方法为 1.17 秒)。我首先为每个离散区域生成一个唯一 ID 的标记数组,计算每个 ID 的大小,屏蔽大小数组以仅关注大小 == 1 块,然后索引原始数组并将 ID 设置为大小 == 1 到 0 :

    def filter_isolated_cells(array, struct):
        """ Return array with completely isolated single cells removed
        :param array: Array with completely isolated single cells
        :param struct: Structure array for generating unique regions
        :return: Array with minimum region size > 1
        """
    
        filtered_array = np.copy(array)
        id_regions, num_ids = ndimage.label(filtered_array, structure=struct)
        id_sizes = np.array(ndimage.sum(array, id_regions, range(num_ids + 1)))
        area_mask = (id_sizes == 1)
        filtered_array[area_mask[id_regions]] = 0
        return filtered_array
    
    # Run function on sample array
    filtered_array = filter_isolated_cells(square, struct=np.ones((3,3)))
    
    # Plot output, with all isolated single cells removed
    plt.imshow(filtered_array, cmap=plt.cm.gray, interpolation='nearest')
    

    结果:

    【讨论】:

      【解决方案2】:

      您可以手动检查邻居并使用矢量化避免循环。

      has_neighbor = np.zeros(square.shape, bool)
      has_neighbor[:, 1:] = np.logical_or(has_neighbor[:, 1:], square[:, :-1] > 0)  # left
      has_neighbor[:, :-1] = np.logical_or(has_neighbor[:, :-1], square[:, 1:] > 0)  # right
      has_neighbor[1:, :] = np.logical_or(has_neighbor[1:, :], square[:-1, :] > 0)  # above
      has_neighbor[:-1, :] = np.logical_or(has_neighbor[:-1, :], square[1:, :] > 0)  # below
      
      square[np.logical_not(has_neighbor)] = 0
      

      这样在正方形上循环是由 numpy 在内部执行的,这比在 python 中循环更有效。这种解决方案有两个缺点:

      1. 如果您的数组非常稀疏,则可能有更有效的方法来检查非零点的邻域。
      2. 如果您的数组非常大,has_neighbor 数组可能会消耗太多内存。在这种情况下,您可以遍历较小尺寸的子数组(python 循环和矢量化之间的权衡)。

      我没有使用 ndimage 的经验,所以可能有更好的解决方案。

      【讨论】:

        【解决方案3】:

        在图像处理中消除孤立像素的典型方法是执行morphological opening,您在scipy.ndimage.morphology.binary_opening 中有现成的实现。不过,这也会影响较大区域的轮廓。

        对于 DIY 解决方案,我会使用 summed area table 来计算每个 3x3 子图像中的项目数,从中减去中心像素的值,然后将结果为零的所有中心点归零。为了正确处理边界,首先用零填充数组:

        sat = np.pad(square, pad_width=1, mode='constant', constant_values=0)
        sat = np.cumsum(np.cumsum(sat, axis=0), axis=1)
        sat = np.pad(sat, ((1, 0), (1, 0)), mode='constant', constant_values=0)
        # These are all the possible overlapping 3x3 windows sums
        sum3x3 = sat[3:, 3:] + sat[:-3, :-3] - sat[3:, :-3] - sat[:-3, 3:]
        # This takes away the central pixel value
        sum3x3 -= square
        # This zeros all the isolated pixels
        square[sum3x3 == 0] = 0
        

        上面的实现是可行的,但是没有特别注意不创建中间数组,所以你可以通过充分重构来节省一些执行时间。

        【讨论】:

          猜你喜欢
          • 2012-05-03
          • 2017-06-06
          • 1970-01-01
          • 2017-09-16
          • 2018-09-26
          • 2013-03-05
          • 1970-01-01
          • 2019-05-24
          • 1970-01-01
          相关资源
          最近更新 更多