function out_im = median_image( in_im, median_order, n_iter )
  
  if nargin < 3
    error( [ mfilename ' needs three input arguments.' ] );
  end
  
  in_class = class( in_im );
  is_double = strcmp( in_class, 'double' );
  if ~is_double
    in_im = double( in_im );
  end
  
  nrows = size( in_im, 1 );
  ncols = size( in_im, 2 );
  nlayers = size( in_im, 3 );
  
  n = 0;
  for b = -median_order:median_order
    for c = -median_order:median_order
      if b^2 + c^2 <= median_order ^ 2
	n = n+1;
      end
    end
  end
  

  out_im = in_im;
  
  for z = 1:n_iter
    
    for a = 1:nlayers
      
      % To fill up the borders ( not very important in our case )
      tmp = median( reshape( out_im( :,:,a ), 1, [] ) );
      tmp = repmat( tmp, [ nrows, ncols, n ] );
      
      i_tmp = 1;
      
      for b = -median_order:median_order

	for c = -median_order:median_order

	  if b^2+c^2 <= median_order^2
	    
	    % Copy the largest possible portion of the original image
	    % with a displacement ( b, c )
	    
	    orig.rowmin = max( 1, 1-b );
	    orig.rowmax = min( nrows, nrows-b );
	    orig.colmin = max( 1, 1-c );
	    orig.colmax = min( ncols, ncols-c );
	    
	    dest.rowmin = orig.rowmin + b;
	    dest.rowmax = orig.rowmax + b;
	    dest.colmin = orig.colmin + c;
	    dest.colmax = orig.colmax + c;
	    
	    tmp( dest.rowmin:dest.rowmax, dest.colmin:dest.colmax, i_tmp ) = out_im( orig.rowmin:orig.rowmax, orig.colmin:orig.colmax, a );
	    
	    i_tmp = i_tmp + 1;

	  end
	  
	end
      end
      
      % Compute median
      
      out_im( :,:,a ) = median( tmp, 3 );
      
    end
    
  end
  
  
  % Return same class
  if ~is_double
    out_im = eval( [ in_class '( out_im );' ] );
  end