1. Introduction
Ce document introduit la transformée de Fourier d’une image, puis la transformée de Fourier discrète (TFD) d’une image échantillonnée. Le calcul de la TFD d’une image avec Python est expliquée. On verra comment représenter le spectre de l’image et comment effectuer un filtrage dans l’espace des fréquences, en multipliant la TFD par une fonction de filtrage.
2. Transformée de Fourier et transformée de Fourier discrète
On considère une image monochrome (niveaux de gris) représentée par une fonction de deux variables réelles, à valeurs complexes, notée u(x,y).
La transformée de Fourier de cette image est la fonction à deux variables réelles et à valeurs complexes définie par :S(fx,fy)=∫-∞∞∫-∞∞u(x,y)exp(-i2π(fxx+fyy))dxdy
On considère à présent un échantillonnage de l’image sur un domaine rectangulaire de dimensions (Lx,Ly) centré en (x,y)=(0,0) et comportant (Nx,Ny) points. La matrice U est définie par :xk=-Lx2+kLxNx0≤k≤Nx-1ym=-Ly2+mLyNy0≤m≤Ny-1Um,k=u(xk,ym)
La matrice U a donc Nx colonnes et Ny lignes.
Considérons une approximation de la transformée de Fourier obtenue par la méthode des rectangles sur le domaine considéré :S(fx,fy)≃LxLyNxNyexpiπ(fxLx+fyLy)∑k=0Nx-1∑m=0Ny-1Um,kexp-i2πfxkLxNxexp-i2πfymLyNy
En pratique, on calcule cette approximation pour les fréquences spatiales suivantes :fx,n=nLx0≤n≤Nx-1fy,l=lLy0≤l≤Ny-1
La transformée de Fourier discrète (TFD) associe à la matrice U une matrice V de mêmes dimensions, définie par :Vn,l=1NxNy∑k=0Nx-1∑m=0Ny-1Um,kexp-i2πknNxexp-i2πmlNy
Remarque : les TFD calculées par les fonctions FFT des logiciels omettent parfois l’inverse de NxNy en facteur.
L’approximation de l’intégrale est :S(fx,n,fy,l)≃LxLyexpiπ(n+l)Vn,l
La fréquence d’échantillonnage est constituée par la paire :fe=NxLx,NyLy
On vérifie facilement la propriété suivante de la TFD:Vn+Nx,l=Vn,lVn,l+Ny=Vn,l
Cette propriété correspond à la périodicité du spectre de l’image échantillonnée. La période du spectre est égale à la fréquence d’échantillonnage. En général, on cherche à calculer une approximation de la transformée de Fourier ne faisant pas apparaitre cet effet de l’échantillonnage. Pour cela, il faut limiter la plus haute fréquence à la moitié de la fréquence d’échantillonnage. Pour simplifier, on suppose que Nx et Ny sont pairs. La TFD donne une approximation de la transformée de Fourier pour les fréquences suivantes :fx,n=nLx0≤n≤Nx2fy,l=lLy0≤l≤Ny2
De plus, la périodicité de la TFD permet d’accéder aux opposées de ces fréquences. En effet :V-n,l=V-n+Nx,lVn,-l=Vn,-l+Ny
Les Nx indices n correspondent aux fréquences suivantes :0,1Lx,2Lx,…,Nx2Lx,-Nx2-1Lx,…,-2Lx,-1Lx
La même relation est valable pour les fréquences de l’axe y.
Pour certains calculs (filtrage, diffraction), il est utile de placer la fréquence nulle au milieu de la matrice, de manière à obtenir des fréquences croissantes en fonction des indices. Notons VC cette matrice centrée. Dans cette matrice, le point de fréquence nulle doit être placé aux indices suivants :(Ny2-1,Nx2-1)
Pour l’axe x, il y a donc Nx2+1 indices correspondant à des fréquences positives (ou nulle) et Nx2-1 indices correspondant à des fréquences strictement négatives.
Pour construire la matrice VC, il faut utiliser la transformation d’indice suivante :nc=Nx2-1+nsi0≤n≤Nx2nc=-Nx2-1+nsiNx2+1≤n≤Nx-1
3. Transformée de Fourier d’une image avec Python
3.a. Principe
Le module imageio permet de lire des fichiers d’image. L’image suivante a 200×100 pixels.
L’image est obtenue sous forme d’un tableau numpy à trois dimensions avec la fonction imageio.imread. On prélève la couche rouge et on convertit le tableau en flottants.
import imageio import numpy imA = imageio.imread("imageA.png") U=numpy.array(imA[:,:,0],dtype=numpy.float64)
import numpy import math from pylab import *
La fonction suivante permet de convertir une matrice en matrice image RGB, en effectuant une normalisation et une correction gamma. La correction gamma permet de simuler la réponse de l’oeil pour les figures de diffraction.
def matriceImage(matrice,gamma,rgb): s = matrice.shape a=1.0/gamma; norm=matrice.max() m = numpy.power(matrice/norm,a) im = numpy.zeros((s[0],s[1],3),dtype=float64) im[:,:,0] = rgb[0]*m im[:,:,1] = rgb[1]*m im[:,:,2] = rgb[2]*m return im
Le module de la TFD décroît très rapidement avec la fréquence. Pour mieux visualiser le spectre, on utilise une échelle logarithmique. La fonction suivante génère une image RGB à partir d’une matrice, en appliquant une fonction logarithme :
def matriceImageLog(matrice,rgb): s = matrice.shape m = numpy.log10(1+matrice) min = m.min() max = m.max() m = (m-min)/(max-min) im = numpy.zeros((s[0],s[1],3),dtype=float64) im[:,:,0] = rgb[0]*m im[:,:,1] = rgb[1]*m im[:,:,2] = rgb[2]*m return im
La fonction suivante affiche le spectre. Les axes sont gradués avec les fréquences spatiales. Les dimensions de l’image sont données en argument :
def plotSpectre(image,Lx,Ly): (Ny,Nx,p) = image.shape fxm = Nx*1.0/(2*Lx) fym = Ny*1.0/(2*Ly) imshow(image,extent=[-fxm,fxm,-fym,fym]) xlabel("fx") ylabel("fy")
from TfdImage import * from numpy.fft import fft2, ifft2, fftshift, ifftshift
La transformée de Fourier discrète est calculée avec la fonction numpy.fft.fft2. Le centrage de la fréquence nulle est obtenue avec la fonction numpy.fft.fftshift. On calcule la puissance, c’est-à-dire le module au carré de la TFD et on trace l’image correspondante en rouge.
V=fft2(U) VC = fftshift(V) P = numpy.power(numpy.absolute(VC),2) img = matriceImage(P,2.0,[1.0,0.0,0.0]) plotSpectre(img,200.0,100.0)
On a choisi les dimensions de l’image égales aux nombres de pixels, ce qui donne une fréquence maximale de 0,5 sur le spectre. La fréquence 1 est la fréquence d’échantillonnage, qui correspond à 1 par pixel. On a ainsi un spectre carré (les pixels sont carrés), bien que la TFD soit une matrice rectangulaire (de mêmes dimensions que l’image de départ). Dans le cas présent, l’image a une longueur Lx deux fois plus grande que sa hauteur Ly, donc la résolution du spectre est deux fois plus grande sur l’axe X.
Voici l’affichage en échelle logarithmique :
img = matriceImageLog(P,[1.0,0.0,0.0]) plotSpectre(img,200.0,100.0)
Pour effectuer un filtrage sur la transformée de Fourier de l’image, on doit définir une fonction de transfert. Voici par exemple un filtrage passe-bas :
def transfert(fx,fy,p): fc=p[0] return 1.0/numpy.power(1+(fx*fx+fy*fy)/(fc*fc),2)
Les fréquences sont relatives à la fréquence d’échantillonnage. Par exemple, une fréquence fx=1/8 signifie une période de 8 pixels. Le filtrage dans le domaine fréquentielle consiste à mutliplier la TFD par une matrice H qui est obtenue (en première approche) par échantillonnage de la fonction de transfert.
La fonction suivante effectue la multiplication d’une matrice par l’échantillonnage d’une fonction de transfert. La fréquence nulle de la matrice est centrée.
def matriceFiltre(matrice,transfert,p): s = matrice.shape Nx = s[1] Ny = s[0] nx = Nx/2 ny = Ny/2 Mat = zeros((Ny,Nx),dtype=numpy.complex) for n in range(Nx): for l in range(Ny): fx = float(n-nx-1)/Nx fy = float(l-ny-1)/Ny Mat[l,n] = matrice[l,n]*transfert(fx,fy,p) return Mat
Voyons l’utilisation de cette fonction pour la fonction de transfert définie plus haut. Initialement, il faut créer une matrice H de même taille que l’image et remplie de 1.
H = numpy.ones(VC.shape,dtype=numpy.complex) H = matriceFiltre(H,transfert,[0.1]) VCF = VC*H PF = numpy.power(numpy.absolute(VCF),2) imPF = matriceImageLog(PF,[1.0,0,0]) plotSpectre(imPF,200,100)
Pour obtenir l’image filtrée, il faut tout d’abord restituer l’ordre initial des fréquences avec la fonction ifftshift puis appliquer la transformée de Fourier discrète inverse avec la fonction ifft2 :
VF = ifftshift(VCF) UF=ifft2(VF) imageF = matriceImage(UF,1.0,[1.0,1.0,1.0]) figure() imshow(imageF)
Comme second exemple de filtre, on définit un filtre passe-haut :
def transfert(fx,fy,p): fc=p[0] return 1.0-1.0/numpy.power(1+(fx*fx+fy*fy)/(fc*fc),2)
H = numpy.ones(VC.shape,dtype=numpy.complex) H = matriceFiltre(H,transfert,[0.4]) VCF = VC*H VF = ifftshift(VCF) UF=ifft2(VF) imageF = matriceImage(UF,1.0,[1.0,1.0,1.0]) figure() imshow(imageF)
3.b. Application : filtrage d’une photographie
imA = imageio.imread("imageC.png") U=numpy.array(imA[:,:,0],dtype=numpy.float64) s = U.shape figure(figsize=(12,6)) imshow(U,cmap='gray')
V=fft2(U) VC = fftshift(V) P = numpy.power(numpy.absolute(VC),2) img = matriceImageLog(P,[1.0,1.0,1.0]) figure(figsize=(6,6)) plotSpectre(img,s[1],s[0])
Cette photographie comporte une périodicité dans la direction X à cause de la répétition des fenêtres. Il y a aussi la répétition des barres des balustrades et des jonctions sur le toit en zinc. La périodicité des fenêtres est de 140 pixels, soit une fréquence de 1/140=0,00714, que l’on discerne sur le spectre sous forme de traits verticaux très rapprochées. Les barres des balustrades sont espacées de 6 pixels, soit une fréquence de 1/6=0.16. Le motif de la frise située sous le premier balcon a une période de 11 pixel, soit une fréquence de 1/11=0,091.
Comme exemple de filtre, considérons un filtre coupe-bande gaussien qui agit sur les fréquences horizontales :
def transfert(fx,fy,p): f0 = p[0] sigma = p[1] ux = numpy.absolute(fx)-f0 return 1.0-numpy.exp(-(ux*ux)/(sigma*sigma))
On se sert de ce filtre pour enlever les traits situés à une fréquence horizontale multiple de 0,091, qui correspondent à une périodicité de 11 pixels.
H = numpy.ones(VC.shape,dtype=numpy.complex) H = matriceFiltre(H,transfert,[0.091,0.01]) H = matriceFiltre(H,transfert,[0.091*2,0.01]) H = matriceFiltre(H,transfert,[0.091*3,0.01]) VCF = VC*H PF = numpy.power(numpy.absolute(VCF),2) imPF = matriceImageLog(PF,[1.0,1.0,1.0]) figure(figsize=(6,6)) plotSpectre(imPF,s[1],s[0])
VF = ifftshift(VCF) UF=ifft2(VF) imageF = matriceImage(UF,1.0,[1.0,1.0,1.0]) figure(figsize=(12,6)) imshow(imageF)
Le filtrage a fait disparaitre la frise située sous le premier balcon. Il y aussi une perte globale de netteté.
Il est intéressant de calculer la réponse impulsionnelle de ce filtre. Pour cela, il faut créer une image noire avec un seul pixel blanc.
impuls = numpy.zeros(s) y0 = s[0]/2 x0 = s[1]/2 impuls[y0,x0] = 1.0 V=fft2(impuls) VC = fftshift(V) VCF = VC*H P = numpy.power(numpy.absolute(VCF),2) img = matriceImage(P,1.0,[1.0,1.0,1.0]) figure(figsize=(6,6)) plotSpectre(img,s[1],s[0])
Ce filtre n’agit que sur les fréquences fx, donc sa réponse impulsionnelle n’a qu’une ligne. On peut donc la tracer avec plot :
reponse = ifft2(ifftshift(VCF)) figure() h = reponse[y0,x0-100:x0+100] plot(h,'.')
À partir de la réponse impulsionnelle, on peut obtenir la réponse fréquentielle dans la direction X, comme on le fait pour un signal à une variable :
import scipy.signal w,g=scipy.signal.freqz(h) figure() plot(w/(2*numpy.pi),numpy.absolute(g)) xlabel("f/fe") ylabel("G") grid()
Lorsqu’on dispose de la réponse impulsionnelle, on peut faire le filtrage par convolution :
from scipy.signal import convolve2d h1 = numpy.zeros((1,h.size)) h1[0,:] = h imgF = convolve2d(U,h1,mode='same') figure(figsize=(12,6)) imshow(imgF,cmap='gray')
Le résultat est similaire à celui obtenu en filtrant dans l’espace des fréquences. On constate néanmoins un effet de bord plus important avec le filtrage par convolution, sous forme de bandes verticales à gauche et à droite de l’image.
L’analyse spectrale de l’image est nécessaire pour concevoir des filtres très sélectifs comme le précédent. Pour la réalisation du filtrage, on a le choix entre un filtrage dans le domaine fréquentiel (appelé aussi filtrage par transformée de Fourier) ou un filtrage dans le domaine spatial (filtrage par convolution). Pour comparer l’efficacité de ces deux méthodes, considérons une image carrée comportant n lignes et n colonnes. La transformée de Fourier rapide pour une ligne prend un temps an ln(n). Pour n lignes, cela fait an2 ln(n). Après application de la TFD sur chaque ligne, on doit appliquer la TFD sur chaque colonne. Le temps total pour la FFT de l’image est donc 2an2 ln(n). Le filtrage dans le domaine fréquentiel consiste à multiplier chaque élément de la TFD, ce qui prend un temps bn2. En tenant compte de la FFT inverse, on obtient finalement pour le filtrage par transformée de Fourier un temps :τ1=n2(4aln(n)+b)
Soit Q la taille de la réponse impulsionnelle, c’est-à-dire le masque de convolution. Pour un filtre ayant une certaine réponse par rapport à la fréquence d’échantillonnage, la taille du masque est fixée. Le nombre d’opérations dans un filtrage par convolution est égal au nombre de pixels de l’image multiplié par la taille du masque :τ2=cQn2
Le rapport des deux temps est donc :τ1τ2=4aln(n)+bcQ(1)
Dans la limite des grandes images (n tendant vers l’infini), le filtrage par convolution est beaucoup plus efficace. Néanmoins, pour une taille d’image fixée, le filtrage par transformée de Fourier peut être plus rapide lorsque Q est grand, c’est-à-dire lorsque la réponse impulsionnelle a une taille du même ordre de grandeur que l’image.
Outre sa plus grande efficacité, le filtrage par convolution possède aussi un autre avantage : le masque de convolution peut varier d’une zone à l’autre de l’image. Autrement dit, on peut définir un filtre avec une réponse impulsionnelle locale. À l’inverse, le filtrage dans le domaine fréquentiel ne peut agir que sur l’image prise globalement.
Pour reprendre l’exemple précédent, nous pouvons appliquer le filtrage par convolution seulement sur une partie de l’image :
partie = U[550:580,0:s[1]] partieF = convolve2d(partie,h1,mode='same') img = U.copy() img[550:580,0:s[1]] = partieF figure(figsize=(12,6)) imshow(img,cmap='gray')