使用 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.

参见 for a complete example, as well as the compilation notes。

Numba：用 @jit decorator 包装上面的函数应该足够了。
PyPy：安装 PyPy + 实验性 numpy branch，可能是另一种值得尝试的选择。尽管如此，您将不得不为所有代码使用 PyPy，这在目前可能是不可能的。

快速实现后，如果需要，您可以使用 multiprocessing 等并行处理不同的图像。或者甚至在外部 for 循环中与 Cython 中的 OpenMP 并行化。

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

剧情是：

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

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

High performance variable blurring in very big images using Python

python

numpy

image-processing

cython

scipy