python绘制八角形_python中的测地距离变换

在python中， scipy.ndimage.morphology 模块中有 distance_transform_edt 函数 . 我将它应用于一个简单的情况，计算与蒙版numpy数组中单个单元格的距离 .

但是，该函数移除数组的掩码，并按预期计算每个单元格的欧几里德距离，其中非空值来自参考单元格，具有空值 .

以下是我在my blog post中给出的一个例子：

%pylab

from scipy.ndimage.morphology import distance_transform_edt

l = 100

x, y = np.indices((l, l))

center1 = (50, 20)

center2 = (28, 24)

center3 = (30, 50)

center4 = (60,48)

radius1, radius2, radius3, radius4 = 15, 12, 19, 12

circle1 = (x - center1[0])**2 + (y - center1[1])**2 < radius1**2

circle2 = (x - center2[0])**2 + (y - center2[1])**2 < radius2**2

circle3 = (x - center3[0])**2 + (y - center3[1])**2 < radius3**2

circle4 = (x - center4[0])**2 + (y - center4[1])**2 < radius4**2

# 3 circles

img = circle1 + circle2 + circle3 + circle4

mask = ~img.astype(bool)

img = img.astype(float)

m = ones_like(img)

m[center1] = 0

#imshow(distance_transform_edt(m), interpolation='nearest')

m = ma.masked_array(distance_transform_edt(m), mask)

imshow(m, interpolation='nearest')

但是我想计算考虑了数组掩码元素的测地距离变换 . 我不想沿着穿过蒙面元素的直线计算欧氏距离 .

我使用Dijkstra算法来获得我想要的结果 . 以下是我提出的实施方案：

def geodesic_distance_transform(m):

mask = m.mask

visit_mask = mask.copy() # mask visited cells

m = m.filled(numpy.inf)

m[m!=0] = numpy.inf

distance_increments = numpy.asarray([sqrt(2), 1., sqrt(2), 1., 1., sqrt(2), 1., sqrt(2)])

connectivity = [(i,j) for i in [-1, 0, 1] for j in [-1, 0, 1] if (not (i == j == 0))]

cc = unravel_index(m.argmin(), m.shape) # current_cell

while (~visit_mask).sum() > 0:

neighbors = [tuple(e) for e in asarray(cc) - connectivity

if not visit_mask[tuple(e)]]

tentative_distance = [distance_increments[i] for i,e in enumerate(asarray(cc) - connectivity)

if not visit_mask[tuple(e)]]

for i,e in enumerate(neighbors):

d = tentative_distance[i] + m[cc]

if d < m[e]:

m[e] = d

visit_mask[cc] = True

m_mask = ma.masked_array(m, visit_mask)

cc = unravel_index(m_mask.argmin(), m.shape)

return m

gdt = geodesic_distance_transform(m)

imshow(gdt, interpolation='nearest')

colorbar()

上面实现的功能运行良好，但对于我开发的需要多次计算测地距离变换的应用来说太慢了 .

以下是欧氏距离变换和测地距离变换的时间基准：

%timeit distance_transform_edt(m)

1000 loops, best of 3: 1.07 ms per loop

%timeit geodesic_distance_transform(m)

1 loops, best of 3: 702 ms per loop

如何获得更快的测地距离变换？