按最近的 "seed" 区域对 Python 数组进行分类？

Question

我有一个生态栖息地栅格，我已将其转换为二维 Python numpy 数组（下面的 example_array）。我还有一个数组，其中包含 "seed" 个具有唯一值的区域（下面的 seed_array），我想用它来对我的栖息地区域进行分类。我想 'grow' 我的种子区域 'into' 我的栖息地区域，以便为栖息地分配最近的种子区域的 ID，如测量 'through' 栖息地区域。 例如：

我最好的方法是使用 ndimage.distance_transform_edt 函数来创建一个数组，该数组描述距离数据集中每个单元格最近的 "seed" 区域，然后将其替换回栖息地数组。然而，这并不是特别有效，因为该函数不测量距离 "through" 我的栖息地区域，例如下面的红色圆圈代表错误分类的单元格：

下面是我的栖息地和种子数据的示例数组，以及我正在寻找的输出类型的示例。我的实际数据集要大得多——超过一百万 habitat/seed 个区域。任何帮助将不胜感激！

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

# Sample study area array
example_array = np.array([[0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0],
                          [0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 1],
                          [0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1],
                          [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0],
                          [1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1],
                          [1, 1, 0, 1, 0, 0, 0, 0, 1, 1, 1, 1],
                          [1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1],
                          [1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0],
                          [1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 0],
                          [1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0],
                          [1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0],
                          [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]])

# Plot example array
plt.imshow(example_array, cmap="spectral", interpolation='nearest')

seed_array = np.array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                       [0, 0, 0, 0, 4, 0, 0, 0, 0, 0, 0, 0],
                       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0],
                       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                       [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                       [0, 1, 1, 1, 0, 0, 2, 2, 0, 0, 0, 0],
                       [0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0],
                       [0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0],
                       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]])

# Plot seeds
plt.imshow(seed_array, cmap="spectral", interpolation='nearest')

desired_output = np.array([[0, 0, 0, 0, 4, 0, 0, 0, 0, 0, 0, 0],
                           [0, 0, 0, 4, 4, 4, 0, 0, 0, 3, 3, 3],
                           [0, 0, 0, 0, 4, 4, 0, 0, 0, 3, 3, 3],
                           [0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 3, 0],
                           [1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 3, 3],
                           [1, 1, 0, 1, 0, 0, 0, 0, 2, 2, 3, 3],
                           [1, 1, 1, 1, 0, 0, 2, 2, 2, 0, 0, 3],
                           [1, 1, 1, 1, 1, 2, 2, 2, 2, 0, 0, 0],
                           [1, 1, 1, 1, 0, 0, 2, 2, 2, 0, 0, 0],
                           [1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0],
                           [1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0],
                           [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]])

# Plot desired output
plt.imshow(desired_output, cmap="spectral", interpolation='nearest')

Answer 1

您可以使用来自 scikits-image 的 watershed segmentation：

距离变换

from scipy import ndimage as nd
distance = nd.distance_transform_edt(example_array)

分水岭分割

from skimage.morphology import watershed, square
result = watershed(-distance, seed_array, mask=example_array, \
                   connectivity=square(3))

结果

subplot(1,2,1)
imshow(-distance, 'spectral', interpolation='none')
subplot(1,2,2)
imshow(result, 'spectral', interpolation='none')

作为另一种变体，按照您最初的方法，您可以使用分水岭来找到与最近种子相连的邻居。正如您在问题中提到的：

计算到种子的距离：

distance = nd.distance_transform_edt(seed_array == 0)

计算距离的分水岭space:

result = watershed(distance, seed_array, mask=example_array, \
                   connectivity=square(3))

绘图结果：

figure(figsize=(9,3))
subplot(1,3,1)
imshow(distance, 'jet', interpolation='none')
subplot(1,3,2)
imshow(np.ma.masked_where(example_array==0, distance), 'jet', interpolation='none')
subplot(1,3,3)
imshow(result, 'spectral', interpolation='none')

进一步讨论：分水岭方法尝试通过图像梯度流动从种子峰增长区域。由于您的图像是二进制的，因此区域将从种子点向各个方向均匀扩展，从而为您提供两个区域之间的点。有关分水岭的更多信息，请参阅 wikipedia.

在第一个例子中，距离变换是在原始图像中计算的，因此区域从种子开始平均扩展，直到达到中间的分裂点。

在第二个示例中，计算从所有像素到任何种子点的距离变换，然后在该 space 中应用分水岭。分水岭基本上会将每个像素分配给它最近的种子，但它会添加连接约束。

注意绘图和分水岭距离图中的符号差异。

注意在距离图中（两个图中的左图），蓝色表示近，红色表示远。

按最近的 "seed" 区域对 Python 数组进行分类？

Classifying Python array by nearest "seed" region?

python

arrays

numpy

scipy

image-segmentation