使用 Python 在非常大的图像中进行高性能变量模糊

Question

我有大量想要模糊处理的大图像（例如 15000x15000 像素）。我需要使用距离函数模糊图像，因此我离图像中的某些区域越远，模糊应该越重。我有一张距离图，描述了给定像素与这些区域的距离。

由于图片数量众多，我必须考虑性能。我看过 NumPY/SciPY，它们有一些很棒的功能，但它们似乎使用固定的内核大小，我需要根据与前面提到的区域的距离来减小或增加内核大小。

如何在 python 中解决这个问题？

---

更新：到目前为止，我的解决方案基于 rth 的回答：

# cython: boundscheck=False
# cython: cdivision=True
# cython: wraparound=False

import numpy as np
cimport numpy as np

def variable_average(int [:, ::1] data, int[:,::1] kernel_size):
    cdef int width, height, i, j, ii, jj
    width = data.shape[1]
    height = data.shape[0]
    cdef double [:, ::1] data_blurred = np.empty([width, height])
    cdef double res
    cdef int sigma, weight

    for i in range(width):
        for j in range(height):
            weight = 0
            res = 0
            sigma =  kernel_size[i, j]
            for ii in range(i - sigma, i + sigma + 1):
                for jj in range(j - sigma, j + sigma + 1):
                    if ii < 0 or ii >= width or jj < 0 or jj >= height:
                        continue
                    res += data[ii, jj]
                    weight += 1
            data_blurred[i, j] = res/weight

    return data_blurred

测试：

data = np.random.randint(256, size=(1024,1024))
kernel = np.random.randint(256, size=(1024,1024)) + 1
result = np.asarray(variable_average(data, kernel))

使用上述设置的方法运行大约需要 186 秒。这是我可以期望最终从方法中挤出来的，还是我可以使用优化来进一步提高性能（仍然使用 Python）？

Answer 1

正如您所指出的，相关的scipy 函数不支持可变大小模糊。您可以使用 for 循环在纯 python 中实现这一点，然后使用 Cython、Numba 或 PyPy 来获得类似 C 的性能。

这是一个低级的python实现，而不是仅使用numpy进行数据存储，

import numpy as np

def variable_blur(data, kernel_size):
    """ Blur with a variable window size
    Parameters:
      - data: 2D ndarray of floats or integers
      - kernel_size: 2D ndarray of integers, same shape as data
    Returns:
      2D ndarray
    """
    data_blurred = np.empty(data.shape)
    Ni, Nj = data.shape
    for i in range(Ni):
        for j in range(Nj):
            res = 0.0
            weight = 0
            sigma =  kernel_size[i, j]
            for ii in range(i - sigma, i+sigma+1):
                for jj in range(j - sigma, j+sigma+1):
                    if ii<0 or ii>=Ni or jj < 0 or jj >= Nj:
                        continue
                    res += data[ii, jj]
                    weight += 1
            data_blurred[i, j] = res/weight
    return data_blurred

data = np.random.rand(50, 20)
kernel_size = 3*np.ones((50, 20), dtype=np.int)
variable_blur(data, kernel_size)

计算具有可变内核大小的像素的算术平均值。就 numpy 而言，这是一个糟糕的实现，从某种意义上说，它不是矢量化的。但是，这样可以方便地移植到其他高性能解决方案：

• Cython：只需静态键入变量，编译即可获得类似 C 的性能，

  def variable_blur(double [:, ::1] data, long [:, ::1] kernel_size):
       cdef double [:, ::1] data_blurred = np.empty(data.shape)
       cdef Py_ssize_t Ni, Nj
       Ni = data.shape[0]
       Nj = data.shape[1]
       for i in range(Ni):
           # [...] etc.
  

  完整示例请参见 this post，以及 compilation notes。

• Numba：用@jit decorator 包装上述函数应该就足够了。

• PyPy：安装 PyPy + 实验性的 numpy branch，可能是另一个值得尝试的替代方案。不过，您将不得不对所有代码使用 PyPy，目前这可能是不可能的。

一旦你有一个快速的实现，你可以使用multiprocessing等来并行处理不同的图像，如果需要的话。甚至在 Cython 中与 OpenMP 并行化外部 for 循环。

Answer 2

我在谷歌搜索时遇到了这个问题，我想我会分享我自己的解决方案，该解决方案主要是矢量化的，不包括像素上的任何 for 循环。您可以通过连续多次运行框模糊来近似高斯模糊。所以我决定使用的方法是对图像进行迭代框模糊，但使用加权函数来改变每个像素的迭代次数。

如果需要较大的模糊半径，迭代次数会呈二次方增长，因此请考虑增加 ksize。

这里是实现

import cv2

def variable_blur(im, sigma, ksize=3):
    """Blur an image with a variable Gaussian kernel.
    
    Parameters
    ----------
    im: numpy array, (h, w)
    
    sigma: numpy array, (h, w)
    
    ksize: int
        The box blur kernel size. Should be an odd number >= 3.
        
    Returns
    -------
    im_blurred: numpy array, (h, w)
    
    """
    variance = box_blur_variance(ksize)
    # Number of times to blur per-pixel
    num_box_blurs = 2 * sigma**2 / variance
    # Number of rounds of blurring
    max_blurs = int(np.ceil(np.max(num_box_blurs))) * 3
    # Approximate blurring a variable number of times
    blur_weight = num_box_blurs / max_blurs

    current_im = im
    for i in range(max_blurs):
        next_im = cv2.blur(current_im, (ksize, ksize))
        current_im = next_im * blur_weight + current_im * (1 - blur_weight)
    return current_im

def box_blur_variance(ksize):
    x = np.arange(ksize) - ksize // 2
    x, y = np.meshgrid(x, x)
    return np.mean(x**2 + y**2)

这是一个例子

im = np.random.rand(300, 300)
sigma = 3

# Variable
x = np.linspace(0, 1, im.shape[1])
y = np.linspace(0, 1, im.shape[0])
x, y = np.meshgrid(x, y)
sigma_arr = sigma * (x + y)
im_variable = variable_blur(im, sigma_arr)

# Gaussian
ksize = sigma * 8 + 1
im_gauss = cv2.GaussianBlur(im, (ksize, ksize), sigma)

# Gaussian replica
sigma_arr = np.full_like(im, sigma)
im_approx = variable_blur(im, sigma_arr)

Blurring results

剧情是：

• 左上角：源图片
• 右上角：可变模糊
• 左下：高斯模糊
• 右下：近似高斯模糊