TFD Image

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 (f_{x}, f_{y}) = \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} u (x, y) exp (- i 2 π (f_{x} x + f_{y} y)) d x d y

On considère à présent un échantillonnage de l'image sur un domaine rectangulaire de dimensions (L_x,L_y) centré en (x,y)=(0,0) et comportant (N_x,N_y) points. La matrice U est définie par :

\begin{matrix} x_{k} = - \frac{L_{x}}{2} + k \frac{L_{x}}{N_{x}} 0 \leq k \leq N_{x} - 1 \\ y_{m} = - \frac{L_{y}}{2} + m \frac{L_{y}}{N_{y}} 0 \leq m \leq N_{y} - 1 \\ U_{m, k} = u (x_{k}, y_{m}) \end{matrix}

La matrice U a donc N_x colonnes et N_y lignes.

Considérons une approximation de la transformée de Fourier obtenue par la méthode des rectangles sur le domaine considéré :

S (f_{x}, f_{y}) ≃ \frac{L_{x} L_{y}}{N_{x} N_{y}} exp (i π (f_{x} L_{x} + f_{y} L_{y})) \sum_{k = 0}^{N_{x} - 1} \sum_{m = 0}^{N_{y} - 1} U_{m, k} exp (- i 2 π f_{x} k \frac{L_{x}}{N_{x}}) exp (- i 2 π f_{y} m \frac{L_{y}}{N_{y}})

En pratique, on calcule cette approximation pour les fréquences spatiales suivantes :

\begin{matrix} f_{x, n} = \frac{n}{L_{x}} 0 \leq n \leq N_{x} - 1 \\ f_{y, l} = \frac{l}{L_{y}} 0 \leq l \leq N_{y} - 1 \end{matrix}

La transformée de Fourier discrète (TFD) associe à la matrice U une matrice V de mêmes dimensions, définie par :

V_{n, l} = \frac{1}{N_{x} N_{y}} \sum_{k = 0}^{N_{x} - 1} \sum_{m = 0}^{N_{y} - 1} U_{m, k} exp (- i 2 π k \frac{n}{N_{x}}) exp (- i 2 π m \frac{l}{N_{y}})

Remarque : les TFD calculées par les fonctions FFT des logiciels omettent parfois l'inverse de N_xN_y en facteur.

L'approximation de l'intégrale est :

S (f_{x, n}, f_{y, l}) ≃ L_{x} L_{y} exp (i π (n + l)) V_{n, l}

La fréquence d'échantillonnage est constituée par la paire :

f_{e} = (\frac{N_{x}}{L_{x}}, \frac{N_{y}}{L_{y}})

On vérifie facilement la propriété suivante de la TFD:

\begin{matrix} V_{n + N_{x}, l} = V_{n, l} \\ V_{n, l + N_{y}} = V_{n, l} \end{matrix}

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 N_x et N_y sont pairs. La TFD donne une approximation de la transformée de Fourier pour les fréquences suivantes :

\begin{matrix} f_{x, n} = \frac{n}{L_{x}} 0 \leq n \leq \frac{N_{x}}{2} \\ f_{y, l} = \frac{l}{L_{y}} 0 \leq l \leq \frac{N_{y}}{2} \end{matrix}

De plus, la périodicité de la TFD permet d'accéder aux opposées de ces fréquences. En effet :

\begin{matrix} V_{- n, l} = V_{- n + N_{x}, l} \\ V_{n, - l} = V_{n, - l + N_{y}} \end{matrix}

Les N_x indices n correspondent aux fréquences suivantes :

0, \frac{1}{L_{x}}, \frac{2}{L_{x}}, \dots, \frac{\frac{N_{x}}{2}}{L_{x}}, - \frac{\frac{N_{x}}{2} - 1}{L_{x}}, \dots, - \frac{2}{L_{x}}, - \frac{1}{L_{x}}

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 :

(\frac{N_{y}}{2} - 1, \frac{N_{x}}{2} - 1)

Pour l'axe x, il y a donc $\frac{N_{x}}{2} + 1$ indices correspondant à des fréquences positives (ou nulle) et $\frac{N_{x}}{2} - 1$ indices correspondant à des fréquences strictement négatives.

Pour construire la matrice VC, il faut utiliser la transformation d'indice suivante :

\begin{matrix} n_{c} = \frac{N_{x}}{2} - 1 + n si 0 \leq n \leq \frac{N_{x}}{2} \\ n_{c} = - \frac{N_{x}}{2} - 1 + n si \frac{N_{x}}{2} + 1 \leq n \leq N_{x} - 1 \end{matrix}

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 200x100 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)

TfdImage.py

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)

plotA.pdf

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 L_x deux fois plus grande que sa hauteur L_y, 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)

plotA1.pdf

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 f_x=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)

plotB.pdf

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)

plotC.pdf

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)

plotD.pdf

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')

plotE.pdf

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])

plotF.pdf

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])

plotG.pdf

VF = ifftshift(VCF)
UF=ifft2(VF)
imageF = matriceImage(UF,1.0,[1.0,1.0,1.0])
figure(figsize=(12,6))
imshow(imageF)

plotH.pdf

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])

plotI.pdf

Ce filtre n'agit que sur les fréquences f_x, 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,'.')

plotJ.pdf

À 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()

plotJ1.pdf

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')

plotK.pdf

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 an² 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 2an² ln(n). Le filtrage dans le domaine fréquentiel consiste à multiplier chaque élément de la TFD, ce qui prend un temps bn². En tenant compte de la FFT inverse, on obtient finalement pour le filtrage par transformée de Fourier un temps :

τ_{1} = n^{2} (4 a ln (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} = c Q n^{2}

Le rapport des deux temps est donc :

\frac{τ_{1}}{τ_{2}} = \frac{4 a ln (n) + b}{c Q} (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')

plotL.pdf