Utilisation de la bibliothèque fftw en traitement d'images

Date de publication : 26/01/2007 , Date de mise à jour : 26/01/2007

Dans ce tutoriel, je présenterai l'utilisation de la bibliothèque fftw dans le cadre du traitement numérique d'images. Je n'entrerai pas dans les détails d'implémentation de cette bibliothèque, et tenterai de rester au plus simple quant à son utilisation.

I. Introduction

Ce cours est la continuité d'un premier cours à caractère purement algorithmique sur l'utilisation des transformées de Fourier en traitement numérique des images. L'index de ce cours est disponible ici. Vous y trouverez le cours purement théorique, un lien vers ce cours, et les sources C++ permettant d'utiliser la bibliothèque fftw.

FFTW est une bibliothèque de fonction écrite en C permettant d'effectuer des calculs de transformée de Fourier discrète en une ou plusieurs dimensions.

Tous les exemples ont été essayé avec Code::blocks sous Windows et avec les fichiers de bibliothèques se trouvant ici (fichier zip)

Cette bibliothèque est disponible sur plusieurs plate-forme (aussi bien sous Windows que sous UNIX). La page officielle se trouve ici. Elle a été développé au MIT par Matteo Frigo et Steven G. Johnson.

Pour être plus précis, cette bibliothèque admet de nombreux avantages :

calcul rapide
utilisable en C et C++ sur toutes les plates-formes
transformation en une ou plusieurs dimensions
taille arbitraire des données (il n'est pas obligatoire d'avoir un signal de taille 2^n)
transformations en cosinus et sinus discret (DCT et DST)
sous license GPL
et d'autres choses... (parallélisation des calculs...)

Attention, certaines anciennes versions de cette bibliothèque nécessite de faire du padding (compléter l'image afin que sa taille soit une puissance de 2^n), mais la dernière ne pose pas de problème.

Nous ne verrons pas toute son utilisation, mais uniquement ce qu'il nous intéresse pour réaliser des transformées de Fourier sur des images.

Nous verrons également différentes parties sur les problèmes que l'on peut rencontrer pour le traitement des images, notamment la visualisation du module d'une image.

Je précise que les sources données ici ne sont pas forcement exactement les mêmes que celles que vous trouverez dans les sources. Ayant utilisé pas mal d'objets sur les images dans mes sources, je n'ai pas voulu les expliciter ici. J'ai donc retranscrit le code pour le langage C.

Vous pouvez avoir l'ensemble des sources C++ (ici)

II. Transformation d'une image

Nous allons voir dans cette section comment calculer la transformée d'une image.

II-A. Les fonctions de FFTW

La bibliothèque fftw fournit principalement deux fonctions permettant d'effectuer les transformées : fftw_plan_dft_2d et fftw_execute. Les prototypes de ces fonctions sont les suivantes :


 fftw_plan fftw_plan_dft_2d(int nx, int ny,
                                fftw_complex *in, fftw_complex *out,
                                int sign, unsigned flags);
                                
 void fftw_execute(const fftw_plan plan);

fftw_plan_dft_2d retourne un plan d'exécution, qui contient des informations : les tailles de l'image, l'image d'entrée, l'image de sortie, le type de transformation (inverse, normal), des options d'exécution.

La fonction prend en entrée :

nx : la taille horizontale de l'image
ny : la taille verticale de l'image
in : l'image d'entrée au format fftw_complex
out : l'image de sortie dans le même format
sign : indique si l'on effectue une transformée normal ou inverse
flags : indique des options d'exécutions

sign peut valoir FFTW_FORWARD ou FFTW_BACKWARD qui correspond en fait au signe dans l'exponentielle (dans la définition de la transformée). FFTW_BACKWARD correspond donc à la transformée inverse et FFTW_FORWARD correspond à la transformée normale.

Nous utiliserons flags avec FFTW_ESTIMATE qui permettra le calcul simple de la transformée. Il y a d'autres options, par exemple FFTW_MEASURE, qui réalisera plusieurs fois la FFT pour optimiser le plan d'exécution.

Les images sont au format fftw_complex*, qui disposent de légers changements suivants la version de la fftw. L'espace s'alloue grâce à la fonction fftw_malloc et se supprime grâce à la fonction fftw_free. Ces fonctions s'utilisent de la même manière que leurs analogues malloc et free, mais disposent de légères différences pour permettre des optimisations.
fftw_complex est un tableau à deux dimensions, la deuxième dimension n'étant que de taille 2. L'une pour les données réelles, l'une pour les données imaginaires.

 
   float* in;
 
   /*
    * CODE
    */
 
   unsigned int i;
   spatial_repr= fftw_malloc(sizeof(fftw_complex)*largeur*hauteur);

   for(i=0;i<largeur*hauteur;i++)
   {
      spatial_repr[i][0] = in[i];
      spatial_repr[i][1] =  0.0f; /*rempli de 0 la partie imaginaire de l'image*/
   }

L'image est une image complexe, ainsi, lors de la transformée normale, on rempliera la partie complexe par des 0 et la partie imaginaires par notre image. in est une image réelle rempli de ligne en ligne. Par exemple, l'image :

1	2	8
4	5	6

Aura comme représentation :


   float* in = {1.0f, 2.0f, 8.0f, 4.0f, 5.0f, 6.0f}

fftw_execute va permettre d'exécuter le plan d'exécution retourné par la fonction précédente.

II-B. Code de la transformée

Avec les éléments que nous avons (et avec quelques autres que je vais vous fournir), nous allons pouvoir calculer la transformée d'une fonction.

Voici, le premier code minimal que nous pouvons obtenir :

 
/* in : image d'entrée
 * reOut : partie réelle de l'image de sortie
 * imOut : partie imaginaire de l'image de sortie
 * largeur : largeur des images d'entrée et de sortie
 * hauteur : idem
 */
void fourierForward(const float* in,
                    float* reOut,
                    float* imOut,
                    unsigned int largeur,
                    unsigned int hauteur)
{
   fftw_complex* spatial_repr;
   fftw_complex* frequency_repr;
   unsigned int i;
   fftw_plan plan;

   spatial_repr= fftw_malloc(sizeof(fftw_complex)*largeur*hauteur);
   frequency_repr= fftw_malloc(sizeof(fftw_complex)*largeur*hauteur);

   /*On remplit la structure qui sera utilisée par fftw*/
   for(i=0;i<largeur*hauteur;i++)
   {
      spatial_repr[i][0] = in[i];
      spatial_repr[i][1] =  0.0f;
   }

   /*on calcule le plan d'exécution*/
   plan=fftw_plan_dft_2d(hauteur, largeur, spatial_repr, frequency_repr, FFTW_FORWARD, FFTW_ESTIMATE);

   /*on effectue la transformée de Fourier*/
   fftw_execute(plan);

   /*on retranscrit le résultat en 2 images, l'une représentant la partie réelle, l'autre
     la partie imaginaire*/
   for(i=0;i<largeur*hauteur;i++)
      {
          reOut[i]=frequency_repr[i][0];
          imOut[i]=frequency_repr[i][1];
      }
   
   /*on détruit les objets*/
   fftw_destroy_plan(plan);
   fftw_free(spatial_repr);
   fftw_free(frequency_repr);
}

Si l'on cherche à obtenir le module (nous verrons après comment le visualiser) de l'image précédente, nous obtenons l'image suivante.

Ceci n'est pas très sympatique. En effet, le point (0,0) de la transformée de Fourier se trouve effectivement au point (0,0) de l'image. On aurait préféré qu'elle se trouve au milieu (ce qui permet d'appliquer les filtres plus simplement). Pour résoudre ce problème, il suffit de faire une translation des points (i,j) vers les points (i+largeur/2) % largeur et (j+longueur/2) % longueur.

Le code devient donc celui-ci :

Code de la transformée de Fourier discrète

 
void fourierForward(const float* in,
                    float* reOut,
                    float* imOut,
                    unsigned int largeur,
                    unsigned int hauteur)
{
   fftw_complex* spatial_repr;
   fftw_complex* frequency_repr;
   unsigned int i;
   unsgiend int j;
   fftw_plan plan;
   int x,y;
   
   spatial_repr= fftw_malloc(sizeof(fftw_complex)*largeur*hauteur);
   frequency_repr= fftw_malloc(sizeof(fftw_complex)*largeur*hauteur);


   for(i=0;i<largeur*hauteur;i++)
   {
      spatial_repr[i][0] = in[i];
      spatial_repr[i][1] =  0.0f;
   }

   /*on calcule le plan d'exécution*/
   plan=fftw_plan_dft_2d(hauteur, largeur, spatial_repr, frequency_repr, FFTW_FORWARD, FFTW_ESTIMATE);

   /*on calcule la transformée*/
   fftw_execute(plan);

  for(j=0;j<hauteur;j++)
      for(i=0;i<largeur;i++)
      {
	        /*on recentre l'image*/
	      x=i;
	      y=j;
	      if (i<largeur/2 && j<hauteur/2){ x=i+largeur/2; y=j+hauteur/2; }
          if (i>=largeur/2 && j<hauteur/2){ x=i-largeur/2; y=j+hauteur/2; }
          if (i<largeur/2 && j>=hauteur/2){ x=i+largeur/2; y=j-hauteur/2; }
          if (i>=largeur/2 && j>=hauteur/2){ x=i-largeur/2; y=j-hauteur/2; }
          reOut[y*largeur+x]=frequency_repr[j*largeur+i][0];
          imOut[y*largeur+x]=frequency_repr[j*largeur+i][1];
      }

   fftw_destroy_plan(plan);
   fftw_free(spatial_repr);
   fftw_free(frequency_repr);

}

Ce qui donne l'image suivante (les basses fréquences étant au milieu, et les hautes fréquences vers les bords).

II-C. Transformée inverse

La transformée inverse est à peu près du même type. On repositionne l'image (que l'on avait centrée avant). Et on applique la fonction fftw_plan_dft_2d avec l'argument FFTW_BACKWARD. Il faut bien se rappeler que cet argument ne fait que changer le signe de l'exponentielle. Il faudra donc bien penser à diviser l'image finale par (largeur * hauteur) pour avoir la bonne formule de la transformée de Fourier inverse.

Code de la transformée de Fourier inverse discrète

 
/* reIn : partie réel de l'image dans l'espace de Fourier
 * imIn : partie imaginaire de l'image dans l'espace de Fourier
 * out : image de sortie
 * largeur : largeur des images d'entrée et de sortie
 */
void fourierBackward(const float* reIn,
                     const float* imIn,
                     float* out,
                     unsigned int largeur,
                     unsigned int hauteur)
{
   fftw_complex* spatial_repr;
   fftw_complex* frequency_repr;
   unsigned int i;
   unsigned int j;
   int x,y;
   fftw_plan plan;

   spatial_repr= fftw_malloc(sizeof(fftw_complex)*largeur*hauteur);
   frequency_repr= fftw_malloc(sizeof(fftw_complex)*largeur*hauteur);

   for(j=0;j<hauteur;j++)
      for(i=0;i<largeur;i++)
      {
          /*on décentre*/
	      x=i;
	      y=j;
	      if (i<largeur/2 && j<hauteur/2){ x=i+largeur/2; y=j+hauteur/2; }
	      if (i>=largeur/2 && j<hauteur/2){ x=i-largeur/2; y=j+hauteur/2; }
	      if (i<largeur/2 && j>=hauteur/2){ x=i+largeur/2; y=j-hauteur/2; }
	      if (i>=largeur/2 && j>=hauteur/2){ x=i-largeur/2; y=j-hauteur/2; }
	      frequency_repr[j*largeur+i][0]=reIn[y*largeur+x];
	      frequency_repr[j*largeur+i][1]=imIn[y*largeur+x];
      }

  plan=fftw_plan_dft_2d(hauteur, largeur, frequency_repr, spatial_repr, FFTW_BACKWARD, FFTW_ESTIMATE);

  fftw_execute(plan);

   /*on retranscrit l'image complexe en image réelle, sans oublier de diviser par largeur*hauteur*/
   for(i=0;i<largeur*hauteur;i++)
   {
      out[i]=spatial_repr[i][0]/(largeur*hauteur);
   }

   fftw_destroy_plan(plan);
   fftw_free(spatial_repr);
   fftw_free(frequency_repr);
}

III. Visualisation du module d'une image

Lorsque l'on passe une image dans l'espace de Fourier, on la représente souvent par son module (même si la phase est un gros support d'information). Mais, si on la calcule, il arrive fréquemment qu'elle atteigne des valeurs hors bornes (par rapport au limite de l'image 0-255). Il y a de nombreuses stratégies, nous allons en voir une en particulier.

III-A. Calcul du module et de la phase

Nous allons tout d'abord voir comment calculer le module d'une image. Ce calcul va en fait se faire assez simplement, en appliquant les formules de bases sur les complexes de passage à coordonnées cartésiennes à coordonnées polaires.

Calcul phase module

 
/*
 * reIn : partie réel de l'image d'entrée
 * imIn : partie imaginaire de l'image d'entrée
 * moduleOut : image représentant le module de l'image
 * phaseOut : image représentant la phase de l'image
 *
 * Précondition : reIn, imIn, moduleOut et phaseOut doivent avoir les mêmes tailles
 */
void fourierI2MP(const float* reIn,
                 const float* imIn,
                 float* moduleOut,
                 float* phaseOut,
                 unsigned int largeur,
                 unsigned int hauteur)
{
 unsigned int i;
 
 for(i=0;i<hauteur*largeur;i++)
   {
      moduleOut[i]=sqrt(reIn[i]*reIn[i] + imIn[i]*imIn[i]); /*on applique la formule de base*/
      if (reIn[i]!=0)
	      phaseOut[i]=atan2(imIn[i],reIn[i]);
   }
}

Au passage, je vais donner le code permettant de réaliser l'opération inverse. C'est à dire de passer de la représentation module/phase en représentation réel/imaginaire.

Calcul phase module

 
void fourierMP2I(const float* moduleIn,
                 const float* phaseIn,
                 float* reOut,
                 float* imOut,
                 unsigned int largeur,
                 unsigned int hauteur)
{
   unsigned int i;
   
   for(i=0;i<largeur*hauteur;i++)
   {
      reOut[i]=moduleIn[i]*cos(phaseIn[i]);
      imOut[i]=moduleIn[i]*sin(phaseIn[i]);
   }
}

III-B. Ajustement du module

Comme nous l'avons dit, les valeurs du module peuvent-être totalement hors bornes. Une manière de faire va consister à calculer le maximum, et à multiplier chaque point par : 255 / maximum. (en considérant que les couleurs vont de 0 à 255).

Cependant, si vous faites directement cela, vous risquez d'avoir un point noir avec un point blanc au milieu. En effet, le point du milieu semble être critique par rapport aux autres et avoir une très forte valeur. Il suffira alors d'exclure un petit cercle de rayon petit au milieu.

Ajustement module

 
static double carre(double a)
{
 return (a*a);
}

.........

 unsigned int i;
 unsigned int j;
 float facteur = 1.0f;
 
 float pourcentImage = 0.01 * (largeur()); /*on exclue un pourcentage du cercle du milieu*/
 
 for(j = 1; j<hauteur; j++)
   for(i = 1; i<largeur; i++)
     if(max<module[i] && (carre((i-largeur/2))+carre((j-hauteur/2)) ) > carre(pourcentImage))
        max = module[i];

 if(max==0) /*pour éviter une division par 0*/
  max = 1;

 for(i = 0; i<largeur*hauteur; i++)
    /*il faut parfois ajuster le facteur*/
    module[i] *= 255 * facteur/ max; /*on calcule la proportionalitée*/

IV. Application de filtres

Nous allons voir ici comment appliquer un filtre déjà calculé (par exemple un flou Gaussien).

IV-A. Code de l'application d'un filtre quelconque

Nous allons appliquer directement un filtre qui se trouve déjà dans l'espace de Fourier. Nous supposons que le filtre et l'image ont la même taille. Le code ne fait en fait, que calculer la multiplication des complexes entre le filtre et l'image en tous les points.

Application d'un filtre

 
void filtreAppliquer(float* reTabImage,
                     float* imTabImage,
                     const float* reTabFiltre,
                     const float* imTabFiltre,
                     unsigned int largeur,
                     unsigned int hauteur)
{
					 
  float a,b,c,d;
  unsigned int i;
  
  for(i=0; i< largeur * hauteur; i++)
  {
    a = reTabImage[i];
    b = imTabImage[i];
    c = reTabFiltre[i];
    d = imTabFiltre[i];

    /*calcul du produit des complexes*/
    reTabImage[i] = a*c - b*d;
    imTabImage[i] = b*c + a*d;

  }

}

IV-B. Filtre gaussien

Je ne donnerai dans ce cours, que l'exemple du filtre Gaussien. D'autres filtres se trouvent dans le code source de ce cours.

Flou gaussien

 
void filtreFlouGaussien(float* reTabImage,
                        float* imTabImage,
                        unsigned int largeur,
                        unsigned int hauteur,
                        float sigma) /*sigma du flou*/
{
 unsigned int i;
 unsigned int j;
 
 float* filtreR = malloc(sizeof(float) * largeur*hauteur);
 float* filtreI = malloc(sizeof(float) * largeur*hauteur):
 
 if(filtreR==NULL || filtreI ==NULL)
 {
   fprintf(stderr, "Probleme d'allocation memoire");
   exit(EXIT_FAILURE);
 }
 
 for(j=0; j<hauteur;j++)
  for(i=0; i<largeur;i++)
  {
    /*on applique la formule du flou gaussien*/
    filtreR[i + largeur*j] =
     exp(- carre(sigma) * (carre(i-largeur/2.0) + carre(j-hauteur/2.0))/(2.0));
    filtreI[i+ largeur *j] = 0.0f;
    /*on a pensé à recentrer en 0, car le point (0,0) correspond au point (largeur/2, hauteur/2)*/
  }

  filtreFourierAppliquer(reTabImage, imTabImage, filtreR, filtreI, largeur, hauteur);

  free(filtreR);
  free(filtreI);
}

Flou gaussien avec sigma = 0.02

Flou gaussien avec sigma = 0.05

Flou gaussien avec sigma = 0.1

Inversion du module de deux images

V. Remerciements

Je tiens à remercier troumad pour les corrections apportées.

Les sources présentées sur cette page sont libres de droits, et vous pouvez les utiliser à votre convenance. Par contre, la page de présentation constitue une oeuvre intellectuelle protégée par les droits d'auteurs. Copyright © 2006 Florent HUMBERT. Aucune reproduction, même partielle, ne peut être faite de ce site et de l'ensemble de son contenu : textes, documents, images, etc sans l'autorisation expresse de l'auteur. Sinon vous encourez selon la loi jusqu'à 3 ans de prison et jusqu'à 300 000 E de dommages et intérêts. Cette page est déposée à la SACD.