Python迭代灰度图像中的连通组件

Question

我有一张灰度图像，像素值在0（黑色）和255（白色）之间。我有一个与灰度图像相同大小的目标矩阵。我需要从灰度图像中的随机像素开始，以深度优先搜索的方式逐个遍历图像中的每个像素，并将其值复制到目标矩阵中的相应位置。显然，我只需对非白色像素执行此操作。我想过可以获取灰度图像的连通分量，并逐个像素遍历每个分量，但我找不到任何合适的连通分量实现。有什么想法吗？

例如，如果我的灰度图像是：

[[255,255,255,255,255,255,255]
[255,255, 0 ,10 ,255,255, 1 ]
[255,30 ,255,255,50 ,255, 9 ]
[51 ,20 ,255,255, 9 ,255,240]
[255,255,80 ,50 ,170,255, 20]
[255,255,255,255,255,255, 0 ]
[255,255,255,255,255,255, 69]]

一种可能的遍历顺序是[0,10,50,9,170,50,80,20,51,30]，然后跟随[1,9,240,20,0,69]，得到[0,10,50,9,170,50,80,20,51,30,1,9,240,20,0,69]。不同对象之间的顺序不重要。

其他可能的遍历顺序包括：[1,9,240,20,0,69,0,10,50,9,170,50,80,20,51,30]或[1,9,240,20,0,69,0,10,50,9,170,50,80,20,30,51]或[1,9,240,20,0,69,10,50,9,170,50,80,20,30,0,51]

等等。

- RaviTej310

---

你能加入一些示例数据吗？我已经读了这个问题5次，但是我没有看出使用两个嵌套的for循环，以及通过复制除值!=255之外的值来排除的问题。或许需要更明确地说明问题。 - benschbob91

添加了一个示例。我不能只复制非255值，因为那不是深度优先搜索。 - RaviTej310

我理解“深度优先”是指在关注每个目录的内容之前，首先沿着层次结构向下遍历目录。那么，在一个扁平的、非分层的图像中，“深度优先”是什么意思呢？ - Mark Setchell

@KartikeySingh修复了这个例子！ - RaviTej310

@MarkSetchell 是的。在这种情况下，我所说的深度优先是指：从一个像素开始，以深度优先的方式遍历相邻的像素（即继续只去任何连接的像素并继续前进，直到一个像素没有更多未访问的连接像素。然后返回到之前的像素并查看它们是否有任何未访问的邻居，并继续这个过程，直到所有非白色像素都被访问）。 - RaviTej310

Answer 1

你可以使用 networkx 库:

from itertools import product, repeat
import numpy as np
import networkx as nx

arr = np.array(
[[255,255,255,255,255,255,255],
 [255,255, 0 ,10 ,255,255, 1 ],
 [255,30 ,255,255,50 ,255, 9 ],
 [51 ,20 ,255,255, 9 ,255,240],
 [255,255,80 ,50 ,170,255, 20],
 [255,255,255,255,255,255, 0 ],
 [255,255,255,255,255,255, 69]])

# generate edges
shift = list(product(*repeat([-1, 0, 1], 2)))
x_max, y_max = arr.shape
edges = []

for x, y in np.ndindex(arr.shape):
    for x_delta, y_delta in shift:
        x_neighb = x + x_delta
        y_neighb = y + y_delta
        if (0 <= x_neighb < x_max) and (0 <= y_neighb < y_max):
            edge = (x, y), (x_neighb, y_neighb)
            edges.append(edge)

# build graph
G = nx.from_edgelist(edges)

# draw graph
pos = {(x, y): (y, x_max-x) for x, y in G.nodes()}
nx.draw(G, with_labels=True, pos=pos, node_color='coral', node_size=1000)

# draw graph with numbers
labels = dict(np.ndenumerate(arr))
node_color = ['coral' if labels[n] == 255 else 'lightgrey' for n in G.nodes()]
nx.draw(G, with_labels=True, pos=pos, labels=labels, node_color=node_color, node_size=1000)

# build subgraph
select = np.argwhere(arr < 255)
G1 = G.subgraph(map(tuple, select))

# draw subgraph
pos = {(x, y): (y, x_max-x) for x, y in G1.nodes()}
labels1 = {n:labels[n] for n in G1.nodes()}
nx.draw(G1, with_labels=True, pos=pos, labels=labels1, node_color='lightgrey', node_size=1000)

# find connected components and DFS trees
for i in nx.connected_components(G1):
    source = next(iter(i))
    idx = nx.dfs_tree(G1, source=source)
    print(arr[tuple(np.array(idx).T)])

输出：

[  0  10  50   9  50  80  20  30  51 170]
[  9   1 240  20   0  69]

- Mykola Zotko

---

看起来非常优雅！ - Mark Setchell

Answer 2

经过大量研究寻找适合的连通组件实现后，我想出了自己的解决方案。为了在性能方面达到最佳效果，我遵循以下规则：

1. 不使用networkx，因为根据基准测试结果，它速度较慢。
2. 尽可能多地使用向量化操作，因为基于Python的迭代速度较慢，根据这个答案。

我正在实现图像的连通组件算法，因为我认为这是这个问题的一个基本部分。

图像连通组件算法

import numpy as np
import numexpr as ne
import pandas as pd
import igraph

def get_coords(arr):
    x, y = np.indices(arr.shape)
    mask = arr != 255
    return  np.array([x[mask], y[mask]]).T

def compare(r1, r2):
    #assuming r1 is a sorted array, returns:
    # 1) locations of r2 items in r1
    # 2) mask array of these locations
    idx = np.searchsorted(r1, r2)
    idx[idx == len(r1)] = 0
    mask = r1[idx] == r2
    return idx, mask

def get_reduction(coords, s):
    d = {'s': s, 'c0': coords[:,0], 'c1': coords[:,1]}
    return ne.evaluate('c0*s+c1', d)

def get_bounds(coords, increment):
    return np.max(coords[1]) + 1 + increment

def get_shift_intersections(coords, shifts):
    # instance that consists of neighbours found for each node [[0,1,2],...]
    s = get_bounds(coords, 10)
    rdim = get_reduction(coords, s)
    shift_mask, shift_idx = [], []
    for sh in shifts:
        sh_rdim = get_reduction(coords + sh, s)
        sh_idx, sh_mask = compare(rdim, sh_rdim)
        shift_idx.append(sh_idx)
        shift_mask.append(sh_mask)
    return np.array(shift_idx).T, np.array(shift_mask).T,

def connected_components(coords, shifts):
    shift_idx, shift_mask = get_shift_intersections(coords, shifts)
    x, y = np.indices((len(shift_idx), len(shift_idx[0])))
    vertices = np.arange(len(coords))
    edges = np.array([x[shift_mask], shift_idx[shift_mask]]).T

    graph = igraph.Graph()
    graph.add_vertices(vertices)
    graph.add_edges(edges)
    graph_tags = graph.clusters().membership
    values = pd.DataFrame(graph_tags).groupby([0]).indices
    return values

coords = get_coords(arr)
shifts=((0,1),(1,0),(1,1),(-1,1))
comps = connected_components(coords, shifts=shifts)

for c in comps:
    print(coords[comps[c]].tolist()) 

结果

[[1, 2], [1, 3], [2, 1], [2, 4], [3, 0], [3, 1], [3, 4], [4, 2], [4, 3], [4, 4]]
[[1, 6], [2, 6], [3, 6], [4, 6], [5, 6], [6, 6]]

说明

算法由以下步骤组成：

• We need to convert image to coordinates of non-white cells. It can be done using function:

  def get_coords(arr):
      x, y = np.indices(arr.shape)
      mask = arr != 255
      return np.array([y[mask], x[mask]]).T
  

  I'll name an outputting array by X for clarity. Here is an output of this array, visually:

  

• Next, we need to consider all the cells of each shift that intersects with X:

  

  In order to do that, we should solve a problem of intersections I posted few days before. I found it quite difficult to do using multidimensional numpy arrays. Thanks to Divakar, he proposes a nice way of dimensionality reduction using numexpr package which fastens operations even more than numpy. I implement it here in this function:

  def get_reduction(coords, s):
      d = {'s': s, 'c0': coords[:,0], 'c1': coords[:,1]}
      return ne.evaluate('c0*s+c1', d)
  

  In order to make it working, we should set a bound s which can be calculated automatically using a function

  def get_bounds(coords, increment):
      return np.max(coords[1]) + 1 + increment
  

  or inputted manually. Since algorithm requires increasing coordinates, pairs of coordinates might be out of bounds, therefore I have used a slight increment here. Finally, as a solution to my post I mentioned here, indexes of coordinates of X (reduced to 1D), that intersects with any other array of coordinates Y (also reduced to 1D) can be accessed via function

  def compare(r1, r2):
      # assuming r1 is a sorted array, returns:
      # 1) locations of r2 items in r1
      # 2) mask array of these locations
      idx = np.searchsorted(r1, r2)
      idx[idx == len(r1)] = 0
      mask = r1[idx] == r2
      return idx, mask
  

• Plugging all the corresponding arrays of shifts. As we can see, abovementioned function outputs two variables: an array of index locations in the main set X and its mask array. A proper indexes can be found using idx[mask] and since this procedure is being applied for each shift, I implemented get_shift_intersections(coords, shifts) method for this case.

• Final: constructing nodes & edges and taking output from igraph. The point here is that igraph performs well only with nodes that are consecutive integers starting from 0. That's why my script was designed to use mask-based access to item locations in X. I'll explain briefly how did I use igraph here:

  • I have calculated coordinate pairs:

      [[1, 2], [1, 3], [1, 6], [2, 1], [2, 4], [2, 6], [3, 0], [3, 1], [3, 4], [3, 6], [4, 2], [4, 3], [4, 4], [4, 6], [5, 6], [6, 6]]
    

  • Then I assigned integers for them:

      [ 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15]
    

  • My edges looks like this:

      [[0, 1], [1, 4], [2, 5], [3, 7], [3, 0], [4, 8], [5, 9], [6, 7], [6, 3], [7, 10], [8, 12], [9, 13], [10, 11], [11, 12], [11, 8], [13, 14], [14, 15]]
    

  • Output of graph.clusters().membership looks like this:

      [0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1]
    

  • And finally, I have used groupby method of Pandas to find indexes of separate groups (I use Pandas here because I found it to be the most efficient way of grouping in Python)

注意事项

下载 igraph 不是一件简单的事情，您可能需要从非官方二进制文件中安装。

安装指南请参考 这里。

- mathfux

---

这么多工作，OP甚至还没有点赞。所以你使用交集来获取那些子图了？ - Mykola Zotko

是的，我已经遍历了每个非白色单元格的移位，并使用初始单元格与移位单元格的交集作为掩码操作。生成的数组长度相等，可以保存为 numpy 矩阵。 - mathfux