5/4/2021 Fourier2D

Transformada de Fourier. Aplicación al

procesamiento de imágenes

Contenidos
Transformada 2D de Fourier
Aplicación: compresión de imágenes
Aplicación: orientación de vehículos

Transformada 2D de Fourier

[Link] 1/23
5/4/2021 Fourier2D

Sea f (x, y) una imagen M × N , con x = 0, 1, 2, … , M − 1 (columnas) y

y = 0, 1, 2, … , N − 1 (filas).

La Transformada Discreta de Fourier (DFT en inglés) de f, que llamaremos

F (m, n), viene dada por

M −1 N −1
1 x y
F (m, n) = ∑ ∑ f (x, y) exp(−2πi( m + n)),
MN M N
x=0 y=0

con m = 0, 1, 2, … , M − 1 y n = 0, 1, 2, … , N − 1.

La región rectangular M × 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

M −1 N −1
m n
f (x, y) = M N ∑ ∑ F (m, n) exp(2πi( x + y)),
M N
m=0 n=0

con x = 0, 1, 2, … , M − 1 y y = 0, 1, 2, … , N − 1 .

Nota. Hay otras formas de normalizar la DFT y su inversa por medio de los factores
1/M N y M N . 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)F (m, n),

donde ¯
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

[Link] 2/23
5/4/2021 Fourier2D

In [1]:

from __future__ import division # hace que se utilice la división en p

unto flotante
import numpy as np
import [Link] as plt
from PIL import Image
from [Link] 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 = [Link]( 0, 2*hW/(2*hW +1), 2*hW+1) # columnas (Anchura)
y = [Link]( 0, 2*hH/(2*hH +1), 2*hH+1) # filas (Altura)

[X,Y] = [Link](x,y)
A = [Link](hFrec*2*[Link]*X)

[Link](A, extent=[0,1,0,1], cmap ='gray');

H,W = [Link](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

[Link] 3/23
5/4/2021 Fourier2D

In [4]:

xx = [Link](0,1,1200)
[Link](xx,[Link](hFrec*2*[Link]*xx))
[Link]()

1
que es una función periódica de frecuencia f = 10.5 Hz y periodo T =
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 = [Link](F)
[Link](P, extent = [-hW,hW,-hH,hH]);

Hacemos zoom (nos acercamos) para ver mejor los valores altos de P a las frecuenccias
±hF req

[Link] 4/23
5/4/2021 Fourier2D

In [7]:

[Link](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]:

[Link](range(-25,25),P[hH,hW-25:hW+25])
[Link]()

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

[Link] 5/23
5/4/2021 Fourier2D

In [9]:

hFrec = 10.5 # Frecuencia horizontal

vFrec = 20.5 # Frecuencia vertical

A1 = [Link](hFrec*2*[Link]*X) + [Link](vFrec*2*[Link]*Y)

[Link]()
[Link](A1, cmap = 'gray', extent=[0,1,0,1]);

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

In [10]:

xx = [Link](0,1,1200)
yy = 0.5*[Link](1200)
[Link](xx,[Link](hFrec*2*[Link]*xx) + [Link](vFrec*2*[Link]*yy))
[Link]()

y si cortamos por el plano x = 0.5

[Link] 6/23
5/4/2021 Fourier2D

In [11]:

yy = [Link](0,1,1200)
xx = 0.5*[Link](1200)
[Link](yy,[Link](hFrec*2*[Link]*xx) + [Link](vFrec*2*[Link]*yy))
[Link]()

In [12]:

F1 = fft2(A1)/(W*H)
F1 = fftshift(F1)
P1 = [Link](F1)

[Link]()
[Link](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

[Link] 7/23
5/4/2021 Fourier2D

In [13]:

[Link](range(-25,25),P1[hH,hW-25:hW+25])
[Link]()

Y con un plano que contiene a OY

In [14]:

[Link](range(-25,25),P1[hH-25:hH+25,hW])
[Link]()

[Link] 8/23
5/4/2021 Fourier2D

In [15]:

hFrec = 10.5 # Frecuencia horizontal

vFrec = 20.5 # Frecuencia vertical

A2 = [Link](hFrec*2*[Link]*X + vFrec*2*[Link]*Y)

[Link]()
[Link](A2, cmap = 'gray',extent=[0,1,0,1]);

F2 = fft2(A2)/(W*H)
F2 = fftshift(F2)
P2 = [Link](F2)

r = 25
[Link]()
[Link](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:

[Link] 9/23
5/4/2021 Fourier2D

In [16]:

I = [Link]("[Link]")
I = [Link]('L') # 'L' para convertir a escala de grises
A3 = [Link](I, dtype = np.float32) # Convertimos la imagen I en
# un array Numpy A3 con elementos tipo flo
at32
H,W = [Link](A3)
[Link](A3, cmap = 'gray');

In [17]:

F3 = fft2(A3)/float(W*H)
F3 = fftshift(F3)
P3 = [Link](F3)

r = 100
mW = int([Link](0.5*W)) # Entero que vale, aproximadamente, la mitad de W
mH = int([Link](0.5*H)) # Entero que vale, aproximadamente, la mitad de H
[Link]()
[Link]([Link](1+P3[mH-r:mH+r,mW-r:mW+r]), extent=[-r,r,-r,r]);

Inicio

[Link] 10/23
5/4/2021 Fourier2D

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:

ST = {0 ≤ m ≤ M − 1, 0 ≤ n ≤ N − 1 : |F (m, n)| ≥ T } .

Llamemos g(T ) al número de elementos del array (m, n) , contenidos en ST . Por una
parte, para T = 0 tenemos g(0) = M N , porque |F (m, n)| son no negativos. Por otra
parte, para T > max |F (m, n)| 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 = [Link](0.05,1,puntos)
g = [Link](puntos, dtype = 'float32')
n = 0

for Tcont in Trango:

g[n] = np.count_nonzero(P3 >= Tcont)
n += 1

[Link](Trango,g);
[Link]('T (valor del umbral de corte)',fontsize=14)
[Link]('Elementos no nulos en la matriz P3',fontsize=14);

[Link] 11/23
5/4/2021 Fourier2D

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
[Link]([Link](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

[Link] 12/23
5/4/2021 Fourier2D

Ejercicio

Tenemos una imagen de tamaño W × 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.

g(T )
Dicho de otra manera, queremos calcular T de forma que se cumpla que = r.
WH

Para conseguirlo, tenemos que calcular la raíz de la función

h(T ) = g(T ) − rW H .

Para escribir el programa [Link]:

Cargar la imagen [Link] 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 [Link]

Relación de compresión = 50.0083937428

[Link] 13/23
5/4/2021 Fourier2D

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 - [Link]())/([Link]() - [Link]()) # normalizamos a [0,1]

I = [Link](np.uint8(A4*255)) # convertimos a imagen PIL

rot_angulo = -65
B = [Link](rot_angulo, expand = True) # rotamos sin recortar

B4 = [Link](B,dtype=np.float32)
sY,sX = [Link](B4)

[Link]()
[Link](B,cmap = 'gray');

[Link] 14/23
5/4/2021 Fourier2D

In [25]:

F4 = fft2(B4)/(sX*sY)
F4 = fftshift(F4)
P4 = [Link](F4)

cX = int([Link](sX/2)) # punto medio de la imagen

cY = int([Link](sY/2))

[Link]()
[Link](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.

[Link] 15/23
5/4/2021 Fourier2D

In [26]:

F4 = fft2([Link]())/(sX*sY)
F4 = fftshift(F4)
P4 = [Link](F4)

[Link]()
[Link](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 = [Link](P4 == [Link](P4))

maxY = (indices[0][0]-cY)/sY
maxX = (indices[1][0]-cX)/sX

alpha=[Link](maxY/maxX)*180/[Link]

print "alpha = ", alpha

alpha = 65.1189251526

Ahora rotamos y recuperamos la imagen original

[Link] 16/23
5/4/2021 Fourier2D

In [28]:

C=[Link](alpha)
[Link]()
[Link](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 [Link]

[Link] 17/23
5/4/2021 Fourier2D

In [30]:

I5 = [Link]("[Link]")

alpha5, C5, P5 = angulo(I5)

print alpha5
[Link]()
[Link](I5);
[Link]()
[Link](C5);
[Link]()
[Link](P5);

[Link] 18/23
5/4/2021 Fourier2D

-30.865815124

[Link] 19/23
5/4/2021 Fourier2D

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?

[Link] 20/23
5/4/2021 Fourier2D

In [31]:

I6 = [Link]("[Link]")
alpha6, C6, P6 = angulo(I6)

print alpha6
[Link]()
[Link](I6);
[Link]()
[Link](C6);
[Link]()
[Link](P6);

[Link] 21/23
5/4/2021 Fourier2D

-90.0
<string>:26: RuntimeWarning: divide by zero encountered in double_sc
alars

Inicio

[Link] 22/23
5/4/2021 Fourier2D

Imágenes
([Link]

[Link] 23/23