Matlab计算图像的局部方差

在图像处理中经常会用到方差和局部方差的概念，这里就给出计算图像局部方法的代码，方便使用。

注：代码是关于导向滤波论文中的源码，大家可以自己下载

function dataVariance = getLocalVariance(srcData, nHeight, nWidth, r)
% 计算局部方差dataVariance
% r为局部窗口大小
N               = boxfilter(ones(nHeight, nWidth), r);
srcDataMean     = boxfilter(srcData, r) ./ N;
srcDataDataMean = boxfilter(srcData.*srcData, r) ./ N;
dataVariance    = srcDataDataMean - srcDataMean .* srcDataMean;
end

function imDst = boxfilter(imSrc, r)

%   BOXFILTER   O(1) time box filtering using cumulative sum
%
%   - Definition imDst(x, y)=sum(sum(imSrc(x-r:x+r,y-r:y+r)));
%   - Running time independent of r; 
%   - Equivalent to the function: colfilt(imSrc, [2*r+1, 2*r+1], 'sliding', @sum);
%   - But much faster.

[hei, wid] = size(imSrc);
imSrc = double(imSrc) %数据类型转换，防止数据超出类型可表示的范围
imDst = zeros(size(imSrc));

%cumulative sum over Y axis
imCum = cumsum(imSrc, 1);
%difference over Y axis
imDst(1:r+1, :) = imCum(1+r:2*r+1, :);
imDst(r+2:hei-r, :) = imCum(2*r+2:hei, :) - imCum(1:hei-2*r-1, :);
imDst(hei-r+1:hei, :) = repmat(imCum(hei, :), [r, 1]) - imCum(hei-2*r:hei-r-1, :);

%cumulative sum over X axis
imCum = cumsum(imDst, 2);
%difference over Y axis
imDst(:, 1:r+1) = imCum(:, 1+r:2*r+1);
imDst(:, r+2:wid-r) = imCum(:, 2*r+2:wid) - imCum(:, 1:wid-2*r-1);
imDst(:, wid-r+1:wid) = repmat(imCum(:, wid), [1, r]) - imCum(:, wid-2*r:wid-r-1);
end