问题描述
我正在尝试用不同的方法在 Python 中计算图像的窗口标准偏差。
首先,我编写了一个“天真的”方法,它在元素上使用双循环,并使用 np.std() 返回窗口值的标准偏差。然后,我使用了建议的实现 here 和 here
注意:我首先将图像 dtype 转换为 int
,因为在第二个和第三个算法中平方图像值会溢出 uint8
,这是导入的 dtype。
使用第二种和第三种算法计算的结果显示了一些奇怪的扇形效果,在云区域尤其明显。 Std dev using dtype int
如果我将输入图像值转换为 float32
,这种效果就会消失:Std dev using dtype float32
我知道计算 here 的方差时存在数值稳定性问题。对于第二种算法,我采用方差绝对值的平方根来解决这个问题。
为什么在使用 int
而不是 float32
的输入数据类型时会出现这种“扇形”,当 int
无论如何在幕后转换为 float
以执行浮点计算?这是我正在使用的测试图像:castle.png
import numpy as np
import matplotlib.pyplot as plt
import cv2
from scipy.ndimage.filters import uniform_filter
from timeit import timeit
img = cv2.imread('castle.png')
img = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY).astype(int) # CHANGE ME
# Naive double loop method
def stddev_naive(img,kernel_size):
out = np.zeros_like(img)
# Applying reflected boundary conditions for kernel padding
ext = int((kernel_size-1)/2)
x_symm = np.concatenate((img[1:1+ext,:][::-1],img,img[-ext-1:-1,:][::-1])) # Reflected boundary condition applied in x
img_bc = np.concatenate((x_symm.transpose()[1:1+ext,x_symm.transpose(),x_symm.transpose()[-ext-1:-1,:][::-1])).transpose() # Reflected boundary condition applied in y
for i in np.arange(img.shape[0]):
for j in np.arange(img.shape[1]):
out[i,j] = np.std(img_bc[i:i+2*ext+1,j:j+2*ext+1])
return out
# https://stackoverflow.com/a/18422519/16078850
def window_stdev(arr,radius):
c1 = uniform_filter(arr,radius*2,mode='constant',origin=-radius)
c2 = uniform_filter(arr*arr,origin=-radius)
return ((c2 - c1*c1)**.5)[:-radius*2+1,:-radius*2+1]
# https://stackoverflow.com/a/36266187
def winVar(img,kernel_size):
wmean,wsqrmean = (cv2.BoxFilter(x,-1,(kernel_size,kernel_size),borderType=cv2.BORDER_REFLECT) for x in (img,img**2))
return np.abs(wsqrmean - wmean*wmean)**0.5
out1 = stddev_naive(img,11)
out2 = window_stdev(img,5)
out3 = winVar(img,11)
# Show outputs
fig,(ax1,ax2,ax3) = plt.subplots(1,3)
ax1.imshow(out1,cmap='gray')
ax1.set_title('Naive method',fontsize=24)
ax2.imshow(out2,cmap='gray')
ax2.set_title('Using scipy.ndimage.uniform_filter',fontsize=24)
ax3.imshow(out3,cmap='gray')
ax3.set_title('Using cv2.BoxFilter',fontsize=24)
plt.suptitle('Using dtype \'int\'',fontsize=30)
plt.show()
解决方法
暂无找到可以解决该程序问题的有效方法,小编努力寻找整理中!
如果你已经找到好的解决方法,欢迎将解决方案带上本链接一起发送给小编。
小编邮箱:dio#foxmail.com (将#修改为@)