Filtre médian

1. Introduction

La médiane d'une liste est l'élément m de cette liste tel qu'il y ait autant d'éléments inférieurs à m que d'éléments supérieurs à m. Par exemple, la médiane de la liste [0,1,5,10,100] est 5.

Un filtre médian est un filtre statistique qui utilise une fenêtre glissante sur le signal (comme un filtre de convolution). La fenêtre comporte un nombre impair de points. Le point central de la fenêtre est remplacé par la médiane des valeurs du signal sur toute la fenêtre.

2. Sélection d'un élément dans une liste

2.a. Définition

L'opération de sélection de l'élément de rang r dans une liste consiste à trouver l'élément de cette liste qui est supérieur ou égal à r-1 éléments, et inférieur tous les autres. Par exemple, l'élément de rang 2 de la liste [5,0,10,1,3] est 1. En général, la liste n'est pas ordonnée. La sélection d'un élément s'apparente donc à une opération de tri.

2.b. Tri rapide

Le tri rapide est un algorithme récursif qui repose sur l'opération de partitionnement d'une liste, qui consiste à choisir un élément de cette liste, que l'on appelera le pivot, et à réarranger la liste pour que tous les éléments inférieurs ou égaux au pivot soient situés avant le pivot et tous les éléments supérieurs au pivot soient situés après le pivot.

L'algorithme de partitionnement procède par échange d'éléments dans la liste. Il parcourt les éléments dans l'orde. Un indice i repère la frontière entre les éléments inférieurs ou égaux au pivot et les éléments supérieurs au pivot. Chaque élément est comparé au pivot et placé du bon côté de la frontière mobile.

Voici une fonction qui effectue le partitionnement en prenant le dernier élément comme pivot. On transmet la liste à traiter et les indices du début et de la fin de la partie de la liste à traiter. La fonction renvoie l'indice de la position finale du pivot dans la sous-liste traitée.

def partition(liste,debut,fin):
    pivot = liste[fin]
    i = debut
    for j in range(debut,fin):
        if liste[j]<=pivot:
            liste[i],liste[j] = liste[j],liste[i]
            i += 1
    liste[i],liste[fin] = liste[fin],liste[i]
    return i

L'algorithme de tri rapide d'une liste consiste à partionner la liste et à trier les deux sous-listes situés de part et d'autre du pivot, par un appel récursif. Voici la fonction de tri rapide :

def tri_rapide(liste,debut,fin):
    if debut<fin:
        i = partition(liste,debut,fin)
        tri_rapide(liste,debut,i-1)
        tri_rapide(liste,i+1,fin)
    return liste

On écrit aussi une fonction d'appel (pour un tableau numpy) :

def trier(liste):
    return tri_rapide(liste.copy(),0,len(liste)-1)

Exemple :

import numpy.random
liste = numpy.random.randint(0,100,size=10)
liste = trier(liste)

print(liste)
--> array([ 6,  8, 13, 23, 61, 69, 85, 96, 99, 99])

2.c. Sélection d'un élément de rang donné

On cherche à sélectionner l'élément de rang r d'une liste. L'algorithme ressemble à celui du tri rapide. On considère une sous-liste allant des indices d à f, dans laquelle on cherche l'élément de rang r. On commence par partitionner cette sous-liste. Soit i l'indice du pivot à la fin du partitionnement. On a bien sûr $d \leq i \leq f$ . Soit L₁ la partie de la sous-liste située avant le pivot; tous ses éléments sont inférieurs ou égaux au pivot. Soit L₂ la partie de la sous-liste située après le pivot; tous ses éléments sont supérieurs au pivot. On note k la longueur de L₁, c'est à dire la position du pivot après le partitionnement, relativement au premier élément de la sous-liste.

Il y a trois cas à considérer :

Si k=r, le pivot est l'élément de rang r cherché.
Si k>r, l'élément cherché se trouve dans la sous-liste L₁, au rang r.
Si k<r, l'élément cherché se trouve dans la sous-liste L₂ au rang r-k.

Pour les deux derniers cas, la fonction de sélection est appelée récursivement, respectivement sur les listes L₁ et L₂.

Voici la fonction de sélection :

def selection(liste,debut,fin,rang): 
    if debut==fin:
        return liste[debut]
    i = partition(liste,debut,fin)
    k = i-debut+1
    if rang==k:
        return liste[i]
    elif rang < k:
        return selection(liste,debut,i-1,rang)
    else:
        return selection(liste,i+1,fin,rang-k)

et la fonction d'appel pour un tableau numpy :

def selectionner(liste,rang):
    return selection(liste.copy(),0,len(liste)-1,rang)

Exemple :

liste = numpy.random.randint(0,100,size=10)
element = selectionner(liste,2)
liste_trie = trier(liste)

print(liste)
--> array([28, 79, 32, 81, 18, 99,  0, 19, 34, 40])

print(element)
--> 18

print(liste_trie)
--> array([ 0, 18, 19, 28, 32, 34, 40, 79, 81, 99])

Enfin voici une fonction qui fournit la médiane d'une liste. Pour une liste comportant un nombre pair d'éléments, on prend la médiane inférieure.

def mediane(liste):
    rang = (len(liste)+1)/2
    return selectionner(liste,rang)

m = mediane(liste)

print(m)
--> 32

3. Filtre médian

3.a. Définition

On se limite ici à un signal à une dimension. Pour le filtrage d'une image, le principe est le même.

Soit une fenêtre de longueur l sur le signal allant des indices d à d+l. La fonction suivante renvoie la médiane du signal sur cette fenêtre :

def median_fenetre(y,debut,longueur):
    liste = y[debut:debut+longueur]
    return mediane(liste)

Pour filtrer un signal, on fait glisser la fenêtre en remplissant une pile avec les valeurs médianes :

def filtre_median(y,longueur):
    m = []
    debut = 0
    while debut+longueur < N:
        m.append(median_fenetre(y,debut,longueur))
        debut += 1
    return m

3.b. Exemple

Le filtre médian est très efficace pour enlever un bruit impulsif, qui se présente sous forme de parasites isolés, qui peuvent avoir une amplitude très grande. Il permet de supprimer le bruit sans altérer les détails de l'image.

On traite comme exemple une sinusoïde, dont quelques points sont affectés d'un bruit parasite.

import random
N=1000
x = numpy.linspace(0,5,N)
y = numpy.sin(2*numpy.pi*x)
P=50
for p in range(P):
    i = random.randrange(N)
    y[i] += (random.random()-0.5)*0.4
 
from matplotlib.pyplot import *
figure()
plot(y)

figA.pdf

Voici le résultat de l'application du filtre médian :

m = filtre_median(y,5)
figure()
plot(m)

figB.pdf

Il est intéressant de comparer ce filtre à un filtre linéaire à réponse impulsionnelle finie, par exemple un filtre gaussien, utilisé couramment pour réduire le bruit.

import math
P_gauss=2
h_gauss = numpy.zeros(2*P_gauss+1)
epsilon=0.01
sigma=P_gauss/math.sqrt(-2.0*math.log(epsilon))
som = 0.0
for k in range(2*P_gauss+1):
    h_gauss[k] = math.exp(-(k-P_gauss)**2/(2*sigma**2))
    som += h_gauss[k]
h_gauss = h_gauss/som

Voici son effet sur le signal précédent :

import scipy.signal
z = scipy.signal.convolve(y,h_gauss,mode='valid')

figure()
plot(z)

figC.pdf

Pour une fenêtre de taille identique, le filtre gaussien est très nettement moins bon.