Introduction au Domaine de Recherche :

Débruitage d'Images Médicales

Jean Feydy

11 octobre 2015

Résumé
Les images médicales ont, quelle que soit la modalité utilisée, vo-
cation à être fortement bruitées et détériorées car obtenues avec des
doses de radiation ou des temps d'acquisition aussi faibles que pos-
sible. En vue de traiter ces données, de nombreuses techniques ont vu
le jour : nous décrirons ici quelques méthodes historiques et populaires
dans l'industrie, avant de traiter plus spéciquement des seuillages en
bases d'ondelettes.

Table des matières

1 Introduction et motivations 1
2 Méthodes linéaires, ltrage 2
2.1 Un exemple de méthode historique : le ltrage adaptatif 2
2.2 La Transformée de Fourier Discrète . . . . . . . . . . . 5

3 Débruitage non-linéaire par seuillage d'une représenta-

tion parcimonieuse 7
3.1 Cahier des charges d'une bonne transformée . . . . . . 7
3.2 Les ondelettes, résultats de Donoho . . . . . . . . . . . 8
3.3 Les ondelettes orientables . . . . . . . . . . . . . . . . 10

4 Ouvertures et conclusion 11

1 Introduction et motivations
En matière d'imagerie médicale, qualité rime souvent avec danger.
En eet, pour acquérir images ou volumes de l'anatomie d'un patient,

1
on expose souvent celui-ci à une irradiation potentiellement nocive.
Les données cliniques seront donc obtenues avec une qualité minimale,
permettant tout juste au médecin d'observer les structures d'intérêt,
sans excès de zèle. Rendre l'image plus lisible, soulager le praticien
ou faciliter un traitement automatisé est donc un enjeu clef pour les
constructeurs : nous proposons ici un très bref apercu des approches
possibles, en nous concentrant sur l'une d'elles, le seuillage en bases
d'ondelettes.
On modélisera les détériorations par un bruit additif s : étant donné
une image propre f , on suppose que le jeu de données bruité est donnée
par :
f = f0 + s, (1)

où s est une variable aléatoire de loi connue. En pratique, la nature

du bruit dépend de la modalité d'imagerie utilisée, comme illustré Fi-
gure 1. Gaussien et blanc (i.i.d.) dans les cas les plus proches de la
physique (phase en IRM, par exemple), il peut être bien plus com-
plexe à analyser si un processus de reconstruction sophistiqué est ap-
pliqué entre les mesures des capteurs et l'achage de l'image nale : la
transformée de Radon utilisée dans les scanners CT est un cas d'école.
On se limitera ici au cas le plus simple à traiter mathématique-
ment : celui du bruit blanc gaussien où les pixels sont perturbés par
des variables aléatoires i.i.d. de loi N (0, σ 2 ).

2 Méthodes linéaires, ltrage

La manière la plus simple de débruiter une image est d'opérer en
chaque point à un moyennage local, les pondérations étant données
par un ltre. Si ce dernier reste constant sur toute l'image, on a alors
simplement aaire à un opérateur de convolution, qui peut être implé-
menté très rapidement en parallèle  de nos jours, sur carte graphique.

2.1 Un exemple de méthode historique : le l-

trage adaptatif
De manière plus sophistiquée, on peut vouloir choisir en chaque
point un ltre adaptée à la structure locale de l'image. Si, dans une
région lisse et isotrope, un ltre outant est un bon choix, conserver
une arête bien nette implique une certaine prudence dans la délimita-
tion du support de ltre : l'idée est illustrée Figure 2.

2
 (a) Fluoroscopie à forte dose. (b) Fluoroscopie à faible dose.

(c) Zoom permettant d'observer le (d) Tranche bruitée d'une acquisi-

motif en grille du bruit sur les uoro- tion par scanner CT. La reconstruc-
scopies - dû à la répartition des cap- tion par inversion d'une transformée
teurs. de Radon induit un bruit en lignes
parallèles dûes aux projections dans
toutes les directions.

Figure 1  Quelques exemples de données médicales. L'objectif de notre

travail est de retrouver une qualité d'image proche de (a) à partir de (b),
acquise avec une plus faible dose en rayons X, et donc moins de risques pour
le patient. Comme on le voit sur (c) et (d), la modélisation de la détérioration
par un bruit blanc gaussien additif (non structuré) n'est dans bien des cas
pratiques qu'une première approximation.

3
 (a) Flot de contrôle.

(b) À l'aide d'un descripteur de structure (tenseur), un ltre adapté est choisi en chaque
position.

Figure 2  Illustrations schématisant l'action d'un pipeline de ltrage adap-

tatif typique.

4
2.2 La Transformée de Fourier Discrète
Pour comprendre ces ltres et autres opérateurs de convolution, il
est de bon ton d'introduire la transformée de Fourier d'une image :
on se souvient en eet que celle-ci diagonalise les convolutions dans
Rd . Attention toutefois : nous travaillons ici dans un cadre discret, et
nos images sont dénies sur un domaine ni [[1, l1 ]] × · · · ×[[1, ld ]] ce qui
impose de prendre quelques précautions.
Examinons donc le problème en dimension 1  la transformée de
Fourier étant séparable, cela n'est guère limitatif. On voit Figure 3
qu'échantillonnage et périodisation sont des opérations duales par trans-
formée de Fourier. Formellement, le résultat clé est que la transformée
de Fourier du peigne de Dirac III1 est III2π , et la dualité entre produit
point-par-point et produit de convolution.
Aussi, travaillant avec des images échantillonnées sur une grille
nie, on périodisera cette dernière par commodité : les transformées
de Fourier de nos images sont alors des images de même taille.

La Transformée de Fourier Rapide Finalement, la transfor-

mée de Fourier discrète peut s'écrire en 1D comme un produit matriciel
entre un opérateur carré harmonique et le vecteur des échantillons sur
une période [[0, l − 1]] de f :

w0·0 w0·1 w0·(l−1)

    
fb(0k) ... f (0)
 fb(1k)   w1·0 w1·1 ... w1·(l−1)  f (1) 
fb =   = ·
    
. . . .. . .
. . . . .

 .   . . . .  . 
fb((l − 1)k) w(l−1)·0 w(l−1)·1 . . . w(l−1)·(l−1) f (l − 1)
(2)
avec k = 2π/l la fréquence fondamentale et w = e−ik son analogue
complexe.
Un algorithme naïf requerrait O(l2 ) opérations, mais on sait depuis
1965 que grâce à sa structure harmonique, l'opérateur (w
n·m )
n,m peut
être factorisé en produits parcimonieux, ce qui mène à un algorithme
diviser pour régner en O(l log(l)) opérations : la clé de l'utilisation
pratique de la transformée de Fourier, qui se généralise aux dimensions
supérieures par séparabilité [2].
En quelques paragraphes, nous avons présenté la transformée de
Fourier discrète (périodique) : il s'agit d'un outil et d'un exemple fon-
damental, le premier changement de représentation auquel seront com-
paré tous les autres.

5
 1 copies Fourier domain
40
1
30

20
0.5
10

0 0
-8l -7l -6l -5l -4l -3l -2l -1l 0l 1l 2l 3l 4l 5l 6l 7l 8l −16 · 2π/l −8 · 2π/l 0 · 2π/l 8 · 2π/l 16 · 2π/l

3 copies Fourier domain

300
1
200

0.5
100

0 0
-8l -7l -6l -5l -4l -3l -2l -1l 0l 1l 2l 3l 4l 5l 6l 7l 8l −16 · 2π/l −8 · 2π/l 0 · 2π/l 8 · 2π/l 16 · 2π/l

5 copies Fourier domain

800
1
600

400
0.5
200

0 0
-8l -7l -6l -5l -4l -3l -2l -1l 0l 1l 2l 3l 4l 5l 6l 7l 8l −16 · 2π/l −8 · 2π/l 0 · 2π/l 8 · 2π/l 16 · 2π/l

7 copies Fourier domain

2,000
1
1,500

1,000
0.5
500

0 0
-8l -7l -6l -5l -4l -3l -2l -1l 0l 1l 2l 3l 4l 5l 6l 7l 8l −16 · 2π/l −8 · 2π/l 0 · 2π/l 8 · 2π/l 16 · 2π/l

Figure 3  De la transformée aux séries de Fourier par périodisation. En pa-

vant le domaine spatial avec des copies d'une fonction à support compact, on
échantillonne sa transformée de Fourier sur le spectre du pavé correspondant.

6
3 Débruitage non-linéaire par seuillage
d'une représentation parcimonieuse
À l'aide de la transformée de Fourier, on peut aisément et complè-
tement décrire toutes les opérations de ltrage linéaires : si f = f0 + s
est une image bruitée, g un ltre, on a

∗ g = fb · gb.
f[ (3)

Le travail de sélection du ltre consiste alors à choisir gb(ω) en

2
fonction du rapport signal-sur-bruit entre les puissances fb0 et s|2 ,
|b
estimées a priori  dans le cas d'un bruit blanc, on a par exemple
s|2 = σ 2 .
|b
Comme illustré Figure 5, il est possible d'aller bien au delà de
ces estimées linéaires qui se limitent en général à de simples ltres
passe-bas. Notre objectif est d'utiliser notre connaissance a priori de
la structure des images propres (présence de contours abrupts délimi-
tant des zones lisses, etc.) pour discriminer le signal original du bruit,
déstructuré.

À la recherche d'une représentation pertinente Comme

on l'a vu plus haut, la transformée de Fourier permet de passer de
la traditionnelle somme de diracs à la représentation en somme d'har-
moniques, complètement délocalisées. Aucune de ces deux visions n'est
vraiment adapté à notre problème : les éléments d'intérêt qui font sens
sur une image (un coin, une arête, un dégradé, etc.) vivent en eet
dans un entre-deux, étant localisés à la fois en espace (position) et
en fréquence (netteté). Pensons, par analogie, à une piste musicale :
un chier échantillonné .wav présente une suite incompréhensible de
pics, tandis que sa transformée de Fourier mélange les instants, dis-
simulant toute l'information temporelle dans des décalages de phase
par-
inintelligibles... La bonne représentation, lisible, traitable, est la
tition, dans laquelle une musique est décrite comme une succession de
notes, aux hauteurs et durées nies, non nulles, déterminées.

3.1 Cahier des charges d'une bonne transfor-

mée
Ces vingt dernières années, un grand sujet de recherche a donc été
la quête d'une bonne représentation, d'un dictionnaire de notes
adapté aux images, naturelles ou médicales. Avant de présenter les
solutions proposées dans la littérature, listons une suite de propriétés
qui doivent être vériées pour permettre un usage pratique.

7
Banque de ltres Tout d'abord, notre représentation  ou trans-
formée  doit être facile à calculer, pour permettre un traitement en
temps réel des données. On se restreint donc aux banques de ltres
linéaires : on choisira avec soin une famille de ltres linéaires (ϕi )i∈I ,
avant de prendre pour représentation d'une image f la suite des ré-
ponses (hf, ϕi i)i∈I calculées par convolution, ou produit dans le do-
maine de Fourier.
Dans le cadre ni où nous travaillons, on peut donc voir notre
opérateur comme une matrice Φ, et la représentation de f comme un
vecteur colonne
f˜ = Φf. (4)

Complétude et tightness Après traitement, on doit pouvoir re-

construire f à partir de f˜. On demande donc que Φ possède un inverse
à gauche, et soit donc de rang maximal  avec possiblement plus de
lignes que de colonnes. Comme on souhaite de plus pouvoir calculer cet
inverse ecacement, on requiert souvent une propriété supplémentaire,
la tightness :
∃a ∈ C, Φ? Φ = aId. (5)

Normalisation et décorrelation De manière complémentaire,

on demande que les ltres ϕi , les lignes de Φ, aient une norme L2 uni-
taire. Cela peut compliquer l'obtention de la tightness, mais assure la
distribution uniforme d'un bruit blanc Gaussien sur toutes les compo-
santes de notre représentation : si f0 est une image corrompue par un
tel bruit blanc s de variance σ2, on aura

f^
0 + s = Φ(f0 + s) (6)

= Φf + Φs (7)

= fe + se, (8)

où se est lui-même un bruit Gaussien de variance σ2. Si, de plus, ses

composantes sont indépendantes, il est légitime de débruiter la re-
présentation uniformément et ligne-par-ligne, le cas le plus simple à
comprendre et implémenter ecacement. En pratique, et avec un opé-
rateur rectangulaire, une indépendance totale des lignes de se = Φs est
impossible à obtenir. On se contente donc d'une décorrelation maxi-
male, obtenue en un certain sens géométrique si Φ est tight [5].

3.2 Les ondelettes, résultats de Donoho

Parcimonie Jusqu'ici, aucune des propriétés requises n'est vrai-
ment contraignante : la transformée de Fourier discrète les vérie

8
toutes en tant que décomposition sur une base orthonormale. Pour
aller plus loin, le véritable critère est celui de la parcimonie : pour des
images f0 prises dans une classe d'intérêt (ici, celle des scans médi-
caux), on souhaite que Φf0 possède aussi peu de coecients non-nuls
que possible. Étant donnée une telle transformée parcimonieuse, on
peut alors utiliser des méthodes simples mais ecaces de débruitage
non-linéaire comme le seuillage doux proposé par Donoho [3].
L'idée sous-jacente est la suivante : supposons que notre image
propre f0 X possède une représentation parcimonieuse Φf0
dénie sur
avec au plus M << |X| coecients non-nuls : modulo une propriété
de conservation d'énergie (vériée dans le cas des tight frames ), ces
coecients auront nécessairement de fortes valeurs. De l'autre côté, et
si la transformée est bien construite, ajouter un bruit déstructuré de
variance σ2 corrompra les composantes de Φf0 uniformément, l'énergie
du bruit étant étalée sur le domaine transformé tout entier. Faire la
distinction entre information structurée et bruit est alors chose aisée,
les grands coecients provenant selon toute probabilité du signal sous-
jacent, tandis que ceux plus modestes proviennent du bruit seul.

Seuillage doux Idéalement, on aimerait projeter la représentation

bruitée Φf = Φ(f0 + s) sur le sous-espace V de dimension M engendré
par les vecteurs de base qui explique la plus grande part de l'énergie,
c'est-à-dire qui minimise l'erreur

E(V ) = kf − pV (f )k22 . (9)

En pratique, ce problème combinatoire est NP-complet, et ne saurait

être résolu ecacement. On se contente donc typiquement de trouver

p(f ) = arg minkf − xk22 + λkxk1 , (10)

x

avec λ un paramètre (le seuil). Sans rentrer dans les détails d'un énoncé
qui ne sont pas utiles en pratique  et nécessaiteraient une longue
introduction des risques associés aux estimateurs , on peut montrer
que l'estimée est obtenue en appliquant coordonnée par coordonnée
un seuillage doux  voir Figure 4 , et qu'il est quasi-optimal parmi
les estimateurs coordonnée-par-coordonnée pour une valeur du seuil
bien choisie  les théorèmes originaux sont dûs à Donoho [3], mais
la théorie s'est bien développé depuis : on consultera avec prot le
chapitre Débruitage de [6].
En utilisant cette méthode, le bruit n'est pas complètement enlevé,
car les forts coecients étaient eux aussi corrompus. Mais l'énergie
du bruit répartie sur ceux-ci n'est plus qu'un simple multiple de M σ2
2
au lieu du |X| σ initial, tout en restant comparé à l'énergie du signal

9
initial kf0 k2 : quelle que soit la méthode de seuillage choisie sur la
transformée, la perte d'énergie résultante ne saurait guère être plus
grande que celle engendrée dans le domaine spatial.

Les ondelettes Construire une transformée parcimonieuse : cela

revient à capturer dans un dictionnaire l'essence de la structure d'une
classe de données. La transformée de Fourier, par exemple, ne convient
pas du tout : les transitions nettes, les bords, engendrent en eet des
coecients importants sur toute la largeur du spectre. Ces dernières
décennies, un eort considérable a donc été fourni pour élaborer une
transformée qui décrive des bords réguliers avec aussi peu de coef-
cients que possible. Le choix de loin le plus populaire est celui des
ondelettes, qui combine simplicité, exibilité et ecacité : son principe
fondamental est de choisir les ltres ϕi comme des versions translatées
et zoomées d'une poignée de ltres générateurs w1 , . . ., wp .
En travaillant avec un ensemble d'échelles S et de positions Υ,
notre banque de ltres peut alors être écrite

p
(ϕi )i∈I = (ws,y )(p,s,y)∈[[1,P ]]×S×Υ (11)

où  
p 1 p x−y
ws,y : x 7→ w (12)
sd/2 s
est la pième ondelette à l'échelle s et position y, de norme L2 égale à
kwp k2 .
Le faible nombre p de formes de base et la structure multi-échelle
permettent de calculer une telle transformée de manière très ecace,
avec de nombreuses propriétés théoriques compilées dans [6].

L'exemple des curvelets Les images présentent des structures

orientées, et leur compréhension est cruciale pour tout processus de
débruitage ecace. Des familles d'ondelettes ont donc été spécialement
concues avec des ltres orientés, capable d'engendrer des décomposi-
tions parcimonieuses à un coup raisonnable. Les curvelets, introduite
par Candès et Donoho dans [1] et illustrées Figure 6 vérient ainsi
l'estimée h i
E kf0 − p(f )k2 = O(|log(σ)|2 σ 4/3 ) (13)

quand la variance σ2 du bruit tend vers 0, pour f une image C2 par

morceaux avec bords C 2.

3.3 Les ondelettes orientables

Plutôt que suivre la piste des curvelets en discrétisant les directions
et pavant le plan de Fourier, on peut aussi s'intéresser aux familles

10
 5

4.5

4
Thresholded coecient
3.5

2.5

1.5
Identity
1 Soft threshold
Hard threshold
0.5
Wiener threshold

0
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
Noisy coecient

Figure 4  Rampes de seuillage, T = 1. Le seuillage doux correspond à

une pénalité pro-parcimonie en kxk1 , les deux autres sont proposées à titre
d'exemple, et utilisées en pratique.

d'ondelettes construites à partir de ltres orientables, c'est à dire telles

que pour toute orientation θ, on puisse synthétiser les ltres tournés
w1 ◦Rθ , . . ., wp ◦Rθ comme combinaisons linéaires des ltres originaux
w1 , . . ., wp . La dénition est illustrée Figure 7, et l'intérêt de la notion
pour le débruitage d'images expliqué Figure 8. On développera en n
de soutenance la construction et l'utilisation de familles d'ondelettes
orientables adaptées aux images d'intérêt, en se basant sur les travaux
de Freeman et Adelson [4], Unser, Van De Ville et Chenouard [7, 8],
abordés au cours de ma thèse de M2.

4 Ouvertures et conclusion
Au cours de cette introduction, nous avons esquissé un très bref
panorama du monde du débruitage d'images médicales, nous concen-
trant sur une des approches les plus prometteuses, celle des ondelettes.
Gardons à l'esprit qu'en matière de débruitage, aucun résultat dénitif
n'est à attendre, et que la comparaison de méthodes a priori fondamen-
talement distinctes reste une part essentielle du travail du chercheur.
Aussi n'avons-nous pas trop insisté sur les théorèmes d'optimalité et
d'estimations asymptotiques qui peuvent être démontrés selon les fa-

11
Figure 5  Illustration schématique d'un algorithme de seuillage dans une
base redondante d'ondelettes séparables. Tirée de [6].

12
Figure 6  Bases de Curvelets. Le plan de Fourier est pavé et le domaine
angulaire, discrétisé. Tiré de [6].

milles, et préféré mettre l'accent sur les principes, contraintes et idées

qui guident le choix d'une transformée.
Si, dans le cadre de cet exposé, nous nous sommes limité au cas
des fonctions à valeurs réelles, soulignons aussi les problèmes posés
par les modalités d'imagerie récentes aux espaces d'arrivées plus exo-
tiques, telles les séquences d'IRM de diusion, aux valeurs tensorielles.
Le monde de l'imagerie médicale est source de nombreux dés, qu'il
importe de relever.

Références
[1] Emmanuel J Candes and David L Donoho. Recovering edges in
ill-posed inverse problems : Optimality of curvelet frames. Annal
of statistics, pages 784842, 2002.
[2] James W Cooley and John W Tukey. An algorithm for the machine
calculation of complex fourier series. Mathematics of computation,
19(90) :297301, 1965.

[3] David L Donoho. Denoising by soft-thresholding. Information

Theory, IEEE Transactions on, 41(3) :613627, 1995.
[4] William T. Freeman and Edward H Adelson. The design and use
of steerable lters. IEEE Transactions on Pattern Analysis & Ma-
chine Intelligence, (9) :891906, 1991.
[5] Jelena Kova£evi¢ and Amina Chebira. An introduction to frames.
Foundations and Trends in Signal Processing, 2(1) :194, 2008.

13
 (a) Ondelettes de Meyer-Riesz (b) Ondelettes de Simoncelli-Riesz
d'ordre N = 3 tournées à un angle d'ordreN = 3 tournées à un angle
θ = π/6. θ = 3π/4.

Figure 7  Orientabilité de deux familles d'ondelettes. Sur la première ligne,

les ltres originaux - similaires à des dérivées directionnelles d'ordre 3. Au
milieu, la matrice des coecents d'orientation : les analogues orientés (à
droite) des ltres cartésiens sont simplement obtenus en sommant ses lignes,
c'est à dire en prenant une combinaison linéaire des ltres de base.

14
(a) Avec des ondelettes séparables standards. De gauche à droite : une frontière propre, son
analogue bruité, et l'estimation débruitée. Un bruit isotrope corrompt chaque orientation
de manière similaire : le patch bruité, représenté par la croix verte, se trouve aléatoire-
ment dans un voisinage du patch original  croix bleue , avec des coecients forts sur
chaque coordonnée. Dans le cadre d'un seuillage d'ondelette, on ne peut guère faire mieux
qu'atténuer indistinctement tous ces coecients. Au nal, l'orientation a été corrompue et
l'atténuation sur chaque coordonnée nous a fait perdre beaucoup d'énergie.

(b) Avec des ondelettes orientables. De gauche à droite : une frontière propre, son analogue
bruité, et l'estimation débruitée. Une estimée extérieure (able) de l'orientation est utilisée
par le passage de la base d'ondelettes cartésienne (en noir) à sa version convenablement
orientée (en bleu). Dans cette dernière, le signal possède une représentation parcimonieuse,
tandis que le bruit isotrope reste toujours aussi épars : le seuillage se comporte de manière
bien plus intéressante, comme une projection sur l'axe privilégié. Au nal, moins d'énergie
structurée est perdue (une seule atténuation au lieu de deux), et l'orientation estimée à été
encouragée dans le résultat nal.

Figure 8  Diagrammes illustrants la manière dont les ondelettes orientables

améliorent la parcimonie tout en permettant de promouvoir une information
a priori (le champs d'orientation, estimé de manière able). Les patchs sont
représentés dans le domaine d'ondelettes par des points dont les coordonnées
sont les réponses de ltres en leurs centres.
15
[6] Stephane Mallat. A wavelet tour of signal processing : the sparse
way. Academic press, 2008.

[7] Michael Unser and Nicolas Chenouard. A unifying parametric fra-

mework for 2d steerable wavelet transforms. SIAM Journal on
Imaging Sciences, 6(1) :102135, 2013.
[8] Michael Unser, Nicolas Chenouard, and Dimitri Van De Ville. Stee-
rable pyramids and tight wavelet frames in l2(rd).Image Proces-
sing, IEEE Transactions on, 20(10) :27052721, 2011.

16