PONTIFICIA UNIVERSIDAD CATÓLICA DEL PERÚ FACULTAD DE CIENCIAS E INGENIERÍA REGISTRO DE IMÁGENES MULTIESPECTRALES AÉREAS UTILIZANDO LA TRANSFORMADA DE WAVELETS Tesis para optar el Título de Ingeniero Electrónico, que presenta el bachiller: Susan Lucy Palacios Salcedo ASESOR: Donato Andrés Flores Espinoza Lima, Julio del 2012 RESUMEN La presente tesis tiene como objetivo la obtención de imágenes aéreas multiespectrales de campos de cultivos a partir de dos imágenes aéreas de una sola banda espectral, este procedimiento es llamado registro de imágenes. En el primer capítulo, se presenta el registro de imágenes como una tarea crucial en el procesamiento de imágenes digitales y particularmente en aplicaciones de percepción remota. La percepción remota es empleada para el monitoreo de cultivos en el proyecto que viene realizando la Pontificia Universidad Católica del Perú y el Centro Internacional de la Papa denominado “Agricultura de Precisión para la Producción de Cultivos de Seguridad Alimentaria y de Agro Exportación”. En el segundo capítulo, se presenta diversos métodos empleados para el registro de imágenes y se propone un método basado en la transformada discreta de wavelets. El algoritmo propuesto es el Dual Tree- Complex Wavelet Transform (DT-CWT), utilizado en el registro de imágenes de diferentes espectros; debido a que presenta ventajas como la multiresolución y su sensibilidad de direccionamiento para resaltar características significativas en las imágenes. En el tercer capítulo, se explica el desarrollo del algoritmo en la herramienta matemática MATLAB. En la primera etapa se obtienen los coeficientes de wavelets, estos deben ayudar a resaltar detalles en las imágenes de la mejor manera posible. En la segunda etapa se muestra la implementación del algoritmo (DT-CWT) con bancos de filtros y finalmente se utiliza la correlación cruzada para poder realizar el registro. Finalmente, se muestra la evaluación del algoritmo en las primeras pruebas se empleo la imagen de Lena y posteriormente se utilizó imágenes de campos de cultivo de olivo y camote. Asimismo, se presentan las conclusiones y recomendaciones del presente trabajo. A Dios, por guiar mi camino. A mis padres y hermano, por enseñarme que todo es posible con perseverancia y creer en mí. Al profesor Andrés Flores, por su confianza y ayuda en todo momento. A mis amigos y familiares, por alentarme a seguir adelante y estar a mí lado. Gracias a todos. ÍNDICE INTRODUCCÍÓN.........................................................................................................1 CAPÍTULO 1: REGISTRO DE IMÁGENES MULTIESPECTRALES..........................2 1.1 Alcances de la Agricultura de Precisión.................................................................2 1.2 Alcances de la Percepción Remota ......................................................................2 1.3 Alcances de Registro de Imágenes ......................................................................3 1.4 Declaración del Problema .....................................................................................3 CAPÍTULO 2: LA TRANSFORMADA DE WAVELET: UN EFICIENTE ENFOQUE PARA EL REGISTRO DE IMÁGENES MULTIESPECTRALES ................................6 2.1 Estado del Arte ......................................................................................................6 2.1.1 Presentación del Asunto de Estudio...............................................................6 2.1.2 Estado de la Investigación..............................................................................7 2.1.3 Síntesis sobre el Asunto de Estudio...............................................................8 2.2 Modelo Teórico ......................................................................................................9 2.2.2 Alcances de la Transformada Discreta de Wavelet (DWT) .......................9 2.2.2.1 La Transformada Discreta de wavelet (DWT)...............................11 2.2.2.2 Algoritmo Dual Tree Complex Wavelet Transform (DT-CWT)…...14 2.2.3 Registro de Imágenes ...............................................................................17 2.3 Metodología de la Investigación.............................................................................18 CAPÍTULO 3: DISEÑO DEL ALGORITMO PARA EL REGISTRO AUTOMÁTICO DE IMÁGENES ...........................................................................................................18 3.1 Objetivos ...............................................................................................................19 3.2 Requerimientos de Diseño ....................................................................................19 3.3 Diseño del Algoritmo para el Registro de Imágenes basado en Wavelets.............20 3.3.1 Diseño del Algoritmo para la obtención de Wavelets ................................20 3.3.1.1 Diseño de los Coeficientes del Banco de Filtros..............................20 3.3.1.2 Diseño del Algoritmo de Banco de Filtro..........................................25 3.3.2 Diseño del Algoritmo para el Registro de Imágenes .................................27 CAPÍTULO 4 :EVALUACIÓN DEL ALGORITMO BASADO EN LA TRANSFORMADA DE WAVELETS ..........................................................................29 4.1 Evaluación de Wavelets ........................................................................................29 4.2 Evaluación del Algoritmo de Registro.....................................................................35 4.3 Resultados de Registro de Imágenes Aéreas de Cultivo ......................................37 4.4 Evaluación del Registro de imágenes con la Transformada de Wavelet Discreta.........................................................................................................................49 CONCLUSIONES.........................................................................................................55 RECOMENDACIONES ................................................................................................56 BIBLIOGRAFÍA ............................................................................................................57 ANEXOS 1 INTRODUCCIÓN La mejora de la calidad en los procesos de la producción agrícola a nivel mundial se debe a la aplicación de tecnología en la agricultura, esto se conoce como agricultura de precisión. Este concepto se basa en la obtención de datos de sensores remotos, sistema de información geográfica, GPS, entre otros; para ayudar a la detección de anomalías en los cultivos y posteriormente buscar medidas correctivas, los cuales contribuirán a la obtención de una producción de calidad. Uno de los centros que viene desarrollando dicho concepto en nuestro país es el Centro Internacional de la Papa (CIP) [1] [2], en donde para obtener una adecuada información de diferentes cultivos se utiliza la percepción remota. Esto se basa en monitorear los cultivos por medio de imágenes multiespectrales, una imagen multiespectral es obtenida por cámaras especiales que detectan determinados espectros ó por la alineación de imágenes provenientes de diferentes sensores espectrales. En el CIP la percepción remota se realiza con dos sensores diferentes, uno trabaja con el espectro infrarrojo y el otro con el espectro rojo; ambos captan imágenes a una determinada altura sobre la superficie utilizando como medio a equipos de aeromodelo. El proceso de alinear espacialmente imágenes de una misma escena utilizado en la percepción remota se denomina registro de imágenes, esta técnica cuyas aplicaciones no sólo tienen importancia en la agricultura de precisión; sino también en otras áreas como la medicina, visión computarizada, entre otras. El registro de imágenes presenta diferentes métodos para sus diferentes aplicaciones, siendo unos más efectivos que otros. El objetivo del presente documento consiste en establecer los diferentes métodos de registro de imágenes y presentar un enfoque basado en la transformada compleja de wavelet, como método eficaz para el registro de imágenes aéreas de diferente espectro. 2 CAPÍTULO 1 REGISTRO DE IMÁGENES MULTIESPECTRALES En este primer capítulo, se mencionan los alcances del registro de imágenes aplicado en la agricultura de precisión. 1.1 Alcances de la Agricultura de Precisión: A partir de los años 70 se desarrolla una nueva filosofía en la agricultura a nivel mundial, se inician los estudios de automatización de máquinas agrícolas, en las décadas siguientes se empieza a realizar el manejo localizado de las prácticas agrícolas con el desarrollo de equipos inteligentes como el GPS. Logrando una aplicación más eficiente de los insumos, reduciendo el impacto ambiental, y por ende reduciendo los costos de la producción de alimentos. La agricultura de precisión utiliza diferentes técnicas orientadas a mejorar la productividad agrícola. Para lograr su objetivo estudia la variabilidad espacial y temporal de las características del suelo y del cultivo; este análisis permite identificar, cuantificar, mapear e incluso georeferenciar la variabilidad; de modo que se determine la cantidad de insumos de acuerdo a la necesidad de los puntos del área de cultivo. Este nuevo concepto de producción agrícola es posible debido a la evolución de tecnologías como el sistema de posicionamiento global (GPS), sistema de información geográfica (SIG), percepción remota, tecnología de dosis variable (sensores, controladores) y el análisis de datos georeferenciados (geoestadística). La aplicación de tecnologías a la agricultura de precisión proporciona la capacidad de recolección de datos, procesamiento e interpretación de la información y la aplicación de los insumos. [3] [4] 1.2 Alcances de la Percepción Remota: La percepción remota es aplicada en la observación de la variabilidad de los cultivos, su idea básica es capturar imágenes aéreas o satelitales mediante sensores con el propósito de recolectar datos e información de áreas de cultivo de manera remota. Los sensores utilizados son dispositivos radiométricos que detectan y registran la radiación electromagnética en una banda determinada del espectro electromagnético, además de generar información para su posterior interpretación. 3 Las imágenes son el resultado de la radiación reflejada de los elementos del área de cultivo por ejemplo un sensor de espectro infrarrojo cercano muestra una imagen caracterizada por la radiación muy reflejada de la vegetación, por otro lado una imagen del sensor de espectro rojo se caracteriza por la radiación muy reflejada de la tierra. La información adquirida es transformada en índices de vegetación como el NDVI (Índice de Vegetación Diferencial Normalizado), a fin de identificar la variabilidad en los cultivos. El NDVI es obtenido a partir de una imagen multiespectral de donde se extraen los parámetros de reflectancia de las bandas espectrales rojo e infrarrojo cercano. [3] 1.3 Alcances de Registro de Imágenes: El registro de imágenes constituye una tarea fundamental en la percepción remota, debido a que las imágenes son adquiridas en diferentes condiciones y estas necesitan ser alineadas. Este proceso consiste en alinear geométricamente dos ó más imágenes adquiridas de diferentes puntos de vista (distintas traslaciones, rotaciones ó escala), en diferentes tiempos (distintas condiciones de iluminación) ó por diferentes tipos de sensores de una misma escena. La percepción remota utilizado en el monitoreo de cultivos, presenta el problema de alinear imágenes obtenidas, usualmente, por diferentes tipos de sensores como el rojo e infrarrojo, pues se requiere construir una imagen multiespectral de donde se pueda procesar la información como el NDVI. 1.4 Declaración del Problema En nuestro país el uso de tecnología en la agricultura aún es limitada, solo se vienen ejecutando proyectos en los cuales se busca demostrar los beneficios de su aplicación. Un proyecto importante que se está ejecutando es “Agricultura de Precisión para la Producción de Cultivos de Seguridad Alimentaria y de Agro Exportación” financiado por el Programa de Ciencia y Tecnología (FINCYT) de la Presidencia del Consejo de Ministros (PCM), aportados por el Banco Interamericano de Desarrollo (BID) y el tesoro público; ejecutado por la PUCP y el Centro Internacional de la Papa (CIP) [1] [2]. Este proyecto tiene como objetivo desarrollar tecnologías de monitoreo de áreas de cultivo mediante el análisis del NDVI utilizando imágenes multiespectrales. La obtención de imágenes multiespectrales es partir de imágenes de un solo espectro, 4 estas últimas son capturadas por dos sensores (Ver figura 1.1) uno de espectro infrarrojo cercano (850nm) y el otro de espectro rojo (650nm), ambas cámaras son previamente calibradas y son ubicadas en aeromodelos como aviones, helicópteros ó globos. La figura 1.2 nos muestra la ubicación de ambas cámaras, como se puede observar estas están separadas 15 mm. Las imágenes obtenidas de ambas cámaras están trasladadas una respecto a otra; esto se debe a la distribución de ambos sensores, y a la captura de las imágenes por las cámaras cuando viajan en los aeromodelos a una altura de la superficie. Entonces la obtención de imágenes multiespectrales constituye un problema de registro de imágenes de diferente espectro y además que estén trasladadas una respecto a la otra. Figura 1.1 Sensores utilizados. 5 En el mercado actual existen equipos que realizan la captura de imágenes en diferentes bandas espectrales, como por ejemplo las cámaras multiespectrales del fabricante Tetracam de serie ADC [5] que ofrece ventajas para el manejo de imágenes de cultivo. Sin embargo; se tiene la desventaja de un elevado costo, en contraste la utilización de dos sensores de diferentes bandas espectrales reducen los costos; pero se presenta el problema de buscar un método adecuado para el registro de las imágenes aéreas de diferentes espectros. 15 mm Sensor de Banda espectral rojo. Sensor de Banda espectral infrarrojo cercano. Figura 1.2 Distancia de separación de ambos sensores. 6 CAPÍTULO 2 LA TRANSFORMADA DE WAVELET: UN EFICIENTE ENFOQUE PARA EL REGISTRO DE IMÁGENES MULTIESPECTRALES En este capítulo se explica algunos métodos utilizados en el registro de imágenes y se propone la utilización de un algoritmo basado en la Transformada Compleja de Wavelets como una técnica eficiente para lograr una alineación correcta. 2.1 Estado del Arte: 2.1.1 Presentación del Asunto de Estudio: Diversos métodos son propuestos para el registro de imágenes en la percepción remota, estos se pueden dividir en métodos basados en la intensidad, en la frecuencia y los basados en características. Los métodos basados en la intensidad se encuentran con un difícil escenario debido a que las imágenes de percepción remota provienen de diferentes bandas espectrales o diferentes sensores, capturando así imágenes con diferentes valores de intensidad en los pixeles de la misma escena, lo que dificulta la comparación de las imágenes. Los métodos basados en las características, tienen un enfoque de extracción de características (esquinas, caminos, construcciones, etc) para realizar la comparación y alinear las imágenes, pero este método se enfrenta a algunos inconvenientes [6] [7]: 1. Sensibilidad al ruido: El ruido del sensor puede afectar significativamente la detección de puntos de control, lo que causa imágenes con altos niveles de ruido y con detalles irrelevantes. 2. Sensibilidad al contraste: Las imágenes de percepción remota están sometidos a diferentes condiciones de contraste tanto globales y locales, causando que las características de las imágenes contengan diferentes valores de intensidad. 3. Inexactitud de localización de los puntos de control: En algunos métodos se asume que la ubicación de los puntos de control son los mismos en las respectivas imágenes pero a veces no es cierto, porque durante el proceso de detección de puntos se puede incurrir en inexactitudes. 7 2.1.2 Estado de la Investigación:  Método basado en la intensidad Un método eficiente basado en la Información mutua es descrito en [8], en donde se presenta un enfoque basado en la evaluación de la entropía y la información mutua, al cual se denomina ENMA (Empirical Entropy Manipulation and Analysis). La alineación propuesta por ENMA se basa en la intensidad pero adicionalmente considera el efecto del cambio de intensidades entre las imágenes que se están registrando. Esta técnica nos dice que un modelo y la imagen están alineados cuando existe una relación coherente entre la intensidad de la imagen y el modelo normal.  Método basado en características En [9] se presenta un enfoque basado en la extracción de características geométricas. El método de Straight Line features Matching on Aerial Imagery propone una estrategia de búsqueda de concordancia de jerarquía de líneas rectas, las características lineales pueden proporcionar más información que las características de los puntos; ya que no son tan propensos a ser afectados por el ruido y es más conveniente ser extraído y ser descrito en comparación con características basadas en regiones. Un método basado en la extracción de regiones como característica es presentado en [10], este enfoque utiliza regiones con límites cerrados para establecer correspondencia entre las imágenes. Para obtener dichas regiones se utiliza la segmentación de imágenes en este paso se pueden identificar varias técnicas. Después de identificar las correspondientes regiones en las imágenes estas son refinadas hasta lograr una similitud tanto como sea posible. Luego son obtenidas las coordenadas de los centros de gravedad de dichas regiones con una exactitud de subpixel, entonces estos centros de gravedad son usados como puntos de control. Con este método podemos registrar imágenes con rotación, traslación y escalas diferentes. Otro método de extracción de características es presentado en [11]. Este algoritmo denominado SIFT (Scale Invariant Feature Transform) tiene como propósito extraer características invariantes a la escala y rotación de la imagen; para ello utiliza la función DoG (Difference-of-Gaussian) para detectar los puntos de interés y un histograma de orientación. En [12] es presentado el algoritmo SURF (Speeded Up Robust Feature), este al igual que SIFT se basa en la extracción de características 8 geométricas, pero ofrece una mayor eficiencia en la velocidad de procesamiento ya que emplea el detector rápido de Hessian y además utiliza la información de la transformada de wavelet Haar.  Método basado en el análisis de Wavelet En [7], [13], [14] y [15] se proponen diferentes métodos para el registro de imágenes basadas en la descomposición de wavelet [16] y [17]. La elección de utilizar wavelet en el registro de imágenes de percepción remota se justifica por su característica de multiresolución. La multiresolución de wavelets permite descomponer las imágenes sin mucha pérdida de información, conservando importantes características de la información original; lo que hace posible la detección de texturas en las imágenes. Adicionalmente la descomposición de la imagen contribuye a la rapidez de ejecución de los algoritmos. El primer método propuesto [13] se basa en la descomposición discreta de wavelet y la correlación. El desarrollo del algoritmo se describe de acuerdo a las etapas de registro. Para ello utiliza “los máximos coeficientes de wavelet”, se llama así a los puntos cuyo valor de intensidades están por encima de un x% de los histogramas. El segundo método [7] introduce ruido invariante, contrasta la invariante detección de puntos de control y puntos de control similares obtenidos en base a características de wavelet complejo. El tercer método [14] se basa en la extracción de regiones para ello utiliza un algoritmo llamado Dual-Tree Complex Wavlelet Tranform (DT-CWT) [18], realiza la segmentación de características de las imágenes, construyendo mapas de regiones que permitirán realizar la alineación. El último método descrito en [15] al igual que el anterior utiliza el DT-CWT, pero adicionalmente utiliza detección de borde, la correlación cruzada así como también emplea la información mutua. 2.1.3 Síntesis sobre el Asunto de Estudio El registro de imágenes trasladadas, rotadas, de escalas y espectros diferentes tienen un componente importante que es encontrar características ó puntos de control en las imágenes los más adecuados posibles para realizar un registro de imágenes lo más exacta posible. Existen variados métodos para desarrollar este proceso, unos se basan en diferentes características de la imágenes como la intensidad de determinados pixeles; formas definidas como líneas rectas o intersecciones de líneas; también 9 regiones definidas dentro de la imagen. Estos métodos de diferentes enfoques, pueden ser unos mejores que otros dependiendo de su aplicación. Estudios realizados en los últimos años sugieren que los métodos basados en transformada de wavelets se presenta como una buena alternativa para el registro de imágenes de diferentes sensores ya que es un método robusto que ayuda a detectar características en las imágenes resaltando los detalles, especialmente la transformada compleja de wavelets. Por tal motivo, el presente documento propone utilizar un método basado en la trasformada de wavelets para el registro de las imágenes de diferentes espectros. 2.2 Modelo Teórico 2.2.2 Alcances de la Transformada Discreta de Wavelet (DWT): A continuación se definirá algunos términos definidos en “Wavelet Transform”, Elements of Wavelets for Engineers and Scientists [17]. Vector espacio: El vector espacio es un conjunto V = {vi} que se caracteriza por realizar las siguientes dos operaciones y además cumplir las propiedades de vectores geométricos. Operaciones:  Dos vectores pertenecientes a V pueden sumarse y formar un tercer vector.  Un escalar ai perteneciente a un campo escalar A se puede multiplicar con un vector para formar un tercer vector. Bases: Las bases de un conjunto V se forman a partir de un subconjunto {v1,v2,….vi} perteneciente a V, donde los vectores vi son linealmente independientes y además cualquier vector diferente de cero perteneciente a V hace al subconjunto dependiente. Base Ortogonal: Las bases ortogonales se forman a partir de un conjunto S={v1,v2…vn}, donde el producto interno de los vectores vi y vk es cero. (2.1) 10 Base Ortonormal: El conjunto ortogonal S forma una base ortonormal si el producto interno consigo misma es igual a 1: , para i=k (2.2) Complemento Ortogonal: Si V0 es un subconjunto del vector espacio V1, se define el complemento ortogonal W0, como el conjunto de todos los vectores de V1 que son ortogonales a cada elemento de V0. (2.3) Además, se puede decir que V1 es la suma directa ( ) de V0 y W0 y se puede escribir como: (2.4) Bases de escala y Wavelet: Las funciones {φjk} son llamadas funciones de escalas y { jk} son llamadas funciones de wavelets. con j, k  (2.5) con j, k  (2.6) En la ecuación (2.5) se muestra la función escala, donde se representa a la función de escala básica φjk(t) por dos factores, el primero 2 j/2 normaliza cada función para que tenga una norma unitaria y el segundo factor φ(2jt-k), donde se encuentra dos términos 2j dilata y k desplaza la función de escala básica. Por otro lado, en la ecuación (2.6) se muestra la función wavelet que al igual que la función escala, se representa por dos factores de dilatación y traslación de la función wavelet básica, a la cual se denomina wavelet madre. 11 Las funciones de escalas básica {φjk} que se utilizan satisfacen la condición de ortogonalidad, tal que las traslaciones discretas {φ(2jt-k)} forman un conjunto ortonormal. Así mismo las funciones de wavelet { jk} también satisfacen las condiciones mencionadas. Más aún se tiene que las funciones de escala forman una base Vj y las funciones wavelet forman otra base W j entonces Vj y Wj satisfacen la condición de complemento ortogonal. Se puede observar que la combinación de funciones de escala y wavelet forman dos bases ortonormales para un siguiente nivel espacial (por ejemplo, V0 + W0 forma la base de V1, también V1+W1 forma la base de W2, etc.). Sin embargo, se puede emplear otras combinaciones como las funciones de escala consigo misma y las funciones de wavelet consigo mismas. 2.2.2.1 Transformada Discreta de Wavelet (DWT): La DWT en el campo de procesamiento de señales tiene varias aplicaciones como la comprensión de video, reconocimiento de objetos, comprensión de comunicaciones vía internet, registro de imágenes y análisis numéricos. [16] La DWT permite descomponer una señal x(t) en términos de wavelets y escalas : Los coeficientes de escala cjk y los coeficientes de wavelets djk son calculados por el producto interno. (2.8) (2.9) Existe un algoritmo eficiente para calcular los coeficientes de DWT a partir de muestras de una señal digital. Este algoritmo se propuso en 1989 por Mallat, el cual se basa en el análisis de multiresolución y para su implementación se utiliza banco de filtros. [19] [20] [21]. La idea que permite relacionar la transformada de wavelets y el banco de filtros se basa en aplicar recursivamente los coeficientes del filtro para obtener la función escala, así se consigue de forma recursiva varias resoluciones como resultado 12 de convolucionar el vector de entrada con el filtro y realizar un cambio de escala. Entonces, toda la información de las funciones de escala y wavelets es obtenida a partir de los coeficientes de los filtros. [16][17] Función de escala: Función de wavelet: Donde h0 y h1 son filtros pasa bajo y pasa alto respectivamente. En la figura 2.1, se puede apreciar un banco de filtros, estos son llamados filtros de análisis, los cuales permite encontrar la transformada discreta de wavelets y en la figura 2.2 se tiene los filtros de síntesis que representa la transformada discreta inversa de wavelet. En la figura 2.1 se observa la descomposición de la señal x(n) en sus coeficientes c y d, además se realiza la operación de sub muestreo para realizar el cambio de escala. Mientras que en la figura 2.2 se observa la reconstrucción de la señal a partir de los coeficientes c y d, además se aplica la operación de sobre muestreo. Los filtros de análisis (h0 y h1) y los filtros de síntesis (H0 y H1) son diseñados de modo que cumplan las condiciones de reconstrucción perfecta (PR) para que y(n)=x(n) ó en el caso más general y(n)=x(n-n0). Estas condiciones pueden ser QMF (Quadrature Mirror Filter) ó CQF (Conjugate Quadrature Filter) [16][17][19]. Adicionalmente, a veces en el procesamiento de señales se requiere que la wavelet este bien localizada en el tiempo, además que los filtros a implementar deben tener respuesta finita (FIR). Figura 2.1 Filtros de Análisis. h0(n) 2 h0(n) 2 h1(n) 2 h1(n) 2 x(n) cj+2(n) dj+2(n) cj+1(n) dj+1(n) 13 Figura 2.2 Filtros de Síntesis. La implementación del algoritmo para imágenes se realiza de la siguiente manera, primero se aplica el algoritmo de una dimensión a las filas y luego a las columnas (ver figura 2.3). A continuación, se describe la obtención de la función escala y tres wavelets:  Primero, la sub banda LL es producto de aplicar la función de filtro pasa bajo φ(x) a la primera dimensión y luego se aplica la función de filtro pasa bajo φ(y) a la segunda dimensión. φLL(x,y)= φ(x) φ(y) (2.12)  El wavelet de sub banda LH es producto de aplicar la función de filtro pasa bajo φ(x) a la primera dimensión y luego se aplica la función de filtro pasa alto (y) a la segunda dimensión. LH(x,y)= φ(x) (y) (2.13)  El wavelet de sub banda HL se obtiene aplicando la función de filtro pasa alto (x) a la primera dimensión y luego se aplica la función de filtro pasa bajo φ(y) a la segunda dimensión HL(x,y)= (x) φ(y) (2.14) dj+1(n) 2 H0(n) 2 H0(n) 2 H1(n) 2 H1(n) cj+2(n) dj+2(n) cj+1(n) y(n) 14  El último wavelet de sub banda HH se obtiene aplicando la función de filtro pasa alto (x) a la primera dimensión y luego se aplica la función de filtro pasa alto (y) a la segunda dimensión. HH(x,y)= (x) (y) (2.15) La característica principal de obtención de los wavelets en diferentes sub bandas es que las sub bandas LH y HL permite resaltar características orientadas verticalmente y horizontalmente; mientras que el wavelet de sub banda HH tiene una orientación de composición +45 y -45. Figura 2.3 Filtros de Análisis para dos dimensiones. [22] 2.2.2.2 Algoritmo Dual-Tree Complex Wavelet Transform (DT-CWT): El algoritmo de DT-CWT descrito en [18] ofrece mejoras significativas como la obtención de selectivas direccionalidades en dos o más dimensiones, sobre el algoritmo de DWT a cambio de la redundancia. Para su implementación también se utiliza bancos de filtro, sin embargo este algoritmo opta por la utilización de wavelets complejos y no reales. La razón de no implementar la DWT con wavelets reales se debe a cuatro problemas principales: el primero es las oscilaciones, esto debido a que desde que los wavelets son funciones pasa bandas tienden a oscilar positivamente y h0(n) h0(n) h1(n) h1(n) Imagen inicial correspondiente al nivel m-1 2 2 h0(n) h1(n) I(x,y) Filas Columnas 2 2 2 2 LL LH HL HH Imagen correspondiente al nivel m Detalles de las imágenes correspondientes al nivel m-1 15 negativamente alrededor de un punto de singularidad complicando la extracción de estos puntos; el segundo es el cambio en la varianza, es decir un pequeño cambio en la señal perturba significativamente la oscilación del coeficiente de wavelet patrón; el tercer problema es el aliasing, esto se produce debido a que los coeficientes de wavelets son calculados por iteración discreta en el tiempo y operaciones de sub muestreo intercaladas utilizando filtros pasa bajos y pasa altos no ideales y el último problema es la falta de selectividad direccional la cual complica el procesamiento de las características geométricas de la imagen como bordes y crestas, estos son descritos claramente en [18]. Este algoritmo inspirado en la representación de Fourier, donde los componentes de seno (parte imaginaria) y coseno (parte real) forman un par de transformada de Hilbert, busca construir una función de wavelets compleja a partir de wavelets reales. (2.16) Donde la parte imaginaria debe ser aproximadamente la transformada de Hilbert de la parte real { .[12] (2.17) Para su implementación se emplea dos bancos de filtros de DWT reales, donde el primer DWT proporciona la transformada de la parte real, mientras que la segunda proporciona la parte imaginaria. El algoritmo busca que los coeficientes de los filtros utilizados no sean imaginarios por ello utiliza la condición de Half Sample Delay, para mayor detalle revisar [18]. La implementación de este algoritmo en imágenes como se mencionó anteriormente permite obtener selectividad en el direccionamiento, a diferencia del DWT donde se tenían 3 direccionamientos horizontal, vertical y diagonal; con este algoritmo podemos resaltar características orientada en seis direcciones ±15, ±45 y ±75. En la figura 2.4 se puede observar los filtros de análisis donde h0 y h1 son el filtro pasa bajo y pasa alto respectivamente, estos proporcionan la transformada de la parte real; por otro lado g0 y g1 son el filtro pasa bajo y pasa alto que proporciona la parte imaginaria de la transformada. 16 Figura 2.4 Filtros de análisis de DT-CWT para 2 D. A continuación se presentará la obtención de los 6 direccionamientos, la demostración se encuentra en [18]: (2.18) ) (2.19) (2.20) (2.21) (2.22) (2.23) El motivo de utilizar este algoritmo es porque al utilizar wavelets complejos soluciona los cuatro problemas mencionados anteriormente. Dos de estos problemas principalmente originan inconvenientes en el análisis de imágenes ya que el problema Filas Columnas h1(n) g0(n) g0(n) g1(n) g1(n) Imagen inicial correspondiente al nivel m-1 2 2 g0(n) g1(n) I(x,y) 2 2 2 2 LLi LHi HLi HHi Imagen correspondiente al nivel m Detalles de las imágenes de la parte imaginaria correspondientes al nivel m-1 h0(n) h0(n) h1(n) 2 2 h0(n) h1(n) 2 2 2 2 LLr LHr HLr HHr Imagen correspondiente al nivel m Detalles de las imágenes de la parte real correspondientes al nivel m-1 17 de oscilaciones dificulta la extracción de singularidades; y la pobre direccionalidad, complica el procesamiento de imágenes con características geométricas. 2.2.3 Registro de Imágenes El registro de imágenes es el proceso de determinar la correspondencia entre puntos de imágenes de una misma escena. El proceso de registro de imágenes consiste en tres pasos: 1. Detección características, en esta etapa se pueden realizar cualquiera de los dos procedimientos que se describen a continuación. El primero es seleccionar puntos de control como líneas, bordes, esquina, etc. de una imagen referencia y de otra imagen fuente; y el segundo busca resaltar los detalles, las formas geométricas, en ambas imágenes. 2. Estimación de parámetros de transformación, en esta etapa se busca los parámetros de transformación adecuados para el registro de la imagen fuente y de referencia. 3. Re muestreo de imágenes, en esta última parte la imagen fuente se transforma utilizando los parámetros hallados en la etapa anterior. El método de registro de imagen se puede explicar de la siguiente forma [15]: Se tiene una imagen fuente If(x,y) y una imagen de referencia Ir(x,y). Ambas imágenes tienen una región en común , entonces se puede decir que existe una transformación geométrica T(x,y) la cual permite que If(x,y)≈ Ir(T(x,y)). La transformación más común entre imágenes es la que presenta diferente rotación ( ), traslación ( ) y escala (s). Considerando las características mencionadas se define a la transformación geométrica como: (2.24) Entonces el problema del registro es encontrar un vector de transformación. (2.25) 18 2.3 Metodología de la Investigación. El procedimiento a seguir para el diseño del algoritmo de registro de imágenes multiespectrales se divide de la siguiente manera: Etapa 1: Investigación preliminar:  Estudio de los diferentes métodos de alineación de imágenes de diferente sensores.  Estudio de la transformada discreta de wavelet.  Estudio de la técnica Dual-Tree Complex Wavelet Transform. Etapa 2: Análisis y diseño del algoritmo  Análisis del algoritmo DT- CWT y determinación de sus requerimientos. Etapa 3: Desarrollo de la solución  Diseño de los coeficientes del banco de filtros de DT-CWT en MATLAB  Diseño del banco de filtros de DT-CWT en MATLAB  Diseño del algoritmo para el registro de las imágenes en MATLAB Etapa 4: Resultados y documentación  Evaluación de los resultados obtenidos con el DT-CWT para determinar la utilidad en el registro de imágenes de cultivo.  Elaboración del documento final. 19 CAPÍTULO 3 DISEÑO DEL ALGORITMO PARA EL REGISTRO AUTOMÁTICO DE IMÁGENES En este capítulo, se indican los objetivos y se explica el diseño del algoritmo. 3.1 Objetivos: • Estudio de la transformada de wavelets para la aplicación del registro de imágenes multiespectrales. • Determinar los puntos de control óptimos para lograr una exitosa alineación de imágenes. • Desarrollar e implementar un banco de filtros utilizando la herramienta de software MATLAB. • Evaluar la pertinencia de los puntos de control para poder determinar su utilidad en el registro de imágenes multiespectrales de cultivo. 3.2 Requerimiento de diseño: Para obtener una imagen multiespectral de las áreas de cultivo, a partir de imágenes aéreas que presentan traslaciones y además son de diferentes espectros (infrarrojo cercano y rojo), se debe emplear el procesamiento de imágenes que permita resaltar características significativas en ambas imágenes, facilitando de esta forma el registro de imágenes. Figura 3.1 IR(x,y) DT-CWT DT-CWT Estimación de parámetros de transformación Alineación de Imágenes Primera etapa Segunda etapa Tercera etapa IFR(x,y) 20 Donde: IFR(x,y): Imagen de espectro infrarrojo IR(x,y): Imagen de espectro rojo En la figura 3.1 se puede observar el diagrama de bloques propuesto para el registro de imágenes, el cual se divide en 3 etapas. En la primera etapa aplicaremos a cada imagen el algoritmo de DT-CWT para poder resaltar los detalles que presenta cada imagen. En la segunda etapa se buscará el vector de traslación para ello se utilizará la correlación cruzada y en la última etapa se obtendrá la imagen multespectral. 3.3 Diseño del Algoritmo para el Registro Imágenes basado en Wavelets: En esta parte del capítulo 3, se describe el desarrollo del algoritmo para el registro de imágenes, para ello se dividirá en dos etapas; el desarrollo del algoritmo basado en wavelets para resaltar las características en ambas imágenes y el desarrollo del algoritmo para el registro. 3.3.1 Diseño del Algoritmo para la obtención de Wavelets: En esta parte se describe el desarrollo del algoritmo para la obtención de wavelets utilizando el algoritmo de DT-CWT; para ello se dividirá en dos etapas: la búsqueda de los coeficientes de los filtros y el diseño del algoritmo de banco de filtros. 3.3.1.1 Diseño de los Coeficientes del banco de Filtros: En [18] se propone tres métodos para la obtención de los coeficientes de los filtros de análisis y síntesis. Para el desarrollo del documento se utilizará el método de Solución del factor Común. Este método es desarrollado en [23], donde se describe el procedimiento de encontrar pares de bases de wavelets que forman aproximadamente un par de transformada de Hilbert. Consideraciones: La solución biortogonal presentada en [23] sugiere que para mejorar la sensibilidad en el direccionamiento se debe emplear wavelets de fase linear, que sigue el enfoque de obtención de wavelets simétricos biortogonales. Los filtros diseñados serán del tipo FIR 21 donde los coeficientes de los filtros de análisis (ver figura 3.2) y de síntesis (ver figura 3.3) deben formar pares de transformada de Hilbert. Figura 3.2 Filtros de Análisis. Figura 3.3 Filtros de Síntesis. Diseño de pares de filtros que cumplan con la transformada de Hilbert: Aquí se utiliza la condición descrita en [18] de Half Sample Delay. Para también se debe cumplir la condición. El dominio de Fourier se tiene: (3.5) (3.6) (3.7) En el dominio Z se tiene: (3.8) Para cumplir la condición (3.8) se utiliza un sistema de filtro pasa todo, tal que satisfaga (3.9) h0(n) 2 g0(n) 2 g1(n) 2 h1(n) 2 (n) 2 (n) 2 (n) 2 (n) 2 22 En la última ecuación se tiene a L como orden del filtro y (3.11) Para encontrar los coeficientes d(n) se utiliza la fórmula de Thiran: (3.12) Donde debe de ser 0.5, para que cumpla la condición de (3.8) Diseñando lo coeficientes de los filtros con el enfoque de simetría de biortogonalidad: Condiciones de biortogonalidad [19][20]: (3.13) Ahora se define la convolución de filtros como en [23] (3.14) (3.15) Se escribe (3.14) y (3.15) de modo que se cumpla con (3.13) , donde n0 es un entero impar. Entonces ph y pg deben cumplir con la propiedad de Halfband. Entonces, los coeficientes de los filtros pasa altos se encuentran de la siguiente forma [23]: (3.16) (3.17) (3.18) (3.19) Aquí se define los coeficientes de los filtros pasa bajo de análisis y síntesis utilizando el sistema de filtro pasa todo. [23] 23 (3.20) (3.21) (3.22) (3.23) En el dominio Z resolviendo (3.20) y (3.21) se obtiene la relación buscada en (3.9), también (3.22) y (3.23) cumplen con esta condición. A continuación las ecuaciones (3.14) y (3.15) son pasadas al dominio Z, obteniéndose [23]: (3.24) Para garantizar la propiedad de momentos de desvanecimiento de los wavelets se define y de la siguiente manera con + debe ser impar [23]: (3.25) (3.26) Remplazando (3.25) y (3.26) en (3.24) se obtiene: (3.27) Donde (3.28) y (3.29) Como ph y pg deben cumplir la propiedad de Halfband [23] (3.30) Despejando r(n) de la ecuación (3.30) se obtiene una secuencia simétrica, luego por factorización espectral se obtiene y ambos son simétricos debido a ello cumplen con la relación: (3.31) (3.32) 24 El desarrollo del algoritmo de encontrar los coeficientes se implementó en MATLAB, para su elaboración se tomó como base el programa sugerido por el autor Selesnick [23]. A continuación se muestra el Pseudo-Código del programa y en el Anexo 1 se muestra el código. Parámetros iniciales: k - Cantidad de ceros en z=-1, de los filtros de análisis K- Cantidad de ceros en z=-1, de los filtros de síntesis L – Orden del sistema de filtro pasa todo Tau- Retardo de las muestras Pseudocódigo 3.1 Coeficientes de Filtros Inicio d - Es el desarrollo de los coeficientes d(n) determinados con la fórmula de Thiran [23] s1 - Es el desarrollo del binomio , en el dominio del tiempo. s2 - Es el desarrollo del producto , en el dominio del tiempo. s - Es el resultado de la convolución de s1 y s2. r - Es obtenido mediante la ecuación que cumple la propiedad de Halfband. - Es una función simétrica resultado de factorizar espectralmente r. - Es una función simétrica resultado de factorizar espectralmente r. - Es el resultado de , en el dominio del tiempo. - Es el resultado de , en el dominio del tiempo. Fin Variables de salida: h0 - coeficientes del filtro pasa bajo de análisis de la parte real. h1 - coeficientes del filtro pasa alto de análisis de la parte real. H0 - coeficientes del filtro pasa bajo de síntesis de la parte real. H1 - coeficientes del filtro pasa alto de síntesis de la parte real. g0 - coeficientes del filtro pasa bajo de análisis de la parte imaginaria. g1 - coeficientes del filtro pasa alto de análisis de la parte imaginaria. G0 - coeficientes del filtro pasa bajo de síntesis de la parte imaginaria. G1 - coeficientes del filtro pasa alto de síntesis de la parte imaginaria. 25 3.3.1.2 Diseño del Algoritmo de Banco de Filtros: El diseño del banco de filtros se realizó hasta el nivel tres. A continuación se presenta el Pseudo-Código del algoritmo en MATLAB , y en el Anexo 2 se encuentra el código. Variables de entrada: I(x,y) – Imagen . n – orden del nivel de filtros. Coeficientes de los filtros de primer nivel. Coeficientes de los filtros para niveles superiores a 1. Pseudocódigo 3.2 Banco de Filtros Inicio M ← Cantidad de filas de I. N ← Cantidad de columnas de I. // En esta parte se obtiene las imágenes de sub bandas LLr y LHr Para i=1:M Im - Es el resultado de convolucionar cada fila de I con el filtro pasa bajo real de análisis Fin Im - Es el resultado de sub muestrear las filas de Im a razón de 2. Para i=1:N ImLLr - Es el resultado de convolucionar cada columna con el filtro pasa bajo real de análisis Fin ImLLr - Es el resultado de sub muestrear las columnas de ImLLr a razón de 2 Para i=1:N ImLHr - Es el resultado de convolucionar cada columna con el filtro pasa alto real de análisis Fin ImLHr - Es el resultado de sub muestrear las columnas de ImLHr a razón de 2. // En esta parte se obtiene las imágenes de subbandas HLr y HHr Para i=1:M Im - Es el resultado de convolucionar cada fila de I con el filtro pasa alto real de análisis Fin Im - Es el resultado de sub muestrear las filas de Im a razón de 2 Para i=N 26 ImHLr - Es el resultado de convolucionar cada columna con el filtro pasa bajo real de análisis Fin ImHLr - Es el resultado de sub muestrear las columnas de ImHLr a razón de 2 Para i=N ImHHr - Es el resultado de convolucionar cada columna con el filtro pasa alto real de análisis Fin ImHHr - Es el resultado de sub muestrear las columnas de ImHHr a razón de 2. // Las imágenes de las sub bandas de la parte imaginaria ImLLi, ImLHi, ImHLi y ImHHi se hallan de la misma manera pero con los filtros imaginarios. //Con todo lo anterior se obtiene las susbbandas del primer nivel. wv1- Es el resultado de operar (ImLHr-ImLHi) / wv2- Es el resultado de operar (ImHLr-ImHLi) / wv3- Es el resultado de operar (ImHHr-ImHHi) / wv4- Es el resultado de operar (ImLHr+ImLHi) / wv5- Es el resultado de operar (ImHLr+ImHLi) / wv6- Es el resultado de operar (ImHHr+ImHHi) / // Para hallar las sub bandas de los siguientes niveles se hallan de la misma forma pero My N cambian ahora serán el tamaño de la sub banda ImLLr para la parte real e ImLLi para la parte imaginaria. Los wavelets (wv1,wv2,wv3,wv4,wv5,wv6) de niveles superiores también se hallan de la misma forma descrita anteriormente Variables de Salida: wv1 – Imagen con características resaltada en la dirección de -75. wv2 - Imagen con características resaltada en la dirección de -15. wv3 - Imagen con características resaltada en la dirección de -45. wv4 - Imagen con características resaltada en la dirección de +75. wv5 - Imagen con características resaltada en la dirección de +15. wv6 - Imagen con características resaltada en la dirección de +45. 27 3.3.2 Diseño del Algoritmo para el Registro de Imágenes: En esta parte se describirá el método utilizado para encontrar el vector de transformación. Este vector solo tendrá parámetros de traslación (x,y); ya que las imágenes aéreas de los sensores presentan características de traslación. En [13] y [15] se sugiere utilizar la correlación cruzada entre las características de wavelets de la imagen roja e infrarroja. La correlación es empleada en el procesamiento de imágenes en la búsqueda de objetos en una imagen, alineación de imágenes, etc. [25]. El método propuesto utiliza la correlación entre las imágenes obtenidas por la descomposición de wavelets para obtener el vector de traslación.  Difn(x,y) , es la descomposición de la imagen de espectro infrarrojo.  Drn(x, y), es la descomposición de la imagen de espectro rojo.  n, representa el nivel de descomposición.  , representa las coordenadas (x,y) del valor máximo de la auto correlación de la imagen E, el cual siempre se encontrará en el centro geométrico de la imagen de correlación .  representa las coordenadas (x,y) del valor máximo de la correlación entre dos imágenes E y F. Luego: (3.33) y además Hallando el vector de transformación: (3.35) (3.34) 28 Como se puede apreciar, los parámetros del vector de transformación se determinan como la diferencia entre las coordenadas del valor máximo de la correlación cruzada entre las imágenes roja e infrarroja y las coordenadas del valor máximo de la auto correlación de una de las imágenes. Esta diferencia es multiplicada por un factor ya que las imágenes al descomponerse reducen sus dimensiones en . Esto debido a las operaciones de sub muestreo de filas y columnas. Las variables de entrada del algoritmo son las imágenes resultado de la descomposición de wavelets (wv1, wv2, wv3, wv4, wv5, wv6) y las variables de salida son los vectores de traslación para cada descomposición. Para ejecutar el algoritmo de registro de imágenes en primer lugar se deberá ejecutar el archivo vect1.m, el cual determinara los vectores de traslación y finalmente se debe ejecutar el archivo alin.m, este mostrará la imagen alineada y guardará la imagen. En el Anexo 3 se encuentra el código de vect1 y alin desarrollados en MATLAB. 29 CAPÍTULO 4 EVALUACIÓN DEL ALGORITMO BASADO EN LA TRANSFORMADA DE WAVELETS Este último capítulo comprende la evaluación y los resultados obtenidos del diseño realizado. 4.1 Evaluación de los Wavelets: En esta primera parte se evalúa el orden de los coeficientes de los filtros de análisis (k), síntesis (K) y además, el orden del filtro pasa todo (L).Para ello se procesará una imagen hasta el nivel 3 y se buscará el orden adecuado de los parámetros con los cuales se resaltan significativamente las características (bordes) de la imagen, para cada nivel se mostrará los seis direccionamientos (wv1, wv2, wv3, wv4, wv5 y wv6) obtenidos de la imagen. Figura 4.2 Imágenes de los seis direccionamientos del nivel1 para K=1, k=2 y L=1. Imagen wv1 Imagen wv2 Imagen wv3 Imagen wv4 Imagen wv5 Imagen wv6 Figura 4.1 Imagen de prueba - Lena 30 Figura 4.3 Imágenes de los seis direccionamientos del nivel 2 para K=1, k=2 y L=1. Figura 4.4 Imágenes de los seis direccionamientos del nivel 3 para K=1, k=2 y L=1 Como se puede apreciar en la Wavelets de nivel 3 ya no se aprecian claramente las características por ello en los siguientes imágenes solo se mostrará los resultados hasta el nivel 2. Imagen wv1 Imagen wv2 Imagen wv3 Imagen wv4 Imagen wv5 Imagen wv6 Imagen wv1 Imagen wv2 Imagen wv3 Imagen wv4 Imagen wv5 Imagen wv6 31 Figura 4.5 Imágenes de los seis direccionamientos del nivel 1 para K=3, k=2 y L=1 Figura 4.6 Imágenes de los seis direccionamientos del nivel 2 para K=3, k=2 y L=1 Imagen wv1 Imagen wv2 Imagen wv3 Imagen wv4 Imagen wv5 Imagen wv6 Imagen wv1 Imagen wv2 Imagen wv3 Imagen wv4 Imagen wv5 Imagen wv6 32 Figura 4.7 Imágenes de los seis direccionamientos del nivel 1 para K=3, k=2 y L=2 Figura 4.8 Imágenes de los seis direccionamientos del nivel 2 para K=3, k=2 y L=2 Imagen wv1 Imagen wv2 Imagen wv3 Imagen wv4 Imagen wv5 Imagen wv6 Imagen wv1 Imagen wv2 Imagen wv3 Imagen wv4 Imagen wv5 Imagen wv6 33 Figura 4.9 Imágenes de los seis direccionamientos del nivel 1 para K=3, k=2 y L=3 Figura 4.10 Imágenes de los seis direccionamientos del nivel 2 para K=3, k=2 y L=3 Imagen wv1 Imagen wv2 Imagen wv3 Imagen wv4 Imagen wv5 Imagen wv6 Imagen wv1 Imagen wv2 Imagen wv3 Imagen wv4 Imagen wv5 Imagen wv6 34 Figura 4.11 Imágenes de los seis direccionamientos del nivel 2 para K=3, k=2 y L=5 Figura 4.12 Imágenes de los seis direccionamientos del nivel 2 para K=3, k=2 y L=5 De las imágenes obtenidas para diferentes valores K, k y L se puede observar que se resaltan significativamente los bordes de la imagen de Lena para los coeficientes de K=1, k=2 y L=1; mientras que para otros valores se resaltan menos los bordes. La búsqueda de resaltar significativamente las características de bordes es debido a que las imágenes de cultivo provenientes de ambos sensores muestran formas Imagen wv1 Imagen wv2 Imagen wv3 Imagen wv4 Imagen wv5 Imagen wv6 Imagen wv1 Imagen wv2 Imagen wv3 Imagen wv4 Imagen wv5 Imagen wv6 35 geométricas y de esta forma se disminuye la complejidad del registro de imágenes. En el Anexo 4 se puede observar los resultados para otros valores de K, k y L. 4.2 Evaluación del Algoritmo de Registro: En esta parte se evalúa la obtención del vector de traslación de las seis sub bandas en diferentes niveles utilizando el algoritmo propuesto; para ello se observará el error como la diferencia entre el vector de traslación estimado con el real. Los coeficientes de los filtros utilizados cumplían las condiciones de K=1, k=2 y L=1, DT-CWT Figura 4.13 Imagen de prueba Lena. Figura 4.14 Imagen de prueba Lena con ruido y desplazada. Esta imagen fue obtenida en el MatLab para ello se adicionó ruido Gaussiano de varianza 0.025. Figura 4.15 Descomposición de wavelets de la Imagen de prueba de Lena. . 36 De acuerdo a la ecuación (3.33) se obtiene el vector de traslación [tx,ty] para cada direccionamiento (wv1, wv2, wv3, wv4, wv5 y wv6) en cada nivel. Tabla 4.1 Nivel 1 Nivel 2 Nivel 3 wv1 [10,20] [5,10] [2 ,5] wv2 [10,20] [5,10] [3,4] wv3 [10,20] [5,10] [3,4] wv4 [10,20] [5,10] [2,5] wv5 [10,20] [5,10] [2,-1] wv6 [10,20] [5,10] [3,6] Traslación real [10,20] [5,10] [2,5] . En la tabla anterior se muestra la información de los vectores de traslación obtenidos utilizando el algoritmo de registro, además se muestra la traslación real la cual es obtenida ubicando de manera visual puntos pertenecientes a ambas imágenes. Tabla 4.2 Error [tx,ty] Nivel 1 Nivel 2 Nivel 3 wv1 [0,0] [0,0] [0,0] wv2 [0,0] [0,0] [1,-1] wv3 [0,0] [0,0] [1,-1] wv4 [0,0] [0,0] [0,0] wv5 [0,0] [0,0] [0,-6] wv6 [0,0] [0,0] [1,1] DT-CWT Figura 4.16 Descomposición de wavelets de la Imagen de prueba de Lena con ruido. . 37 El error es la diferencia de los parámetros del vector de traslación estimado y real. De los resultados de error mostrados en la tabla 4.2, se puede apreciar que para niveles 1 y 2 en sus distintos direccionamientos se tiene un error de cero. Mientras que el error obtenido en el nivel 3 varia, esto se debe a que en el nivel 3 las imágenes son de dimensiones pequeñas perdiéndose información en común entre las imagen sin ruido y con ruido. 4.3 Resultados de registro de imágenes aéreas de cultivo: Se utilizaron 20 pares de imágenes aéreas de cultivo para evaluar el algoritmo. A continuación se mostrará las etapas (ver figura 3.1) de alineación de las imágenes de cultivo para el par de imágenes de la figura 4.17 y 4.18 y se presentará el resultado para dos pares más; estas imágenes fueron proporcionadas por el CIP. En el Anexo 5 se encuentran los resultados para diecisiete pares de imágenes más. Primer grupo de imágenes: Figura 4.17 Imagen de espectro infrarrojo cercano de cultivo de olivo N0 0147. 38 En la primera etapa se obtiene para cada una de las imágenes seis descomposiciones de wavelets cada una de estas presenta detalles resaltados en seis direcciones diferentes, como se puede observar en las figuras 4.19 (Imagen infrarroja) y 4.20 (Imagen Roja. Imagen wvi1 Imagen wvi2 Imagen wvi3 Imagen wvi4 Imagen wvi5 Figura 4.19 Imagen wvi6 Figura 4.18 Imagen de espectro roja de cultivo de olivo N0 0148. 39 En la segunda etapa se obtiene el vector de traslación para ello se utiliza la correlación, aplicada a cada par de imágenes (wvi y wvr) , en este caso aplicamos al par wvi1 y wvr1 cuando obtenemos la correlación de ambas imágenes en 2 dimensiones ubicamos la coordenada del valor máximo, como se puede observar en la figura 4.21. Luego realizamos la auto correlación de una imagen del par, en este caso se aplicará la auto correlación a wv1. Aquí también ubicamos la coordenada del valor máximo el cual siempre se ubicará en el centro geométrico de la correlación; esto se puede apreciar en la figura 4.22. Después, mediante la diferencia entre las coordenadas podemos encontrar el vector de traslación. wv1=[tx,ty]=(312-307, 293-231)=(5,62) Vector de traslación. La correlación es aplicada a cada par wvi y wvr; por ello obtenemos 6 vectores estos pueden ser diferentes o iguales, pero lo importante de ello es que podemos obtener el vector de traslación correcto dentro de los seis. Imagen wvr1 Imagen wvr2 Imagen wvr3 Imagen wvr4 Imagen wvr5 Figura 4.20 Imagen wvr6 40 Figura 4.21 Figura 4.22 41 Tabla 4.3 Nivel 1 Nivel 2 Nivel 3 wv1 [5, 62] [3, 31] [2, 0] wv2 [5, 62] [3, 31] [1, 16] wv3 [5, 62] [3, 31] [1, 16] wv4 [6, 63] [3, 31] [1, 0] wv5 [7, 64] [3, 32] [0, 0] wv6 [7, 64] [3, 32] [1, 15] Traslación real [6,63] [3,33] [1,14] En la Tabla 4.3 se muestra la información de los vectores de traslación obtenidos utilizando el algoritmo de registro, además se muestra la traslación real la cual es obtenida ubicando de manera visual puntos pertenecientes a ambas imágenes. Tabla 4.4 Error [tx,ty] Nivel 1 Nivel 2 Nivel 3 wv1 [-1,-1] [0,-2] [1,-14] wv2 [-1,-1] [0,-2] [0,2] wv3 [-1,-1] [0,-2] [0,2] wv4 [0,0] [0,-2] [0,-14] wv5 [1,1] [0,-1] [-1,-14] wv6 [1,1] [0,-1] [0,-1] El error es la diferencia de los parámetros del vector de traslación estimado y real. De los resultados de la tabla 4.4 se observa que para el cuarto direccionamientos del nivel 1 se tiene un error de cero, ello quiere decir que las características de las imágenes estaban significativamente en este direccionamientos; sin embargo para los demás niveles aún se tiene el error ±1 ó ±2 pixeles y en algunos casos se tiene un error de 14 pixeles. 42 Resultado Final Del resultado final se observa que el registro de imágenes tiene éxito sin embargo; se pierde parte de la información en ambas imágenes ya que solo se conserva el área común entre ellas. Figura 4.23 Imagen registrada, obtenida con un vector de traslación de [6, 63] de nivel 1. 43 Segundo grupo de imágenes: Figura 4.24 Imagen de espectro infrarrojo cercano de cultivo N0 0087 Figura 4.25 Imagen de espectro rojo de cultivo N0 0088 44 Tabla 4.5 Nivel 1 Nivel 2 Nivel 3 wv1 [23, -32] [11, -16] [0, 0] wv2 [188, 88] [94, 44] [0, 3] wv3 [-22, -12] [12, -16] [-5, -13] wv4 [20, -31] [12, -16] [-3, 0] wv5 [19, -8] [12, 4] [0, 2] wv6 [75, 1] [-42,7] [-20, 13] Traslación real [21,-29] [11,-14] [5,-6] En la Tabla 4.5 se muestra la información de los vectores de traslación obtenidos utilizando el algoritmo de registro, además se muestra la traslación real la cual es obtenida ubicando de manera visual puntos pertenecientes a ambas imágenes. Tabla 4.6 Error [tx,ty] Nivel 1 Nivel 2 Nivel 3 wv1 [2,-3] [0,-2] [-5,6] wv2 [167,117] [83,58] [-5,9] wv3 [-43,17] [1,-2] [0,-7] wv4 [-1,-2] [1,-2] [-8,6] wv5 [-2,21] [1,18] [-5,8] wv6 [54,30] [-53,21] [-25,19] El error es la diferencia de los parámetros del vector de traslación estimado y real. De los resultados de la tabla 4.6 se observa que para el segundo direccionamiento (wv2) en los tres niveles se tiene un error grande esto debido a que no se ha resaltado significativamente las características en esa dirección en ambas imágenes. Mientras que en otros niveles se tiene errores pequeños. 45 Resultado Final De la figura 4.26 se puede apreciar igual que el caso anterior perdida de información, pero la alineación es casi correcta. Para realizar la alineación se tomo el vector traslación [23,-32]. Figura 4.26 Imagen registrada. Obtenida con un vector de traslación de [23,-32] de nivel 1 46 Tercer grupo de Imágenes: Figura 4.28 Imagen de espectro rojo cercano de cultivo N0 3172 Figura 4.27 Imagen de espectro infrarrojo de cultivo N0 3171 47 Tabla 4.7 Nivel 1 Nivel 2 Nivel 3 wv1 [-1,34] [0,16] [9,0] wv2 [-1,34] [0,16] [0,8] wv3 [-2,35] [0,16] [0,8] wv4 [33,21] [0,16] [0,0] wv5 [-2,35] [7,19] [0,5] wv6 [1,31] [3,11] [17,-3] Traslación real [0,8] [0,4] [0,2] En la Tabla 4.7 se muestra la información de los vectores de traslación obtenidos utilizando el algoritmo de registro, además se muestra la traslación real la cual es obtenida ubicando de manera visual puntos pertenecientes a ambas imágenes. Tabla 4.8 Error [tx,ty] Nivel 1 Nivel 2 Nivel 3 wv1 [-1,26] [0,12] [9,0] wv2 [-1,26] [0,12] [0,6] wv3 [-2,27] [0,12] [0,6] wv4 [33,13] [0,12] [0,-2] wv5 [-2,27] [7,15] [0,3] wv6 [1,23] [3,7] [17,-5] El error es la diferencia de los parámetros del vector de traslación estimado y real. De los resultados de la tabla 4.8 se observa que existe un elevado error esto debido a que no se ha resaltado las características significativamente ya que no hay caminos u otras formas que se puedan definir entre las imágenes. 48 Resultado Final Figura 4.29 Imagen registrada. Obtenida con un vector de traslación de [1,31] de nivel 1. Figura 4.30 Imagen registrada. Obtenida con un vector de traslación de [3,11] de nivel 2. 49 De las figuras mostradas se puede apreciar que no se tiene un correcto registro sin embargo, para niveles superiores a 1, el error del vector traslación disminuye. 4.4 Evaluación del Registro de Imágenes con la Transformada de Wavelet Discreta: En esta parte se evalúa el registro de imágenes para tres pares de imágenes utilizando la transformada discreta de wavelets estándar, la wavelet utilizada es la de Daubechies número 3 (W. Db3). Se mostrará el vector traslación para las tres sub bandas: horizontal, vertical y diagonal; además del error entre el vector estimado y el real. En el Anexo 6 se encuentran los resultados para diecisiete pares de imágenes y también el algoritmo utilizado. Imágenes N0 0147 y N0 0148 La tabla 4.9 muestra los vectores de traslaciones de las tres sub bandas que se obtienen horizontal (wh1), vertical (wv2) y diagonal (wd3). Tabla 4.9 Nivel 1 wh1 [5, 63] wv1 [7, 64] wd1 [32, 8] Traslación real [6,63] Figura 4.31 Imagen registrada. Obtenida con un vector de traslación de [0,5] de nivel 3. 50 . Tabla 4.10 Error [tx,ty] Nivel 1 wh1 [-1, 0] wv2 [1, 1] wd3 [26, -55] La tabla 4.10 muestra la diferencia de los parámetros del vector de traslación estimado y real. De acuerdo a los resultados obtenidos de la tabla anterior se observa que las sub bandas de características horizontales y verticales presenta menos error que la sub banda de características diagonales. Resultado del registro El registro de imagen realizado con el vector de traslación de la sub banda de características horizontales muestra un resultado de un registro correcto. Figura 4.32 Imagen registrada. Obtenida con un vector de traslación de [5,63] de la sub banda de características horizontales. . 51 Imágenes N0 0087 y N0 0088 La tabla 4.11 muestra los vectores de traslaciones de las tres sub bandas que se obtienen horizontal (wh1), vertical (wv2) y diagonal (wd3). Tabla 4.12 Error [tx,ty] Nivel 1 wh1 [-55,25] wv2 [-28, -83] wd3 [-148, 15] La tabla 4.12 muestra la diferencia de los parámetros del vector de traslación estimado y real. Resultado del registro Tabla 4.11 Nivel 1 wh1 [-34,- 4] wv2 [-7, -112] wd3 [-127,- 14] Traslación real [21,-29] Figura 4.33 Imagen registrada. Obtenida con un vector de traslación de [-34,-4] de la sub banda de características horizontales. 52 El registro de imagen de la figura anterior muestra un resultado de una alineación incorrecta. Como se puede apreciar en la tabla 4.12 se muestran valores altos de error en las tres sub bandas. Imágenes N0 3171 y N0 3172 Tabla 4.13 Nivel 1 wh1 [-1,8] wv2 [14,22] wd3 [27,9] Traslación real [0,8] La tabla 4.13 muestra los vectores de traslaciones de las tres sub bandas que se obtienen horizontal (wh1), vertical (wv2) y diagonal (wd3). Tabla 4.14 Error [tx,ty] Nivel 1 wh1 [-1,0] wv2 [14,14] wd3 [27,-1] La tabla 4.14 muestra la diferencia de los parámetros del vector de traslación estimado y real. 53 Resultado del registro La figura anterior tiene como resultado un registro correcto. La tabla 4.15 es un resumen de los vectores de traslaciones determinados utilizando el algoritmo de DT-CWT y el algoritmo con wavelets estándar; así también se presenta el error entre el vector estimado y el vector real para los 20 pares de imágenes utilizados. En la tabla se observa que para los pares de imágenes número 1, 8, 11, 12, 13, 14, 16, 17, 18 y 19 el registro utilizando el DT-CWT determinó correctamente los vectores de traslación; también se puede ver que los vectores de las imágenes 2, 5, 6, 7, 9, 10, 15 y 20 presentan un error pequeño menor de 3 pixeles acercándose al valor real. Por otra parte, el registro utilizando wavelets estándar determino correctamente los vectores de traslación de los pares de imagen número 11, 17, 18 y 19; además para las imágenes 1, 3, 5, 6, 12, 13, 14, 15 y 20 se tienen vectores de traslación con un error pequeño menor a 3 pixeles. Figura 4.34 Imagen registrada. Obtenida con un vector de traslación de [-1,8] de la sub banda de características horizontales. 54 Tabla 4.15 Pares de Imágenes Vector de Traslacio n Real [tx,ty] Vector de traslación ( DT-CWT) [tx,ty] Vector de traslación ( W. Db3) [tx,ty] Error de traslación (DT- CWT) [tx,ty] Error de traslación (W. Db3) [tx,ty] 1 0147 y 0148 [6,63] [6,63] [5,63] [0,0] [-1,0] 2 0087 y 0088 [21,-29] [20,-31] [-34,4] [1,-2] [-55,33] 3 3171 y 3172 [0,8] [1,31] [-1,8] [1,23] [-1,0] 4 0241 y 0242 [21,27] [21,57] [43,39] [0,30] [22,12] 5 2419 y 2420 [3,10] [2,12] [3,8] [-1,2] [0,-2] 6 3131 y 3132 [0,11] [-2,13] [-1,12] [-2,2] [-1,1] 7 0034 y 0033 [22,15] [22,12] [0,0] [0,-3] [-22,-15] 8 0100 y 0099 [-10,5] [-10,5] [0,0] [0,0] [10,-5] 9 0238 y 0237 [7,-13] [6,-13] [0,0] [-1,0] [-7,13] 10 0412 y 0413 [-6,9] [-7,10] [0,0] [-1,1] [6,-9] 11 0017 y 0016 [-30,16] [-30,16] [-30,16] [0,0] [0,0] 12 0149 y 0150 [-20,12] [-20,12] [-20,10] [0,0] [0,-2] 13 0131 y 0130 [-21,11] [-21,11] [-20,10] [0,0] [1,-1] 14 0135 y 0134 [-20,11] [-20,11] [-20,10] [0,0] [0,-1] 15 0167 y 0166 [-19,12] [-20,12] [-20,11] [-1,0] [-1,-1] 16 0211 y 0210 [-21,13] [-21,13] [2,-3] [0,0] [23,-16] 17 3109 y 3110 [-14,13] [-14,13] [-14,13] [0,0] [0,0] 18 3687 y 3686 [-30,18] [-30,18] [-30,18] [0,0] [0,0] 19 3743 y 3742 [-29,21] [-29,21] [-29,21] [0,0] [0,0] 20 3779 y 3778 [-26,17] [-28,16] [-27,15] [-2,-1] [-1,-2] 55 CONCLUSIONES  La presente tesis propone un algoritmo basado en Dual tree complex wavelet transform (DT-CWT) para el registro de imágenes aéreas de cultivos, obtenidas de diferentes sensores. Este algoritmo ayuda a que el proceso de obtención de imágenes multiespectrales de los cultivos sea un proceso sencillo.  La característica del DT-CWT de resaltar detalles en seis selectivas direcciones en la descomposición de wavelets contribuye a determinar el vector de traslación utilizando métodos sencillos como la correlación. Además, la obtención de seis vectores de traslación en un nivel de descomposición y la característica de multiresolución, que nos permite descomponer en más niveles, disminuye la probabilidad de cometer inexactitudes al realizar el registro.  Se puede observar en las imágenes de prueba que las características se resaltan significativamente para coeficientes con órdenes de K=1, k=2 y L=1. Mientras que en imágenes con otros órdenes las características se resaltan pero no lo suficiente.  Se realizaron pruebas con 20 pares imágenes aéreas de cultivo (Tabla 4.15) de la cuales el 50% tuvo éxito en el registro, otro 40% presento un error pequeño y solo un 10% tiene un error elevado utilizando el método de DT-CWT. En contraste, las pruebas realizadas con wavelets estándar se obtuvo un 20% de éxito en el registro, un 45% presento un error pequeño y otro 35% tiene un error elevado. Comprobando de esta forma la eficacia de la utilización del registro utilizando el DT-CWT.  En las imágenes utilizadas en las pruebas se puede ver que en las imágenes donde no se observa detalles definidos como formas geométricas, se presentan un alto error; sin embargo, el error se reduce para los niveles siguientes comprobando que la multiresolución contribuye en realizar un registro correcto. 56 RECOMENDACIONES  Para poder realizar el registro exitosamente de imágenes que presentan características de rotación se recomienda emplear otro método aparte de la correlación para asegurar su eficacia, como por ejemplo el de información mutua, propuesto en “A multi- Modal Automatic Image Registration Technique Based on Complex Wavelet” [15].  Para obtener descomposiciones de wavelets con direccionamientos más definidos se recomienda implementar el banco de filtro con el método propuesto en "The Dual-Tree Complex Wavelet Transform “ [18] denominado Oriented 2- Dual Tree CWT .  La obtención de los coeficientes del banco de filtros también puede ser desarrollado utilizando los otros métodos descritos en en "The Dual-Tree Complex Wavelet Transform “ [18], los cuales son Linear Phase Biorthogonal Solution y q- Shift Solution.  Como se sabe la utilización del algoritmo de correlación implementada en su definición matemática disminuye la eficiencia computacional para ello se debe utilizar el método de la correlación utilizando la Transformada Rápida de Fourier. 57 BIBLIOGRAFÍA [1] P. Chavez, C. Yarlequé, O. Piro, A. Posadas, V. Mares, H. Loayza, C. Chuquillanqui, P. Zorosgastúa, J, Flexas y R. Quiroz, “Applying Multifractal Analysis to Remotely Sensed Data for Assessing PYVV Infection in Potato (Solanum tuberosum L.) Crops”, Remote Sensing ,2010. [en línea]. Disponible en: [Consultado el 01/11/2010] [2] (16 Febrero, 2010). “Utilizan aeronaves a control remoto para predicción de cosechas y tomas fotográficas de cultivos”, ANDINA: Agencia peruana de noticias [en línea]. Disponible en: [Consultado el 07/10/2010]. [3] R. Bongiovanni, E. Chartuni, S. Best y A. Roel, Agricultura de Precisión: Integrando conocimientos para una agricultura moderna y sustentable. PROCISUR//IICA [en línea]. Disponible en: [Consultado el 21/10/2010]. [4] P. Kreimer, Las TICS en la Agricultura de Precisión.Madrid.CEDITEC,2003. [5] TETRACAM INC. (2007).“TETRACAM INC: ADC Series” [en línea]. Disponible en: < http://www.tetracam.com/adc_air.html > [Consultado 29/09/2010] [6] A. Wong y D. Clausi, “ARRSI: Automatic registration of remote sensing images,” IEEE Trans. Geosci. Remote Sens., vol. 45, no 5.Part II, pp1483-1493, 2007. [7] A. Wong y D. Clausi, “Automatic Registration of Inter-Band and Inter-sensor Images using Robust Complex Wavelet Feature Representations,” In:Proc. of 5th IAPR Workshop on Pattern Recognition in Remote Sensing (PRRS 2008). 58 [8] P. Viola y W. M. Wells III, “Alignment by Maximization of Mutual Information,” International Journal of Computer Vision, 1997. [9] F. Zhongliang y S. Zhiqun, “An algorithm of Straight line Features matching on aerial Imagery,” The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Beijing ,Vol. XXXVII. Part B3b, 2008. [10] A. Goshtasby, G. C. Stockman Y C. V. Page, “A Region-Based approach to Digital Image registration with Subpixel Accuracy, ” IEEE Transactions on Geosciencie and Remote sensing, Vol. GE-24, Nº 3,1986. [11] D. G. Lowe, “Distinctive Image Features from Scale-Invariant Keypoints,” International Journal of Computer Vision, Computer Science Department, University of British Columbia, 2004. [12 ] H. Bay, T. Tuytelaars, L. V. Gool, “SURF: Speeded Up Robust Features,” European Conference on Computer Vision, pages 404-417,2006. [13] J. Le Moigne,W. Cambell, y R. Cromp, “An Automated Parallel Image Registration Technique Based on the Correlation of Wavelet Features,” IEEE Trans. Geoscience and Remote Sensing, Vol.40, No.8, pp.1849-1864, 2002. [14] J. J. Lewis , R. J. O’Callaghan, S. G. Nikolov, D. R. Bull ,C. N. Canagarajah,” Region-Based Image Fusion Using Complex Wavelets,” In Proc. Seventh Int.Conf on Information Fusion,pp.555-562, 2004. [15] M. Ghantous, S. Ghosh, M. Bayoumi, “A multi-modal automatic image resgistration technique based on Complex Wavelet, ” International Conference on Image Processing (ICIP),16th IEEE International Conf, El cairo,pp 173-176, 2009. [16] M. Weeks, “The Wavelet Transform”, Digital Signal Processing using Matlab and Wavelets , Massachusetts, Infinity Science Press, 2007. 59 [17] D. F. Mix y K. J. Olejniczak, Elements of Wavelets for Engineers and Scientists, New Jersey, John Wiley & Sons Inc.,2003. [18] I. W. Selesnick, R. G. BAraniuk, y N. Kingsbury, “The Dual-Tree Complex Wavelet Transform”,IEEE Signal Processing Magazine,Vol.22,No. 6, pp 123-151, 2005. [19] S.Mallat, A Wavelet tour of signal processing, 2da Edition, Cambridge, Academic Press, 1999. [20] M. Antonini, M.Barlaud, P.Mathieu and I.Daubechies, ”Image Coding Using Wavelet Transform”,IEEE Transactions on Image Processing, Vol 1 , No 2, 1992. [21] G.Hong, “Image Fusion, Image Registration, and Radiometric normalization for High Resolution Image Processing”, Ph.D. dissertation, Department of Geodesy and Geomatics Engineering, Technical Report No. 247, University of New Brunswick, Fredericton, New Brunswick, Canada, 2007. [22] Khalid Sayood, Morgan Kaufmann Publishers, 2nd edition, Burlington, United States of America, 2000. [23] I.W. Selesnick, “The design of approximate Hilbert transform pairs of wavelet bases,” IEEE Trans. Signal Processing, vol. 50, no. 5, pp. 1144–1152, 2002. [24] A. Ardeshir Goshtasby, Introduction 2-D and 3-D Image Registration For Medical, Remote Sensing, and Industrial Applications, John Wiley & Sons, Inc. 2005. [25] J. C. Russ, The Image Processing Handbook, Taylor & Francis Group, 5ta. Edition, 2007.