回答此问题可获得 20 贡献值,回答如果被采纳可获得 50 分。
<p>在python中,<code>scipy.ndimage.morphology</code>模块中有<code>distance_transform_edt</code>函数。我把它应用到一个简单的例子中,计算出一个蒙面核阵列中单个细胞的距离。</p>
<p>但是,该函数删除数组的掩码,并按预期计算每个具有非空值的单元格与具有空值的引用单元格之间的欧几里德距离。</p>
<p>下面是我在<a href="http://bloggb.fr/python/2015/01/26/geodesic-distance-transform-in-python.html" rel="noreferrer">my blog post</a>中给出的一个示例:</p>
<pre><code>%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')
</code></pre>
<p><img src="https://i.stack.imgur.com/NALKQ.png" alt="Euclidean distance transform"/></p>
<p>不过,我想计算测地距离变换,考虑到数组的屏蔽元素。我不想计算穿过屏蔽元素的直线上的欧几里德距离。</p>
<p>我用Dijkstra算法得到了我想要的结果。以下是我提议的实施方案:</p>
<pre><code>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()
</code></pre>
<p><img src="https://i.stack.imgur.com/QjhgG.png" alt="enter image description here"/></p>
<p>上面实现的函数运行良好,但对于我开发的需要多次计算测地距离变换的应用程序来说太慢了。</p>
<p>以下是欧氏距离变换和测地距离变换的时间基准:</p>
<pre><code>%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
</code></pre>
<p>如何获得更快的测地距离变换?</p>