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 icitransformée de Fourier. 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és avec Code::blocks sous Windows et avec les fichiers de bibliothèques se trouvant ici (fichier zip)FFTW
Cette bibliothèque est disponible sur plusieurs plates-formes (aussi bien sous Windows que sous UNIX). La page officielle se trouve ici. Elle a été développée 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 licence GPL ;
- et d'autres choses… (parallélisation des calculs…).
Attention, certaines anciennes versions de cette bibliothèque nécessitent 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 forcément 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 normale 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 suivant 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 remplie 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 sympathique. 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 :
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.
/* reIn : partie réelle 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 aux limites 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 base sur les complexes de passage à coordonnées cartésiennes à coordonnées polaires.
/*
* reIn : partie réelle 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éelle/imaginaire.
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.
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 proportionnalité*/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.
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.
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);
}
V. Remerciements▲
Je tiens à remercier troumad pour les corrections apportées.