<p>稍微快一点(大约10倍)的实现,可以获得与您的<code>geodesic_distance_transform</code>相同的结果:</p>
<pre><code>def getMissingMask(slab):
nan_mask=numpy.where(numpy.isnan(slab),1,0)
if not hasattr(slab,'mask'):
mask_mask=numpy.zeros(slab.shape)
else:
if slab.mask.size==1 and slab.mask==False:
mask_mask=numpy.zeros(slab.shape)
else:
mask_mask=numpy.where(slab.mask,1,0)
mask=numpy.where(mask_mask+nan_mask>0,1,0)
return mask
def geodesic(img,seed):
seedy,seedx=seed
mask=getMissingMask(img)
#----Call distance_transform_edt if no missing----
if mask.sum()==0:
slab=numpy.ones(img.shape)
slab[seedy,seedx]=0
return distance_transform_edt(slab)
target=(1-mask).sum()
dist=numpy.ones(img.shape)*numpy.inf
dist[seedy,seedx]=0
def expandDir(img,direction):
if direction=='n':
l1=img[0,:]
img=numpy.roll(img,1,axis=0)
img[0,:]==l1
elif direction=='s':
l1=img[-1,:]
img=numpy.roll(img,-1,axis=0)
img[-1,:]==l1
elif direction=='e':
l1=img[:,0]
img=numpy.roll(img,1,axis=1)
img[:,0]=l1
elif direction=='w':
l1=img[:,-1]
img=numpy.roll(img,-1,axis=1)
img[:,-1]==l1
elif direction=='ne':
img=expandDir(img,'n')
img=expandDir(img,'e')
elif direction=='nw':
img=expandDir(img,'n')
img=expandDir(img,'w')
elif direction=='sw':
img=expandDir(img,'s')
img=expandDir(img,'w')
elif direction=='se':
img=expandDir(img,'s')
img=expandDir(img,'e')
return img
def expandIter(img):
sqrt2=numpy.sqrt(2)
tmps=[]
for dirii,dd in zip(['n','s','e','w','ne','nw','sw','se'],\
[1,]*4+[sqrt2,]*4):
tmpii=expandDir(img,dirii)+dd
tmpii=numpy.minimum(tmpii,img)
tmps.append(tmpii)
img=reduce(lambda x,y:numpy.minimum(x,y),tmps)
return img
#----------------Iteratively expand----------------
dist_old=dist
while True:
expand=expandIter(dist)
dist=numpy.where(mask,dist,expand)
nc=dist.size-len(numpy.where(dist==numpy.inf)[0])
if nc>=target or numpy.all(dist_old==dist):
break
dist_old=dist
return dist
</code></pre>
<p>还要注意,如果蒙版形成超过1个连接区域(例如,添加另一个不接触其他区域的圆),则函数将陷入一个无休止的循环。</p>
<p><strong>更新</strong>:</p>
<p>我发现了一个快速扫描方法<a href="https://github.com/znah/notebooks/blob/master/geodesic.ipynb" rel="nofollow noreferrer">in this notebook</a>的Cython实现,它可以以可能相当的速度实现与<code>scikit-fmm</code>相同的结果。我们只需要输入一个二进制标记矩阵(1作为可行点,否则为<code>inf</code>)作为<code>GDT()</code>函数的代价。</p>