本想再用滑动窗口计算写代码计算局部方差,网上查了下,matlab提供stdfilt函数来计算局部方差:
matlab内的解释是:
stdfilt
图像的局部方差
语法
J = stdfilt(I)
J = stdfilt(I, NHOOD)
Description
J = stdfilt(I) returns thearray J, where each output pixel contains the standarddeviation of the 3-by-3 neighborhood around the corresponding pixelin the input image I. I canhave any dimension. The output image J is the samesize as the input image I.
For pixels on the borders of I, stdfilt usessymmetric padding. In symmetric padding, the values of padding pixelsare a mirror reflection of the border pixels in I.
J = stdfilt(I, NHOOD) calculatesthe local standard deviation of the input image I,where you specify the neighborhood in NHOOD. NHOOD isa multidimensional array of zeros and ones where the nonzero elementsspecify the neighbors. NHOOD's size must be oddin each dimension.
By default, stdfilt uses the neighborhood ones(3). stdfilt determinesthe center element of the neighborhood by floor((size(NHOOD)+ 1)/2).
Class Support
I can be logical or numeric and must bereal and nonsparse. NHOOD can be logical or numericand must contain zeros and/or ones. J is of classdouble.
Notes
To specify neighborhoods of various shapes, such as a disk,use the strel function to create a structuringelement object and then use the getnhood functionto extract the neighborhood from the structuring element object.
Examples
I = imread('circuit.tif');
J = stdfilt(I);
imshow(I);
figure, imshow(J,[]);
matlab就是如此强大。
针对我们的问题,我们可以这样调用:
I = imread('E:\\110.jpg');
I = im2double(rgb2gray(I));%转化为灰度图像
J= stdfilt(I) ;%计算局部方差
imshow(J)%显示图像
dlmwrite('E:\\a.txt',J,'delimiter',' ','newline','pc'); %当然可以保存数据
但更有意思的是,我们可以在Python上模拟实现:
import cv2
import numpy as np
img = cv2.imread('fruits.jpg', True)
img = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY)
img = img / 255.0
# c = imfilter(I,h,'symmetric');
h = np.ones((3,3))
n = h.sum()
n1 = n - 1
c1 = cv2.filter2D(img**2, -1, h/n1, borderType=cv2.BORDER_REFLECT)
c2 = cv2.filter2D(img, -1, h, borderType=cv2.BORDER_REFLECT)**2 / (n*n1)
J = np.sqrt( np.maximum(c1-c2,0) )
cv2.imshow('stdfilt', J)
cv2.waitKey(0)
cv2.destroyWindow('stdfilt')