计算集贸市场的MCP几何

Question

我试图使用几何学 find_costs函数计算市场集。它一直在出色地计算最低成本的路线，但我不想找到最近的来源的旅行成本，我想计算最近的来源的指数。

示例代码

import numpy as np
import skimage.graph as graph
import copy

img = np.array([[1,1,2,2],[2,1,1,3],[3,2,1,2],[2,2,2,1]])
mcp = graph.MCP_Geometric(img)

destinations = [[0,0],[3,3]]
costs, traceback = mcp.find_costs(destinations)
print(costs)

[[0.         1.         2.5        4.5       ]
 [1.5        1.41421356 2.41421356 4.        ]
 [4.         2.91421356 1.41421356 1.5       ]
 [5.5        3.5        1.5        0.        ]]

这如预期的工作，并创造了一个不错的旅行成本光栅。但是，我希望(对于每个单元)知道哪个目的地是最近的。我找到的最佳解决方案是分别运行每个目的地，然后通过最小计算将它们组合起来。它起作用，但速度慢，而且一直没有大规模发挥作用。

all_c = []
for dest in destinations:
    costs, traceback = mcp.find_costs([dest])
    all_c.append(copy.deepcopy(costs))

res = np.dstack(all_c)
res_min = np.amin(res, axis=2)
output = np.zeros([res_min.shape[0], res_min.shape[1]])
for idx in range(0, res.shape[2]):
    cur_data = res[:,:,idx]
    cur_val = (cur_data == res_min).astype(np.byte) * idx
    output = output + cur_val
output = output.astype(np.byte)

print(output)
array([[0, 0, 0, 0],
       [0, 0, 1, 1],
       [0, 1, 1, 1],
       [1, 1, 1, 1]], dtype=int8)

我一直在研究MCP_Geometric和MCP_Flexible函数的重载问题，但是我找不到提供目的地索引的任何信息。

希望这能提供足够的信息来复制和理解我想要做的事情，谢谢！

Answer 1

好吧，这是一段旅程，但弄清楚这很有趣。我不清楚它的速度会有多快，但我认为在很多目的地和舒适的RAM图像中，它应该是相当快的。

关键是traceback返回值，它某种程度上告诉您邻居索引要到达最近的目的地。因此，只要找到一点路径，你就能找到目的地。能快点吗？事实证明，它可以，与一些NumPy索引争论，scipy.sparse矩阵，和connected_components从scipy.sparse.csgraph！

让我们从相同的成本数组开始，这两个目的地：

import numpy as np

image = np.array(
    [[1, 1, 2, 2],
     [2, 1, 1, 3],
     [3, 2, 1, 2],
     [2, 2, 2, 1]]
)
destinations = [[0, 0], [3, 3]]

然后我们绘制图表，得到成本和追溯：

from skimage import graph

mcp = graph.MCP_Geometric(image)
costs, traceback = mcp.find_costs(destinations)
print(traceback)

给予：

[[-1  4  4  4]
 [ 6  7  7  1]
 [ 6  6  0  1]
 [ 3  3  3 -1]]

现在，我必须查找回溯是什么的文档：

与costs数组相同的形状；该数组包含其先前索引中对任何给定索引的偏移量。将偏移索引索引到offsets属性中，这是一个n-d偏移量数组。在2-d情况下，如果偏移量[tracebackx，y]是(-1，-1)，这意味着x，y在到某个起始位置的最小成本路径中的前身是x+1，y+1。注意，如果offset_index是-1，则不考虑给定的索引。

由于某种原因，我的mcp对象没有一个偏移属性--可能是Cython？不知道，稍后再研究--但是搜索源代码会告诉我，偏移是用函数定义的。因此，我做了一件坏事，从那个私有模块导入，这样我就可以声称什么是正确的--偏移列表，它从traceback中的数字转换为图像坐标中的偏移：

from skimage.graph import _mcp

offsets = _mcp.make_offsets(2, True)
print(offsets)

这意味着：

[array([-1, -1]),
 array([-1,  0]),
 array([-1,  1]),
 array([ 0, -1]),
 array([0, 1]),
 array([ 1, -1]),
 array([1, 0]),
 array([1, 1])]

现在，与偏移量有关的最后一件事是:您将注意到，目标在回溯中标记为"-1"，这与偏移量数组的最后一个元素不对应。因此，我们附加np.array([0, 0])，然后traceback中的每个值都对应于一个真正的偏移量。在目的地的情况下，你有一个自我优势，但这是好的。

offsets.append(np.array([0, 0]))
offsets_arr = np.array(offsets)  # shape (9, 2)

现在，我们可以通过偏移量、像素坐标和像素ids构建一个图形。首先，我们使用np.indices为图像中的每个像素获取一个索引：

indices = np.indices(traceback.shape)
print(indices.shape)

给予：

(2, 4, 4)

为了获得一个数组，该数组对于每个像素都具有它的邻居的偏移量，我们使用花式数组索引

offset_to_neighbor = offsets_arr[traceback]
print(offset_to_neighbor.shape)

这意味着：

(4, 4, 2)

追溯索引和numpy索引之间的轴是不同的，但是任何小的转换都不能修复：

neighbor_index = indices - offset_to_neighbor.transpose((2, 0, 1))

最后，我们希望处理整数像素in，以便创建一个所有像素的图形，而不是坐标。为此，我们使用np.ravel_multi_index。

ids = np.arange(traceback.size).reshape(image.shape)
neighbor_ids = np.ravel_multi_index(
    tuple(neighbor_index), traceback.shape
)

这为每个像素提供了唯一的ID，然后为每个像素提供了唯一的“到达目的地的步骤”：

print(ids)
print(neighbor_ids)

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

然后，我们可以使用SciPy稀疏矩阵将其转化为一个图。我们不关心这个图的权重，所以我们只对边使用1的值。

from scipy import sparse

g = sparse.coo_matrix((
    np.ones(traceback.size),
    (ids.flat, neighbor_ids.flat),
    shape=(ids.size, ids.size),
)).tocsr()

(这使用稀疏COOrdinate矩阵的(值、(行、列)或(数据，(i，j))输入格式)。

最后，我们使用连通元件获取图--距离每个目标最近的像素组。该函数返回组件的数量，并将“像素id”映射到组件：

n, components = sparse.csgraph.connected_components(g)
basins = components.reshape(image.shape)
print(basins)

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

(请注意，这个结果与您的结果略有不同，因为有关像素的成本与目标0和1相同，因此标记哪个像素是任意的。)

print(costs)

[[0.         1.         2.5        4.5       ]
 [1.5        1.41421356 2.41421356 4.        ]
 [4.         2.91421356 1.41421356 1.5       ]
 [5.5        3.5        1.5        0.        ]]

希望这能有所帮助！