16. Procesado óptico de la información

Autor:Luis Miguel Sánchez Brea
Revisor:José Luis Vilas Prieto

Se entiende por procesado de imágenes a cualquier transformación que pueda ser realizada en la imagen ya sea esta digital u óptica. El procedimiento para cualquiera de los dos tipos de imágenes es el mismo, y en gran parte del procesado de imágenes implica pasar al espacio de frecuencias es decir, realizar la transformada de Fourier. Se puede probar como un haz colimado focaliza en el foco de la lente, obteniéndose en dicho punto la transformada de Fourier del objeto. La transformada de Fourier de la imagen, contiene la misma información pero en el espacio de frecuencias, ubicándose en los bajos órdenes de difracción las bajas frecuencias los cuales se corresponden con la imagen general del objeto. Así mismo todos aquellos puntos alejados del origen, esto es los ordenas más altos de difracción representan aquellos pequeños detalles de la imagen, por tanto si se eliminan se pierde definición pero se mantiene la estructura de la imagen. El objetivo de los filtros de frecuencias consiste en obtener la transformada de Fourier de una imagen y tapar la región de la transformada con las frecuencias a eliminar, finalmente se realizará la transformada inversa para recuperar la imagen filtrada.

Herramientas utilizadas en este tutorial:

  • scipy: manipulación de arrays de datos.
  • matplotlib.pyplot: generación de dibujos similar a Matlab.
  • fuentesXY: Definición de haces de luz.
  • camposXY: Diversas funciones de tratamiento de imágenes e inicialización.

Archivos necesarios para el procesado de imágenes:

Contenidos de los archivos:

  • procesadoOptico.py:

    • procesar(). Realiza el procesamiento de la imagen, desde la definición del campo inicial, el plano de Fourier, la aplicación del filtro y la obtención del campo filtrado.
    • dibujar(). Representa el proceso de filtrado, intensidad del objeto, su fase, el plano de Fourier, el filtro de intensidad, el filtro de fase y el campo filtrado.
    • dibujarCampoFiltrado (). Representa el campo filtrado.

16.1. Funcionamiento de la clase procesadoOptico

A continuación se muestran una serie de funciones necesarias para el procesado óptico.

Al igual que en un experimento físico se dispone de una fuente, un objeto y un filtro, en simulación es preciso definir estos ‘parámetros’ del experimento para realizar un procesamiento óptico de la imagen. Así se inicia la clase con la función __init__ especificando tales condiciones. Esta función dispone de un if, de modo que en ausencia de cualquiera de estos tres elementos no se realice filtrado alguno.

	def __init__(self, fuente, objeto, filtro):
		"""Se inicia un nuevo experimento"""
		#Se definen los parametros de describen la fuente, el objeto y el filtro
		self.fuente = fuente
		self.objeto = objeto
		self.filtro = filtro

		#Si alguno de ellos falla no puede realizarse procesamiento alguno
		if fuente == None or objeto == None or filtro == None:
			self.planoFourier = None
			self.campoFiltrado = None
		else:
			self.procesar()

Así mismo se define la función de procesamiento, la cual calcula el campo producto del objeto con la iluminación. Puesto que se trabaja en el plano focal y en este se obtiene la transformada de Fourier del objeto bajo un haz colimado, se calcula la transformada de Fourier. Finalmente se aplica el filtro realizando el producto de la transformada de Fourier por la función filtro. Finalmente, se deshace la transformada, para obtener el campo filtrado.

	def procesar(self):
		#Campo inicial
		campoInicial = self.fuente * self.objeto

		#Transformada de Fourier del campo inicial
		planoFourier = campoInicial.fft(quitar0 = False, shift = True)

		#Aplicacion del filtro
		planoFourierFiltrado = planoFourier * self.filtro

		#Transformada de Fourier para obtener el campo filtrado
		campoFiltrado = planoFourierFiltrado.fft(quitar0 = True, shift = False)

		self.planoFourier = planoFourier
		self.campoFiltrado = campoFiltrado

Resta realizar una representación de los campos a fin de comprar el incidente y el filtrado, esto se lo hace la función dibujar. Así se genera un array de imágenes donde se ubicarán el objeto, su fase, la transformada de Fourier, el filtro de intensidad, el filtro de fase y finalmente el campo filtrado. Así se selecciona la componente del array donde representar la imagen (objeto, fase, transformada de Fourier...). Finalmente se dibuja con plt, el resto de comandos determinan el título, los ejes, entre otros.

	def dibujar(self, nombreArchivo = ''):
		#Se crea uan figura
		id = plt.figure(figsize = (9, 6), facecolor = 'w', edgecolor = 'k')

		#objeto intensidad		
		#Seleccion de la region dela figura donde se pretende representar
		plt.subplot(2, 3, 1)
		#Intensidad
		dibujo = abs(self.objeto.u) ** 2
		#Representacion
		so=self.objeto
		h1 = plt.imshow(dibujo, interpolation = 'bilinear', aspect = 'auto',
		 	 origin = 'lower', extent = [so.x.min(), so.x.max(), so.y.min(),\
													 so.y.max()])
		#Titulo
		plt.title('objeto intensidad', fontsize = 18)
		#Seleccion del mapa de color
		h1.set_cmap("gist_heat") #RdBu
		#Los ejes no se representan
		plt.axis('off')

		#objeto fase		
		plt.subplot(2, 3, 2)
		dibujo = sp.angle(self.objeto.u)
		#Representacion
		h1 = plt.imshow(dibujo, interpolation = 'bilinear', aspect = 'auto',
		 	 origin = 'lower', extent = [so.x.min(), so.x.max(), so.y.min(),\
													 so.y.max()])
		#Titulo
		plt.title('objeto fase', fontsize = 18)
		#Seleccion del mapa de color
		h1.set_cmap("RdBu") #RdBu
		#Los ejes no se representan
		plt.axis('off')

		#filtro	intensidad	
		plt.subplot(2, 3, 4)
		#Intensidad
		dibujo = abs(self.filtro.u) ** 2
		#Representacion
		sf=self.filtro
		h1 = plt.imshow(dibujo, interpolation = 'bilinear', aspect = 'auto',\
								origin = 'lower',  \
			extent = [sf.x.min(), sf.x.max(), sf.y.min(), sf.y.max()])
		#Titulo
		plt.title('filtro intensidad', fontsize = 18)
		#Seleccion del mapa de color
		h1.set_cmap("gist_heat") #RdBu
		#Los ejes no se representan
		plt.axis('off')

		#filtro fase		
		plt.subplot(2, 3, 5)
		dibujo = sp.angle(self.filtro.u)
		#Representacion
		h1 = plt.imshow(dibujo, interpolation = 'bilinear', aspect = 'auto',
		 	 origin = 'lower', extent = [sf.x.min(), sf.x.max(), sf.y.min(),\
													 sf.y.max()])
		#Titulo
		plt.title('filtro fase', fontsize = 18)
		#Seleccion del mapa de color
		h1.set_cmap("RdBu") #RdBu
		#Los ejes no se representan
		plt.axis('off')

		#planoFourier		
		plt.subplot(2, 3, 3)
		dibujo = sp.log(abs(self.planoFourier.u) ** 2 + 1)
		#Representacion
		pf=self.planoFourier
		h1 = plt.imshow(dibujo, interpolation = 'bilinear', aspect = 'auto',
		 	 origin = 'lower', extent = [pf.x.min(), pf.x.max(), pf.y.min(),\
													 pf.y.max()])
		plt.title('plano Fourier', fontsize = 18)
		#Seleccion del mapa de color
		h1.set_cmap("gist_heat") #RdBu
		#Los ejes no se representan
		plt.axis('off')

		#campoFiltrado		
		plt.subplot(2, 3, 6)
		dibujo = sp.fliplr(sp.flipud(abs(self.campoFiltrado.u) ** 2))
		#Representacion
		cf=self.campoFiltrado
		h1 = plt.imshow(dibujo, interpolation = 'bilinear', aspect = 'auto',
		 	 origin = 'lower', extent = [cf.x.min(), cf.x.max(), cf.y.min(),\
													 cf.y.max()])
		#Titulo
		plt.title('campo Filtrado', fontsize = 18)
		#Seleccion del mapa de color
		h1.set_cmap("gist_heat") #RdBu
		#Los ejes no se representan
		plt.axis('off')

		if not nombreArchivo == '':
			plt.savefig(nombreArchivo, dpi = 100, bbox_inches = 'tight',\
									pad_inches = 0.1) 



	

16.2. Filtro pasa baja

Un filtro pasa es aquel que atenúa las altas frecuencias y permitiendo únicamente el paso de las bajas frecuencias. De este modo, un filtro pasa baja conserva los detalles más generales de la imagen, eliminando los pequeños detalles. Esto se logra en el plano de Fourier diafragmando los mayores órdenes de difracción, por ser estos los correspondientes a altas frecuencias. En las imágenes de plano Fourier y filtro intensidad se puede apreciar este hecho. La abertura del filtro de intensidad se ubica en el plano de Fourier, atravesando el filtro las bajas frecuencias.

#!/usr/bin/env python
# -*- coding: utf-8 -*-
#------------------------------------
# Autor:	Luis Miguel Sanchez Brea
# Fecha		2013/10/24 (version 1.0)
# Licencia:	GPL
#-------------------------------------
"""Simulaciones relacionadas con la leccion de filtrado optico"""

#Importacion de paquetes y librerias
import locale
locale.setlocale(locale.LC_NUMERIC, 'C')
import sys
sys.path.append('../')
import scipy as sp
from procesadoOptico import *
from fuentesXY import fuenteXY
from mascarasXY import mascaraXY

#Tamano
tamano = 1000 * um
x0 = sp.linspace(-tamano/2 * um, tamano/2 * um, 512)
y0 = sp.linspace(-tamano/2 * um, tamano/2 * um, 512)

#Definicion de la fuente
fuente1 = fuenteXY(x=x0, y=y0)
fuente1.onda_plana(A=1, theta=0*grados, phi=0*grados)

def test_filtroPasaBaja(): #amplitud
	#Definicion del objeto
	objeto = mascaraXY(x=x0, y=y0)
	objeto.imagen(nombre='lena.png', normalizar=True, canal=0, tamanoImagen=False, invertir=False, angulo=0)
	#Definicion del filtro
	filtro = mascaraXY(x=x0, y=y0)
	filtro.circulo(r0=(0*um,0*um), radius=(20*um,20*um), angulo=0*grados)
	#Campo filtrado, procesado
	p1=procesadoOptico(fuente1, objeto, filtro)
	p1.dibujar2(idfigs=(1,0,2,3,0,4))
	p1.dibujarCampoFiltrado()

if __name__ == '__main__':
	test_filtroPasaBaja()
	plt.show()

(Source code)

16.3. Eliminación de ruido

Una aplicación interesante del filtro pasa baja es la eliminación de ruido en una imagen. Para comprobarlo, tomamos la foto de lena.png y añadimos un ruido aleatorio.

#!/usr/bin/env python
# -*- coding: utf-8 -*-
#------------------------------------
# Autor:	Luis Miguel Sanchez Brea
# Fecha		2013/10/24 (version 1.0)
# Licencia:	GPL
#-------------------------------------
"""Simulaciones relacionadas con la leccion de filtrado optico"""

#Carga de paquetes y librerias
import locale
locale.setlocale(locale.LC_NUMERIC, 'C')
import sys
sys.path.append('../')
import scipy as sp
from procesadoOptico import *
from fuentesXY import fuenteXY
from mascarasXY import mascaraXY

tamano = 1000 * um
x0 = sp.linspace(-tamano/2 * um, tamano/2 * um, 512)
y0 = sp.linspace(-tamano/2 * um, tamano/2 * um, 512)

def test_filtroEliminarRuido():
	#Definicion de la fuente
	fuente1 = fuenteXY(x=x0, y=y0)
	fuente1.onda_plana(A=1, theta=0*grados, phi=0*grados)
	#Definicion del objeto
	objeto = mascaraXY(x=x0, y=y0)
	objeto.imagen(nombre='lena.png', normalizar=True, canal=0, tamanoImagen=False, invertir=False, angulo=0)
	objeto.u=objeto.u+sp.rand(len(x0),len(y0))
	#Definicion del filtro
	filtro = mascaraXY(x=x0, y=y0)
	filtro.circulo(r0=(0*um,0*um), radius=(50*um,50*um), angulo=0*grados)
	#Campo filtrado, procesado de la imagen
	p1=procesadoOptico(fuente1, objeto, filtro)
	p1.dibujar2(idfigs=(1,0,2,3,0,4))
	p1.dibujarCampoFiltrado()

test_filtroEliminarRuido()
plt.show()

(Source code)

16.4. Filtro pasa alta

Este tipo de filtro es el opuesto al filtro pasa baja, conserva las altas frecuencias y atenúa las bajas. De este modo, los detalles más gruesos de la imagen ser pierden, pero los pequeños se mantienen. Por ejemplo un conjunto de franjas equiespaciadas muy separadas entre sí, se atenuarían respecto de otro conjunto de franjas poco espaciadas entre sí. La consecución de este tipo de filtro se logra tapando los órdenes más bajos de difracción en el plano de Fourier, atravesando únicamente el filtro, los más altos. Nuevamente, véase la imagen, como el plano de Fourier el filtro de intensidad tapa dichos órdenes, de modo que el campo filtrado muestra los pequeños detalles de la imagen (altas frecuencias).

#!/usr/bin/env python
# -*- coding: utf-8 -*-
#------------------------------------
# Autor:	Luis Miguel Sanchez Brea
# Fecha		2013/10/24 (version 1.0)
# Licencia:	GPL
#-------------------------------------
"""Simulaciones relacionadas con la leccion de filtrado optico"""

#Carga de librerias y paquetes
import sys
sys.path.append('../')
from procesadoOptico import *
from fuentesXY import fuenteXY
from mascarasXY import mascaraXY

#Tamano
tamano = 1000 * um
x0 = sp.linspace(-tamano/2 * um, tamano/2 * um, 512)
y0 = sp.linspace(-tamano/2 * um, tamano/2 * um, 512)

def test_filtroPasaAlta():
	#Definicion de la fuente
	fuente1 = fuenteXY(x=x0, y=y0)
	fuente1.onda_plana(A=1, theta=0*grados, phi=0*grados)
	#Definicion del Objeto
	objeto = mascaraXY(x=x0, y=y0)
	objeto.imagen(nombre='lena.png', normalizar=True, canal=0, tamanoImagen=False, invertir=False, angulo=0)
	#Definicion del filtro
	filtro = mascaraXY(x=x0, y=y0)
	filtro.circulo(r0=(0*um,0*um), radius=(50*um,50*um), angulo=0*grados)
	filtro.inversaAmplitud()
	#Campo filtrado, procesado de la imagen
	p1=procesadoOptico(fuente1, objeto, filtro)
	p1.dibujar2(idfigs=(1,0,2,3,0,4))
	p1.dibujarCampoFiltrado()


test_filtroPasaAlta()
plt.show()

(Source code)

veamos otro ejemplo de detección de bordes, donde éstos están más definidos.

#!/usr/bin/env python
# -*- coding: utf-8 -*-
#------------------------------------
# Autor:	Luis Miguel Sanchez Brea
# Fecha		2013/10/24 (version 1.0)
# Licencia:	GPL
#-------------------------------------
"""Simulaciones relacionadas con la leccion de filtrado optico"""

#Carga de librerias y paquetes
import sys
sys.path.append('../')
from procesadoOptico import *
from fuentesXY import fuenteXY
from mascarasXY import mascaraXY

tamano = 1000 * um
x0 = sp.linspace(-tamano/2 * um, tamano/2 * um, 512)
y0 = sp.linspace(-tamano/2 * um, tamano/2 * um, 512)

def test_filtroPasaAlta2():
	#Definicion de la fuente
	fuente1 = fuenteXY(x=x0, y=y0)
	fuente1.onda_plana(A=1, theta=0*grados, phi=0*grados)
	#Definicion del bjeto
	objeto = mascaraXY(x=x0, y=y0)
	objeto.cuadrado(r0=(0,0), size=(250*um,250*um), angulo=45*grados)
	#Definicion del filtro
	filtro = mascaraXY(x=x0, y=y0)
	filtro.circulo(r0=(0*um,0*um), radius=(40*um,40*um), angulo=0*grados)
	filtro.inversaAmplitud()
	#Campo filtrado, procesamiento de la imagen
	p1=procesadoOptico(fuente1, objeto, filtro)
	p1.dibujar2(idfigs=(1,0,2,3,0,4))
	p1.dibujarCampoFiltrado()

test_filtroPasaAlta2()
plt.show()

(Source code)

16.5. Filtro pasa banda

Los filtros pasa banda se ubican entre los filtros pasa alta y los filtros pasa baja, permitiendo únicamente el paso de frecuencias dentro de una banda espectral, filtrando las altas y las bajas frecuencias. Para ello, se bloquean en el plano de Fourier los bajos y altos órdenes de difracción, así el filtro de intensidad tendrá forma de corona circular, siendo los puntos interiores a la corona, aquellos que permiten en paso de luz.

#!/usr/bin/env python
# -*- coding: utf-8 -*-
#------------------------------------
# Autor:	Luis Miguel Sanchez Brea
# Fecha		2013/10/24 (version 1.0)
# Licencia:	GPL
#-------------------------------------
"""Simulaciones relacionadas con la leccion de filtrado optico"""

#Carga de librerias y paquetes
import sys
sys.path.append('../')
from procesadoOptico import *
from fuentesXY import fuenteXY
from mascarasXY import mascaraXY

tamano = 1000 * um
x0 = sp.linspace(-tamano/2 * um, tamano/2 * um, 512)
y0 = sp.linspace(-tamano/2 * um, tamano/2 * um, 512)

def test_filtroPasaBanda(): #amplitud
	#Definicion de la fuente
	fuente1 = fuenteXY(x=x0, y=y0)
	fuente1.onda_plana(A=1, theta=0*grados, phi=0*grados)
	#Definicion del objeto
	objeto = mascaraXY(x=x0, y=y0)
	objeto.imagen(nombre='lena.png', normalizar=True, canal=0, tamanoImagen=False, invertir=False, angulo=0)
	#Definicion del filtro
	filtro = mascaraXY(x=x0, y=y0)
	filtro.anillo(r0=(0*um,0*um), radius1=(50*um,50*um), radius2=(80*um,80*um), angulo=0*grados)
	#Campo filtrado, procesado de la imagen
	p1=procesadoOptico(fuente1, objeto, filtro)
	p1.dibujar2(idfigs=(1,0,2,3,0,4))
	p1.dibujarCampoFiltrado()

test_filtroPasaBanda()
plt.show()

(Source code)

16.6. Filtro de contraste de fase

Un filtro de contraste de fase potencia las diferencias de fase entre las distintas regiones del objeto. Para ello se ilumina el objeto con luz colimada y se introduce una lámina/anillo desfasador posterior al sistema de iluminación, el objetivo de este anillo es retardar la fase de la luz.

#!/usr/bin/env python
# -*- coding: utf-8 -*-
#------------------------------------
# Autor:	Luis Miguel Sanchez Brea
# Fecha		2013/12/10 (version 1.5) - dibujar 2
# Licencia:	GPL
#-------------------------------------
"""Simulaciones relacionadas con la leccion de filtrado optico"""

#Carga de librerias y paquetes
import sys
sys.path.append('../')
from procesadoOptico import *
from fuentesXY import fuenteXY
from mascarasXY import mascaraXY

def filtroContrasteFase():
	#Definicion de la fuente
	fuente1 = fuenteXY(x=x0, y=y0)
	fuente1.onda_plana(A=1, theta=0*grados, phi=0*grados)
	#Definicion del objeto
	objeto = mascaraXY(x=x0, y=y0)
	objeto.imagen(nombre='lena.png', normalizar=True, canal=0, tamanoImagen=False, invertir=False, angulo=0)
	objeto.setFase(q=1, fase_min=0, fase_max=sp.pi)
	#Definicion del filtro
	filtro = mascaraXY(x=x0, y=y0)
	filtro.circulo(r0=(0*um,0*um), radius=(4*um,4*um), angulo=0*grados)
	filtro.setFase(q=1, fase_min=0, fase_max=sp.pi)
	#Campo filtrado y procesado de la imagen
	p1=procesadoOptico(fuente1, objeto, filtro)
	p1.dibujar2(idfigs=(0,1,2,0,3,4))
	p1.dibujarCampoFiltrado()
filtroContrasteFase()
plt.show()

(Source code)