function J = adaptiveMedianFilter(I, S_max)
% Adaptive median filter (two-stage algorithm from the lecture)
%   I     : noisy image (double)
%   S_max : maximum allowed window size
[m, n] = size(I);
J = zeros(m, n);
pad = floor(S_max/2);
P = padarray(I, [pad pad], 'symmetric');
for r = 1:m
    for c = 1:n
        ip = r + pad; jp = c + pad;
        z_xy = P(ip, jp);
        windowSize = 3;
        done = false;
        while ~done
            h = floor(windowSize/2);
            local = P(ip-h:ip+h, jp-h:jp+h);
            z_min = min(local(:));
            z_max = max(local(:));
            z_med = median(local(:));
            % Stage A
            A1 = z_med - z_min;
            A2 = z_med - z_max;
            if (A1 > 0) && (A2 < 0)
                % Stage B
                B1 = z_xy - z_min;
                B2 = z_xy - z_max;
                if (B1 > 0) && (B2 < 0)
                    J(r, c) = z_xy;
                else
                    J(r, c) = z_med;
                end
                done = true;
            else
                windowSize = windowSize + 2;
                if windowSize > S_max
                    J(r, c) = z_med;
                    done = true;
                end
            end
        end
    end
end
end