clear all;
% close all;

OriginalImage = imread('cameraman.tif');
f = imnoise(OriginalImage, 'salt & pepper', 0.25);  %加椒盐噪声，噪声密度＝0.25.
[M,N] = size(f);
figure;
imshow(OriginalImage(4:M-3,4:N-3), [0 255]);
title('NoisyImage');

f = double(f);

result1 = medfilt2(f,[7,7]);
figure;
imshow(result1(4:M-3,4:N-3), [0 255]);
title('中值滤波');
clear result1;

%自适应中值滤波
result2 = zeros(M,N);
Smax = 7;
for i=4:M-3,
    for j=4:N-3,
        S = 3;
        Zxy = f(i,j);
        i_min = i-1; i_max = i+1; j_min = j-1;  j_max = j+1;
        while S<=Smax,
            tmp = f(i_min:i_max,j_min:j_max);
            tmp = tmp(:);    % 把二维矩阵编程列向量
            Zmed = median(tmp);
            Zmin = min(tmp);
            Zmax = max(tmp);
            A1 = Zmed - Zmin;
            A2 = Zmed - Zmax;
            if (A1>0)&(A2<0),   %B层
                B1 = Zxy - Zmin;
                B2 = Zxy - Zmax;
                if (B1>0)&(B2<0),
                    result2(i,j) = Zxy;
                else
                    result2(i,j) = Zmed;
                end
                S=Smax+2;  %为了跳出while循环
            else
                S = S+2;
                if S<=Smax,
                    i_min = i_min-1; i_max = i_max+1; j_min = j_min-1;  j_max = j_max+1;
                else
                    result2(i,j) = Zxy;   %Zmed
                    S=Smax+2;  %为了跳出while循环
                end
            end
        end  % end of (while S<=Smax,)
    end  % end of (for j=4:N-3,)
end  % end of (for i=4:M-3,)
figure;
imshow(result2(4:M-3,4:N-3), [0 255]);
title('自适应中值滤波');