python - 在 numpy/scipy 中优化旋转掩码的实现

Question

这是我第一次尝试在 numpy 中使用 strides，与不同过滤器上的简单迭代相比，它确实提高了速度，但它仍然很慢(感觉至少有一两件事是完全多余的或效率低下)。

所以我的问题是:是否有更好的方法来执行此操作或调整我的代码以使其速度显着加快？

该算法对每个像素执行 9 种不同过滤器的局部评估，并选择标准偏差最小的过滤器(我尝试实现 Nagau 和 Matsuyma(1980 年)“复杂区域照片的结构分析”，如图像分析书)。结果是平滑和边缘锐化的图像(如果你问我的话，非常酷!)

import numpy as np
from scipy import ndimage
from numpy.lib import stride_tricks

def get_rotating_kernels():

    kernels = list()

    protokernel = np.arange(9).reshape(3,  3)

    for k in xrange(9):

        ax1, ax2 = np.where(protokernel==k)
        kernel = np.zeros((5,5), dtype=bool)
        kernel[ax1: ax1+3, ax2: ax2+3] = 1
        kernels.append(kernel)

    return kernels


def get_rotation_smooth(im, **kwargs):

    kernels = np.array([k.ravel() for k in get_rotating_kernels()],
                dtype=bool)

    def rotation_matrix(section):

        multi_s = stride_tricks.as_strided(section, shape=(9,25),
            strides=(0, section.itemsize))

        rot_filters = multi_s[kernels].reshape(9,9)

        return rot_filters[rot_filters.std(1).argmin(),:].mean()

    return ndimage.filters.generic_filter(im, rotation_matrix, size=5, **kwargs)

from scipy import lena
im = lena()
im2 = get_rotation_smooth(im)

(只是评论，get_rotating_kernel 并没有真正优化，因为几乎没有时间花在那里)

在我的上网本上，花了 126 秒，Lena 毕竟是一个很小的图像。

编辑:

我得到的建议是将 rot_filters.std(1) 更改为 rot_filters.var(1) 以保存相当多的平方根，并且它削减了一些东西5s 的顺序。

Answer 1

我相信您将很难使用 Python + scipy 进行显着优化。但是，我能够通过使用 as_strided 直接生成 rot_filters 而不是通过 bool 索引进行小的改进。这是基于一个非常简单的 n 维 windows 函数。 (在我意识到 scipy 中存在二维卷积函数之前，我写它是为了解决 this problem。)下面的代码在我的机器上提供了 10% 的适度加速；有关其工作原理的说明，请参见下文:

import numpy as np
from scipy import ndimage
from numpy.lib import stride_tricks

# pass in `as_strided` as a default arg to save a global lookup
def rotation_matrix2(section, _as_strided=stride_tricks.as_strided):
    section = section.reshape(5, 5)  # sqrt(section.size), sqrt(section.size)
    windows_shape = (3, 3, 3, 3)     # 5 - 3 + 1, 5 - 3 + 1, 3, 3
    windows_strides = section.strides + section.strides
    windows = _as_strided(section, windows_shape, windows_strides)
    rot_filters = windows.reshape(9, 9)
    return rot_filters[rot_filters.std(1).argmin(),:].mean()

def get_rotation_smooth(im, _rm=rotation_matrix2, **kwargs):
    return ndimage.filters.generic_filter(im, _rm, size=5, **kwargs)

if __name__ == '__main__':
    import matplotlib.pyplot as plt
    from scipy.misc import lena
    im = lena()
    im2 = get_rotation_smooth(im)
    #plt.gray()      # Uncomment these lines for
    #plt.imshow(im2) # demo purposes.
    #plt.show()

上面的函数 rotation_matrix2 等同于下面两个函数(加在一起实际上比你原来的函数慢一点，因为 windows 更通用)。这与您的原始代码完全相同——将 9 个 3x3 窗口创建为 5x5 数组，然后将它们重新整形为 9x9 数组以进行处理。

def windows(a, w, _as_strided=stride_tricks.as_strided):
    windows_shape = tuple(sa - sw + 1 for sa, sw in zip(a.shape, w))
    windows_shape += w
    windows_strides = a.strides + a.strides
    return _as_strided(a, windows_shape, windows_strides)

def rotation_matrix1(section, _windows=windows):
    rot_filters = windows(section.reshape(5, 5), (3, 3)).reshape(9, 9)
    return rot_filters[rot_filters.std(1).argmin(),:].mean()

windows 适用于任何维度的数组，只要窗口具有相同数量的维度。以下是其工作原理的分割:

    windows_shape = tuple(sa - sw + 1 for sa, sw in zip(a.shape, w))

我们可以将 windows 数组视为 n 维数组的 n 维数组。外部 n-d 阵列的形状由较大阵列内窗口的自由度决定；在每个维度中，窗口可以占据的位置数等于较大数组的长度减去窗口的长度加一。在这种情况下，我们有一个 3x3 窗口到 5x5 数组，因此外部二维数组是一个 3x3 数组。

    windows_shape += w

内部 n 维数组的形状与窗口本身的形状相同。在我们的例子中，这又是一个 3x3 数组。

现在大踏步前进。我们必须为外部 n-d 数组和内部 n-d 数组定义步幅。但事实证明它们是一样的!毕竟，窗口在较大数组中移动的方式与单个索引在数组中移动的方式相同，对吧？

    windows_strides = a.strides + a.strides

现在我们有了创建窗口所需的所有信息:

    return _as_strided(a, windows_shape, windows_strides)