function M = perform_median_filtering(M,k)% perform_median_filtering - perform moving average median%%   M = perform_median_filtering(M,k);%%   k is the half width of the window (detult k=1).%%   This filtering method is very efficient to remove impulsive or salt and%   peper noise.%%   Copyright (c) 2007 Gabriel Peyreif nargin<2    k = 1;endw = 2*k+1;n = size(M,1);options.sampling = 'uniform';H = compute_patch_library(M,k,options);H = reshape(H, [w*w n*n]);H = median(H); M = reshape(H, n,n);function [H,X,Y] = compute_patch_library(M,w,options)% [H,X,Y] = compute_patch_library(M,w,options);%%   M is the texture%   w is the half-width of the patch, so that each patch%       has size (2*w+1,2*w+1,s) where s is the number of colors.%%   H(:,:,:,i) is the ith patch (can be a color patch).%   X(i) is the x location of H(:,:,:,i) in the image M.%   Y(i) is the y location of H(:,:,:,i) in the image M.%%   options.sampling can be set to 'random' or 'uniform'%   If options.sampling=='random' then you can define%       options.nbr as the number of random patches and/or%       option.locations_x and option.locations_y.%%   If w is a vector of integer size, then the algorithm %   compute a set of library, one for each size, %   and use the same location X,Y for each patch.%   H{k} is the k-th library.%%   You can define options.wmax to avoid too big patch.%   In this case, if w>wmax, the texture M will be filtered, and the patches%   will be subsampled to match the size.%%   Copyright (c) 2006 Gabriel Peyr