Brevet FR2766598A1 — publié le 19990129

Procédé de filtrage d'images numériques par convolution avec un noyau déterminé par la méthode de Savitzky et Golay La présente invention concerne un procédé de filtrage d'images numériques par convolution avec un noyau réalisant une approximation locale polynomiale de l'image. Plus précisément, cette approximation est solution de la méthode des moindres carrés. L'invention trouve notamment une application dans le domaine des images numériques issues de microscopes électroniques à haute résolution. 1111 ' 11111 De façon connue, une image numérique est constituée d'un tableau à deux dimensions de points ou pixels (pour picture element en anglais). On peut distinguer deux sortes bien différentes d'images numériques: les images de synthèse et les images issues d'un capteur (microscope électronique, scanner, caméra, etc.). L'invention est plus particulièrement relative aux images de la seconde espèce. En effet, les capteurs n'étant pas parfaits, les images qui en sont issues présentent un bruit dont la nature dépend du capteur. Dans le cas du microscope électronique susmentionné, une des composantes du bruit est dite fractale. Il s'agit d'un des types de bruits les plus difficiles à supprimer par les techniques connues de traitement d'images. D'une façon générale, la majorité des techniques récentes de filtrage d'images numériques se basent sur une approche fréquentielle du filtrage. De façon toujours générale, ces techniques consistent à effectuer une transformée de type transformée de Fourier, puis à réaliser un traitement dans le domaine fréquentiel, typiquement un filtrage passe- bande (par exemple une apodisation), et enfin a effectuer une transformée inverse afin de revenir dans l'espace initial. Classiquement, le traitement réalisé dans le domaine fréquentiel consiste à séparer le spectre fréquentiel de l'image de celui du bruit afin d'éliminer ce dernier. Ces techniques sont actuellement très utilisées, mais elles présentent cependant de nombreux inconvénients: Premièrement, même dans les mises en oeuvre de type FFT (pour Fast Fourier Transform en anglais), une transformée de Fourier reste une opération coûteuse en temps de calcul. En effet, la complexité algorithmique qui donne l'évolution du nombre d'opérations en fonction du nombre N des données d'entrées, est égale a NxLog2(N). Il en va de même pour d'autres transformées mathématiques comme, par exemple, la transformée en cosinus ou la FHT (pour Fast Hartley Transform en anglais). De plus, les approches de type FFT nécessitent que le signal' se présente sous la forme d'une matrice carrée dont les dimensions sont des puissances de 2, ce qui est bien sûr gênant pour le traitement d'images qui peuvent a priori être de dimensions quelconques. Par ailleurs, un des effets pernicieux d'un filtrage fréquentiel basé sur une transformée de Fourier est de décaler le signal traité. En effet, par définition, apres une transformée de Fourier, on a une somme de fonctions de type sinus. Or le calcul d'un sinus sur un ordinateur implique de nombreux arrondis qui font que, d'une façon générale, une méthode basée sur une transformée de type Fourier induit une imprécision du résultat. Enfin, les techniques basées sur la FFT n'obtiennent pas de bonnes performances avec un bruit de type bruit fractal, car celui-ci n'est pas localisé dans une bande de fréquence précise qui pourrait être supprimée par un filtre de type passe-bande. De plus, les techniques d'apodisation mentionnées précédemment sont difficiles à automatiser dans le cas général. La convolution d'une image numérique avec une matrice de relativement petite taille, appelée noyau, est une autre technique connue de l'état de l'art. Dans la suite, nous allons adopter les notations suivantes. L'image avant filtrage est notée 1, et chaque point de celle-ci est désigné par l(x, y) o x et y sont les coordonnées du point. De la même façon, I'image après filtrage est notée 1', et chaque point de celle-ci est noté l'(x, y). On note N le noyau de la convolution. On note 2a+1 et 2b+ 1 ses dimensions. La plupart du temps, les noyaux de convolution sont carrés et, par conséquent, a=b. Les points de l'image 1' peuvent s'exprimer par l'équation suivante: 1 a b I'(xY)=- ='I(x + i, y+ j). N(i, j) i=-a j=-b dans laquelle K est un coefficient normalisateur dont le but est de donner a l'image 1' la même dynamique (c'est-à-dire le même intervalle de valeurs utilisées) que celle de I. Par exemple, pour une image I comprenant 256 niveaux de gris (c'est-à-dire pour laquelle la couleur est codée sur un octet), K doit être choisi de sorte que les valeurs de 1' se répartissent au mieux dans l'intervalle [0; 255]. L'effet de la convolution d'une image avec un noyau dépend de ce dernier. De nombreux noyaux existent, chacun ayant des propriétés et des performances très différentes. Parmi les plus connus, on peut par exemple citer les noyaux de Prewitt et de Sobel dont le but est le renforcement des contours. On peut indiquer aussi un noyau de la forme suivante: dont le but est d'effectuer un filtrage simple d'une image numérique. La convolution correspond alors à remplacer la valeur de chaque pixel par la valeur moyenne calculée sur l'ensemble des voisins de ce pixel. Ce procédé est néanmoins assez lourd dans la mesure o, toutes les valeurs de la matrice étant différentes de 0, le nombre d'opérations à effectuer pour calculer la valeur d'un pixel est important. De plus, un simple calcul de moyenne ne peut pas donner de résultats satisfaisant dans la mesure o le signal et le bruit ne sont pas différenciés lors de ce traitement. Pour avoir plus de renseignements sur les techniques de convolution d'une image numérique par un noyau, le lecteur peut se référer à l'un des multiples ouvrages traitant de l'analyse d'images numériques. On peut par exemple citer l'ouvrage Handbook of Image Processing Operators de Reinhard Klette et Piero Zamperoni aux éditions John Wiley and Sons. Le but de la présente invention est de proposer un procédé de filtrage d'une image numérique, qui ne présente pas les inconvénients des techniques fréquentielles ni de la convolution avec un noyau comme celui mentionné précédemment. À cet effet, I'invention a pour objet un procédé de filtrage d'une image numérique, consistant à effectuer une convolution entre cette image et un noyau de convolution. Elle se caractérise en ce que les coefficients dudit noyau sont déterminés, sur un ou plusieurs alignements de ce noyau, par la méthode de Savitzky et Golay, et nuls partout ailleurs. Un second but de l'invention est de proposer un procédé de dérivation d'une image numérique, à un degré quelconque, et présentant les mêmes avantages. L'invention telle que revendiquée ainsi que d'autres de ses avantages apparaîtront de façon plus claire dans la description suivante, faite en relation avec les dessins annexés. La figure 1 illustre schématiquement la méthode de Savitzky et Golay en en fixant les notations utilisées. La figure 2 montre une mise en oeuvre d'utilisation des résultats de la méthode de Savitzky et Golay pour construire un noyau selon l'invention. Les figures 3 et 4 donnent deux autres mises en oeuvre non exhaustives de dispositions des coefficients calculés pour former des noyaux de convolution selon l'invention. L'article Smoothing and Differentiation of Data by Simplified Least Squares Procedures de Abraham Savitzky et Marcel J. E. Golay, paru dans Analytical Chemistry, volume 36, numéro 8, juillet 1964 se proposait de résoudre le problème technique du filtrage d'un ensemble de données de mesure. Ces Ii valeurs de mesure se présentent sous la forme d'une suite ordonnée de valeurs, chacune de ces valeurs ayant une composante correspondant à une grandeur physique mesurée, et une composante de bruit. Pour cela, la méthode décrite consiste à convoluer l'ensemble des données de mesure par un vecteur dont les valeurs sont calculées de façon à approximer localement l'ensemble des données de mesure par un polynôme de degré choisi. L'approximation est une minimisation de l'erreur quadratique moyenne entre les valeurs des données de mesure et les valeurs du polynôme d'approximation. La figure 1 représente un signal s dont on connaît des échantillons si pris à intervalles réguliers. On veut approximer ce signal s par un polynôme p de degré n sur un intervalle [-m; m] contenant 2m+1 valeurs de s. On peut exprimer le polynôme p de la façon suivante: n p(i) = bnk ik k=O Selon cette notation b.k est le coefficient du terme de degré k pour un polynôme p de degré n. La variable i est une variable discrète qui varie entre -m et +m. Minimiser l'erreur quadratique moyenne entre le signal s et le polynôme p revient à chercher le minimum de l'expression suivante: i=-M dans laquelle si est l'échantillon du signal s, repéré par l'indice i. Afin de trouver le minimum de cette expression, on recherche quand s'annulent les dérivées de l'expression précédente par rapport aux différents coefficients bnk du polynôme d'approximation i=+m 2 1=1 a.z (s,-p.i)) O,: =0 Vk E[0;n] Obnk D'o: i=+m k=n E =+bmikI s) sjir = 0, avec rE[0; n] i=-M k=o ou encore: k=n i=+m i=+m bnl ik+r= si ir, avec rE[0O; n] [I] k=O i=-m i=-m À partir de cette équation, on obtient, pour chaque valeur du degré n du polynôme, un système d'équations déterminant bnk en fonction des échantillons si du signal. On arrive finalement à une expression de la forme: i + m br= oi si [Il] i=-m On considère alors que l'ensemble des valeurs ci forme le vecteur avec lequel est convolué le signal. Ce vecteur forme le résultat de la méthode de Savitzky et Golay. Il dépend de trois paramètres: le degré n du polynôme, la taille m de l'intervalle et le coefficient k. Ce coefficient k correspond au degré de dérivation. Ainsi, lorsque k=0, le vecteur formé par les ai permet de filtrer le signal. Lorsque k=l, le vecteur formé par les ai permet de filtrer le signal et de le dériver lors de la même convolution. Lorsque k=2, le résultat obtenu par convolution du signal et du vecteur formé par les a( est un signal filtré correspondant à la dérivée seconde du signal de départ. Le choix du paramètre k dépend donc de l'effet recherché par l'utilisateur de la méthode. Afin de rendre plus claire la méthode de Savitzky et Golay, on peut prendre un exemple particulier et fixer n=3, m=4 (soit un intervalle comprenant 9 échantillons) et k=3 (on cherche donc un vecteur permettant d'obtenir la dérivée troisième du signal). L'équation [I] donne: k=3 i=+4 i=+4 Yb3k +' = si ir, avec rE[O; 3] k=O i=-4,= - Sr+k Fr ce qui donne, de façon plus extensive, le système à 4 équations suivant: b30S0 + b31S1 + b32S2 + b33S3 = F0 b30S1 + b3S2 + b32S3 + b33S4 =F1 b30S2 + b31S3 + b32S4 + b33S5 = F2 b30S3 + b31S4 + b32S5 + b33S6 = F3 Dans ce système les coefficients Sr+ k sont connus et facilement calculables. On peut d'ailleurs remarquer que Sr+k=0 lorsque r+k est un nombre impair. Le système de ci-dessus peut donc se simplifier en: b S +b Sb2 -F. (équation 1) b3o0 + b32S2 - O b31S2 + b33S4 F (équation 2) b30S3 + b32S4 = F2 (équation 3) F, (équation 3) b31S4+ b33S6 =F3 (équation 4) On cherche alors à exprimer le coefficient b33 en fonction des termes Fr et Sr+k, ce qui s'obtient par les équations 2 et 4. On a alors: b33 = F3S2 -F1S4 S2S6 - S2 Finalement, en exprimant et calculant les termes Sr+k et F, on obtient l'équation suivante: b33 =-14- s + 7 -s3 + 13 s2 + 9 s +0s+ 9-s +13-s + 7-s3-14s4 a(-4 0-3 (., a,- ( a 1 a2 a3 a4 Sur cette équation, on identifie les valeurs des ai comme étant les coefficients devant les si, comme cela est défini par l'équation [Il]. Le vecteur que l'on retient pour la convolution est donc: [-14; 7; 13; 9; 0; 9; 13; 7; -14] Le procédé selon l'invention consiste à convoluer une image numérique avec un noyau dont les coefficients sont déterminés sur certains alignements par la méthode de Savitzky et Golay. La figure 2 indique la façon dont on peut utiliser le résultat de la méthode de Savitzky et Golay pour la construction d'un noyau de convolution. La référence V indique un vecteur obtenu par la méthode de Savitzky et Golay. Ce vecteur est formé de 2m+1 coefficients et peut s'écrire [a m... ao...am]. Ces coefficients sont alors disposés selon un ou plusieurs alignements dans une matrice carrée de taille 2m+1. Sur la figure 2, ces coefficients sont disposés selon deux alignements de cette matrice N (c'est-à-dire du noyau de convolution) suivant la ligne AH passant par le centre du noyau N, et suivant la colonne Av passant par le centre du noyau N. Les zones 1, 2, 3 et 4 contiennent des coefficients nuls. La figure 3 montre un deuxième exemple de construction d'un noyau à partir des coefficients calculés par la méthode de Savitzky et Golay. En plus des alignements AH et Av, les coefficients issus de la méthode de Savitzky et Golay sont disposés selon les deux diagonales AD1 et AD2 passant par le centre du noyau (c'est-à-dire les bissectrices du noyau N). Les zones 5, 6, 7, 8, 9, 10, 11 et 12 contiennent uniquement des coefficients nuls. La figure 4 montre un troisième exemple de disposition: les coefficients calculés par la méthode de Savitzky et Golay sont disposés selon les deux diagonales AD1 et AD2. Les zones 13, 14, 15 et 16 contiennent uniquement des coefficients nuls. Ces trois exemples de mise en oeuvre de l'invention ne sont bien sûr pas exhaustifs, et il est possible, pour l'homme du métier, d'imaginer d'autres dispositions des coefficients pour former un noyau de convolution. On remarque que la méthode de Savitzky et Golay donne un vecteur de valeurs ai symétrique par rapport à son centre (c'est-à-dire que pour tout ic[- m; m], a i=ca.). En raison de la méthode utilisée pour la construction des noyaux, ceux-ci héritent de cette propriété et sont, par conséquent, invariants par rotation. Autrement dit, le résultat du filtrage d'une image par le procédé selon l'invention ne dépend pas de l'orientation de celle-ci. On remarque aussi que les coefficients non nuls ne sont présents que sur certains alignements du noyau. La proportion de ces valeurs non nulles diminue lorsque la taille du noyau augmente. On a donc une majorité de 0 (c'est- à-dire une matrice creuse). Par conséquent, la convolution avec un noyau selon l'invention est très rapide, comparée avec une convolution avec un autre noyau. De plus, les valeurs des coefficients non nuls sont des valeurs entières. Or on sait que sur un ordinateur classique, les opérations sur des valeurs entières sont beaucoup plus rapides que les opérations sur des valeurs réelles. Selon l'article de A. Savitzky et M. Golay susmentionné, il est possible d'obtenir des coefficients qui au lieu de simplement filtrer (i.e. retirer le bruit) de l'image, permet d'obtenir une image dérivée de l'image de départ. Ainsi, les coefficients ai obtenus en calculant le coefficient blk du polynôme p permettent de calculer une image dérivée k fois et approximée par un polynôme de degré n. L'ensemble des coefficients obtenu en calculant b.0 permet donc d'effectuer un simple filtrage de l'image de départ. De la même façon, l'ensemble des coefficients obtenu en calculant b.1 permet d'effectuer un filtrage et une dérivation de l'image de départ. Par ailleurs, pour un même niveau de dérivation (c'est-à-dire pour une même valeur de k), on peut obtenir une pluralité de noyaux différents, en fonction du degré n du polynôme d'approximation. En reprenant les valeurs des coefficients ai qui avaient été calculés dans PI l'exemple de mise en oeuvre de la méthode de Savitzky et Golay, donné précédemment, on peut obtenir les noyaux N1, N2 et N3 suivants: O O O O -14 0 0 0 O0 St,'1 0 0 0 0 7 0 0 0 0 rI1'" 0 0 0 0 13 0 0 0 0 O O O 0 9 0 0 0 0 N = -14 7 13 9 0 9 13 7 -14 O O O 0 9 0 0 0 0 O O O 0 13 0 0 0 0 O O O 0 7 0 0 0 0 O O O O -14 0 0 0 0 -14 0 0 0 -14 0 0 0 -14 0 7 0 0 7 0 0 7 0 0 0 13 0 13 0 13 0 0 O O 0 9 9 9 0 0 0 N,= -14 7 13 9 0 9 13 7 -14 O O 0 9 9 9 0 0 0 O 0 13 0 13 0 13 0 0 0 7 0 0 7 0 0 7 0 -14 0 0 0 -14 0 0 0 -14 -14 0 0 0 0 0 0 0 -14 0 7 0 0 0 0 0 7 0 O 0 13 0 0 0 13 0 0 O O 0 9 0 9 0 0 0 Nô= 0 0 0 0 0 0 0 0 0 O 0 0 9 0 9 0 0 0 O 0 13 0 0 0 13 0 0 0 7 0 0 0 0 0 7 0 -14 0 0 0 0 0 0 0 -14 Comme tous les procédés de convolution, le procédé selon l'invention permet de favoriser les phénomènes périodiques présents dans l'image ayant une périodicité proche de la taille 2m+1 du noyau de convolution. Par exemple, dans le cas d'une image satellite d'une surface d'eau (océan, lac...), si l'on veut observer la houle, il convient de choisir la taille m du noyau en fonction de la distance entre deux crêtes de vagues successives et du rapport d'échelle (c'est-à-dire de la relation qui lie la distance entre deux pixels et la distance entre les points géographiques qu'ils représentent). REVENDICATIONS 1) Procédé de filtrage d'image numérique, consistant à effectuer une convolution entre ladite image et un noyau de convolution, caractérisé en ce que les coefficients dudit noyau sont déterminés, sur un ou plusieurs alignements dudit noyau, par la méthode de Savitzky et Golay, et nuls partout ailleurs. 2) Procédé selon la revendication précédente, caractérisé en ce que lesdits alignements comprennent la ligne dudit noyau qui passe par le centre dudit noyau, et la colonne dudit noyau qui passe par le centre dudit noyau 3) Procédé selon l'une des revendications précédentes, caractérisé en ce que lesdits alignements comprennent les deux diagonales dudit noyau qui passent par le centre dudit noyau.