FR 3059128 A1 20180525 FR 1661425 A 20161124 1 PROCÉDÉ DE CARACTERISATION DE L'ANISOTROPIE DE LA TEXTURE D'UNE IMAGE NUMÉRIQUE [on] L'invention concerne un procédé de caractérisation de l'anisotropie de la texture d'une image numérique. L'invention concerne également un procédé pour classer des images numériques selon l'anisotropie de leur texture. L'invention concerne enfin un support d'enregistrement d'informations et un calculateur électronique pour mettre en oeuvre ces procédés. [2] La demande W02016/042269A1 décrit un procédé qui permet d'estimer l'exposant H de Hurst de la texture d'une image et des termes [3, qui varient en fonction des caractéristiques de la texture de cette image dans une direction particulière correspondant à un angle a,. Ce procédé fonctionne très bien pour identifier l'anisotropie d'une image. [3] Il a également été proposé de construire un indice A qui caractérise 15 l'anisotropie de la texture de l'image à partir des termes 13,. Cet indice est appelé « indice d'anisotropie ». Par exemple, le calcul d'un indice A d'anisotropie à partir des termes [3, est décrit dans les articles suivants : - F.J.P Richard : « Analysis of anisotropic brownian textures and application to lesion detection in mammograms », Procedia Environmental Sciences 27 (2015) 16 - 20, 20 2015, et - F.J.P Richard : « Anisotropy indices for the characterization of brownian textures and their application for breast images », 2016. [004] Dans ces articles, l'indice A est fonction de la moyenne des termes [3,. Pour un angle a, donné, le terme [3, varie en fonction des caractéristiques de la texture 25 dans cette direction donnée mais aussi en fonction de l'exposant H de Hurst. On rappelle que l'exposant H de Hurst est une caractéristique globale de la texture qui est indépendante de l'orientation de l'image. Ainsi, quand une variation du terme [3, est observée, il n'est pas possible de savoir simplement si cette variation est due à une modification de l'anisotropie de l'image ou à une modification de la texture dans 30 son ensemble et donc de l'exposant H de Hurst. Dès lors, cet indice A varie aussi bien en fonction l'anisotropie de la texture de l'image qu'en fonction de l'exposant H de Hurst de la texture de cette image. [0051 L'invention vise à proposer un procédé de caractérisation de l'anisotropie d'une image à l'aide d'un indice d'anisotropie qui varie en fonction de l'anisotropie de 35 la texture tout en étant beaucoup moins sensible aux variations de l'exposant H de Hurst de cette même texture. Elle a donc pour objet un tel procédé comportant : a) l'acquisition d'une image numérique formée de pixels, chaque pixel étant associé à une intensité lumineuse et à une position dans l'espace Zd, où d est un entier naturel supérieur ou égal à deux ; 3059128 2 b) la transformation automatique de l'image acquise pour obtenir une image transformée I,,k, la transformation comportant l'application d'une modification -ro 4. 3059128 16 [0051] Les K-incréments de l'image T,,k(I), notés Vi,k, sont alors calculés grâce au filtre v comme suit : V [ = V[12] Z[111 T k 19] 5 p40 d où - V,,k[m] est un K-incrément calculé sur l'image T,,k(I) pour le pixel de position m, avec m un vecteur appartenant à un ensemble E qui sera défini dans ce qui suit ; - le produit T,,k.p correspond à l'application de la modification T,,k au pixel de position p 10 de l'image I et exprime les coordonnées dans Z', après application de la modification Tj,k, du pixel qui avait initialement la position p dans l'image I, - v[p] est la valeur du filtre v pour la valeur de p. [0052] Pour chaque image T,,k(I), le calcul du K-incrément est réalisé uniquement sur les pixels de l'image T,,k(I) dont les positions appartiennent à un ensemble E. 15 L'ensemble E contient seulement des positions : - qui appartiennent à l'image I, et - qui, quelle que soit la modification Tj,k appliquée, occupent une position qui existe déjà dans l'image I après application de cette modification Tj,k. [0053] En outre, pour toute position m appartenant à E, les pixels de position 20 « m - T,,k.p » occupent une position contenue dans l'image I. [0054] On note nE le nombre des positions qui appartiennent à l'ensemble E. [0055] Ainsi, le calcul des variations quadratiques est réalisé uniquement sur les points de l'image transformée pour lequel aucune interpolation n'est nécessaire. On peut ainsi avoir recours à des rotations j suivant n'importe quel angle, au contraire de 25 ce qui se passe avec les projections. En effet, si une projection est réalisée suivant une direction par exemple diagonale de l'image, des points projetés ont une position qui n'appartient pas à l'ensemble [0,NI]d . Autrement dit, ils ne font plus partie de la matrice. Il faut donc avoir recours à une interpolation pour déterminer la valeur d'intensité associée à des points qui eux appartiennent à cette matrice. Cela introduit 30 une approximation et donc une erreur. Avec les modifications -ro 4. [0062] Dans cet exemple, en dimension d = 2, le vecteur L a deux composantes, L1 et L2. Pour éliminer une tendance de degré polynomial Po=2, il faut choisir L1 et L2 tels que ILI soit égal à quatre. On choisit de préférence L1=4 et L2=0. En effet, en choisissant des valeurs pour les coordonnées du vecteur « L » suffisamment 20 éloignées les unes des autres, le filtre présente une plus grande sensibilité directionnelle. Ainsi, il réagira de façon plus marquée, donc filtrera plus efficacement, des variations qui sont orientées selon une direction particulière. Au contraire, un filtre pour lequel on choisirait L1=2 et L2=2 serait moins sensible à un signal directionnel et présenterait une efficacité moindre. 25 [0063] Dans ce mode de réalisation, le filtre v est défini par la relation suivante : (1) 1-1! si le vecteur p appartient à [0,L]x{0} et sinon v[p] = 0, où L1 appartient à 101\101. Ici, L1=4. Dans ces conditions, l'ordre du noyau v est égal à K = L1-1. [0064] Lors de cette étape 22, pour chaque valeur différente de j et k, le calculateur 14 réalise successivement les opérations suivantes : 30 - application de la modification T,,k à l'image 2 pour obtenir l'image T,,k(I), puis 3059128 18 - application du filtre v à l'image T,,k(I) pour obtenir l'image transformée [0065] Ensuite, lors d'une étape 24, pour chaque image Ij,k, on calcule la p-variation Wj,k associée à cette image Ij,k. La notion de p-variations est bien connue de l'homme du métier dans le domaine des statistiques et des probabilités. Ici, elle est calculée de 5 la façon suivante : [/n 10 [0066] Dans l'équation ci-dessus, on a utilisé le symbole « q » à la place du symbole « p », classiquement utilisé dans cette équation, pour éviter toute confusion avec le symbole « p » utilisé dans cette description pour désigner la position d'un pixel. Dans cet exemple, on utilise une forme particulière des p-variations : les « variations quadratiques » ou « 2-variations », pour lesquelles q=2. Ainsi, on calcule la variation 15 quadratique Wj,k de l'image Ij,k de la façon suivante, à partir des K-incréments calculés après filtrage lors de l'étape 22 : W,, 1 V nE meE 20 [0067] Ces variations quadratiques Wj,k contiennent des informations importantes pour l'identification de l'anisotropie. Pour extraire ces informations, on procède comme suit. [0068] Lors d'une étape 26, une analyse de covariance comportant une régression 25 statistique est effectuée sur toutes les variations Wj,k calculées pour chacune des images Ij,k afin d'estimer : - la valeur de l'exposant de Hurst H de l'image I et, - un terme [3,. [0069] La régression statistique est définie par la relation suivante : 30 log(IW,,k I) = log(lu,k I2)*H + [3, + e,,k , où : - lu,k12 est la norme euclidienne au carré du vecteur u,k; - H est l'exposant de Hurst de l'image I ; - p, est une quantité qui ne dépend pas du changement d'échelle k ; ici ce paramètre est analogue à un paramètre dit d'intercept de la régression, sauf qu'il dépend des 35 rotations j. - £,,k est un terme d'erreur de la régression dont les propriétés statistiques sont prédéterminées et fixées par l'utilisateur. Par exemple, les termes d'erreur Ej,k sont des variables aléatoires gaussiennes corrélées entre elles. 3059128 19 [0070] On rappelle que l'exposant de Hurst H est une grandeur physique indépendante des rotations de l'image. [0071] On obtient ainsi un nombre n, de termes [3,, n, étant le nombre de rotations différentes appliquées à l'image I. 5 [0072] Par exemple, si l'on s'est contenté de faire les deux rotations j1 et j2 précédemment décrites, on effectue la régression sur la base de toutes les variations quadratiques calculées pour j1 et pour j2. On obtient ainsi deux termes, [3,1 et [3,2. [0073] À ce stade, lors d'une étape 28, le calculateur 14 estime les coefficients scalaires Tm d'une fonction paire T(0) appelée fonction de topothésie asymptotique. La 10 fonction T(0) est continue sur l'intervalle [0 ; 2n]. La fonction T(0) est la fonction qui minimise le critère C suivant : (2) * où : - le symbole « T*F(0) » désigne la fonction obtenue en réalisant le produit de convolution circulaire entre les fonctions T(0) et F(0), 15 - F(0) est la fonction définie par la relation suivante : r(o) b/(PO) p où : - V est la transformée de Fourier discrète du noyau v, - H est l'exposant de Hurst de l'image acquise, - p est la variable d'intégration, [0074] Les coefficients scalaires Tm de le fonction T(0) sont définis par la relation 20 suivante : t(())=1-rnifrn(6) n=0 où : - M est un nombre entier supérieur à un, - Tm sont les coefficients scalaires de la fonction T(0), - fm(0) sont les fonctions d'une base des fonctions n-périodiques définies sur 25 l'intervalle [0 ; 2n]. [0075] Une fonction n-périodique est une fonction périodique de période n. 3059128 20 [0076] Dans ce mode de réalisation, la base de fonction n-périodique utilisée est une base de Fourier. Par conséquent, la fonction T(e) est ici définie par la relation suivante : rn,.1 où To, Ti,m et T2,m sont les coefficients scalaires de la fonction T(e), 5 [0077] Le nombre M est une constante prédéfinie, par exemple par l'utilisateur. Généralement, ce nombre M est inférieur au nombre n, d'angles a, différents. Typiquement, ce nombre M est également supérieur ou égale à 2 ou 4. Typiquement, le nombre M est choisi de manière à ce que le nombre de coefficients scalaires de la fonction T(e) soit compris dans l'intervalle [0,35n, ; 0,75n,] ou dans l'intervalle [0,45n, ; 10 0,55n,]. [0078] Dans ce mode de réalisation et dans le contexte précisé ci-dessus, il a été possible d'établir une relation linéaire entre une approximation T* des coefficients de la fonction T(e) et les termes [3, estimés. Cette relation est la suivante : (3) -C=(LT L-Fk R)i1LTI où : 15 - T* est le vecteur (To*, T2,1*, T1,2*, T2,2*, , T2,M-1*, T1,M*, T2,M*)r, les coefficients To*, T1,1*, T2,1*, T1,2*, T2,2*, T1,M-1*, T2,M-1*, T1,M*, T2,M* étant les estimations, respectivement, des coefficients To, T1,1, T2,1, T1,2, - T 2,2, - - - , T1,M-1, T2,M-1, T1,M, T2,M, - L est une matrice de dimension n, x (2M+1) dont la k-ième colonne est définie par la relation suivante : (00},M1cos(ak),11[1]sin(clk),...,p ]cos(Mak),µ[M]sin(Mak)) où il [0], cl, [1], Il[M] sont les coefficients de la transformée de Fourier discrète de la fonction p(8) suivante : tt 0)= cos \ 0 , où H est l'exposant de Hurst de l'image acquise, 20 - À est un paramètre prédéterminé, - R est une matrice diagonale de dimension (2M+1) x (2M+1) dont les coefficients sur la diagonale sont, dans l'ordre : 0, 2, 2, 5, 5, ..., (1+M2), (1+M2), - le symbole «T » désigne l'opération transposée, - a est le vecteur ([31, [32, ..., [3n,_1, [3n,)T. 3059128 21 [0079] Les inventeurs ont également établi que la valeur optimale du paramètre À est égale ou très proche d'une valeur À". La valeur À" est obtenue à l'aide de la relation suivante : (4) K trace(V (p)) (1+ M 2) 012 où : 5 - K est égal à v+/v_, et v_ étant, respectivement, la plus grande et la plus petite valeur propre de la matrice LTL, - trace(X) est la fonction qui retourne la somme des coefficients diagonaux d'une matrice carrée X, - V(a) est la matrice de covariance de l'estimateur du vecteur Z, 10 - la12 est la norme euclidienne au carré du vecteur 5. [00801 Ainsi, dans ce mode de réalisation, lors de l'étape 28, le calculateur 14 calcule la valeur À" à l'aide de la relation ci-dessus. Ensuite, il choisit la valeur du paramètre À proche de la valeur À". Par exemple, il choisit, par exemple aléatoirement, la valeur du paramètre À dans l'intervalle [0 ; 1,32C] ou [0 ; 1,1À*]. Le plus souvent, la 15 valeur du paramètre À est choisie dans l'intervalle [0,7À*; 1,3À1 ou [0,9À*; 1,14 Ici, le paramètre À est systématiquement choisi égal à la valeur À". Enfin, le calculateur estime les coefficients To, TZ,m, T2,m à l'aide de la relation (3). [00811 La relation (4) a pu être établie en recherchant l'expression numérique des coefficients To*, T2,m* qui minimise non pas directement le critère C mais un 20 critère pénalisé CA. Ce critère pénalisé CA est par exemple le suivant : (5) CA= C- [0082] Le terme qui vient se soustraire au critère C dans la relation (4) est connue sous l'expression « pénalité ». [0083] De plus, dans le cas où le filtre v est celui défini par la relation (1), alors le produit de convolution circulaire T*F(0) peut s'écrire sous la forme suivante : (6) T F = * {, ') 25 où : - p(e) est défini par la relation suivante : 3059128 22 (-)Hco )J - y est une constante indépendante de la valeur 0, - le symbole « * » est le produit de convolution circulaire de sorte que T*p(0) désigne la fonction résultant du produit de convolution circulaire entre les fonctions T(0) et p(e). 5 [0084] La constante y est définie par la relation suivante : où les symboles utilisés dans la relation ci-dessus on déjà été définis précédemment et 1...12 est la norme euclidienne au carré. [0085] Dès lors, le critère approché CA s'écrit sous la forme suivante : CA 2+ T I? T où les symboles L, T, a, À et R ont précédemment été définis et 1...12 est la norme 10 euclidienne au carré. Le vecteur T* est la solution numérique du critère CA. [0086] Pour une valeur donnée de l'angle 0 comprise entre 0 et n, la valeur de la fonction T(0) dépend des caractéristiques de la texture de l'image dans la direction 0. Cette valeur est indépendante de la valeur de l'exposant H de Hurst. Dès lors, la dispersion statistique des valeurs de la fonction T(0) pour 0 variant entre 0 et n est 15 représentative de l'anisotropie de la texture de l'image. Par exemple, la dispersion statistique de la fonction T(0) est représentée par un indice A fonction de la somme des écarts suivants : T - Mt I Lp, pour 0 variant entre 0 et n, où : - T est la fonction T(0), - NA, est une valeur moyenne des valeurs de la fonction T(0) pour 0 variant de 0 à n ou 20 une approximation de cette moyenne, et - 1...1 Lp est la norme L1 si Lp est égal à 1, la norme L2 si Lp est égal à 2 et ainsi de suite, Lp étant strictement supérieur à zéro. [0087] Ainsi, la dispersion statistique de la fonction T(0) peut être la variance ou l'écart type des valeurs de cette fonction sur [0 ; n]. Par exemple, ici, Lp est choisi 25 égal à 2 et l'indice A calculé est égal à la racine carrée de la somme définie ci-dessus. Dès lors, plus la valeur de l'indice A est grande, plus l'anisotropie de l'image est grande. Dans ces conditions, l'indice A est défini par la relation suivante : 3059128 23 (7) [00881 Ici, lors d'une étape 30, le calculateur 14 calcule l'indice A à l'aide de la formule suivante qui correspond à la relation (7) : [0089] L'ordinogramme de la figure 4 décrit un exemple d'application du procédé ci- dessus pour classer automatiquement des images 2 les unes par rapport aux autres 5 en fonction de leur texture. Cet exemple est donné dans le cas particulier où les images sont des photographies de feuilles de papier prises à l'aide d'un microscope et sous une illumination rasante. De telles prises de vues sont illustrées sur les figures 5A à 5F. Les figures 5A à 5C représentent des images de feuilles de papier brillant de trois fabricants différents. Le figures 5D et 5E représentent des images de 10 feuilles de papier satiné. La figure 5F représente une image d'une feuille de papier mate. Les bases de données contenant ces images sont décrites en détails dans les articles suivants : - R. Johnson, P. Messier, W. Sethares, et al : « Pursuing automated classification of historic photographic papers from raking light images », J. AM. Inst. Conserv. 15 53(3):159-170, 2014, et - P. Messier, R. Johnson, H. Wilhelm, W. Sethares, A. Klein, et Al : « Automated surface texture classification of inkjet and photographic media », In NIP & Digital Fabrication Conference, pages 85-91, Society for Imaging Science and Technology, 2013. 20 [0090] Lors d'une étape 40, une pluralité d'images numériques 2 sont automatiquement acquises. Parmi les images acquises, certaines correspondent à du papier brillant, d'autres à du papier satiné et d'autres à du papier mate. [0091] Lors d'une étape 42, pour chacune d'entre elles, l'indice A d'anisotropie et l'exposant H de Hurst sont calculés en mettant en oeuvre le procédé de la figure 3. 25 Dans ce cas particulier, le procédé de la figure 3 a été mis en oeuvre avec M=23 et n, = 96. [0092] Lors d'une étape 44, les images acquises sont classées automatiquement les unes par rapport aux autres en fonction de leur indice A et de leur exposant H calculés lors de l'étape 42. Cette classification est par exemple réalisée au moyen 3059128 24 d'un classificateur (« classifier » en langue anglaise) tel qu'un algorithme de classement à base de réseaux de neurones ou une machine à vecteurs de support (« support vector machine » en langue anglaise). [0093] Le graphe de la figure 6 représente pour chaque image un point de 5 coordonnées (H, A), où H et A sont, respectivement, l'exposant de Hurst et l'indice d'anisotropie calculés pour cette image. Ici, la fonction T(0) a été normalisé. Sur ce graphe, les points représentés par des croix, des ronds et des losanges correspondent à des images d'un papier, respectivement, brillant, satiné et mate. Ce graphe montre que la combinaison de l'exposant H de Hurst et de l'indice A 10 d'anisotropie permet de distinguer efficacement les différents types de papier entre eux. Ici, tous les papiers brillants, satinés et mates sont dans des zones très différentes. Ces zones sont entourées sur le graphe de la figure 6. De plus, à l'intérieur d'une même zone, un amas de points tous regroupés dans une zone encore plus étroite correspond souvent à un fabricant particulier ou à un tirage 15 particulier. Ainsi, la classification peut non seulement distinguer les différents types de papier mais aussi, pour un même type de papier, des fabricants différents ou des tirages différents. [0094] Variantes de l'image : [0095] Les pixels de l'image 2 peuvent présenter d'autres valeurs d'intensités. La 20 valeur d'intensité de chaque pixel peut être une valeur réelle. Ou bien elle peut être supérieure à 256. Par exemple, l'image 2 est encodée en couleurs. Dans ce cas, l'image en couleurs est séparée en une pluralité d'images monochromes correspondant chacune à des canaux colorimétriques qui composent l'image en couleurs. Le procédé est alors appliqué séparément pour chacune de ces images 25 monochromes. [0096] Limage 2 peut présenter une forme non carrée. Par exemple, en dimension d = 2, l'image présente une forme rectangulaire ou même trapézoïdale. Lorsque l'image ne présente pas une forme régulière, les notions de direction « horizontale » et « verticale » sont remplacées par des directions de références adaptées à la 30 géométrie de l'image. Par exemple, dans le cas d'une image de forme triangulaire, on prendra comme référence la base et la hauteur du triangle. [0097] La dimension d des images peut être supérieure à deux. Par exemple, l'image 2 peut être un hypercube de dimension d. [0098] Limage 2 peut être autre chose qu'un cliché de mammographie ou d'une 35 feuille de papier. Par exemple, il peut s'agir d'un cliché d'un tissu osseux. L'anisotropie de la texture de l'image renseigne alors sur la présence de pathologies osseuses, comme l'ostéoporose. D'autres domaines d'application plus larges peuvent être envisagés, comme d'autres types de tissus biologiques, des images aériennes ou satellitaires, des images géologiques, ou des clichés de matériaux. De façon 40 générale, la méthode s'applique à n'importe quel type d'image irrégulière et texturée 3059128 25 tel qu'une image obtenue à partir d'un appareil électronique quelconque de prise d'images. [0099] Variante du procédé d'estimation des termes (3, et H : [00loc)] D'autres modifications T,,k peuvent être utilisées. Par exemple, en dimension 5 d = 3, les modifications Tbk réalisent une rotation j autour d'un axe de rotation donné et un changement d'échelle k selon une direction donnée de l'image. Par exemple, les modifications suivantes peuvent être utilisées : k= YI, - - rcosf (1.1 -sin(a. cos (a» 0 Yk sin 0 T ,k= yk cos 0 ) 0 Yk 0 sin (ai, cos( Yk 0 T = 0 cos -sin , 0 sin cos 10 Les modifications Tj,k ci-dessus réalisent une rotation autour d'un axe et un changement d'échelle dans une direction non-parallèle à cet axe. [00101] Les valeurs de l'angle a, peuvent être différentes. De préférence, on choisit des valeurs de l'angle a, qui ne nécessitent pas d'interpolations. Toutefois, il est aussi possible de choisir des valeurs de l'angle a, qui nécessitent une interpolations des 15 pixels de l'image transformée pour retrouver les valeurs associées à chaque position p comprise dans l'ensemble E. [00102] En variante, la rotation et le changement d'échelle ne sont pas appliqués en même temps. [00103] D'autres filtres v peuvent être utilisés pour calculer les K-incréments. Par 20 exemple, pour calculer les K-incréments on peut également utiliser tout filtre dont le noyau de convolution est défini comme étant égal au produit de convolution : - d'un noyau vl de convolution quelconque, et - d'un noyau v2 de convolution égal au noyau v précédemment décrit. 3059128 26 Dans le cas particulier où le noyau v1 est une matrice identité, on retrouve le filtre précédemment décrit. A l'inverse, choisir un noyau v1 différent de la matrice identité permet de construire un grand nombre de filtres différents des filtres v précédemment décrits mais qui conviennent tous pour calculer les K-incréments. 5 [00104] Le filtrage peut être implémenté différemment lors de l'étape 22. En particulier, la transformation et le filtrage ne sont pas forcément appliqués simultanément, mais dans des formules séparées. [00105] En variante, toutes les transformations T,,k sont d'abord appliquées sur l'image I, puis, dans un second temps, les filtres sont appliqués sur chacune des images 10 T,,k(I). [00106] La valeur de K peut être différente. En particulier, avec le filtre v choisi dans l'exemple, lorsque l'image I ne présente pas de tendance, c'est-à-dire que M = 0, alors on prendra préférentiellement ILI = 2 ou, si d>4, ILI = 1+c1/4. [00107] En variante, plusieurs filtres y, sont appliqués sur chaque image T,,k(I) lors de 15 l'étape 22. On note n, le nombre de filtres y, différents appliqués sur une image donnée T,,k(I). Dans ce cas, on note I,,,,k l'image transformée obtenue en appliquant le filtre y, sur l'image T,,k(I) et V,,,,k[m] le K-incrément de cette image à la position m dans cette image, où « i » est un indice qui identifie de façon unique le filtre y, appliqué. Cet indice « i » est ici distinct de l'indice « i » utilisé précédemment comme variable 20 muette notamment en référence à la dérivée partielle du polynôme Qv(z). De fait, on peut alors calculer plusieurs variations quadratiques pour chaque image T,,k(I), une pour chaque filtre y, appliqué à cette image T,,k(I). On note ainsi W,,,,k la variation quadratique calculée pour l'image 1,4k. polos] À cet effet, l'étape 22 comporte une opération de sélection des filtres v,, par 25 exemple parmi une bibliothèque de filtres prédéfinie. [00109] Les variations quadratiques W,,,,k sont alors calculées lors de l'étape 24 a l'aide de la relation suivante : 30 ponol Lors de l'étape 26, la régression est alors effectuée de la manière suivante en tenant compte des ni filtres appliqués : I) = log(14 I2)*H + [3ij + eilk, où : - pij est le terme [3; associé au filtre vi ; - Ei,j,k est le terme d'erreur e,,k associé au filtre vi. 35 [00111] On obtient ainsi un nombre nb de termes [3;,,, où nb = n,*ni, n, étant le nombre de rotations différentes appliquées à l'image I. Lors de l'étape 28, le procédé de la figure 3 est mis en oeuvre pour chaque filtre vi. Ainsi, on obtient les coefficients scalaires de i fonctions Ti(e). L'indice A d'anisotropie est alors calculé à partir des coefficients approximés pour chacune des ces fonctions Ti(e). Par exemple, dans un 40 mode de réalisation simplifié, un indice A; d'anisotropie est calculé comme décrit 3059128 27 précédemment pour chacune des fonctions T,(0). Ensuite, l'indice A calculé est la moyenne de ces indices A,. [00112] Lorsque plusieurs filtres y, sont utilisés, le nombre de filtres y, appliqués peut varier d'une image T,,k(I) à l'autre, à condition toutefois qu'à un filtre i correspondent 5 au moins deux rotations j et, pour chacune de ces rotations j, au moins deux changements d'échelle k. [00113] Variantes de l'estimation de T* : [00114] La pénalité utilisée dans le critère CA peut être différente. Tant que la pénalité est une fonction différentiable, alors il est possible de déterminer une relation 10 linéaire, telle que la relation (3), qui exprime directement l'estimation des coefficients Tm en fonction des termes [3,. En particulier, il est possible de trouver une telle relation linéaire quel que soit le filtre v et la base de fonctions fm(0) n-périodiques utilisés. Le filtre v peut donc être différent de celui défini par la relation (1) et la base utilisée peut aussi être différente de la base de Fourier. Lorsque le filtre v est différent de celui 15 défini par la relation (1) ou lorsque la base est différente de la base de Fourier, la relation linéaire est différente de celle définie par la relation (3). [00115] La pénalité utilisée dans le critère CA peut aussi être une fonction non différentiable. Dans ce cas, il peut être difficile voire impossible d'établir une relation linéaire entre l'estimation des coefficients Tm et les termes [3,. Par exemple, la pénalité 20 peut utiliser une norme L1 de la fonction T(e) qui est non différentiable. Dans ce cas, d'autres méthodes sont possibles pour approximer les coefficients de la fonction T(e) qui minimise ce critère pénalisé. Par exemple, l'estimation des coefficients To, T1,m, T2,m qui minimisent le critère CA sont estimées en exécutant un algorithme connu de minimisation d'un tel critère tel que l'algorithme ISTA (« Iterative Shrinkage- 25 Thresholding Algorithm ») ou FISTA («Fast Iterative Shrinkage-Thresholding Algorithm »). [00116] La variante décrite ici permet d'estimer les valeurs des coefficients To, TZ,m, T2,m qui minimisent le critère CA sans disposer pour cela d'une relation numérique linéaire, comme la relation (3), qui permet d'obtenir directement une estimation de ces 30 coefficients To, TZ,m, Tzrn à partir des valeurs des termes [3,. [00117] Le mode de réalisation précédent a été décrit dans le cas particulier où la fonction T(e) est décomposée sur une base de Fourier de fonction n-périodique. Toutefois, le procédé décrit précédemment fonctionne aussi si la fonction T(e) est décomposée sur n'importe quelle autre base de fonction fm n-périodique où chaque 35 fonction fm est définie sur l'intervalle [0 ; 2n]. En particulier, il n'est pas nécessaire que la base des fonctions fm n-périodique soit une base orthogonale. Ainsi, dans la base des fonctions fm, la fonction T(e) est définie la par la relation suivante : 3059128 28 T (01- f.(0) où les Tm sont les coefficients scalaires de la fonction T(0). Par exemple, en variante, les fonctions fm sont des fonctions n-périodiques constantes par morceaux sur [0 ; n]. Une fonction constante par morceaux est une fonction qui prend des valeurs constantes sur plusieurs sous-intervalles immédiatement successifs et compris entre 5 [0 ; n]. [001181 La méthode de minimisation du critère C ou du critère CA au moyen d'algorithmes connus de minimisation d'un tel critère peut être mise en oeuvre quelle que soit la forme de la fonction fm retenue. [00119] En variante, le nombre M peut être supérieur ou égal au nombre n,. 10 [00120] Variantes du calcul de l'indice A : [00121] En variante, l'indice A est calculée uniquement pour les angles 0 égaux aux angles a, et non pas pour toutes les valeurs de 0 comprises entre 0 et n. Dans ce cas, par exemple, l'indice A est uniquement fonction de la somme des écarts suivants : 15 [00122] La classification peut être réalisée différemment lors de l'étape 42. Par exemple, l'ordre de classement des images peut être choisi différemment. REVENDICATIONS1. Procédé de caractérisation de l'anisotropie de la texture d'une image numérique, comportant : a) l'acquisition (20) d'une image numérique formée de pixels, chaque pixel étant associé à une intensité lumineuse et à une position dans l'espace Z', où d est un entier naturel supérieur ou égal à deux ; b) la transformation automatique (22) de l'image acquise pour obtenir une image transformée I,,k, la transformation comportant l'application d'une modification -ro 2. Procédé selon la revendication 1, dans lequel le noyau v utilisé lors de l'étape b) 15 est égal au produit de convolution d'un noyau de convolution quelconque avec le noyau défini de la façon suivante : 41)1= (-1) p Ll xn 1-=-1" p,k(1, p/). si le vecteur p appartient à [0,L]d et v[p] = 0 sinon, où les termes CPL désignent des coefficients binomiaux, la constante K étant alors égale à K = ILI -1. 20 3. Procédé selon la revendication 2, dans lequel : - les fonctions fm sont les fonctions de la base de Fourier et la fonction T(0) est définie par la relation suivante : T(H)=T0 où To, Ti,m et T2,m sont les coefficients scalaires de la fonction T(0), 3059128 32 - le noyau v est définie par la relation suivante : v [Pl=(-1)'' L1! -191)! - lors de l'étape e) (28), l'estimation de ces coefficients est calculée à l'aide de la relation suivante : T*=( LT L +À, R r1,T13 où : 5 - T* est le vecteur (To*, T2,1*, T1,2*, T2,2*, , T2,M-1*, T1,M*, T2,M*)r, les coefficients To*, T1,1*, T2,1*, T1,2*, T2,2*, - - - , T1,M-1*, T2,M-1*, T1,M*, T2,M* étant les estimations, respectivement, des coefficients To, T1,1, T2,1, T1,2, - T 2,2, - - - , T1,M-1, T2,M-1, T1,M, T2,M, - L est la matrice de dimension n, x (2M+1) dont la k-ième colonne est définie par la relation suivante : 41[0],M1-_ cos (ak) [ cos (M ak),Làs[M]sin(Mak)) où 0, [0], ct[1],..., Çi [144] sont les coefficients de la transformée de Fourier discrète de la fonction p(8) suivante : cos( oi , où H est l'exposant de Hurst de l'image acquise, 10 - À est un paramètre prédéterminé, - R est une matrice diagonale de dimension (2M+1) x (2M+1) dont les coefficients sur la diagonale sont, dans l'ordre : 0, 2, 2, 5, 5, ..., (1+M2), (1+M2), - le symbole «T » désigne l'opération transposée, - a est le vecteur ([31, [32, ..., [3n,)T. 15 4. Procédé selon la revendication 3, dans lequel le procédé comporte : - le calcul (28) d'une valeur À* définie par la relation suivante : K trace(Ve (1+M 2) 1012 où : - K est égal à v+/v_, v+ et v_ étant, respectivement, la plus grande et la plus 20 petite valeur propre de la matrice LTL, - trace(X) est la fonction qui retourne la somme des coefficients diagonaux d'une matrice carrée X, - V((3) est la matrice de covariance du vecteur 5, - 1512 est la norme euclidienne au carré du vecteur a, et 3059128 33 - le choix (28) automatique du paramètre À dans l'intervalle [0 ; 1,3À*]. 5. Procédé selon l'une quelconque des revendications précédentes, dans lequel l'indice d'anisotropie est calculé (30) à partir de la somme des écarts suivants : T(Ci M, Lp 5 où : - NA, est une estimation d'une valeur moyenne de la fonction T(0) pour 0 variant entre 0 et n, - 1...14 est la norme L1 si Lp est égal à 1, la norme L2 si Lp est égal à deux et ainsi de suite, Lp étant strictement supérieur à zéro. 10 6. Procédé selon la revendication 5, dans lequel l'indice A d'anisotropie est calculé (30) à l'aide de la relation suivante : A= 7. Procédé selon l'une quelconque des revendications précédentes, dans lequel : 15 - le vecteur u,k est un vecteur de Z2\{(0,0)}, - la modification T,,k présente : - pour d = 2, la forme matricielle suivante : -cos( sin sin i) cos u - et, pour d = 3, l'une des formes matricielles suivantes ou une composition de ces formes matricielles : cos (cc -sin( r. sin (". 1 cos 0 0 Yk 3059128 l'image transformée étant obtenue en multipliant les coordonnées de la position de chaque pixel par la matrice Tj,k, et - le calcul des K-incréments est réalisé à partir des seuls pixels de l'image qui occupent une position m appartenant à un ensemble E, cet ensemble E comportant 5 seulement des positions m qui existent déjà dans l'image I et qui, quelle que soit la modification Tj,k, après application de cette modification T,,k, occupent une position qui existe aussi déjà dans l'image I et pour lesquels la position « m- T,,k.p » occupe une position qui existe aussi déjà dans l'image I ; l'image transformée I,,k obtenue à l'issue du calcul de ce K-incrément comportant uniquement des pixels dont les positions 10 appartiennent à l'ensemble E. 8. Procédé selon la revendication 7, dans lequel les p-variations calculées lors de l'étape c) sont des variations quadratiques calculées selon la formule suivante : 15 ity où q=2 et nE le nombre des positions qui appartiennent à l'ensemble E. 9. Procédé de classement automatique d'images numériques en fonction de 20 l'anisotropie de leur texture, ce procédé comportant l'acquisition (40) d'une pluralité d'images formées chacune d'une pluralité de pixels ; caractérisé en ce qu'il comporte : - le calcul automatique (42), pour chacune des images acquises, d'un indice d'anisotropie respectif au moyen d'un procédé conforme à l'une quelconque des 25 revendications précédentes, et - le classement (44) des images numériques acquises, à l'aide d'un classificateur automatique, en fonction de l'indice d'anisotropie calculé pour chacune desdites images. 34 T j,k k Yk 0 T = kJ 0 cos fa -sin 0 cos Y 3059128 10. Support (16) d'enregistrement d'informations, caractérisé en ce qu'il comporte des instructions pour la réalisation d'un procédé conforme à l'une quelconque des revendications précédentes, lorsque ces instructions sont exécutées par un 5 calculateur électronique. 11. Calculateur électronique (14) pour la mise en oeuvre de l'une quelconque des revendications 1 à 9, ce calculateur étant programmé pour exécuter les étapes suivantes : 10 a) l'acquisition (20) d'une image numérique formée de pixels, chaque pixel étant associé à une intensité lumineuse et à une position dans l'espace Zd, où d est un entier naturel supérieur ou égal à deux ; b) la transformation automatique (22) de l'image acquise pour obtenir une image transformée I,,k, la transformation comportant l'application d'une modification T,,k de 15 l'image qui fait tourner d'un angle a, chaque pixel de l'image acquise d'une position à une autre autour d'un point ou d'un axe et qui agrandit ou réduit l'image d'un facteur yk, où a,= arg(u,k) et yk = lu,k12, u,k étant un vecteur qui caractérise complètement la modification T,,k, les indices j et k identifiant respectivement et de façon unique l'angle a, et le facteur yk , 20 puis, pour chaque image sur laquelle une modification T,,k a été appliquée, le calcul d'un K-incrément V,,k[m] pour chaque pixel de position m de l'image transformée, ce K-incrément étant calculé par application d'un noyau de convolution v au moyen de la formule suivante : 25 V, [i = v[p]- Z[m où : - le produit T,,k.p correspond à l'application de la modification T,,k au pixel qui avait initialement la position p dans l'image I ; 30 - le noyau de convolution v réalise un filtrage linéaire et possède un polynôme caractéristique Qv(z) et un support fini [0,L]d, v[p] étant la valeur du noyau v de convolution pour la position p, le polynôme caractéristique Qv(z) étant défini par la formule suivante : 35 Vz E -9,( ) = v[P]z' pl id et satisfaisant la condition suivante : 3059128 36 VaG 0,K tel que alors où - L est un vecteur acquis de [0, Nid qui paramétrise le noyau v, - N est un vecteur appartenant à %d qui code la taille de l'image et dont les composantes sont des entiers naturels strictement positifs ; 5 - la constante K est un entier naturel non nul acquis; - z un vecteur de composantes z1, z2, zd ; - zP désigne le monôme z1P1* z2P2*.-.*zdPd ; - apiQv/aziai...azdad est la dérivée partielle du polynôme Qv(z) par rapport aux composantes du vecteur z, le symbole az,a' indiquant une différentiation du polynôme 10 Qv(z) d'ordre a, par rapport à la variable z,, où z, désigne la i-ième composante du vecteur z et a, la i-ième composante du vecteur a, i étant un indice entier supérieur ou égal à 0 et inférieur ou égal à d ; l'étape b) étant exécutée avec n, angles a, différents et, pour chaque angle avec au moins deux facteurs yk différents, n, étant un entier supérieur ou égal à deux de 15 manière à obtenir au moins quatre images transformées I,,k différentes ; c) pour chaque image transformée I,,k différente, le calcul (24) d'une p-variation W,,k de cette image transformée à partir desdits K-incréments calculés ; d) l'estimation (26) des termes , de la régression statistique suivante : log(IW,,k I) = log(lu,k I2)*H + [3 p , + Ebk , OÙ 20 - H est l'exposant de Hurst de l'image acquise ; - Ej,k est un terme d'erreur de la régression dont les propriétés statistiques sont prédéterminées ; caractérisé en ce que le calculateur est également programmé pour exécuter les étapes suivantes : 25 e) l'estimation (28) des coefficients scalaires Tm d'une fonction paire T(0) définie sur [0 ; 2n] qui minimise le critère C suivant : = - T * où : - p, sont les termes estimés lors de l'étape d), 3059128 37 - T(0) est la fonction définie par la relation suivante pour tout angle 0 appartenant à [0 ; 2n] : r(e)= (0 m=a où : - M est un nombre entier supérieur à un acquis et constant, 5 - Tm sont les coefficients scalaires de la fonction T(0), - fm(0) sont les fonctions d'une base des fonctions n-périodiques définies sur l'intervalle [0 ; 2n], - F(0) est la fonction définie par la relation suivante : r(o où : - v est la transformée de Fourier discrète du noyau v, 10 - H est l'exposant de Hurst de l'image acquise, - p est la variable d'intégration, - le symbole « * » désigne le produit de convolution circulaire entre les fonctions T(0) et F(0), f) puis, le calcul (30), en fonction de l'estimation des coefficients scalaires Tm, d'un 15 indice d'anisotropie qui caractérise l'anisotropie de l'image, cette indice variant de façon monotone en fonction de la dispersion statistique des valeurs de la fonction T(0) pour 0 variant entre 0 et n.