function MMF = MMF(image,N);

%Multilevel Median filter

%reference is http://video.ee.ntu.edu.tw/~pony/fuzzy.html

% returns an image of half the size of image filtered by a multilevel
% median filter. 2N+1 is the size of the window considered when taking
% the median



[row,col] = size(image);
rep = replicateBounds(image, N, N);

Z1 = zeros([row col 2*N+1]);
Z2 = zeros([row col 2*N+1]);
Z3 = zeros([row col 2*N+1]);
Z4 = zeros([row col 2*N+1]);

s = [row col];
for i = -N:N
   Z1(:,:,i+N+1) = getShift(rep, 0, i, s);
   Z2(:,:,i+N+1) = getShift(rep, i, 0, s);
   Z3(:,:,i+N+1) = getShift(rep, i, i, s);
   Z4(:,:,i+N+1) = getShift(rep, i, -i, s);
end

z = zeros([row col 4]);
z(:,:,1) = median(Z1,3);
z(:,:,2) = median(Z2,3);
z(:,:,3) = median(Z3,3);
z(:,:,4) = median(Z4,3);

y = zeros([row col 3]);
y(:,:,1) = min(z,[],3);
y(:,:,2) = max(z,[],3);
y(:,:,3) = image;

MMF = median(y, 3);

