高效计算Python图像的方差

Question

我正在一个项目中工作，需要获取图像的方差。
目前我有2种方法（两种方法都可以工作，但速度非常慢）：

1. 逐个计算每个像素点的方差：

这是使用numpy编写的代码，varianceMatrix是输出结果

varianceMatrix = np.zeros(im.shape,np.uint8)
w = 1              # the radius of pixels neighbors 
ny = len(im)
nx = len(im[0])


for i in range(w,nx-w):
    for j in range(w,ny-w):

        sampleframe = im[j-w:j+w, i-w:i+w]
        variance    = np.var(sampleframe)
        varianceMatrix[j][i] = int(variance)

return varianceMatrix   

2. 使用现有的 scipy 函数：

这是 scipy 函数：

from scipy import ndimage

varianceMatrix = ndimage.generic_filter(im, np.var, size = 3)

scipy函数更快一些，但并不是很快。我正在寻找更好的计算方差的替代方法。

有任何想法吗？？？

- est.tenorio

Answer 1

这里是使用OpenCV的快速解决方案：

import cv2

def winVar(img, wlen):
  wmean, wsqrmean = (cv2.boxFilter(x, -1, (wlen, wlen),
    borderType=cv2.BORDER_REFLECT) for x in (img, img*img))
  return wsqrmean - wmean*wmean

在我的电脑上，以以下示例为例，winVar() 比 ndimage.generic_filter() 快2915倍，并且比 sliding_img_var() 快10.8倍（请参见pv的答案）。

In [66]: img = np.random.randint(0, 256, (500,500)).astype(np.float)

In [67]: %timeit winVar(img, 3)
100 loops, best of 3: 1.76 ms per loop

In [68]: %timeit ndimage.generic_filter(img, np.var, size=3)
1 loops, best of 3: 5.13 s per loop

In [69]: %timeit sliding_img_var(img, 1)
100 loops, best of 3: 19 ms per loop

结果与 ndimage.generic_filter() 相匹配:

In [70]: np.allclose(winVar(img, 3), ndimage.generic_filter(img, np.var, size=3))
Out[70]: True

- Ulrich Stern

---

可以使用ndimage.uniform_filter()代替cv2.boxFilter()，请参见此答案以获取类似问题的解决方案。对于我在这里使用的示例，OpenCV版本快了4.1倍。 - Ulrich Stern

要计算标准差（即将winVar()转换为winStd()），只需更改为return np.sqrt(wsqrmean - wmean*wmean)。 - Ulrich Stern

Answer 2

您可以使用一种广为人知的滑动窗口步幅技巧来加速计算。它在不复制数据的情况下，在数组末尾添加了两个“虚拟维度”，然后对它们进行方差计算。

请注意，在您的代码中，im[j-w:j+w, ..]遍历的索引是j-w,j-w+1,...,j+w-1，最后一个索引是排除在外的，这可能不是您想要的。此外，方差大于uint8范围，因此您最终会得到整数环绕。

import numpy as np
import time
np.random.seed(1234)

img = (np.random.rand(200, 200)*256).astype(np.uint8)

def sliding_window(a, window, axis=-1):
    shape = list(a.shape) + [window]
    shape[axis] -= window - 1
    if shape[axis] < 0:
        raise ValueError("Array too small")
    strides = a.strides + (a.strides[axis],)
    return np.lib.stride_tricks.as_strided(a, shape=shape, strides=strides)

def sliding_img_var(img, window):
    if window <= 0:
        raise ValueError("invalid window size")
    buf = sliding_window(img, 2*window, 0)
    buf = sliding_window(buf, 2*window, 1)

    out = np.zeros(img.shape, dtype=np.float32)
    np.var(buf[:-1,:-1], axis=(-1,-2), out=out[window:-window,window:-window])
    return out

def looping_img_var(im, w):
    nx, ny = img.shape
    varianceMatrix = np.zeros(im.shape, np.float32)
    for i in range(w,nx-w):
        for j in range(w,ny-w):
            sampleframe = im[j-w:j+w, i-w:i+w]
            variance    = np.var(sampleframe)
            varianceMatrix[j][i] = variance
    return varianceMatrix

np.set_printoptions(linewidth=1000, edgeitems=5)
start = time.time()
print(sliding_img_var(img, 1))
time_sliding = time.time() - start
start = time.time()
print(looping_img_var(img, 1))
time_looping = time.time() - start
print("duration: sliding: {0} s, looping: {1} s".format(time_sliding, time_looping))

- pv.

---

这是我机器上输出的最后一行，显示加速比：持续时间：滑动：0.00510311126709秒，循环：0.955919027328秒。 - Ulrich Stern

你如何处理二维数组中的NaN值？ - webapp

Answer 3

如果使用ndimage.generic_filter方法不够快，您可以在Cython中编写自己优化的方差计算实现。

- rth