加快 Python 中图像距离变换的计算

Question

我想在不使用 scipy 包 distance_trnsform_edt() 的情况下以最快的方式找到二进制图像的距离变换。图像是 256 x 256。我不想使用 scipy 的原因是因为在 tensorflow 中使用它很困难。每次我想使用这个包时，我都需要开始一个新的会话，这需要很多时间。所以我想做一个只使用numpy的自定义函数。

我的方法如下：找到图像中所有 1 和所有 0 的坐标。找到每个零像素 (a) 和一个像素 (b) 之间的欧几里德距离，然后每个 (a) 位置的值是到 a (b) 像素的最小距离。我为每个 0 像素执行此操作。生成的图像具有与原始二值图相同的尺寸。我的尝试如下所示。

我尝试尽可能快地做到这一点，不使用循环，只使用矢量化。但是我的功能仍然不能像 scipy 包那样快。当我对代码进行计时时，看起来分配给变量“a”的时间最长。但是不知道有没有办法加快速度。

如果有人对解决这个距离变换问题的不同算法有任何其他建议，或者可以指导我使用 python 中的其他实现，将不胜感激。

def get_dst_transform_img(og): #og is a numpy array of original image
    ones_loc = np.where(og == 1)
    ones = np.asarray(ones_loc).T # coords of all ones in og
    zeros_loc = np.where(og == 0)
    zeros = np.asarray(zeros_loc).T # coords of all zeros in og

    a = -2 * np.dot(zeros, ones.T) 
    b = np.sum(np.square(ones), axis=1) 
    c = np.sum(np.square(zeros), axis=1)[:,np.newaxis]
    dists = a + b + c
    dists = np.sqrt(dists.min(axis=1)) # min dist of each zero pixel to one pixel
    x = og.shape[0]
    y = og.shape[1]
    dist_transform = np.zeros((x,y))
    dist_transform[zeros[:,0], zeros[:,1]] = dists 

    plt.figure()
    plt.imshow(dist_transform)

Answer 1

OP 中的实现是距离变换的蛮力方法。该算法是 O(n²)，因为它计算每个背景像素到每个前景像素的距离。此外，由于它是矢量化的方式，它需要大量内存。在我的计算机上，它无法计算 256x256 图像的距离变换而不会颠簸。文献中描述了许多其他算法，下面我将讨论两种 O(n) 算法。

注意：通常，距离变换是针对对象像素（值 1）计算到最近的背景像素（值 0）。 OP 中的代码正好相反，所以我在下面粘贴的代码遵循 OP 的约定，而不是更常见的约定。

---

IMO 最容易实现的是倒角距离算法。这是一种递归算法，对图像进行两次遍历：一次从左到右和从上到下，一次从右到左和从下到上。在每次传递中，传播为先前像素计算的距离。该算法可以使用邻居之间的整数距离或浮点距离来实现。当然，后者会产生较小的错误。但是在这两种情况下，可以通过增加在此传播中查询的邻居数量来显着减少错误。该算法较旧，但 G. Borgefors 对其进行了分析并提出了合适的邻居距离 (G. Borgefors, Distance Transformations in Digital Images, Computer Vision, Graphics, and Image Processing 34:344-371, 1986)。

这是一个使用 3-4 距离的实现（到边缘连接的邻居的距离是 3，到顶点连接的邻居的距离是 4）：

def chamfer_distance(img):
   w, h = img.shape
   dt = np.zeros((w,h), np.uint32)
   # Forward pass
   x = 0
   y = 0
   if img[x,y] == 0:
      dt[x,y] = 65535 # some large value
   for x in range(1, w):
      if img[x,y] == 0:
         dt[x,y] = 3 + dt[x-1,y]
   for y in range(1, h):
      x = 0
      if img[x,y] == 0:
         dt[x,y] = min(3 + dt[x,y-1], 4 + dt[x+1,y-1])
      for x in range(1, w-1):
         if img[x,y] == 0:
            dt[x,y] = min(4 + dt[x-1,y-1], 3 + dt[x,y-1], 4 + dt[x+1,y-1], 3 + dt[x-1,y])
      x = w-1
      if img[x,y] == 0:
         dt[x,y] = min(4 + dt[x-1,y-1], 3 + dt[x,y-1], 3 + dt[x-1,y])
   # Backward pass
   for x in range(w-2, -1, -1):
      y = h-1
      if img[x,y] == 0:
         dt[x,y] = min(dt[x,y], 3 + dt[x+1,y])
   for y in range(h-2, -1, -1):
      x = w-1
      if img[x,y] == 0:
         dt[x,y] = min(dt[x,y], 3 + dt[x,y+1], 4 + dt[x-1,y+1])
      for x in range(1, w-1):
         if img[x,y] == 0:
            dt[x,y] = min(dt[x,y], 4 + dt[x+1,y+1], 3 + dt[x,y+1], 4 + dt[x-1,y+1], 3 + dt[x+1,y])
      x = 0
      if img[x,y] == 0:
         dt[x,y] = min(dt[x,y], 4 + dt[x+1,y+1], 3 + dt[x,y+1], 3 + dt[x+1,y])
   return dt

请注意，这里的很多复杂情况是为了避免索引超出范围，但仍要计算一直到图像边缘的距离。如果我们简单地跳过图像边框周围的像素，代码就会变得简单得多。

因为它是一种递归算法，所以不可能对其实现进行向量化。 Python 代码效率不会很高。但是用 C 或类似语言编程会产生一个非常快速的算法，它可以很好地近似欧几里得距离。

OpenCV's cv.distanceTransform 实现了这个算法。

---

另一种非常有效的算法计算距离变换的平方。平方距离是可分离的（即可以为每个轴独立计算并添加）。这导致了一种易于并行化的算法。对于每个图像行，该算法执行前向和后向传递。对于结果中的每一列，该算法然后进行另一次向前和向后传递。这个过程会导致精确的欧几里得距离变换。

该算法由 R. van den Boomgaard 在his Ph.D. thesis in 1992 中首次提出。不幸的是，这没有引起注意。该算法随后再次由 A. Meijster、J.B.T.M. 提出。鲁尔丁克和 W.H. Hesselink (A General Algorithm for Computing Distance Transforms in Linear Time, Mathematical Morphology and its Applications to Image and Signal Processing, pp 331-340, 2002)，再次由 P. Felzenszwalb 和 D. Huttenlocher (Distance transforms of sampled functions, Technical report, Cornell University, 2004) 撰写。

这是已知的最有效的算法，部分原因是它是唯一一种可以轻松高效地并行化的算法（每个图像行的计算，然后是每个图像列的计算，独立于其他行/列）。

很遗憾，我没有任何 Python 代码可供分享，但您可以在网上找到实现。例如OpenCV's cv.distanceTransform 实现了这个算法，DIPlib's dip.EuclideanDistanceTransform 也实现了。