Transformada de Fourier. Aplicación al procesamiento de imágenes

Contenidos

Transformada 2D de Fourier

Sea $f(x,y)$ una imagen $M\times N$, con $x=0,1,2,\ldots,M-1$ (columnas) y $y=0,1,2,\ldots,N-1$ (filas).

La Transformada Discreta de Fourier (DFT en inglés) de $f$, que llamaremos $F(m,n)$, viene dada por

$$ F(m,n)= \frac{1}{MN}\sum_{x=0}^{M-1}\sum_{y=0}^{N-1} f(x,y) \exp(-2\pi i(\frac{x}{M}m+\frac{y}{N}n)),$$

con $m=0,1,2,\ldots,M-1$ y $n=0,1,2,\ldots,N-1$.

La región rectangular $M\times N$ definida para $ (m,n)$ se llama dominio de frecuencias, y los valores de $F(m,n)$ se llaman coeficientes de Fourier.

La Transformada Inversa Discreta de Fourier viene dada por

$$ f(x,y)=MN\sum_{m=0}^{M-1}\sum_{n=0}^{N-1} F(m,n) \exp(2\pi i(\frac{m}{M}x+\frac{n}{N}y)),$$

con $x=0,1,2,\ldots,M-1$ y $y=0,1,2,\ldots,N-1$.

Nota. Hay otras formas de normalizar la DFT y su inversa por medio de los factores $1/MN$ y $MN$. El que hemos escogido es el que se usa habitualmente en Python.

Incluso si $f(x,y)$ es real, su transformada $F(m,n)$ es, en general, una función compleja. Para visualizarla utilizaremos el espectro de potencias,

$$P(m,n)= F(m,n) \bar F(m,n),$$

donde $\bar F$ es el conjugado del complejo $F$.

Propiedad: Si la imagen es periódica, es decir, si toma los mismos valores en los bordes superiores en inferiores y en los bordes derecho e izquierdo, entonce, la DFT es periódica también y el espectro de potencia es simétrico respecto del origen: $$ P(m,n)=P(-m,-n). $$

Si la imagen no es periódica, entonces $F$ puede tener discontinuidades. Esto se conoce como el fenómeno de Gibbs.

Gracias a la propiedad mencionada, podemos representar la DFT en un cuadrado centrado en $ (0,0) $, lo cual es más conveniente para su visualización.

La DFT y su inversa se obtienen en la práctica usando el algoritmo de la Transformada Rápida de Fourier (Fast Fourier Transform, FFT). Las funciones correspondientes de Numpy para la transformada directa y la inversa son fft2 y ifft2.

Para calcular el espectro de potencia usamos la función Numpy abs que nos da el módulo del complejo. Para mover el origen de la transformada al centro del rectángulo de frecuencias usamos fftshift. Finalmente, si queremos ver mejor el resultado, usaremos una escala logarítmica.

Ejemplo

In [1]:
from __future__ import division             # hace que se utilice la división en punto flotante
import numpy as np                          
import matplotlib.pyplot as plt             
from PIL import Image                      
from numpy.fft import fft2, fftshift, ifft2

Muestra los gráficos insertados.

In [2]:
%matplotlib inline

Empezamos creando una imagen periódica de tamaño $(601,1201)$. El periodo, $10.5$, aparece en la dirección horizontal. En vertical, la imagen es constante:

In [3]:
hW, hH = 600, 300    
hFrec = 10.5     

# Crea una malla en el cuadrado de dimensiones [0,1)x[0,1)
x = np.linspace( 0, 2*hW/(2*hW +1), 2*hW+1)     # columnas (Anchura)
y = np.linspace( 0, 2*hH/(2*hH +1), 2*hH+1)     # filas    (Altura)

[X,Y] = np.meshgrid(x,y)
A = np.sin(hFrec*2*np.pi*X)

plt.imshow(A, extent=[0,1,0,1], cmap ='gray');
H,W = np.shape(A)   # Dimensiones de la imagen A

Si entendemos que estamos representando la función $A = f(X,Y)$ como una superficie, un corte por un plano paralelo al plano $OXZ$ de dicha superficie sería

In [4]:
xx = np.linspace(0,1,1200)
plt.plot(xx,np.sin(hFrec*2*np.pi*xx))
plt.show()

que es una función periódica de frecuencia $f = 10.5$ Hz y periodo $T=\dfrac{1}{f}$

In [5]:
print "T = ", 1./hFrec

T =  0.0952380952381

El paso siguiente es calcular la DFT centrada en el origen y mostrar la figura del espectro de potencia (su raíz cuadrada).

In [6]:
F = fft2(A)/(W*H)                          
F = fftshift(F)
P = np.abs(F)                            
plt.imshow(P, extent = [-hW,hW,-hH,hH]);

Hacemos zoom (nos acercamos) para ver mejor los valores altos de $P$ a las frecuenccias $\pm hFreq$

In [7]:
plt.imshow(P[hH-25:hH+25,hW-25:hW+25], extent=[-25,25,-25,25]);

Si cortamos la superficie representada por un plano vertical que contenga el eje $OX$ tenemos

In [8]:
plt.plot(range(-25,25),P[hH,hW-25:hW+25])
plt.show()

Otros ejemplos: Con la malla anterior, creamos otras imágenes.

In [9]:
hFrec = 10.5   # Frecuencia horizontal
vFrec = 20.5   # Frecuencia vertical

A1 = np.sin(hFrec*2*np.pi*X) + np.sin(vFrec*2*np.pi*Y)

plt.figure()
plt.imshow(A1, cmap = 'gray', extent=[0,1,0,1]);

Si cortamos esta superficie por el plano $y=0.5$ tenemos la curva

In [10]:
xx = np.linspace(0,1,1200)
yy = 0.5*np.ones(1200)
plt.plot(xx,np.sin(hFrec*2*np.pi*xx) + np.sin(vFrec*2*np.pi*yy))
plt.show()

y si cortamos por el plano $x =0.5$

In [11]:
yy = np.linspace(0,1,1200)
xx = 0.5*np.ones(1200)
plt.plot(yy,np.sin(hFrec*2*np.pi*xx) + np.sin(vFrec*2*np.pi*yy))
plt.show()
In [12]:
F1 = fft2(A1)/(W*H)                          
F1 = fftshift(F1)
P1 = np.abs(F1)

plt.figure()
plt.imshow(P1[hH-25:hH+25,hW-25:hW+25], extent=[-25,25,-25,25]);

Si cortamos la superficie con un plano perpendicular al plano $OXY$ que contiene al eje $OX$

In [13]:
plt.plot(range(-25,25),P1[hH,hW-25:hW+25])
plt.show()

Y con un plano que contiene a $OY$

In [14]:
plt.plot(range(-25,25),P1[hH-25:hH+25,hW])
plt.show()
In [15]:
hFrec = 10.5   # Frecuencia horizontal
vFrec = 20.5   # Frecuencia vertical

A2 = np.sin(hFrec*2*np.pi*X + vFrec*2*np.pi*Y)

plt.figure()
plt.imshow(A2, cmap = 'gray',extent=[0,1,0,1]);

F2 = fft2(A2)/(W*H)                          
F2 = fftshift(F2)
P2 = np.abs(F2)

r = 25
plt.figure()
plt.imshow(P2[hH-r:hH+r,hW-r:hW+r], extent=[-r,r,-r,r]);

Finalmente vamos a usar una imagen real, con una frecuencia alta a lo largo del eje horizontal, lo que se puede ver porque el máximo relativo se alínea con el eje horizontal:

In [16]:
I = Image.open("canes.jpg")
I = I.convert('L')                      # 'L' para convertir a escala de grises
A3 = np.asarray(I, dtype = np.float32)  # Convertimos la imagen I en
                                        # un array Numpy A3 con elementos tipo float32
H,W = np.shape(A3)
plt.imshow(A3, cmap = 'gray');
In [17]:
F3 = fft2(A3)/float(W*H)                          
F3 = fftshift(F3)
P3 = np.abs(F3)

r = 100
mW = int(np.fix(0.5*W))   # Entero que vale, aproximadamente, la mitad de W
mH = int(np.fix(0.5*H))   # Entero que vale, aproximadamente, la mitad de H
plt.figure()
plt.imshow(np.log(1+P3[mH-r:mH+r,mW-r:mW+r]), extent=[-r,r,-r,r]);

Inicio

Aplicación: compresión de imágenes

Una imagen se recupera a partir de los coeficientes $F$ usando la Transformada Inversa de Fourier. Como hemos visto, en general, no todos los coeficientes tienen el mismo valor. Para usar sólo los más importantes, podemos seleccionar aquellos cuyo módulo sea mayor que uno dado llamado umbral $T$. Es decir, hacemos cero todos los coeficientes no contenidos en el conjunto siguiente:

$$ S_T=\left\{0\leq m\leq M-1, 0\leq n \leq N-1 : \left| F(m,n)\right| \geq T\right\} .$$

Llamemos $g(T)$ al número de elementos del array $(m,n)$, contenidos en $S_T$. Por una parte, para $T=0$ tenemos $g(0)= MN$, porque $\left| F(m,n)\right| $ son no negativos. Por otra parte, para $T> \max \left| F(m,n)\right| $ tenemos $g(T)=0$. Observar que $g$ es una función de $T$ que es decreciente. Construiremos un plot de $g$ usando la función count_nonzero, que da el número de elementos no nulos de una matriz.

In [18]:
puntos = 100
Trango = np.linspace(0.05,1,puntos)
g = np.zeros(puntos, dtype = 'float32')
n = 0

for Tcont in Trango:
    g[n] = np.count_nonzero(P3 >= Tcont)
    n += 1   
    
plt.plot(Trango,g);    
plt.xlabel('T (valor del umbral de corte)',fontsize=14)
plt.title('Elementos no nulos en la matriz P3',fontsize=14);

Ahora aplicamos el método del valor umbral, es decir, mandamos a cero los coeficientes menores que un $T$ dado. ¿Cómo escogemos $T$?

  • $T$ debe ser lo bastante grande para eliminar algún coeficiente,
  • $T$ debe ser lo bastante pequeño para quedarse con algún coeficiente

Tomemos, por ejemplo, $T=0.1$ y veamos el resultado:

In [19]:
T = 0.1
c = F3 * (P3 >= T)
fM = ifft2(c)*W*H
plt.imshow(np.abs(fM), cmap = 'gray');

Finalmente, calcularemos el ratio entre el número de coeficientes de la Transformada de Fourier distintos de cero y aquellos que superan el umbral contenidos en la matriz $c$. Así obtenemos una relación de compresión:

In [20]:
out1 = np.count_nonzero(F3)
out2 = np.count_nonzero(c)
print 'Número inicial de elementos distintos de cero = ',out1
print 'Número final de elementos distintos de cero   = ',out2
print 'Relación de compresión                        = ',out1/out2

Número inicial de elementos distintos de cero =  1228800
Número final de elementos distintos de cero   =  5099
Relación de compresión                        =  240.988429104

Ejercicio

Tenemos una imagen de tamaño $W\times H$, y queremos comprimirla de acuerdo con un ratio dado $r$ manteniendo los coeficientes de la transformada de Fourier más significativos, es decir, queremos encontrar el umbral, $T$, tal que el ratio de compresión sea $r$.

Dicho de otra manera, queremos calcular $T$ de forma que se cumpla que $ \displaystyle\frac{g(T)}{WH}=r $.

Para conseguirlo, tenemos que calcular la raíz de la función $$h(T)=g(T) -rWH.$$

Para escribir el programa compresion.py:

  • Cargar la imagen cameraman.tif and definir $r=1/50$.
  • Calcular su transformada de Fourier, $F$, junto con el máximo y el mínimo de su módulo, $maxF,minF$.
  • La función $ h(T) $ se podría definir como en los ejemplos anteriores, usando np.count_nonzero. Para hacerse una idea de cómo es $h$, dibujar esta función para el intervalo $[minF,0.6]$.
  • Como $h(T)$ es decreciente, podemos buscar el $T$ para el que cambia el signo de $h$ (y por lo tanto, dond $h$ se aproxima a cero) comprobando sus signo desde $T=minF$ hasta $T=maxF$ a intervalos regulares. Una vez que hemos determinado gráficamente el valor de $T$ cerca del cula cambia el signo de $h$, podemos definir la matriz $c$ obtenida a partir de la matriz $F$ aplicando el umbral y la imagen comprimida fM como hicimos en este apartado. Dibujar las imágenes original y comprimida.
  • Comprobar que el ratio de compresión está próximo a $r$ (out1/out2).
  • Adaptar cualquiera de los programas del método de Bisección para calcular una aproximación de la única raíz de $h$. Por cierto, ¿por qué es única?. Comprobar que las dos soluciones obtenidas son parecidas.
  • Finalmente, transformar el programa en función con argumentos de entradas la imagen y la relación de compresión y con argumentos de salida, el umbral T y una figura de la imagen original y la comprimida.
In [21]:
%run compresion.py

Relación de compresión =  50.0083937428

Inicio

Aplicación: orientación de vehículos

En esta aplicación, tenemos un vehículo que queremos orientar de forma que siga líneas verticales.

Primero, tomamos la imagen original $A$ y la rotamos:

In [24]:
A4 = (A - A.min())/(A.max() - A.min())    # normalizamos a [0,1]
I = Image.fromarray(np.uint8(A4*255))     # convertimos a imagen PIL

rot_angulo = -65
B = I.rotate(rot_angulo, expand = True)   # rotamos sin recortar

B4 = np.asarray(B,dtype=np.float32)
sY,sX = np.shape(B4)

plt.figure()
plt.imshow(B,cmap = 'gray');
In [25]:
F4 = fft2(B4)/(sX*sY)
F4 = fftshift(F4)
P4 = np.abs(F4)

cX = int(np.fix(sX/2))                    # punto medio de la imagen
cY = int(np.fix(sY/2))

plt.figure()
plt.imshow(P4[cY-25:cY+25,cX-25:cX+25], extent=[-25,25,-25,25]);

Como hay una región grande de valor constante (negro), el valor del espectro de potencia en el origen es alto, haciendo que no resalten los valores de las frecuencias importantes. Eliminamos este efecto restando el valor medio de la imagen.

In [26]:
F4 = fft2(B4-B4.mean())/(sX*sY)
F4 = fftshift(F4)
P4 = np.abs(F4)

plt.figure()
plt.imshow(P4[cY-25:cY+25,cX-25:cX+25], extent=[-25,25,-25,25]);

Calculamos las coordenadas de los coeficientes de más alta frecuencia y calculamos el ángulo

In [27]:
indices = np.where(P4 == np.max(P4))
maxY = (indices[0][0]-cY)/sY 
maxX = (indices[1][0]-cX)/sX 

alpha=np.arctan(maxY/maxX)*180/np.pi

print "alpha = ", alpha

alpha =  65.1189251526

Ahora rotamos y recuperamos la imagen original

In [28]:
C=B.rotate(alpha)
plt.figure()
plt.imshow(C,cmap = 'gray');

Ejercicio

Crear una función que tenga como:

  • Entrada: una imagen (con un comportamiento periódico)
  • Salida: el ángulo para enderezar la imagen, la imagen enderezada, y el espectro de potencia

Aplicarlo a las siguientes imágenes:

In [29]:
%run angulo.py
In [30]:
I5 = Image.open("surcos1.jpg")

alpha5, C5, P5 = angulo(I5)

print alpha5
plt.figure()
plt.imshow(I5);
plt.figure()
plt.imshow(C5);
plt.figure()
plt.imshow(P5);

-30.865815124

En este ejemplo, el algoritmo funcionó bien. Fijarse en que los valores altos de P5 están lejos del origen.

Sin embargo, en el ejemplo siguiente, el algoritmo no funciona. ¿Por qué? ¿Cómo podría arreglarse?

In [31]:
I6 = Image.open("surcos2.jpg")
alpha6, C6, P6 = angulo(I6)

print alpha6
plt.figure()
plt.imshow(I6);
plt.figure()
plt.imshow(C6);
plt.figure()
plt.imshow(P6);

-90.0

<string>:26: RuntimeWarning: divide by zero encountered in double_scalars

Inicio

Imágenes