用于计算市场的 MCP 几何学
MCP Geometrics for calculating marketsheds
我正在尝试使用 skimage.MCP_geometric 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 的函数，但我找不到任何提供目标索引信息的东西。
希望提供足够的信息来复制和理解我想做什么，谢谢！
好吧，这有点费劲，但弄明白很有趣。我不清楚它的速度有多快，但我认为在许多目的地和舒适的 RAM 图像的情况下它应该非常快。
关键是 traceback return 值，有点像告诉您到达最近目的地的邻居索引。因此，通过一些寻路，您应该能够找到该目的地。可以这么快吗？事实证明它可以，通过一些 NumPy 索引争论，scipy.sparse matrices, and connected_components from 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]]
现在，我不得不查找文档以了解回溯是什么：
Same shape as the costs array; this array contains the offset to
any given index from its predecessor index. The offset indices
index into the offsets attribute, which is a array of n-d
offsets. In the 2-d case, if offsets[traceback[x, y]] is (-1, -1),
that means that the predecessor of [x, y] in the minimum cost path
to some start position is [x+1, y+1]. Note that if the
offset_index is -1, then the given index was not considered.
出于某种原因，我的 mcp 对象没有偏移量属性 — 可能是 Cython 继承错误？不知道，稍后会调查 - 但搜索 source code shows me that offsets is defined with the skimage.graph._mcp.make_offsets function。所以我做了一件坏事并从那个私有模块导入，所以我可以声明我的权利 - 偏移量列表，它从 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)
现在，我们可以根据偏移量、像素坐标和像素 ID 构建图形。首先，我们使用 np.indices 为图像中的每个像素获取索引：
indices = np.indices(traceback.shape)
print(indices.shape)
给出：
(2, 4, 4)
要获得一个数组，每个像素都有与其相邻像素的偏移量，我们使用 fancy array indexing:
offset_to_neighbor = offsets_arr[traceback]
print(offset_to_neighbor.shape)
给出：
(4, 4, 2)
traceback 和 numpy 索引之间的轴不同，但一点点换位都解决不了：
neighbor_index = indices - offset_to_neighbor.transpose((2, 0, 1))
最后，我们要处理整数像素 ID，以便创建所有像素而不是坐标的图形。为此，我们使用 np.ravel_multi_index.
ids = np.arange(traceback.size).reshape(image.shape)
neighbor_ids = np.ravel_multi_index(
tuple(neighbor_index), traceback.shape
)
这为每个像素提供了一个唯一的 ID，然后为每个像素提供了一个唯一的 "step towards the destination"：
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()
（这使用 sparse COOrdinate matrices 的（值，（行，列））或（数据，（i，j））输入格式。）
最后，我们使用 connected components 来获取图表 — 离每个目的地最近的像素组。函数 return 的组件数和 "pixel 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. ]]
希望对您有所帮助！
我正在尝试使用 skimage.MCP_geometric 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 的函数，但我找不到任何提供目标索引信息的东西。
希望提供足够的信息来复制和理解我想做什么，谢谢！
好吧，这有点费劲，但弄明白很有趣。我不清楚它的速度有多快，但我认为在许多目的地和舒适的 RAM 图像的情况下它应该非常快。
关键是 traceback return 值，有点像告诉您到达最近目的地的邻居索引。因此，通过一些寻路，您应该能够找到该目的地。可以这么快吗？事实证明它可以，通过一些 NumPy 索引争论，scipy.sparse matrices, and connected_components from 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]]
现在，我不得不查找文档以了解回溯是什么：
Same shape as the
costsarray; this array contains the offset to any given index from its predecessor index. The offset indices index into theoffsetsattribute, which is a array of n-d offsets. In the 2-d case, if offsets[traceback[x, y]] is (-1, -1), that means that the predecessor of [x, y] in the minimum cost path to some start position is [x+1, y+1]. Note that if the offset_index is -1, then the given index was not considered.
出于某种原因，我的 mcp 对象没有偏移量属性 — 可能是 Cython 继承错误？不知道，稍后会调查 - 但搜索 source code shows me that offsets is defined with the skimage.graph._mcp.make_offsets function。所以我做了一件坏事并从那个私有模块导入，所以我可以声明我的权利 - 偏移量列表，它从 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)
现在，我们可以根据偏移量、像素坐标和像素 ID 构建图形。首先，我们使用 np.indices 为图像中的每个像素获取索引：
indices = np.indices(traceback.shape)
print(indices.shape)
给出：
(2, 4, 4)
要获得一个数组，每个像素都有与其相邻像素的偏移量，我们使用 fancy array indexing:
offset_to_neighbor = offsets_arr[traceback]
print(offset_to_neighbor.shape)
给出：
(4, 4, 2)
traceback 和 numpy 索引之间的轴不同，但一点点换位都解决不了：
neighbor_index = indices - offset_to_neighbor.transpose((2, 0, 1))
最后，我们要处理整数像素 ID，以便创建所有像素而不是坐标的图形。为此，我们使用 np.ravel_multi_index.
ids = np.arange(traceback.size).reshape(image.shape)
neighbor_ids = np.ravel_multi_index(
tuple(neighbor_index), traceback.shape
)
这为每个像素提供了一个唯一的 ID，然后为每个像素提供了一个唯一的 "step towards the destination"：
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()
（这使用 sparse COOrdinate matrices 的（值，（行，列））或（数据，（i，j））输入格式。）
最后，我们使用 connected components 来获取图表 — 离每个目的地最近的像素组。函数 return 的组件数和 "pixel 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. ]]
希望对您有所帮助！