Laboratoire 5

Dans ce laboratoire, nous allons à nouveau exploiter la transformée de Fourier discrète, mais en l’appliquant cette fois à l’analyse et à la manipulation d’images.

La transformée de Fourier discrète à deux dimensions

Dans le laboratoire 4, nous avons étudié la transformée de Fourier discrète d’un ensemble d’échantillons \(x[0]\), \(x[1]\), …, \(x[N - 1]\). Cette transformée peut être généralisée de façon à pouvoir l’appliquer à un ensemble de données de dimension quelconque. Pour traiter des images, nous allons utiliser la transformée de Fourier discrète à 2 dimensions, définie de la façon suivante. La transformée des échantillons

\[\begin{split}\begin{array}{cccc} x[0, 0] & x[0, 1] & \cdots & x[0, N_2 - 1] \\ x[1, 0] & x[1, 1] & \cdots & x[1, N_2 - 1] \\ \vdots & \vdots & \ddots & \vdots \\ x[N_1 - 1, 0] & x[N_1 - 1, 1] & \cdots & x[N_1 - 1, N_2 - 1] \end{array}\end{split}\]

est donnée par

\[X[k_1, k_2] = \sum_{n_1 = 0}^{N_1 - 1} \sum_{n_2 = 0}^{N_2 - 1} x[n_1, n_2]\, e^{-2i\pi\left(\frac{k_1 n_1}{N_1}+\frac{k_2 n_2}{N_2} \right)}\]

pour \(k_1\) et \(k_2\) entiers.

Principales propriétés

Les propriétés suivantes généralisent celles de la transformée à une dimension:

1. Pour tous \(k_1\) et \(k_2\) entiers, on a \(X[k_1 + N_1, k_2] = X[k_1, k_2 + N_2] = X[k_1, k_2]\). En d’autres termes, les valeurs de la transformée sont périodiques de période \(N_1\) selon la première dimension, et de période \(N_2\) selon la seconde.

   Par conséquent, pour connaître entièrement la transformée de Fourier discrète à deux dimensions des échantillons \(x[0, 0]\), \(x[0, 1]\), …, \(x[N_1 - 1, N_2 - 1]\), il est suffisant de se limiter au tableau

   \[\begin{split}\begin{array}{cccc} X[0, 0] & X[0, 1] & \cdots & X[0, N_2 - 1] \\ X[1, 0] & X[1, 1] & \cdots & X[1, N_2 - 1] \\ \vdots & \vdots & \ddots & \vdots \\ X[N_1 - 1, 0] & X[N_1 - 1, 1] & \cdots & X[N_1 - 1, N_2 - 1]. \end{array}\end{split}\]

2. Si tous les échantillons \(x[0, 0]\), \(x[0, 1]\), …, \(x[N_1 - 1, N_2 - 1]\) sont réels, alors on a \(X[k_1, k_2] = (X[N_1 - k_1, N_2 - k_2])^*\), où \(z^*\) dénote le conjugué de \(z\).

   Cela signifie que si l’on interprète \(X[0, 0]\), \(X[0, 1]\), …, \(X[N_1 - 1, N_2 - 1]\) comme les pixels d’une image, alors celle-ci présentera une symétrie centrale.

3. La transformée de Fourier discrète inverse à deux dimensions, permettant de retrouver \(x[0, 0]\), \(x[0, 1]\), …, \(x[N_1 - 1, N_2 - 1]\) à partir de \(X[0, 0]\), \(X[0, 1]\), …, \(X[N_1 - 1, N_2 - 1]\) est donnée par:

   \[x[n_1, n_2] = \frac{1}{N_1 N_2}\sum_{k_1 = 0}^{N_1 - 1} \sum_{k_2 = 0}^{N_2 - 1} X[k_1, k_2]\, e^{2i\pi\left(\frac{k_1 n_1}{N_1}+\frac{k_2 n_2}{N_2} \right)}.\]

De cette dernière propriété, on déduit que, à une constante de proportionnalité près, \(X[k_1, k_2]\) est le coefficient d’une composante périodique \(e^{2i\pi\left(\frac{k_1 n_1}{N_1}+\frac{k_2 n_2}{N_2}\right)}\) présente dans le signal échantillonné. Cette composante possède donc une fréquence égale à \(\displaystyle\frac{k_1}{N_1}\) selon la première dimension, et \(\displaystyle\frac{k_2}{N_2}\) selon la seconde. La superposition de toutes ces composantes permet de reconstruire le signal de départ.

Calcul

Pour calculer la transformée de Fourier discrète à deux dimensions d’un tableau d’échantillons \(x[0, 0]\), \(x[0, 1]\), …, \(x[N_1 - 1, N_2 - 1]\), le plus simple consiste à exécuter l’algorithme suivant, qui exploite une propriété de séparabilité de cette transformée:

1. Pour \(k = 0, 1, \ldots, N_1-1\), remplacer la ligne

   \(x[k, 0],\ x[k, 1],\, \ldots,\, x[k, N_2 - 1]\)

   du tableau par la transformée de Fourier discrète (à une dimension) de cette séquence de valeurs.

2. Pour \(k = 0, 1, \ldots, N_2-1\), remplacer la colonne

   \(x[0, k],\ x[1, k],\, \ldots,\, x[N_1 - 1, k]\)

   du tableau par sa transformée de Fourier discrète (à une dimension).

La même stratégie permet également de calculer la transformée de Fourier discrète inverse d’un tableau, en appliquant successivement la transformée inverse aux lignes puis aux colonnes de celui-ci.

Implémentation

Vous pourriez implémenter la procédure précédente sous la forme d’une fonction acceptant un tableau bidimensionnel d’échantillons, et retournant un tableau de la même forme. Cependant, afin de pouvoir facilement réutiliser la fonction fft() que vous avez développée dans le laboratoire 4, il est plus commode de représenter les données à deux dimensions sous la forme d’un seul vecteur à une dimension.

En d’autres termes, plutôt que de placer les échantillons

\[\begin{split}\begin{array}{cccc} x[0, 0] & x[0, 1] & \cdots & x[0, N_2 - 1] \\ x[1, 0] & x[1, 1] & \cdots & x[1, N_2 - 1] \\ \vdots & \vdots & \ddots & \vdots \\ x[N_1 - 1, 0] & x[N_1 - 1, 1] & \cdots & x[N_1 - 1, N_2 - 1] \end{array}\end{split}\]

dans une matrice à \(N_1\) lignes et \(N_2\) colonnes, nous allons les représenter sous la forme d’une seule liste

\[[ x[0, 0], x[0, 1], \ldots, x[0, N_2 - 1], x[1, 0], x[1, 1], \ldots, x[1, N_2 - 1], x[2, 0], \ldots, x[N_1 - 1, N_2 - 1] ]\]

de taille \(N_1\,N_2\).

Vous êtes maintenant prêt à implémenter la transformée de Fourier discrète à deux dimensions. (L’idée est bien sûr de réutiliser au maximum les fonctions déjà obtenues dans le laboratoire 4.)

Marche à suivre:

1. Créer une fonction fft_2d() prenant comme arguments:

   • une liste valeurs de nombre complexes dont on souhaite calculer la transformée de Fourier discrète à deux dimensions.

   • deux entiers n1 et n2 représentant respectivement le nombre de lignes et le nombre de colonnes du tableau d’échantillons (La liste valeurs dont donc contenir n1 * n2 éléments.)

     Pour pouvoir réutiliser la fonction fft() du laboratoire 4, nous imposerons que n1 et n2 soient des puissances de 2. (Dans le cas contraire, la fonction fft_2d() doit retourner None.)

   Cette fonction doit retourner la transformée de Fourier discrète à deux dimensions du tableau d’échantillons, sous la même forme que valeurs (c’est-à-dire, une liste de n1 * n2 nombres complexes).

2. Tester votre fonction. Pour le tableau d’échantillons

   \[\begin{split}\begin{array}{cccc} 1 + i & -2 & 5 - 4i & 7 + 2i\\ 1 + i & -4 & 2 + i & -2i, \end{array}\end{split}\]

   vous devriez obtenir comme résultat

   \[\begin{split}\begin{array}{cccc} 10 - i & -5 + 18i & 8 - i & -5 -8i \\ 12 - i & -7 + 10i & -6 - 9i & 1. \end{array}\end{split}\]

3. Rédiger une fonction inv_fft_2d() implémentant la transformée de Fourier discrète inverse à deux dimensions, avec la même interface et les mêmes modalités d’utilisation que fft_2d().

4. Tester cette fonction inv_fft_2d() en vérifiant que si vous l’appliquez au tableau obtenu au point 2, vous retrouvez bien les échantillons de départ (aux imprécisions numériques près).

Traitement d’images

Nous allons maintenant programmer des mécanismes permettant d’appliquer à des images la transformée de Fourier discrète à deux dimensions.

Chargement et conversion en niveaux de gris

La première étape a pour objectif d’obtenir une fonction capable de charger une image depuis un fichier, et de la transformer en un tableau d’échantillons. Le plus simple consiste à utiliser les primitives de Pygame pour lire une image dans un format standard (en particulier JPEG ou PNG), et de convertir ensuite cette image en niveaux de gris. Les outils permettant d’effectuer ces opérations sont les suivants:

• Ce fragment de code charge une image depuis le fichier nommé nom_fichier, en calcule les dimensions, et en extrait un tableau de pixels:

  image = pygame.image.load(nom_fichier).convert(24)
  largeur, hauteur = image.get_size()
  pixels = pygame.PixelArray(image)
  

  Note: Ces instructions requièrent que le module Pygame ait été préalablement chargé et initialisé.

• A partir du tableau pixels généré par les instructions précédentes, on peut obtenir les composantes rouge (r), verte (g) et bleue (b) du pixel situé aux coordonnées (x, y) grâce au code suivant:

  v = pixels[x][y]
  r = v >> 16
  g = (v >> 8) & 0xff
  b = v & 0xff
  

  (Chacune de ces composantes possèdera alors une valeur comprise entre 0 et 255.)

• Pour l’espace de couleurs sRGB, La luminance \(l\) d’un pixel de couleur \((r, g, b)\) est donnée par la formule

  \(l = 0.2126 r + 0.7152 g + 0.0722 b\).

Vous pouvez donc à présent définir une fonction charger_image() acceptant comme argument un nom de fichier. Après avoir chargé et converti l’image correspondante, cette fonction doit retourner un triplet (echantillons, largeur, hauteur), où:

• largeur et hauteur donnent les dimensions de l’image,
• echantillons est une liste de taille largeur * hauteur contenant la luminance des pixels de l’image, avec une valeur dans l’intervalle [0, 255].

Sauvegarde

Afin de pouvoir tester votre code, vous pouvez programmer l’opération symétrique, c’est-à-dire une fonction sauver_image() prenant comme arguments

• un triplet (echantillons, largeur, hauteur) de la même forme que ceux retournés par charger_image(),
• une chaîne de caractères nom_fichier fournissant le nom du fichier à créer.

Cette fonction doit sauvegarder dans le fichier l’image correspondant au tableau d’échantillons, en niveaux de gris. Les mécanismes suivants vous seront utiles:

• Ces instructions créent un tableau de pixels de dimension (largeur, hauteur):

  image = pygame.Surface((largeur, hauteur))
  pixels = pygame.PixelArray(image)
  

• Pour attribuer la luminance v au pixel de coordonnées (x, y), on peut exécuter:

  pixels[x][y] = (v, v, v).
  

  Notes:

  • Cela revient à attribuer la valeur v aux composantes rouge, verte et bleue de ce pixel.
  • Pour cette opération, il faut veiller à ce que la valeur de v reste bien comprise entre 0 et 255.

• Après avoir attribué une valeur à chaque pixel du tableau, l’instruction suivante sauve celui-ci sous la forme d’une image dans le fichier nom_fichier:

  pygame.image.save(pixels.make_surface(), nom_fichier)
  

  Note: Le format de l’image est déterminé sur base de l’extension figurant dans le nom du fichier, par exemple “.png”.

N’oubliez pas de tester soigneusement ce code. Vous pouvez le faire en chargeant une image quelconque grâce à charger_image(), en sauvegardant le tableau d’échantillons obtenu vers un autre fichier via sauver_image(), et en examinant le résultat grâce à un outil de visualisation d’images (par exemple, eog ou gimp).

Transformée de Fourier d’une image

Vous avez maintenant tous les outils en main pour écrire un programme qui calcule la transformée de Fourier discrète à deux dimensions d’une image. Ce programme peut soit afficher le résultat obtenu à l’écran, soit le sauvegarder dans des fichiers. Nous allons opter pour la seconde solution, afin de pouvoir plus tard modifier la transformée de l’image à l’aide d’outils externes, et étudier l’impact de ces modifications sur l’image de départ.

Sauvegarder la transformée de Fourier discrète à deux dimensions d’une image vers un fichier présente cependant une difficulté: Cette transformée prend la forme d’un tableau de nombres complexes, et non de valeurs réelles que l’on pourrait directement traduire en luminance de pixels. La solution que nous allons mettre en oeuvre est la suivante:

• On sépare chaque valeur \(X[k_1, k_2]\) de la transformée en son module \(|X[k_1, k_2]|\) et sa phase \(\arg(X[k_1, k_2])\).

• On sauvegarde séparément le module et la phase des valeurs de la transformée dans deux images distinctes, en niveaux de gris:

  • Dans la première, la luminance \(\ell\) de chaque pixel est déterminée d’après le au module \(m\) de la valeur correspondante.

    Afin d’obtenir une luminance comprise dans l’intervalle \([0, 255]\) sans perdre trop d’information, on suit la règle simple suivante:

    • Si \(m = 0\), alors \(\ell = 0\).

    • Si \(m > 0\), alors

      \(\displaystyle\ell = \frac{255 \log_2 m}{7.99 + \log_2 N_1 + \log_2 N_2}\),

      où \(N_1\) et \(N_2\) sont respectivement le nombre de lignes et le nombre de colonnes de l’image.

  • Dans la seconde, la luminance \(\ell\) de chaque pixel représente la phase \(\phi \in [-\pi, \pi]\) de la valeur correspondante, selon la loi \(\ell = 255 \frac{\phi + \pi}{2\pi}\). (En d’autres termes, cette luminance interpole linéairement la phase entre \(\ell = 0\) qui représente \(\phi = -\pi\), et \(\ell = 255\) qui correspond à \(\phi = \pi\).)

• On sauvegarde aussi une troisième image destinée à nous permettre de visualiser la transformée. Cette image est similaire à celle dans laquelle on sauve le module, à ceci près qu’au lieu d’utiliser une loi logarithmique, on calcule la luminance \(\ell\) à partir du module \(m\) grâce à l’équation

  \(\displaystyle\ell = \frac{m}{\sqrt{N_1\,N_2}}\),

  où \(N_1\) et \(N_2\) sont respectivement le nombre de lignes et le nombre de colonnes de l’image.

• Par convention, la valeur \(X[0, 0]\) correspondra au centre des trois images (souvenez-vous que la transformée de Fourier discrète à deux dimensions est périodique, donc l’image peut être décalée comme on le souhaite).

Votre objectif consiste maintenant à développer un programme prog-13.py qui:

• accepte deux arguments nom_image et prefixe. Le premier fournit le nom du fichier contenant l’image à analyser. Le second donne le préfixe qui sera utilisé pour générer le nom des trois fichiers créés par le programme.
• charge l’image depuis nom_image et en calcule la transformée de Fourier discrète à deux dimensions.
• sauvegarde cette transformée dans trois fichiers: Celui représentant le module aura pour nom préfixe suivi par “-mod.png”, celui contenant la phase prefixe suivi par “-phase.png”, celui destiné à la visualisation prefixe suivi par “-visu.png”

Remarques:

• Etant donné que nous n’avons implémenté le calcul de la transformée de Fourier que pour des ensembles d’échantillons dont la taille est une puissance de deux, le programme doit retourner une erreur si une des deux dimensions de l’image fournie n’est pas une puissance de deux. Si vous souhaitez essayer le programme avec une image donnée, il faudra donc préalablement la découper aux dimensions appropriées (par exemple, en utilisant l’outil crop de gimp).
• Implémenter la sauvegarde de la transformée dans une fonction séparée est une bonne idée. Bien sûr, il est pertinent d’exploiter dans cette fonction celle que vous avez déjà programmée pour sauvegarder une image.
• En Python:
  • L’opérateur + permet de concaténer des chaînes de caractères.
  • Les expressions abs(z) et cmath.phase(z) permettent de récupérer respectivement le module et l’argument (c’est-à-dire, la phase) du nombre complexe z.
  • La fonction math.log2() permet de calculer le logarithme en base 2 d’un nombre.

Premiers essais

Pour commencer, vous pouvez essayer votre programme en lui fournissant l’image suivante:

Exemple de signal périodique (1).

Si tout est correct, vous devriez obtenir l’image suivante de la transformée de Fourier discrète à deux dimensions de ce signal (dans le fichier “…-visu.png”):

Transformée de l’exemple précédent.

Remarquez que cette image n’est pas entièrement noire: Elle contient trois points blancs correspondant à

• la composante continue du signal (au centre de l’image),

• la composante périodique du signal et son conjugué.

  Dans notre exemple, l’image d’origine possède un motif qui se répète suivant l’axe vertical, avec une période égale à 32 pixels. L’image ayant une hauteur de 512 pixels, cela correspond à une fréquence de 512/32 = 16. Sur la transformée de Fourier, les pixels correspondant à cette composante périodique se situent donc sur l’axe vertical, à 16 pixels du centre de l’image.

Ensuite, vous pouvez expérimenter sur base des images suivantes, dans lesquelles la fréquence et/ou la direction de la composante périodique diffèrent:

Exemple de signal périodique (2).

Exemple de signal périodique (3).

Exemple de signal périodique (4).

Qu’obtenez-vous comme transformée de ces images? Etes-vous capable de mettre en correspondance ce que vous observez sur les transformées avec les composantes présentes dans les images?

Analyse d’images plus complexes

Les images des exemples précédents ne contenaient qu’une seule composante périodique. En général, les images en contiennent plusieurs, et l’intérêt de la transformée de Fourier discrète à deux dimensions est précisément d’être capable de les séparer.

Par exemple, l’image suivante correspond à la superposition de deux composantes périodiques:

Superposition de deux composantes périodiques (5).

Calculez-en la transformée. Qu’observez-vous? Les deux composantes sont-elles correctement isolées?

Avant de passer à des images plus réalistes, remarquons que celles que nous avons considérées jusqu’ici sont particulières: Les périodes de leurs composantes périodiques divisent exactement les dimensions horizontale et verticale de l’image. Cela se voit facilement en pavant le plan avec des copies d’une de ces images: Un tel pavage ne présente aucun raccord visible entre une copie de l’image et ses voisines.

Dans un cas plus général, lorsque la période des composantes périodique ne divise pas exactement les dimensions de l’image, un phénomène de fuite de spectre apparaît, comme c’était déjà le cas pour la transformée de Fourier discrète à une dimension. Vous pouvez par exemple observer ce phénomène en calculant la transformée des images suivantes:

Exemple de signal périodique (6).

Exemple de signal périodique (7).

Exemple de signal périodique (8).

Exemple de signal périodique (9).

Afin de pouvoir tester votre programme, voici la visualisation de la transformée que vous devriez obtenir pour cette dernière image:

Transformée de l’exemple 9.

Comme on le voit sur ce dernier exemple, malgré le phénomène de fuite spectrale, la transformée de Fourier discrète à deux dimensions reste capable de bien mettre en évidence les composantes périodiques présentes dans l’image. Cela fonctionne aussi avec des images moins artificielles! Essayez ainsi de calculer la transformée des deux images suivantes, et d’en interpréter les résultats:

Exemple d’image (1).

Exemple d’image (2).

Si votre programme fonctionne correctement, placez-en une copie dans le répertoire centralisé des laboratoires, sous le suffixe prog-13.py.

Transformation d’images

Notre dernier objectif consistera à apporter des modifications à une image en en calculant d’abord la transformée, puis en faisant subir certaines opérations à celle-ci, et en reconstruisant enfin l’image grâce à une transformée inverse.

Outil de transformation inverse

Pour ce faire, vous allez développer un programme prog-14.py réalisant l’opération symétrique de prog-13.py. Ce programme doit:

• accepter deux arguments prefixe et nom_sauvegarde. Le premier permet de déterminer le nom des deux fichiers à partir desquels la transformée doit être lue (contenant le module et la phase des valeurs), et le second donne le nom du fichier vers lequel l’image reconstruite doit être écrite.
• charger la transformée de Fourier discrète à deux dimensions dont le module est contenu dans le fichier prefixe suivi par “-mod.png”, et la phase dans prefixe suivi par “-phase.png”. (Si ces fichiers contiennent des images de dimensions différentes, ou non égales à des puissances de deux, le programme doit signaler une erreur.)
• calculer la transformée inverse de façon à reconstruire une image.
• sauvegarder cette image dans le fichier nom_sauvegarde.

Remarques: Attention, n’oubliez pas que:

• Dans le fichier représentant le module, chaque valeur \(m\) de ce module a été transformée par la loi logarithmique suivante:

  • Si \(m = 0\), alors \(\ell = 0\).

  • Si \(m > 0\), alors

    \(\displaystyle\ell = \frac{255 \log_2 m}{7.99 + \log_2 N_1 + \log_2 N_2}\),

    où \(N_1\) et \(N_2\) sont respectivement le nombre de lignes et le nombre de colonnes de l’image.

  Pour reconstruire \(m\) à partir de \(\ell\), il faut donc appliquer la transformation réciproque.

• Dans les fichiers de module et de phase, le centre de l’image correspond à l’élément \(X[0, 0]\) de la transformée. Il ne faut donc pas oublier d’effectuer le décalage correspondant.

Pour tester ce programme, vous pouvez partir de n’importe quelle image déjà traitée, par exemple celle contenant le dessin d’un cube, et essayer de la reconstruire à partir de sa transformée de Fourier.

Cela fonctionne-t-il? Si oui, déposez votre programme dans le répertoire centralisé, avec le suffixe prog-14.py.

Filtre passe-bas

En utilisant prog-13.py, calculez la transformée d’une image de votre choix, par exemple celle du cube. Vous devriez donc obtenir trois fichiers cube-mod.png, cube-phase.png et cube-visu.png.

Dans un éditeur d’images externe, par exemple gimp, ouvrez cube-mod.png et modifiez cette image de façon à n’en garder que la partie centrale (le reste de l’image étant peint en noir). L’objectif est d’arriver à une image similaire à celle-ci (dans cet exemple, le disque central possède un rayon de 32 pixels):

Transformée du cube modifiée par un filtre passe-bas.

Pour rappel, les valeurs de la transformée représentes les composantes fréquentielles de l’image, cette fréquence devenant de plus en plus grande lorsqu’on s’éloigne du centre. Notre modification a donc eu pour effet d’éliminer toutes les composantes au delà d’une certaine fréquence.

Si vous reconstruisez l’image grâce à prog-14.py (en gardant le fichier de phase inchangé), vous devriez retrouver ceci:

Cube reconstruit.

Comme vous le voyez, l’élimination des hautes fréquences a introduit du flou dans l’image, en atténuant les transitions brusques de luminance. L’image contient cependant des artifices, causés par la coupure de fréquences très abrupte de notre filtre: Les fréquences situées à l’intérieur du disque utilisé comme masque sont préservées, et celles à l’extérieur sont éliminées, sans transition.

Pour atténuer cet effet, nous pouvons utiliser un masque aux bords plus progressifs, par exemple en utilisant l’option de feathering de la sélection dans gimp. En effectuant cette opération sur une distance de 40 pixels, la transformée ressemble à la suivante:

Transformée du cube modifiée par un filtre passe-bas progressif.

La reconstruction de l’image donne alors le résultat suivant

Cube reconstruit.

qui représente un flou plus réaliste.

Filtre passe-haut

De nombreuses autres manipulations d’images sont possibles! Comme dernier exemple, nous effectuons l’opération contraire à la précédente, en ne gardant que les hautes fréquences de l’image.

Note: Il faut cependant veiller à ne pas modifier la valeur de l’élément \(X[0, 0]\) de la transformée, situé au centre de l’image, afin de préserver la composante continue du signal.

Sur base de la transformée

Transformée du cube modifiée par un filtre passe-haut.

(remarquez le pixel central), on reconstruit l’image

Cube reconstruit.

dans lequel seuls les bords des lignes ont cette fois été préservés.

Si vous avez réussi à reproduire ces expériences avec votre propre implémentation, bravo! La transformée de Fourier discrète permet bien d’autres opérations; n’hésitez pas à explorer et à expérimenter par vous-même, ainsi que d’appliquer les opérations que nous avons étudiées à d’autres types d’images.