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

Question

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

我的方法如下：找到图像中所有 1 和所有零的协调。找到每个零像素 (a) 和一个像素 (b) 之间的欧几里德距离，然后每个 (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 的约定，而不是更常见的约定。

在我看来，最容易实现的是倒角距离算法。这是一种递归算法，对图像进行两次遍历：一次从左到右、从上到下，一次从右到左、从下到上。在每次传递中，都会传播为先前像素计算的距离。该算法可以使用邻居之间的整数距离或浮点距离来实现。当然，后者产生的误差较小。但在这两种情况下，通过增加传播中查询的邻居数量可以显着减少错误。该算法较旧，但 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 或类似语言编程将产生一种非常快速的算法，该算法产生对欧几里得距离的相当好的近似值。

OpenCVcv.distanceTransform实现了这个算法。

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

该算法最初由 R. van den Boomgaard 在他的博士论文中提出。1992年的论文。不幸的是，这没有引起注意。该算法随后再次由 A. Meijster、JBTM Roerdink 和 WH Hesselink 提出（A General Algorithm forComputing Distance Transforms in Linear Time, Mathematical Morphology and its applications to Image and Signal Processing, pp 331-340, 2002），并再次由 A. Meijster、JBTM Roerdink 和 WH Hesselink 提出P. Felzenszwalb 和 D. Huttenlocher（采样函数的距离变换，技术报告，康奈尔大学，2004 年）。

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

不幸的是，我没有任何 Python 代码可供分享，但您可以在线找到实现。例如，OpenCVcv.distanceTransform实现了该算法，DIPlibdip.EuclideanDistanceTransform也实现了该算法。