kuwahara_filter滤波器----MATLAB

%% kuwahara_filter滤波器
clc,clear,close all  % 清理命令区、清理工作区、关闭显示图形
warning off       % 消除警告
feature jit off      % 加速代码运行
[filename ,pathname]=...
    uigetfile({'*.bmp';'*.tif';'*.jpg';},'选择图片'); %选择图片路径
str=[pathname filename]; % 合成路径+文件名
im = imread(str);        % 读图
% 转化为灰度图像
if size(im,3)==1
    im = im;
else
    im = rgb2gray(im);
end
im = imnoise(im,'gaussian',0,1e-3); % 原图像 + 白噪声

figure,
subplot(121),imshow(im);title('原始图像')
colormap(jet)  % 颜色
shading interp % 消隐
im1 = kuwahara(im,5);
subplot(122),imshow(im1);title('kuwahara滤波图像')
colormap(jet)  % 颜色
shading interp % 消隐
function [Y,Xpad] = kuwahara(X,WINSZ)
% kuwahara_filter滤波器
% 对图像边界轮廓有较强的保护作用
% kuwahara nonlinear edge-preserving filtering
% 函数输入：
%          X：二维图像矩阵
%          WINSZ： window size
% 函数输出：
%          Y：滤波图像
%          Xpad：点扩展矩阵
% The Kuwahara filter 4块，(最直观的如下所示，5x5 pixels). 
%
%    ( a  a  ab   b  b)
%    ( a  a  ab   b  b)
%    (ac ac abcd bd bd)
%    ( c  c  cd   d  d)
%    ( c  c  cd   d  d)

if nargin < 2
    WINSZ = 5;  % 默认值
end

if ~isa(X,'double')
    X = im2double(X);  % double类型
end

Pad = floor(WINSZ/2);  % 向下取整

Xpad = padarray(X,[Pad, Pad],'replicate');  % 点扩展
% A = [ 1 2; 
%       3 4 ];
% B = padarray(A,[3 2],'replicate','post')
% B =
%      1     2     2     2
%      3     4     4     4
%      3     4     4     4
%      3     4     4     4
%      3     4     4     4
[padRows,padCols] = size(Xpad);                    % 求维数
Y = zeros(size(X));

nRowIters = length((Pad+1):(padRows-Pad));
count = 1;
for i = (Pad+1):(padRows-Pad)
    for j = (Pad+1):(padCols-Pad)
        % 分成每一个小窗，Q1~Q4
        w = Xpad((i-Pad):(i+Pad),(j-Pad):(j+Pad));
        wnw = w(1:(Pad+1),1:(Pad+1));
        wne = w(1:(Pad+1),(Pad+1):WINSZ);
        wsw = w((Pad+1):WINSZ,1:(Pad+1));
        wse = w((Pad+1):WINSZ,(Pad+1):WINSZ);
        % 计算方差
        s = var([wnw(:),wne(:),wsw(:),wse(:)]);
        m = mean([wnw(:),wne(:),wsw(:),wse(:)]);
        [y,k] = min(s);  
        % 计算均值
        Y(i,j) = m(k);
    end

end