Generalization of the median filter

annotation

This article describes a unique filter, an article about which appeared in 1990: A. Maslov, V. Sergeev. Identification of a linear distorting system using rank signal processing // Computer Optics. - M., 1990. - Vol.6. - P.97-102. This algorithm is called the “Algorithm of rank processing” and in fact is a generalization of the median filter.
The use of this filter is justified in two cases - to suppress noise and to reduce blurring.

Figure 1 - the original image, 2 - blurred and noisy salt.

Algorithm of signal processing rank - extreme filter

The algorithm consists in processing the image by a local window with the recording of the processing result in a new image:

Suppose we are at the point of the image with coordinates (I, J) - the current count
A local neighborhood of size NxN is considered around the current reference point.
The variational series is constructed from the values of points in the local neighborhood, which we denote by p . The size of this series is N * N.
In the resulting image, the current sample takes on the value according to the following rule:

where k is an algorithm parameter, N is an odd number, Im1 is the original image, Im2 is the resulting image.
If k = (N ^ 2 + 1) / 2 - that is, the center of the variation series - this filter becomes a known median filter. In the future, this parameter will be called indent .

Extreme Filter Properties

The properties of this filter are very useful in practice, as the filter allows not only to compensate for the noise but to eliminate (partially) the effects of blurring the image. The limiting case of this filter with k = (N ^ 2 + 1) / 2 we have a median filter, which only eliminates noise, but does not touch the boundaries, and if the image is blurred, then the blur will remain.
')
When k <(N ^ 2 + 1) / 2, the noise is slightly worse filtered, but the sharpness of the image increases, and at k = 0, the noise is not filtered at all, but is blurred out in the strongest way.

In order not to bother readers with an elementary implementation of this algorithm, I will give Matlab code here.
The testing script will read the image, lubricate it, and add noise. Then the image is restored to this extreme filters.

I1 = imread( 'coins.png' ); h = ones(3,3) / 9; I2 = imfilter(I1,h) ; I3 = imnoise(I2, 'salt & pepper' ,0.02); I4 = im2_rang_filter (I3, 1, 2); figure; imagesc(I1); colormap gray; figure; imagesc(I3); colormap gray; figure; imagesc(I4) colormap gray; * This source code was highlighted with Source Code Highlighter .

I1 = imread( 'coins.png' ); h = ones(3,3) / 9; I2 = imfilter(I1,h) ; I3 = imnoise(I2, 'salt & pepper' ,0.02); I4 = im2_rang_filter (I3, 1, 2); figure; imagesc(I1); colormap gray; figure; imagesc(I3); colormap gray; figure; imagesc(I4) colormap gray; * This source code was highlighted with Source Code Highlighter .
I1 = imread( 'coins.png' ); h = ones(3,3) / 9; I2 = imfilter(I1,h) ; I3 = imnoise(I2, 'salt & pepper' ,0.02); I4 = im2_rang_filter (I3, 1, 2); figure; imagesc(I1); colormap gray; figure; imagesc(I3); colormap gray; figure; imagesc(I4) colormap gray; * This source code was highlighted with Source Code Highlighter .
I1 = imread( 'coins.png' ); h = ones(3,3) / 9; I2 = imfilter(I1,h) ; I3 = imnoise(I2, 'salt & pepper' ,0.02); I4 = im2_rang_filter (I3, 1, 2); figure; imagesc(I1); colormap gray; figure; imagesc(I3); colormap gray; figure; imagesc(I4) colormap gray; * This source code was highlighted with Source Code Highlighter .
I1 = imread( 'coins.png' ); h = ones(3,3) / 9; I2 = imfilter(I1,h) ; I3 = imnoise(I2, 'salt & pepper' ,0.02); I4 = im2_rang_filter (I3, 1, 2); figure; imagesc(I1); colormap gray; figure; imagesc(I3); colormap gray; figure; imagesc(I4) colormap gray; * This source code was highlighted with Source Code Highlighter .
I1 = imread( 'coins.png' ); h = ones(3,3) / 9; I2 = imfilter(I1,h) ; I3 = imnoise(I2, 'salt & pepper' ,0.02); I4 = im2_rang_filter (I3, 1, 2); figure; imagesc(I1); colormap gray; figure; imagesc(I3); colormap gray; figure; imagesc(I4) colormap gray; * This source code was highlighted with Source Code Highlighter .
I1 = imread( 'coins.png' ); h = ones(3,3) / 9; I2 = imfilter(I1,h) ; I3 = imnoise(I2, 'salt & pepper' ,0.02); I4 = im2_rang_filter (I3, 1, 2); figure; imagesc(I1); colormap gray; figure; imagesc(I3); colormap gray; figure; imagesc(I4) colormap gray; * This source code was highlighted with Source Code Highlighter .
I1 = imread( 'coins.png' ); h = ones(3,3) / 9; I2 = imfilter(I1,h) ; I3 = imnoise(I2, 'salt & pepper' ,0.02); I4 = im2_rang_filter (I3, 1, 2); figure; imagesc(I1); colormap gray; figure; imagesc(I3); colormap gray; figure; imagesc(I4) colormap gray; * This source code was highlighted with Source Code Highlighter .
I1 = imread( 'coins.png' ); h = ones(3,3) / 9; I2 = imfilter(I1,h) ; I3 = imnoise(I2, 'salt & pepper' ,0.02); I4 = im2_rang_filter (I3, 1, 2); figure; imagesc(I1); colormap gray; figure; imagesc(I3); colormap gray; figure; imagesc(I4) colormap gray; * This source code was highlighted with Source Code Highlighter .
I1 = imread( 'coins.png' ); h = ones(3,3) / 9; I2 = imfilter(I1,h) ; I3 = imnoise(I2, 'salt & pepper' ,0.02); I4 = im2_rang_filter (I3, 1, 2); figure; imagesc(I1); colormap gray; figure; imagesc(I3); colormap gray; figure; imagesc(I4) colormap gray; * This source code was highlighted with Source Code Highlighter .
I1 = imread( 'coins.png' ); h = ones(3,3) / 9; I2 = imfilter(I1,h) ; I3 = imnoise(I2, 'salt & pepper' ,0.02); I4 = im2_rang_filter (I3, 1, 2); figure; imagesc(I1); colormap gray; figure; imagesc(I3); colormap gray; figure; imagesc(I4) colormap gray; * This source code was highlighted with Source Code Highlighter .
I1 = imread( 'coins.png' ); h = ones(3,3) / 9; I2 = imfilter(I1,h) ; I3 = imnoise(I2, 'salt & pepper' ,0.02); I4 = im2_rang_filter (I3, 1, 2); figure; imagesc(I1); colormap gray; figure; imagesc(I3); colormap gray; figure; imagesc(I4) colormap gray; * This source code was highlighted with Source Code Highlighter .
I1 = imread( 'coins.png' ); h = ones(3,3) / 9; I2 = imfilter(I1,h) ; I3 = imnoise(I2, 'salt & pepper' ,0.02); I4 = im2_rang_filter (I3, 1, 2); figure; imagesc(I1); colormap gray; figure; imagesc(I3); colormap gray; figure; imagesc(I4) colormap gray; * This source code was highlighted with Source Code Highlighter .

The code for the im2_rang_filter filter function is :

function outImage= im2_rang_filter (aImage, aHalfWindowSize, aOtsup) [ver,hor] = size(aImage); wsize = (aHalfWindowSize*2+1)^2; result = zeros(ver,hor); for i = aHalfWindowSize+1 : (ver - aHalfWindowSize) for j = aHalfWindowSize+1 : (hor - aHalfWindowSize) wind = aImage((i-aHalfWindowSize) : (i + aHalfWindowSize), (j-aHalfWindowSize) : (j + aHalfWindowSize)); vec = reshape(wind,1,wsize); vec = sort(vec); if (abs(vec(aOtsup+1) - aImage(i,j)) < abs(vec(wsize - aOtsup) - aImage(i,j)) ) result(i,j) = vec(aOtsup+1); else result(i,j) = vec(wsize - aOtsup); end ; end ; end ; outImage = result; * This source code was highlighted with Source Code Highlighter .

function outImage= im2_rang_filter (aImage, aHalfWindowSize, aOtsup) [ver,hor] = size(aImage); wsize = (aHalfWindowSize*2+1)^2; result = zeros(ver,hor); for i = aHalfWindowSize+1 : (ver - aHalfWindowSize) for j = aHalfWindowSize+1 : (hor - aHalfWindowSize) wind = aImage((i-aHalfWindowSize) : (i + aHalfWindowSize), (j-aHalfWindowSize) : (j + aHalfWindowSize)); vec = reshape(wind,1,wsize); vec = sort(vec); if (abs(vec(aOtsup+1) - aImage(i,j)) < abs(vec(wsize - aOtsup) - aImage(i,j)) ) result(i,j) = vec(aOtsup+1); else result(i,j) = vec(wsize - aOtsup); end ; end ; end ; outImage = result; * This source code was highlighted with Source Code Highlighter .
function outImage= im2_rang_filter (aImage, aHalfWindowSize, aOtsup) [ver,hor] = size(aImage); wsize = (aHalfWindowSize*2+1)^2; result = zeros(ver,hor); for i = aHalfWindowSize+1 : (ver - aHalfWindowSize) for j = aHalfWindowSize+1 : (hor - aHalfWindowSize) wind = aImage((i-aHalfWindowSize) : (i + aHalfWindowSize), (j-aHalfWindowSize) : (j + aHalfWindowSize)); vec = reshape(wind,1,wsize); vec = sort(vec); if (abs(vec(aOtsup+1) - aImage(i,j)) < abs(vec(wsize - aOtsup) - aImage(i,j)) ) result(i,j) = vec(aOtsup+1); else result(i,j) = vec(wsize - aOtsup); end ; end ; end ; outImage = result; * This source code was highlighted with Source Code Highlighter .
function outImage= im2_rang_filter (aImage, aHalfWindowSize, aOtsup) [ver,hor] = size(aImage); wsize = (aHalfWindowSize*2+1)^2; result = zeros(ver,hor); for i = aHalfWindowSize+1 : (ver - aHalfWindowSize) for j = aHalfWindowSize+1 : (hor - aHalfWindowSize) wind = aImage((i-aHalfWindowSize) : (i + aHalfWindowSize), (j-aHalfWindowSize) : (j + aHalfWindowSize)); vec = reshape(wind,1,wsize); vec = sort(vec); if (abs(vec(aOtsup+1) - aImage(i,j)) < abs(vec(wsize - aOtsup) - aImage(i,j)) ) result(i,j) = vec(aOtsup+1); else result(i,j) = vec(wsize - aOtsup); end ; end ; end ; outImage = result; * This source code was highlighted with Source Code Highlighter .
function outImage= im2_rang_filter (aImage, aHalfWindowSize, aOtsup) [ver,hor] = size(aImage); wsize = (aHalfWindowSize*2+1)^2; result = zeros(ver,hor); for i = aHalfWindowSize+1 : (ver - aHalfWindowSize) for j = aHalfWindowSize+1 : (hor - aHalfWindowSize) wind = aImage((i-aHalfWindowSize) : (i + aHalfWindowSize), (j-aHalfWindowSize) : (j + aHalfWindowSize)); vec = reshape(wind,1,wsize); vec = sort(vec); if (abs(vec(aOtsup+1) - aImage(i,j)) < abs(vec(wsize - aOtsup) - aImage(i,j)) ) result(i,j) = vec(aOtsup+1); else result(i,j) = vec(wsize - aOtsup); end ; end ; end ; outImage = result; * This source code was highlighted with Source Code Highlighter .
function outImage= im2_rang_filter (aImage, aHalfWindowSize, aOtsup) [ver,hor] = size(aImage); wsize = (aHalfWindowSize*2+1)^2; result = zeros(ver,hor); for i = aHalfWindowSize+1 : (ver - aHalfWindowSize) for j = aHalfWindowSize+1 : (hor - aHalfWindowSize) wind = aImage((i-aHalfWindowSize) : (i + aHalfWindowSize), (j-aHalfWindowSize) : (j + aHalfWindowSize)); vec = reshape(wind,1,wsize); vec = sort(vec); if (abs(vec(aOtsup+1) - aImage(i,j)) < abs(vec(wsize - aOtsup) - aImage(i,j)) ) result(i,j) = vec(aOtsup+1); else result(i,j) = vec(wsize - aOtsup); end ; end ; end ; outImage = result; * This source code was highlighted with Source Code Highlighter .
function outImage= im2_rang_filter (aImage, aHalfWindowSize, aOtsup) [ver,hor] = size(aImage); wsize = (aHalfWindowSize*2+1)^2; result = zeros(ver,hor); for i = aHalfWindowSize+1 : (ver - aHalfWindowSize) for j = aHalfWindowSize+1 : (hor - aHalfWindowSize) wind = aImage((i-aHalfWindowSize) : (i + aHalfWindowSize), (j-aHalfWindowSize) : (j + aHalfWindowSize)); vec = reshape(wind,1,wsize); vec = sort(vec); if (abs(vec(aOtsup+1) - aImage(i,j)) < abs(vec(wsize - aOtsup) - aImage(i,j)) ) result(i,j) = vec(aOtsup+1); else result(i,j) = vec(wsize - aOtsup); end ; end ; end ; outImage = result; * This source code was highlighted with Source Code Highlighter .
function outImage= im2_rang_filter (aImage, aHalfWindowSize, aOtsup) [ver,hor] = size(aImage); wsize = (aHalfWindowSize*2+1)^2; result = zeros(ver,hor); for i = aHalfWindowSize+1 : (ver - aHalfWindowSize) for j = aHalfWindowSize+1 : (hor - aHalfWindowSize) wind = aImage((i-aHalfWindowSize) : (i + aHalfWindowSize), (j-aHalfWindowSize) : (j + aHalfWindowSize)); vec = reshape(wind,1,wsize); vec = sort(vec); if (abs(vec(aOtsup+1) - aImage(i,j)) < abs(vec(wsize - aOtsup) - aImage(i,j)) ) result(i,j) = vec(aOtsup+1); else result(i,j) = vec(wsize - aOtsup); end ; end ; end ; outImage = result; * This source code was highlighted with Source Code Highlighter .
function outImage= im2_rang_filter (aImage, aHalfWindowSize, aOtsup) [ver,hor] = size(aImage); wsize = (aHalfWindowSize*2+1)^2; result = zeros(ver,hor); for i = aHalfWindowSize+1 : (ver - aHalfWindowSize) for j = aHalfWindowSize+1 : (hor - aHalfWindowSize) wind = aImage((i-aHalfWindowSize) : (i + aHalfWindowSize), (j-aHalfWindowSize) : (j + aHalfWindowSize)); vec = reshape(wind,1,wsize); vec = sort(vec); if (abs(vec(aOtsup+1) - aImage(i,j)) < abs(vec(wsize - aOtsup) - aImage(i,j)) ) result(i,j) = vec(aOtsup+1); else result(i,j) = vec(wsize - aOtsup); end ; end ; end ; outImage = result; * This source code was highlighted with Source Code Highlighter .
function outImage= im2_rang_filter (aImage, aHalfWindowSize, aOtsup) [ver,hor] = size(aImage); wsize = (aHalfWindowSize*2+1)^2; result = zeros(ver,hor); for i = aHalfWindowSize+1 : (ver - aHalfWindowSize) for j = aHalfWindowSize+1 : (hor - aHalfWindowSize) wind = aImage((i-aHalfWindowSize) : (i + aHalfWindowSize), (j-aHalfWindowSize) : (j + aHalfWindowSize)); vec = reshape(wind,1,wsize); vec = sort(vec); if (abs(vec(aOtsup+1) - aImage(i,j)) < abs(vec(wsize - aOtsup) - aImage(i,j)) ) result(i,j) = vec(aOtsup+1); else result(i,j) = vec(wsize - aOtsup); end ; end ; end ; outImage = result; * This source code was highlighted with Source Code Highlighter .
function outImage= im2_rang_filter (aImage, aHalfWindowSize, aOtsup) [ver,hor] = size(aImage); wsize = (aHalfWindowSize*2+1)^2; result = zeros(ver,hor); for i = aHalfWindowSize+1 : (ver - aHalfWindowSize) for j = aHalfWindowSize+1 : (hor - aHalfWindowSize) wind = aImage((i-aHalfWindowSize) : (i + aHalfWindowSize), (j-aHalfWindowSize) : (j + aHalfWindowSize)); vec = reshape(wind,1,wsize); vec = sort(vec); if (abs(vec(aOtsup+1) - aImage(i,j)) < abs(vec(wsize - aOtsup) - aImage(i,j)) ) result(i,j) = vec(aOtsup+1); else result(i,j) = vec(wsize - aOtsup); end ; end ; end ; outImage = result; * This source code was highlighted with Source Code Highlighter .
function outImage= im2_rang_filter (aImage, aHalfWindowSize, aOtsup) [ver,hor] = size(aImage); wsize = (aHalfWindowSize*2+1)^2; result = zeros(ver,hor); for i = aHalfWindowSize+1 : (ver - aHalfWindowSize) for j = aHalfWindowSize+1 : (hor - aHalfWindowSize) wind = aImage((i-aHalfWindowSize) : (i + aHalfWindowSize), (j-aHalfWindowSize) : (j + aHalfWindowSize)); vec = reshape(wind,1,wsize); vec = sort(vec); if (abs(vec(aOtsup+1) - aImage(i,j)) < abs(vec(wsize - aOtsup) - aImage(i,j)) ) result(i,j) = vec(aOtsup+1); else result(i,j) = vec(wsize - aOtsup); end ; end ; end ; outImage = result; * This source code was highlighted with Source Code Highlighter .
function outImage= im2_rang_filter (aImage, aHalfWindowSize, aOtsup) [ver,hor] = size(aImage); wsize = (aHalfWindowSize*2+1)^2; result = zeros(ver,hor); for i = aHalfWindowSize+1 : (ver - aHalfWindowSize) for j = aHalfWindowSize+1 : (hor - aHalfWindowSize) wind = aImage((i-aHalfWindowSize) : (i + aHalfWindowSize), (j-aHalfWindowSize) : (j + aHalfWindowSize)); vec = reshape(wind,1,wsize); vec = sort(vec); if (abs(vec(aOtsup+1) - aImage(i,j)) < abs(vec(wsize - aOtsup) - aImage(i,j)) ) result(i,j) = vec(aOtsup+1); else result(i,j) = vec(wsize - aOtsup); end ; end ; end ; outImage = result; * This source code was highlighted with Source Code Highlighter .
function outImage= im2_rang_filter (aImage, aHalfWindowSize, aOtsup) [ver,hor] = size(aImage); wsize = (aHalfWindowSize*2+1)^2; result = zeros(ver,hor); for i = aHalfWindowSize+1 : (ver - aHalfWindowSize) for j = aHalfWindowSize+1 : (hor - aHalfWindowSize) wind = aImage((i-aHalfWindowSize) : (i + aHalfWindowSize), (j-aHalfWindowSize) : (j + aHalfWindowSize)); vec = reshape(wind,1,wsize); vec = sort(vec); if (abs(vec(aOtsup+1) - aImage(i,j)) < abs(vec(wsize - aOtsup) - aImage(i,j)) ) result(i,j) = vec(aOtsup+1); else result(i,j) = vec(wsize - aOtsup); end ; end ; end ; outImage = result; * This source code was highlighted with Source Code Highlighter .
function outImage= im2_rang_filter (aImage, aHalfWindowSize, aOtsup) [ver,hor] = size(aImage); wsize = (aHalfWindowSize*2+1)^2; result = zeros(ver,hor); for i = aHalfWindowSize+1 : (ver - aHalfWindowSize) for j = aHalfWindowSize+1 : (hor - aHalfWindowSize) wind = aImage((i-aHalfWindowSize) : (i + aHalfWindowSize), (j-aHalfWindowSize) : (j + aHalfWindowSize)); vec = reshape(wind,1,wsize); vec = sort(vec); if (abs(vec(aOtsup+1) - aImage(i,j)) < abs(vec(wsize - aOtsup) - aImage(i,j)) ) result(i,j) = vec(aOtsup+1); else result(i,j) = vec(wsize - aOtsup); end ; end ; end ; outImage = result; * This source code was highlighted with Source Code Highlighter .
function outImage= im2_rang_filter (aImage, aHalfWindowSize, aOtsup) [ver,hor] = size(aImage); wsize = (aHalfWindowSize*2+1)^2; result = zeros(ver,hor); for i = aHalfWindowSize+1 : (ver - aHalfWindowSize) for j = aHalfWindowSize+1 : (hor - aHalfWindowSize) wind = aImage((i-aHalfWindowSize) : (i + aHalfWindowSize), (j-aHalfWindowSize) : (j + aHalfWindowSize)); vec = reshape(wind,1,wsize); vec = sort(vec); if (abs(vec(aOtsup+1) - aImage(i,j)) < abs(vec(wsize - aOtsup) - aImage(i,j)) ) result(i,j) = vec(aOtsup+1); else result(i,j) = vec(wsize - aOtsup); end ; end ; end ; outImage = result; * This source code was highlighted with Source Code Highlighter .
function outImage= im2_rang_filter (aImage, aHalfWindowSize, aOtsup) [ver,hor] = size(aImage); wsize = (aHalfWindowSize*2+1)^2; result = zeros(ver,hor); for i = aHalfWindowSize+1 : (ver - aHalfWindowSize) for j = aHalfWindowSize+1 : (hor - aHalfWindowSize) wind = aImage((i-aHalfWindowSize) : (i + aHalfWindowSize), (j-aHalfWindowSize) : (j + aHalfWindowSize)); vec = reshape(wind,1,wsize); vec = sort(vec); if (abs(vec(aOtsup+1) - aImage(i,j)) < abs(vec(wsize - aOtsup) - aImage(i,j)) ) result(i,j) = vec(aOtsup+1); else result(i,j) = vec(wsize - aOtsup); end ; end ; end ; outImage = result; * This source code was highlighted with Source Code Highlighter .
function outImage= im2_rang_filter (aImage, aHalfWindowSize, aOtsup) [ver,hor] = size(aImage); wsize = (aHalfWindowSize*2+1)^2; result = zeros(ver,hor); for i = aHalfWindowSize+1 : (ver - aHalfWindowSize) for j = aHalfWindowSize+1 : (hor - aHalfWindowSize) wind = aImage((i-aHalfWindowSize) : (i + aHalfWindowSize), (j-aHalfWindowSize) : (j + aHalfWindowSize)); vec = reshape(wind,1,wsize); vec = sort(vec); if (abs(vec(aOtsup+1) - aImage(i,j)) < abs(vec(wsize - aOtsup) - aImage(i,j)) ) result(i,j) = vec(aOtsup+1); else result(i,j) = vec(wsize - aOtsup); end ; end ; end ; outImage = result; * This source code was highlighted with Source Code Highlighter .

Experiments

Below are the results of filtering the window 3 by 3 with the variation of the parameter k = 1,2,5. The first drawing is the original image, the second drawing is distorted, then one after the other are filtered with k = 1, 2, 5:
Original: