Estudio y desarrollo de métodos de zoom para imágenes de gammagrafía obtenidas mediante la aplicación de meta-yodobenzilguanidina

Estudio y desarrollo de métodos de zoom para imágenes de gammagrafía obtenidas mediante la aplicación de meta-yodobenzilguanidina Autor: Juan Cartes Font Director: Samuel Morillas Gómez Ingeniería Informática Escuela Técnica Superior de Ingeniería Informática Universitat Politècnica de València Valencia, 28 de septiembre de 2011 2 Resumen El presente proyecto trata de profundizar en las técnicas actuales para proporcionar un zoom digital aplicable a todo tipo de imágenes digitales, especialmente a imágenes médicas como los son las gammagrafías. Este proyecto se encuentra incluido en otro más amplio que tiene por objeto apoyar a la toma de decisiones por parte de los médicos en el proceso de diagnóstico de neuroblastomas. En esta memoria se detallan una serie de técnicas interpolatorias no lineales y se proponen posibles implementaciones, con el objeto de elegir aquélla que mejores resultados proporcione y que pueda ser incorporada al proyecto marco, como una funcionalidad más. En la sección de experimentos, se realizarán diversas pruebas con cada una de las técnicas que serán evaluadas mediante el uso de ciertos indicadores para imágenes digitales, de manera que quede justificada alguna de ellas frente otras para un determinado tipo de imágenes. PALABRAS CLAVE: Interpolación, Neuroblastoma, zoom, imagen médica, gammagrafía. 3 4 Índice general 1. Introducción 15 1.1. Elneuroblastoma............................ 15 1.1.1. Cuadroclínico.......................... 16 1.1.2. Diagnóstico ........................... 16 1.2. Contexto................................. 17 1.3. Objetivos ................................ 17 1.3.1. Filtrado en el dominio de la frecuencia . . . . . . . . . . . . 18 1.3.2. Filtrado en el dominio del espacio o realce . . . . . . . . . . 18 1.3.3. Zoom ............................... 19 1.4. Estructura de la memoria . . . . . . . . . . . . . . . . . . . . . . . 20 2. Imagen médica 21 2.1. LasimágenesDICOM ......................... 21 2.1.1. Historia ............................. 22 2.1.2. Estructura............................ 22 3. Técnicas de interpolación no lineales 31 3.1. Introducción............................... 31 3.2. Multirresolución de Harten . . . . . . . . . . . . . . . . . . . . . . . 33 3.3. Métodos de interpolación no lineales . . . . . . . . . . . . . . . . . 38 3.3.1. Interpolación ENO . . . . . . . . . . . . . . . . . . . . . . . 40 3.3.2. Interpolación ENO Subcell Resolution ............. 47 3.3.3. Interpolación WENO . . . . . . . . . . . . . . . . . . . . . . 54 3.3.4. Interpolación Racional . . . . . . . . . . . . . . . . . . . . . 61 3.3.5. Interpolación PPH . . . . . . . . . . . . . . . . . . . . . . . 65 4. Experimentos y resultados 71 4.1. Indicadores medibles para la evaluación de la calidad . . . . . . . . 71 4.1.1. Introducción........................... 71 4.1.2. Error Cuadrático Medio . . . . . . . . . . . . . . . . . . . . 71 4.1.3. Relación señal a ruido de pico (PSNR) . . . . . . . . . . . . 72 5 6ÍNDICE GENERAL 4.1.4. Correlación cruzada normalizada . . . . . . . . . . . . . . . 73 4.1.5. Diferencia media . . . . . . . . . . . . . . . . . . . . . . . . 73 4.1.6. Diferencia Máxima . . . . . . . . . . . . . . . . . . . . . . . 73 4.1.7. Error Absoluto Medio . . . . . . . . . . . . . . . . . . . . . 73 4.2. Experimentos realizados . . . . . . . . . . . . . . . . . . . . . . . . 74 4.2.1. Introducción........................... 74 4.2.2. Descripción y procedimientos . . . . . . . . . . . . . . . . . 74 4.3. Resultados................................ 79 4.3.1. Resultados con la imagen Dicom1 ............... 81 4.3.2. Resultados con la imagen Dicom2 ............... 84 4.3.3. Resultados con la imagen Lena ................ 87 4.3.4. Resultados con la imagen Geo ................. 90 4.3.5. Resultados con la imagen Tac1 ................ 93 4.3.6. Resultados con la imagen Tac2 ................ 96 4.3.7. Resultados con la imagen Pet-Tc1 .............. 99 4.3.8. Resultados con la imagen Pet-Tc2 ..............102 4.3.9. Resultados con la imagen Pet1 ................105 4.3.10. Resultados con la imagen Pet2 ................108 4.3.11. Análisis de los resultados obtenidos . . . . . . . . . . . . . . 121 5. Conclusiones 129 i. 131 Índice de figuras 2.1. CapasDICOM.............................. 23 3.1. Ilustración del concepto «píramide de resolución.».......... 34 3.2. Representación de tres niveles de resolución, Xk,Xk+1 2yXk+1 respectivamente. Se parte de un nivel de resolución kcon el conjunto de puntos xk i, yk jJk i,j=0. Mediante la técnica del producto tensor se obtienen los valores interpolados para las nuevas filas, es decir los valores para nxk+1 2 i, yk jo0≤i≤Jk+1,0≤j≤Jk correspondientes al nivel de resolución k+1 2. Por último, se obtienen los valores para las columnas alcanzando el nivel de resolución k+1, es decir, con el conjunto de puntos xk+1 i, yk+1 jJk+1 i,j=0. En la figura se aprecian los píxeles originales (J), los detalles interpolados verticales (⊗), los interpolados horizontales (⊕) y los mixtos (?). .......................... 35 4.1. Imágenes anterior (a) y posterior (b) de un paciente, denominadas respectivamente Dicom1 yDicom2. Imágenes fotográficas Lena.jpg (a) y Geo.jpg (b)............................. 77 4.2. Resto de imágenes utilizadas en los experimentos. Lena.jpg (a), Tac1.jpg (b), Tac2.jpg (c), Pet-Tc1.jpg (d), Pet-Ct2.jpg (e), Pet1.jpg (f) y Pet2.jpg (g). ........................... 78 4.3. Imágenes Dicom1 con un nivel de zoom (x2) acotadas por una ventana de 175 ×151 píxeles. Reconstrucción con ENO no jerárquico (a), con ENO jerárquico (b), con ENO-SR (c), con WENO (d), con Racional (e), y con la técnica PPH (f)................. 82 4.4. Imágenes Dicom1 con un nivel de 2(x4) acotadas por una ventana de 175×151 píxeles. Reconstrucción con ENO no jerárquico (a), con ENO jérárquico (b), con ENO-SR (c), con WENO (d), con Racional (e), y con la técnica PPH (f)...................... 83 7 8ÍNDICE DE FIGURAS 4.5. Imágenes Dicom2 con un nivel de zoom (x2) acotadas por una ventana de 174 ×151 píxeles. Reconstrucción con ENO no jerárquico (a), con ENO jérárquico (b), con ENO-SR (c), con WENO (d), con Racional (e), y con la técnica PPH (f)................. 85 4.6. Imágenes Dicom2 con un nivel 2 (x4) de zoom acotadas por una ventana de 174×151 píxeles. Reconstrucción con ENO no jerárquico (a), con ENO jérárquico (b), con ENO-SR (c), con WENO (d), con Racional (e), y con la técnica PPH (f)................. 86 4.7. Imágenes Lena.jpg con un nivel de zoom (x2) acotadas por una ventana de 256×256 píxeles. Reconstrucción con ENO no jerárquico (a), con ENO jérárquico (b), con ENO-SR (c), con WENO (d), con Racional (e), y con la técnica PPH (f)................. 88 4.8. Imágenes Lena.jpg con un factor de zoom 2 (x4) acotadas por una ventana de 255×255 píxeles. Reconstrucción con ENO no jerárquico (a), con ENO jérárquico (b), con ENO-SR (c), con WENO (d), con Racional (e), y con la técnica PPH (f)................. 89 4.9. Imágenes Geo.jpg con un nivel de zoom (x2) acotadas por una ventana de 255 ×255 píxeles. Reconstrucción con ENO no jerárquico (a), con ENO jérárquico (b), con ENO-SR (c), con WENO (d), con Racional (e), y con la técnica PPH (f)................. 91 4.10. Imágenes Geo.jpg con un factor de zoom 2 (x4) acotadas por una ventana de 255×255 píxeles. Reconstrucción con ENO no jerárquico (a), con ENO jérárquico (b), con ENO-SR (c), con WENO (d), con Racional (e), y con la técnica PPH (f)................. 92 4.11. Imágenes Tac1.jpg con un 1 nivel de zoom (x2) acotadas por una ventana de 512×542 píxeles. Reconstrucción con ENO no jerárquico (a), con ENO jérárquico (b), con ENO-SR (c), con WENO (d), con Racional (e), y con la técnica PPH (f)................. 94 4.12. Imágenes Tac1.jpg con un factor de zoom 2 (x4) acotadas por una ventana de 512×542 píxeles. Reconstrucción con ENO no jerárquico (a), con ENO jérárquico (b), con ENO-SR (c), con WENO (d), con Racional (e), y con la técnica PPH (f)................. 95 4.13. Imágenes Tac2.jpg con un 1 nivel de zoom (x2) acotadas por una ventana de 512×512 píxeles. Reconstrucción con ENO no jerárquico (a), con ENO jérárquico (b), con ENO-SR (c), con WENO (d), con Racional (e), y con la técnica PPH (f)................. 97 4.14. Imágenes Tac2.jpg con un factor de zoom 2 (x4) acotadas por una ventana de 512×512 píxeles. Reconstrucción con ENO no jerárquico (a), con ENO jérárquico (b), con ENO-SR (c), con WENO (d), con Racional (e), y con la técnica PPH (f)................. 98 ÍNDICE DE FIGURAS 9 4.15. Imágenes Pet1.jpg con un nivel de zoom (x2) acotadas por una ventana de 961×512 píxeles. Reconstrucción con ENO no jerárquico (a), con ENO jérárquico (b), con ENO-SR (c), con WENO (d), con Racional (e), y con la técnica PPH (f).................100 4.16. Imágenes Pet-Tc1.jpg con un factor de zoom 2 (x4) acotadas por una ventana de 961 ×512 píxeles. Reconstrucción con ENO no jerárquico (a), con ENO jérárquico (b), con ENO-SR (c), con WENO (d), con Racional (e), y con la técnica PPH (f)............101 4.17. Imágenes Pet-Ct2.jpg con un nivel de zoom (x2) acotadas por una ventana de 516×542 píxeles. Reconstrucción con ENO no jerárquico (a), con ENO jérárquico (b), con ENO-SR (c), con WENO (d), con Racional (e), y con la técnica PPH (f).................103 4.18. Imágenes Pet-Ct2.jpg con un factor de zoom 2 (x4). Reconstrucción con ENO no jerárquico (a), con ENO jérárquico (b), con ENO-SR (c), con WENO (d), con Racional (e), y con la técnica PPH (f). . . 104 4.19. Imágenes Pet1.jpg con un nivel de zoom (x2) acotadas por una ventana de 255×255 píxeles. Reconstrucción con ENO no jerárquico (a), con ENO jérárquico (b), con ENO-SR (c), con WENO (d), con Racional (e), y con la técnica PPH (f).................106 4.20. Imágenes Pet1.jpg con un factor de zoom 2 (x4) acotadas por una ventana de 255×255 píxeles. Reconstrucción con ENO no jerárquico (a), con ENO jérárquico (b), con ENO-SR (c), con WENO (d), con Racional (e), y con la técnica PPH (f).................107 4.21. Imágenes Pet2.jpg con un nivel de zoom (x2) acotadas por una ventana de 400×256 píxeles. Reconstrucción con ENO no jerárquico (a), con ENO jérárquico (b), con ENO-SR (c), con WENO (d), con Racional (e), y con la técnica PPH (f).................109 4.22. Imágenes Pet2.jpg con un factor de zoom 2 (x4) acotadas por una ventana de 400×256 píxeles. Reconstrucción con ENO no jerárquico (a), con ENO jérárquico (b), con ENO-SR (c), con WENO (d), con Racional (e), y con la técnica PPH (f).................110 4.23. Gráfico que muestra los valores del PSNR obtenidos para cada técnica y para cada nivel, aplicados sobre la imagen Dicom1. . . . . . . 111 4.24. Gráfico que muestra los valores del PSNR obtenidos para cada técnica y para cada nivel, aplicados sobre la imagen Dicom2. . . . . . . 111 4.25. Gráfico que muestra los valores del PSNR obtenidos para cada técnica y para cada nivel, aplicados sobre la imagen Lena.jpg. . . . . . 112 4.26. Gráfico que muestra los valores del PSNR obtenidos para cada técnica y para cada nivel, aplicados sobre la imagen Geo.jpg. . . . . . 112 16 CAPÍTULO 1. INTRODUCCIÓN ya se ha metastatizado hacia otros órganos. Debido a la temprana aparición, muchos de los estudios realizados se han centrado en encontrar factores paternos relacionados con la concepción y la gestación. Estos factores incluyen, entre otros, la exposición a productos químicos en industrias específicas, el tabaquismo, el consumo de licor, el uso de medicamentos o fármacos durante el embarazo, además de otros factores relativos al nacimiento; sin embargo los resultados no han sido concluyentes.[1] 1.1.1. Cuadro clínico El neuroblastoma puede ubicarse en cualquier punto a lo largo del sistema simpático y ello provoca que los signos posibles aparezcan en función de la ubicación del mismo. Entre muchos otros síntomas típicos, se hacen frecuentemente patentes dolores en los huesos, ojos protuberantes, hipertensión arterial y otros signos paraneoplásticos, fiebre elevada, anorexia, así como la existencia de masa palpable en el abdomen, el cuello o el tórax. 1.1.2. Diagnóstico El diagnóstico ha de ser confirmado por un patólogo quirúrgico, teniendo en cuenta la presentación clínica, los hallazgos microscópicos y otras pruebas que pudieran haber sido llevadas a cabo en el laboratorio. Este tipo de pruebas son variables, y entre otras, podemos destacar las siguientes: Análisis de orina: En un 90% de los casos en los que se presenta este tipo de cáncer, las células tumorales producen niveles elevados de ciertas hormonas. El cuerpo las convierte en ácidos (llamados HVA y VMA) que se excretan en la orina. Si los resultados son elevados, esta prueba puede ser una manera fácil de seguir la enfermedad y la reacción de un niño al tratamiento. También pueden ser fácilmente medidas después de que la terapia termine para determinar si la enfermedad esta menguando. El problema de esta técnica es que no todos los pacientes han de tener necesariamente altos niveles de concentración de HVA y VMA. Análisis de sangre: Mediante un conteo sanguíneo completo o (CBC) se puede revisar si las cuentas de sangre están bajas, efecto de el crecimiento del tumor en la médula. Asimismo mediante un panel de química sanguínea puede ser supervisada la función de los riñones y el hígado; y buscar la presencia de algunas sustancias, a modo de indicadores, que puedan aumentar como resultado del crecimiento de un posible tumor. 1.2. CONTEXTO 17 Biopsia: Para realizar un diagnóstico de neuroblastoma, se requiere una muestra del tumor del paciente mediante cirugía. El patólogo examina esta parte del tumor extirpado y determina sus características, en aras de facilitar un tratamiento más acertado. Neuroimágenes: Actualmente los estudios de diagnóstico mediante imágenes se realizan para obtener imágenes del interior del cuerpo del paciente con el objetivo de determinar la exacta ubicación de los diferentes tumores. La exploración por gammagrafía con MIBG (meta-yodobenzilguanidina), es la sustancia idónea que actúa como marcador. El MIBG es absorbido por las neuronas simpáticas en un funcionamiento análogo al del neurotransmisor norepinefrina. Las exploraciones con MIBG se han convertido en una delicada y precisa técnica de buscar la propagación de neuroblastomas. Se pueden utilizar dosis más altas de yodo radiactivo para concentrarlo en las células tumorales, lo que da una forma de radioterapia muy localizada, que puede llegar a matar el cáncer. Se trata de un nuevo enfoque al tratamiento que es utilizado cada vez mas en niños en estados avanzados, o en aquéllos que el neuroblastoma ha reincidido después de un tratamiento convencional. 1.2. Contexto Este estudio se encuentra enmarcado dentro de otro proyecto mucho más amplio, que tiene por objeto el diseño de un software que ayude a los médicos en el proceso de toma de decisiones en lo que al diagnóstico, seguimiento y al tratamiento de los pacientes afectados por el neuroblastoma se refiere. Este proyecto marco es fruto de la colaboración entre el Hospital La Fe de Valencia y la Universitat Politècnica de València, y pretende responder a las necesidades que el personal facultativo ha requerido en los sucesivos encuentros mantenidos con el Departamento de Matemática Aplicada para el establecimiento de los requisitos. Así pues, los resultados que en esta memoria se recogen sobre la elección de una técnica interpolatoria para llevar a cabo un zoom digital, serán tratados y analizados a la hora de incluirlos en el mencionado software como una funcionalidad más. 1.3. Objetivos El objetivo del presente estudio, es el de estudiar las técnicas de zoom aplicadas a imágenes digitales. El zooming digital está enmarcado dentro del procesamiento digital de imágenes, que es el conjunto de todas aquellas técnicas aplicadas a las 18 CAPÍTULO 1. INTRODUCCIÓN imágenes digitales, y destinadas a mejorar la calidad de las imágenes o a facilitar la búsqueda de información en ellas. De entre las diversas técnicas de procesamiento digital, se pueden destacar el filtrado, el realce y el zoom. 1.3.1. Filtrado en el dominio de la frecuencia El filtrado digital de imágenes se basa en la operación de convolución entre una imagen y una función filtro. El cambio de dominio de la imagen, del espacio de descripción al frecuencial, permite sustituir las convoluciones por productos, con ventajas para el proceso de cálculo. Este tipo de filtrado permite mayor flexibilidad ya que hace posible seleccionar no solo la dirección de filtrado, sinó también los intervalos de frecuencia que han de ser eliminados. El filtrado en el dominio de la frecuencia es sencillo, poderoso y flexible. A grandes rasgos se trata de aplicar una determinada máscara o función de filtrado sobre una función, en este caso sobre una imagen en el dominio de la frecuencia. Dependiendo del tipo de filtro empleado, se eliminaran unas frecuencias u otras, alterando el espacio frecuencial de la imagen de destino. Como ejemplos de altas frecuencias se pueden citar los bordes, las líneas así como el ruido en ciertas imágenes. En contraposición, las bajas frecuencias son producidas por los cambios graduales de brillo en la imagen. 1.3.2. Filtrado en el dominio del espacio o realce El realce, como parte integrante del procesamiento digital de imágenes, comprende una serie de operaciones que tienen por objeto mejorar la calidad de las imágenes. Estas operaciones permiten realzar las características de brillo y contraste de una imagen, reducir su contenido de ruido, o agudizar ciertos detalles que se puedan presentar en ella. Estas mencionadas operaciones que componen la técnica del realce, pueden ser divididas en dos grupos, según su tipo de procesamiento. Por un lado están aquellas operaciones de procesamiento puntual o de «píxel por píxel»; y por otro, aquellas de procesamiento por grupo de píxeles, o también llamadas «sobre vecindades». Se parte de dos imágenes disponibles; una imagen de entrada, cuyos datos serán procesados y, una imagen de salida, que será el resultado de el realce. El primer tipo de operaciones, tiende a mejorar el contraste tonal de la imagen, es decir, mejoran la diferencia entre los valores más oscuros y los más claros que se visualizan en un monitor. Este procesamiento altera los niveles de gris de los píxeles de una imagen. En la imagen de entrada, cada píxel es modificado por un 1.3. OBJETIVOS 19 nuevo valor mediante una serie de operaciones matemáticas o relaciones lógicas. El valor resultante es colocado en la imagen de salida ocupando la misma posición que poseía en la imagen de entrada. De ahí que reciba el nombre de «píxel por píxel», ya que la transformación sucede a nivel individual y los píxeles en posiciones vecinas no tienen ningún tipo de influencia. El segundo tipo, las llamadas operaciones de procesamiento por vecindades, mejoran el contraste espacial de la imagen, esto es, la diferencia entre el valor digital de brillo de un determinado píxel y la de sus vecinos. El objetivo es suavizar o reforzar estos contrastes espaciales de manera que los valores de brillo de cada píxel se asemejen o se distancien (en términos de brillo) más o menos respecto de sus vecinos. Como se ve, este tipo de operaciones opera sobre un conjunto de píxeles de la imagen de entrada, para producir el valor de un solo píxel en la imagen de salida; mediante la valiosa aportación de sus vecinos. 1.3.3. Zoom El zoom digital es un método para disminuir el ángulo de visión de una imagen digital. Se logra recortando una imagen con el mismo radio de aspecto que la original, e interpolando el resultado. En contraposición al denominado «zoom clásico», el zoom digital puede lograr cualquier aumento aunque este es directamente proporcional a la pérdida de calidad. Las técnicas de ampliación o zooming, emplean polinomios para averiguar el valor de los «nuevos» píxeles que carecen de valor asignado, al redimensionar la matriz original en un cierto factor. En la elección de estos polinomios es donde se pone de manifiesto la linealidad ono linealidad de las técnicas de ampliación. De esta manera, un algoritmo de zooming lineal es aquél que se basa en algún polinomio lineal de interpolación, por ejemplo el de Lagrange, para obtener el valor de los píxeles desconocidos a priori. Este tipo de técnicas siempre se aplica de igual manera, y no tiene en cuenta las particularidades que pueda tener una determinada imagen. Análogamente, un algoritmo zooming no lineal hace uso de una técnica no lineal de interpolación, es decir una técnica que tenga en cuenta las discontinuidades de la imagen que se está tratando. Es por ello que éstos últimos obtienen, en principio, una mayor calidad de imagen, ya que la obtención del valor de un determinado píxel puede variar, dependiendo de la imagen de la que se parta. Es por ello que en este estudio se ha trabajado sobre éstas últimas técnicas de interpolación. 20 CAPÍTULO 1. INTRODUCCIÓN 1.4. Estructura de la memoria La presente memoria ha sido estructurada en cinco capítulos con el objetivo de facilitar tanto la necesidad y motivaciones por las que este proyecto ha surgido, como la aproximación a las técnicas planteadas para proporcionar soluciones al problema planteado. En la primera sección se exponen las motivaciones que desencadenan la realización de este proyecto. También son descritas las condiciones contextuales en las que este se halla enmarcado, así como los objetivos y la estructura que regirá la memoria de este proyecto. En el segundo capítulo, se profundiza en el contexto de la imagen médica y se relata brevemente la situación que a lo largo de los años ha requerido la creación de estándares en este contexto. En este bloque se exponen también algunas características técnicas de las imágenes de gammagrafía, profundizando en aquellos aspectos que conciernen más directamente al objeto de estudio. A continuación, se encuentra un capítulo dedicado a las técnicas interpolatorias no lineales, y al marco teórico que subyace en las técnicas de zooming sobre imágenes digitales, es decir, a la Multirresolución de Harten. Aunque son explicados de forma teórica, se pueden encontrar detallados los algoritmos de las citadas técnicas en sus respectivas secciones. El cuarto capítulo comprende la documentación relativa a los experimentos llevados a cabo de forma práctica, para evaluar de manera directa el impacto de las técnicas anteriormente descritas, y tienen por objeto obtener indicadores medibles, a fin de cuantificar el beneficio que éstas puedan aportar. Por último, se encuentra la sección de conclusiones donde se exponen las ideas que tras la realización de este proyecto han surgido, y otras que pudieran surgir pero que exceden de los objetivos previamente fijados en este trabajo. A este último capítulo le sigue un anexo, donde se puede encontrar documentación relativa implementación que se ha realizado de las distintas técnicas propuestas. Capítulo 2 Imagen médica Recibe el nombre de imagen médica el conjunto de «técnicas y procesos usados para crear imágenes del cuerpo humano, o partes de él, con propósitos clínicos o para la ciencia médica»1. En el campo de la investigación científica, la imagen médica constituye una subdisciplina de la ingeniería biomédica, la física médica o la medicina, dependiendo del contexto de estudio. Este contexto es muy amplio y comprende actividades como la investigación el desarrollo en el área de instrumentación, adquisición de imágenes, el modelado y la cuantificación son normalmente reservadas para la ingeniería biomédica, física médica y ciencias de la computación; la investigación en la aplicación e interpretación de las imágenes médicas se reserva normalmente a la radiología y las subdisciplinas médicas relevantes en la enfermedad médica o área de la ciencia médica bajo investigación. 2.1. Las imágenes DICOM DICOM (Digital Imaging and COmmunication in Medicine) es un estándar reconocido mundialmente para el intercambio de imágenes médicas para el almacenamiento, manipulación, impresión y transmisión de imágenes médicas. Nació como un acuerdo entre la ACR2(American College of Radiology) y la NEMA3 (National Electrical Manufacturers Association) ante la necesidad inminente de interconectar distintos aparatos de adquisición de imagen radiológica, ya que en aquel momento cada equipo de adquisición contaba hasta entonces con su propio protocolo propietario. 1http://es.wikipedia.org/wiki/Imagen_m %C3 %A9dica 2http://www.rheumatology.org/ 3http://www.nema.org/ 21 22 CAPÍTULO 2. IMAGEN MÉDICA 2.1.1. Historia En 1983, el ACR y la NEMA formaron un comité cuya misión era diseñar y desarrollar una interfaz entre el equipamiento existente y cualquier otro dispositivo que el usuario quisiera conectar. Además de las especificaciones para la conexión del hardware, el estándar sería desarrollado para permitir además la inclusión de un diccionario de los elementos de datos necesarios para la interpretación y la manipulación de imágenes. Debido a todo ello, en 1985 surgió la primera versión del estándar y tres años después se lanzó lanzaría la segunda. El principal problema de esta nueva versión era que los usuarios requerían una interfaz entre los distintos dispositivos, y una red, el protocolo de la cual no poseía la robustez necesaria para soportar las comunicaciones necesarias. Este problema propició el posterior rediseño del proceso en su totalidad, dando lugar la tercera versión del estándar, el DICOM 3.0; cuya división en capas podemos ver en la figura (2.1). Con la aparición de los ordenadores y la tecnología de la imagen digital (TAC, Radiología Digital, PET, SPECT,.. . ) fueron desarrollados diversos sistemas con la intención de integrar el historial clínico del paciente y las diferentes pruebas que se le hubieran desarrollado para contribuir a un diagnóstico más aproximado. Estos desarrollos desembocaron en lo que hoy se conoce como PACS (Picture Archiving and Communication Systems), sistemas informáticos que aportan nuevos modos de trabajo a la radiología diagnóstica. Tienen por objetivo final el de permitir el funcionamiento de un servicio de radiología integrando las imágenes y la información clínica. Constan de un sistema central de gestión y archivo, y de diferentes sistemas de adquisición, visualización y archivo de imágenes, unidos por redes de comunicaciones. El problema de interconexión entre éstos equipos de naturaleza heterogénea quedaba solventado así gracias a la tercera versión del estándar DICOM. 2.1.2. Estructura El formato de un fichero DICOM es muy complejo, debido a la gran cantidad de campos que se especifican en la cabecera, así como los diferentes tipos de cabecera que permite, y la multitud de formatos en los que puede estar grabada la imagen. El fichero DICOM se puede dividir en: 1. Un preámbulo y prefijo identificativo del fichero. 2. Una meta-cabecera. 3. Una cabecera. 2.1. LAS IMÁGENES DICOM 23 Figura 2.1: Capas DICOM. 4. La imagen propiamente dicha (un elemento más de la cabecera según el punto de vista de la cabecera). Preámbulo El estándar DICOM especifica que un fichero en formato DICOM ha de comenzar necesariamente con un preámbulo. Éste tiene un tamaño fijo de 128 bytes, y su uso es dependiente de la implementación. Tampoco está especificada la manera en la que los datos han de ser estructurados, delegando ésta decisión a de los encargados del diseño de la implementación. En caso de que no se haga uso de él, debe estar presente con todos sus bytes puestos al valor 00h. Prefijo Se designa prefijo identificativo a aquel conjunto de datos que sigue al preámbulo. Este prefijo consta de cuatro bytes que contienen la cadena de caracteres DICOM. Esta cadena debe estar codificada siempre con las letras en mayúscula, y usando el conjunto de caracteres especificados en la ISO 8859G0. El propósito de dicho prefijo es permitir a las implementaciones diferenciar si un fichero está o no en formato DICOM. 24 CAPÍTULO 2. IMAGEN MÉDICA Elementos de datos El resto de elementos (cabecera y meta-cabecera) consisten en una serie de campos con toda la información sobre la imagen, incluyendo a ésta. En estos campos se encuentra información de muy distinta naturaleza; aunque los más interesantes y valor añadido poseen, desde el punto de vista técnico, son aquellos que contienen información para el procesado y la visualización de la imagen. Al conjunto de la información codificada sobre un campo se le conoce con el nombre de Elemento de Datos o Data Element. A continuación se expondrá cómo se codifican estos Elementos de Datos, paso previo para la descripción posterior de la cabecera y la meta-cabecera. Un Elemento de Datos está definido por los siguientes campos: Etiqueta del Elemento de Datos (Data Element Tag): Su misión es la de identificar cada elemento de datos de forma unívoca. Una etiqueta está constituida por un Número de Grupo (Group Number) y un Número de Elemento (Element Number). En la documentación del estándar están las descripciones de todos los Elementos de Datos, ordenados según ésta etiqueta. Asimismo se explica el propósito de cada uno de ellos y su requerida obligatoriedad o no. Suelen ser representados como un vector de dos dimensiones, en cuya primera dimensión se encuentra el Número de Grupo y en la segunda el Número de Elemento, en hexadecimal, mediante cuatro dígitos. Representación del Valor (Value Representation, VR): Indica la forma en que se codifica el valor del elemento. Este campo no siempre está codificado en un Elemento de Datos, sinó que depende de la sintaxis de transferencia. Longitud del Valor (Value Length): Entero que se corresponde con la longitud del campo Valor. Valor (Value): es el valor del elemento de datos codificado según el campo VR y con la longitud que indica el campo Longitud del Valor. I. Campos Todos los campos definidos por DICOM se encuentran listados en una base de datos que se encuentra en el documento número seis del estándar, y se la conoce como Registro de los Elementos de Datos DICOM (Registry of DICOM Data Elements). Cada elemento está indexado por su etiqueta (Número de Elemento yNúmero de Grupo), y para cada uno de ellos se halla especificado: Nombre: Nombre del elemento y pequeña descripción de su función. 2.1. LAS IMÁGENES DICOM 25 VR: Representación del valor de cada elemento. VM: Cantidad de valores del mismo tipo que puede contener el campo Valor del elemento de datos. También, en caso de que el un elemento de datos esté obsoleto y haya sido retirado en una versión actual del estándar, poseerá el identificador RET. DICOM establece la obligatoriedad de cada uno de sus campos mediante una clasificación basada en tipos. Tipo 1: Este tipo es de inclusión obligatoria. La longitud del campo no puede ser cero, y debe tener un valor válido. Tipo 1C: Tipo de inclusión obligatoria siempre que se den ciertas condiciones. Si éstas tienen lugar, el elemento es, a todos los efectos, perteneciente al grupo 1. Tipo 2: Este tipo también es de inclusión obligatoria, con la salvedad que puede tener una longitud de campo igual a cero y sin campo Valor. Ésto último solo tiene lugar bajo varios supuestos específicos. Tipo 2C: Análogamente al Tipo1, existe un Tipo 2C, que equivale al Tipo 2 solo bajo la existencia de ciertas condiciones. Tipo 3: Este tipo de campos es opcional y carente de las restricciones de los otros tipos. II. Representación del valor El estándar DICOM define una serie de VR con diferentes características, con la intención de que el campo Valor de cada Elemento de Datos esté codificado correctamente según aquello que represente. El listado con las descripciones resumidas puede ser encontrado en las tablas (2.1) y (2.2). 32 CAPÍTULO 3. TÉCNICAS DE INTERPOLACIÓN NO LINEALES El efecto escalera se debe a una ineficiente eliminación de las repeticiones del espectro que introduce el insertador de ceros. El rizado se corresponde con oscilaciones en la amplitud de la imagen interpolada que no están presentes en la imagen original. La razón de ello es que la función interpolante no decrece de forma monótona a medida que |t| crece, sino que presetna oscilaciones que, en general, decrecen a medida que nos alejamos del origen. La multiresolución de Harten es una herramienta para el procesamiento de imágenes. El objetivo de ésta técnica es establecer un marco para llevar a cabo las transformaciones entre distintos niveles de multirresolución, utilizando una serie de operadores. Estos operadores están íntimamente relacionados con la reconstrucción y la discretización de la función objeto de estudio, y permiten conectar diferentes niveles discretos de resolución con un espacio funcional adecuado, el cual es dependiente de las aplicaciones. Es en el operador reconstrucción el que adquiere mas importancia en nuestro caso, porque será el que implemente alguna técnica de interpolación no lineal que determinará la calidad de la aproximación empleada. En una primera aproximación, podríamos pensar que dichas singularidades se verían resueltas aumentando el orden la función, pero, en caso de implementarlo, se puede ver que la discontinuidad acaba por afectar a un mayor numero de conjuntos de puntos, denominados stencils, aumentando una zona de la función o de la imagen en este caso, donde la calidad no es óptima. El punto crucial radica en la selección de los nodos adecuados para construir el interpolante, de modo que el conjunto de puntos elegido no contenga ninguna singularidad. Un primer acercamiento a la solución del problema implicaría la utilización del algoritmo ENO, en el cual el polinomio interpolado se construye tomando información de las zonas donde la función interpolada es suave. De esta manera, si las singularidades de la función están lo suficientemente aisladas, es posible reducir la zona donde la aproximación se ve degradada, al intervalo que lo contiene. Si se conoce la localización exacta de la singularidad, se puede acotar la pérdida de exactitud a un entorno alrededor de la singularidad mediante el algoritmo ENO-SR. Seguidamente, la técnica interpolatoria WENO constituye a priori una mejora notable de la técnica ENO. Consiste en construir la función interpolante mediante combinaciones convexas de todas las aproximaciones obtenidas a partir de stencils que contienen el intervalo a interpolar, de modo que en la combinación se priman las aproximaciones de aquellos puntos de zonas suaves y si la función es suave en todos ellos, se obtiene una aproximación de orden óptimo. A continuación se presentarán las técnicas Racional y PPH. La primera se considera como una modificación de la técnica WENO con una selección particular 3.2. MULTIRRESOLUCIÓN DE HARTEN 33 de los pesos; mientras que la segunda, se detalla como una interpolación con idénticos resultados que la interpolación lineal en regiones suaves, y con resultados aceptables en regiones en las cuales se halla alguna singularidad presente. 3.2. Multirresolución de Harten La multirresolución de Harten es una herramienta muy eficaz para el procesamiento de imágenes. El objetivo de ésta técnica es obtener una reordenación multiescala de la información contenida en un conjunto de datos discretos; y el resultado puede ser interpretado como una aproximación de la información inicial en un nivel de resolución menor, más unos detalles que en principio nos permiten recuperar datos iniciales. Formalmente partimos de un espacio Vk, en el que kindica el nivel de resolución, y de una función fperteneciente a dicho espacio. Un mayor valor de k indica un mayor nivel de resolución. En el caso que nos ocupa, la reconstrucción mediante valores puntuales, podemos considerar que los datos discretos son valores puntuales en una malla dada. La multirresolución se apoya en los operadores decimación y predicción que permiten la transición entre dos niveles consecutivos de resolución. Ambos se definen como sigue: Decimación: Proporciona información discreta a un nivel de resolución k−1, a partir de un nivel de resolución k. Formalmente se denota por Dk−1 k:Vk→ Vk+1 k. Predicción: Es el operador que porporciona una aproximación discreta del nivel ka partir de la información contenida en un nivel k+ 1 y al que además no se le exige que sea lineal. Siguiendo la notación, se denotará como Pk k−1:Vk−1→Vk. Los datos discretos se obtienen a partir de la discretización de una función f, para lo cual existen distintos tipos de operadores. Dependiendo del operador discretización utilizado, la secuencia de datos fkes diferente. El objetivo del enfoque propuesto por Harten es la construcción de esquemas multirresolución adaptados a cada proceso de discretización. Esto se consigue definiendo un operador reconstrucción apropiado. Estos dos últimos operadores, son los elementos a partir de los cuales se construyen los operadores de decimación y predicción del esquema de multirresolución. Para entender la terminología empleada por los distintos operadores convendría consultar el ejemplo 3.2. 34 CAPÍTULO 3. TÉCNICAS DE INTERPOLACIÓN NO LINEALES Figura 3.1: Ilustración del concepto «píramide de resolución.» Para definir estos operadores formalmente, consideremos Fun espacio de funciones: F ⊂ {f|f: Ω ⊂Rm−→ R} Se define el operador discretización,Dk, como aquél operador que asigna a cada elemento de este espacio, f∈ F, una secuencia de fkde datos discretos perteneciente al espacio Vk. De este modo el operador discretización: Dk:F → Vk=Dk(f) que ha de ser lineal y sobreyectivo, y que a cada f∈ F le asocia: fk=Dk(f) La reconstrucción ofrece la equivalencia en sentido inverso, tomando una secuencia de datos discretos para reconstruir, a partir de la información proporcionada por dichos datos, la función de la cual provienen: Rk:Vk→ F A este operador no se le exige que sea lineal, ésta es la principal novedad introducida por Harten. Los operadores decimación y reconstrucción deben verificar una condición de consistencia, la cual pretende asegurar que la reconstrucción de un conjunto discreto de datos contenga exactamente la misma información que el conjunto inicial de datos. DkRkvk,∀vk∈Vk,es decir DkRk=IVk.(3.1) 3.2. MULTIRRESOLUCIÓN DE HARTEN 35 Xk=    JJJJJ JJJJJ JJJJJ JJJJJ     Xk+1 2=    J⊗J⊗J⊗J⊗J J⊗J⊗J⊗J⊗J J⊗J⊗J⊗J⊗J J⊗J⊗J⊗J⊗J     Xk+1 =              J⊗J⊗J⊗J⊗J ⊕?⊕?⊕?⊕?⊕ J⊗J⊗J⊗J⊗J ⊕?⊕?⊕?⊕?⊕ J⊗J⊗J⊗J⊗J ⊕?⊕?⊕?⊕?⊕ J⊗J⊗J⊗J⊗J ⊕?⊕?⊕?⊕?⊕ J⊗J⊗J⊗J⊗J               Figura 3.2: Representación de tres niveles de resolución, Xk,Xk+1 2yXk+1 respectivamente. Se parte de un nivel de resolución kcon el conjunto de puntos xk i, yk jJk i,j=0. Mediante la técnica del producto tensor se obtienen los valores interpolados para las nuevas filas, es decir los valores para nxk+1 2 i, yk jo0≤i≤Jk+1,0≤j≤Jk correspondientes al nivel de resolución k+1 2. Por último, se obtienen los valores para las columnas alcanzando el nivel de resolución k+ 1, es decir, con el conjunto de puntos xk+1 i, yk+1 jJk+1 i,j=0. En la figura se aprecian los píxeles originales (J), los detalles interpolados verticales (⊗), los interpolados horizontales (⊕) y los mixtos (?). 36 CAPÍTULO 3. TÉCNICAS DE INTERPOLACIÓN NO LINEALES Dada una secuencia de operadores discretización {Dk}, y de reconstrucción {Rk}se definen los operadores decimación y predicción de la siguiente manera: Dk−1 k=Dk−1Rk, Pk k−1=DkRk−1. Como se observa, se satisface la restricción de consistencia establecida en (3.1). Se dice que una sucesión de operadores discretización, Dk, es anidada si cumple la siguiente expresión Aunque a priori parezca que el operador decimación dependa de la elección del operador reconstrucción, diremos que una operación de discretización es ennidada si verifica: Dkf= 0 ⇒ Dk−1f= 0,∀f∈ F (3.2) Esta propiedad significa que la información contenida en los datos a un cierto nivel de resolución kno será nunca mayor que la información contenida en un nivel de resolución superior. En caso de que se cumpla esta propiedad, se tendrá la garantía de que el operador decimación será independiente del operador reconstrucción. A modo de demostración, si consideramos dos secuencias de operadores resconstrucción, DkyD0 kque verifican la ecuación (3.1), se tiene: Dk−1(Rkvk−R0 kvk) = DkRkvk−DkR0 kvk=vk−vk= 0,∀vk∈Vk.(3.3) Verificando que ambos operadores son independientes, como sigue: Dk−1(Rkvk−R0 kvk) = 0 ⇒ Dk−1Rkvk=Dk−1R0 kvk,∀vk∈Vk.(3.4) A partir de las definiciones (3.2), (3.3) y (3.4) se deduce la relación de consistencia para los operadores decimación y predicción, análogamente a la ecuación (3.1). Si decimamos la información obtenida a partir de la predicción realizada sobre una información con resolución dada por Vk−1, obtenemos exactamente la misma información de partida, sin haber introducido ningún elemento nuevo. Dk−1 kPk k−1=Dk−1RkDkRk−1=Dk−1Rk−1=IVk−1(3.5) Si denotamos por vka aquella información discreta en un nivel de resolución k, al aplicarle el operador decimación sobre ella, obtenemos vk−1, es decir, la información contenida en el nivel de resolución k−1: vk−1=Dk−1 kvk 3.2. MULTIRRESOLUCIÓN DE HARTEN 37 Dado que Pk k−1Dk−1 kvkconstituye una aproximación a vk, el error queda definido como sigue: ek=vk−Pk k−1Dk−1 kvk= (Ik V−Pk k−1Dk−1 k)vk=Qkvk∈Vk. De esta manera, conocido vk−1Dk−1 kvk∈Vkyekpuede ser recuperado vk, conteniendo la misma información tanto el conjunto vkcomo el conjunto vk−1, ek, es decir: vk≡vk−1, ek(3.6) haciendo obvia la siguiente relación vk=Pk k−1vk+ek. El problema es que siguiendo este procedimiento se tiene información redundante, pues si si suponemos Vkes un espacio de dimensión finita, dimV k=Nk, resulta que vk−1, ekconsta de Nk−1+Nkelementos, aun conteniendo vk−1, ek yvkla misma información. Esta información redundante puede ser eliminada, como sigue: Dk−1 kek=Dk−1 k(Ik V−Pk k−1Dk−1 k)vk =Dk−1 kvk−Dk−1 kPk k−1Dk−1 kvk =Dk−1 kvk−Dk−1 kvk= 0. es decir, ek∈N(Dk−1 k) = vk∈Vk:Dk−1 kvk= 0cuya dimensión es dimN(Dk−1 k) = dimV k−dimVk−1=Nk−Nk−1. Sea µk iel conjunto definido por los elementos que generan el espacio N(Dk−1 k). Entonces el error ekse define como ek=Pdk iµk i. Si definimos Gkcomo el operador que a cada elemento de ek∈N(Dk−1 k)asocia un elemento del conjunto de coeficientes dk icorrespondientes a la base µk i; y sea Ek el operador que dada una serie de coeficientes dk iles asocie Pidk iµk i, se establece la equivalencia siguiente: vk≡nvk−1,dko(3.7) donde ahora ambos conjuntos tienen igual cantidad de elementos, pues el número de elementos de fk−1,dkserá igual a dimFk−1+dimN(Dk−1 k) = Nk−1+ (Nk− Nk−1) = Nk=dimV k. Destacar que mediante las siguientes expresiones queda definida la equivalencia entre fkyfk−1: vk−1 = Dk−1 kvk, dk=Gk(I−Pk k−1Dk−1 k)vk, y el paso contrario, mediante vk=Pk k−1vk−1+Ekdk relación extraída directamente a partir de la equivalencia ek=Ekdk. 38 CAPÍTULO 3. TÉCNICAS DE INTERPOLACIÓN NO LINEALES Mediante la equivalencia anterior (3.6) se obtiene la descomposición multiescala de vk. Por ejemplo si consideramos que los datos originales parten de un nivel L de resolución, se tiene: vk≡v0, dl, . . . , d1 vLdL −→ vL−1dL−1 −→ vL−2... −→ Y los correspondientes algoritmos para obtener la transformación multiescala y su paso inverso, son los siguientes: Algoritmo de transformación directa vL→MvL=v0, d1, . . . , dL=   Hacer k=L, . . . , 1 vk−1 = Dk−1 kfk dk=Gk(vk−Pk k−1vk−1) Algoritmo de transformación inversa MvL→M−1MvL=Hacer k=L, . . . , 1 vk=Pk k−1vk−1+Ekdk Llegados a este punto, es obvio que el paso crucial en la construcción de un esquema de multirresolución es la definición de un operador reconstrucción apropiado para la discretización que se esté considerando. De ello dependerán tanto la calidad final de la imagen como el coste computacional total del proceso de ampliación. Habitualmente se utilizan dos tipos de reconstrucción en la multirresolución de Harten, y son la discretización por valores puntuales y la discretización por medias en celda. A continuación se exponen una serie de algoritmos aptos para ser implementados como operador reconstrucción, todos ellos a partir de valores puntuales debido a la natureleza de los datos de entrada (i.e. un conjunto de píxeles). 3.3. Métodos de interpolación no lineales Los algoritmos de zoom que se exponen a continuación han sido definidos para secuencias de datos dos dimensionales. La estrategia llevada a cabo por las diversas técnicas es la de producto tensor la cual se describe a continuación. Sea fun array bidimensional definido como f= (f0 i,j)J0 (i,j)=0 al que denotamos como A, y donde A=A0. La estrategia seguida en las sucesivas técnicas consiste en aplicar el proceso 3.3. MÉTODOS DE INTERPOLACIÓN NO LINEALES 39 de zoom primero sobre las filas y a continuación sobre las columnas, de manera independiente, en contraposición a aquellos algoritmos que actúan directamente de manera bidimensional sobres los datos, esto es, seleccionando unos stencil de más de una dimensión. 40 CAPÍTULO 3. TÉCNICAS DE INTERPOLACIÓN NO LINEALES 3.3.1. Interpolación ENO La interpolación ENO (Essentially Non Oscillatory) tiene como objetivo construir trozos o partes de polinomios usando, en la medida de lo posible, datos pertenecientes a regiones suaves de una función. El punto clave de esta técnica interpolatoria es el proceso por el cual se selecciona el stencil que se intenta elegir dentro de una región suave de una función dada, f(x), esto es, fes infinitamente diferenciable en todos sus órdenes. Este proceso de selección trabaja de la manera siguiente: para cada intervalo Ii= [xj−1, xj], se consideran todos los posibles conjuntos con r >= 2 puntos, incluyendo los puntos xk j−1,xk j. Después de seleccionar el stencil según alguno de los dos métodos de selección que a continuación veremos, el stencil ENO, queda de la siguiente manera: SENO =nxk sj−1, . . . , xk sj+r−1o siendo r+ 1 el orden de interpolación. Para la selección de dicho existen dos estrategias. Ambas producen asintóticamente conjuntos de puntos de interpolación que se mueven lejos de la discontinuidad. Consecuentemente, el orden de aproximación del operador de predicción ENO sigue siendo r+ 1 siempre que sea posible evitar dichas discontinuidades. Algoritmos para la selección del stencil 1. Selección jerárquica Básicamente consiste en, partiendo de los extremos del intervalo, ir añadiendo progresivamente puntos a derecha o izquierda del mismo, comparando las diferencias divididas correspondientes a los conjuntos formados por los extremos del intervalo, más los puntos añadidos, y escogiendo aquélla de menor valor absoluto. 2. Selección no jerárquica Esta selección, por contra, considera las diferencias divididas de mayor orden correspondientes a todos los stencils posibles y calcula el mínimo entre todos los valores absolutos de dichas diferencias. Tanto si empleamos el primero como el segundo, los nodos xi−1,xi+1 pertenecen al stencil SENO. En caso de que ftenga alguna discontinuidad de salto en xd∈Ii, y sean Sun stencil que no cruza dicha discontinuidad y S∗un stencil conteniendo a los nodos xi−1yxi, ambos con s+ 1 nodos. Tenemos entonces: f[S] = O(1); f[S∗]O1 h2 Y si la discontinuidad pertenece a la primera derivada, obtenemos f[S] = O(1); f[S∗] = O1 hs−1 3.3. MÉTODOS DE INTERPOLACIÓN NO LINEALES 41 Algoritmo 1 Selección jerárquica del stencil ENO Entrada: ~ f,r ~ fVector de datos. rOrden de la reconstrucción. Salida: S S-Stencil seleccionado para formar el polinomio interpolador. 1: for i= 1, . . . , J do 2: s0=i 3: for l= 0, . . . , r −2do 4: si |f[xsl−2, . . . , x[sl+l+1] |<|f[xsl−1, . . . , xsl+l+1]|then 5: sl+1 =sl+ 1 6: fin si 7: fin for 8: si=sr−1 9: fin for Algoritmo 2 Selección no jerárquica del stencil ENO Entrada: ~ f,r ~ fVector de datos. rOrden de la reconstrucción. Salida: S S-Stencil seleccionado para formar el polinomio interpolador. 1: for i= 1, . . . , J do 2: elegir sique verifique 3: for l= 0, . . . , r −2do 4: |f[xsi−1, . . . , xsi+r−1]|< min {| f[xl−1, . . . , xl+r−1]|, i −r+ 1 ≤l≤i} 5: fin for 6: si=sr−1 7: fin for 48 CAPÍTULO 3. TÉCNICAS DE INTERPOLACIÓN NO LINEALES Por tanto, una vez hallado ξ, podemos definir la función interpolante de f en el intervalo Iide modo que a la izquierda de ξcoincida con qi−1y a la derecha de ξtome el valor de qi+1, esto es: ISR(x) =    ql(x), x ∈Il, l 6=i qi−1(x), x ∈[xi−1, ξ] qi+1(x), x ∈[ξ, xi−1] Así definido, el error de interpolación es I(x) = f(x) + O(hr+1)excepto en una pequeña región alrededor de xd. El algoritmo quedara entonces como sigue:                                        Hacer k= 1, . . . , L fk 2j=fk−1 j fk 2j−1=                                ENO          5fk−1 j−3+ 15fk−1 j−2+ 5fk−1 j−1+fk−1 j,si S1 j −fk−1 j−2+9fk−1 j−1+9fk−1 j−fk−1 j+1 16 ,si S2 j 5fk−1 j−1+ 15fk−1 j+ 5fk−1 j+1 +fk−1 +2 ,si S3 j ENO −SR    −5fk−1 j−4+21fk−1 j−3−35fk−1 j−2+35fk−1 j−1 16 si qj−1 35fk−1 j−35fk−1 j+1 +21fk−1 j+2 −5fk−1 j+3 16 ,si qj+1 Algoritmo Seguidamente se presenta el pseudocódigo correspondiente a ENO Subcell Resolution, estudiado en la sección anterior. 3.3. MÉTODOS DE INTERPOLACIÓN NO LINEALES 49 Algoritmo 7 Interpolación ENO-SR Entrada: niv,orden,im,met_particion niv - Niveles de zoom. orden - Orden del polinomio interpolatorio. im - Imagen original sobre la cual aplicar el zoom. Salida: b bImagen interpolada. 1: [n m] = size(a) 2: //Bucle para los niveles de zoom. 3: for k= 1 hasta niv do 4: n_filas = 2 ∗n−1 5: n_columnas = 2 ∗m−1 6: b=zeros(n_filas, n_columnas) 7: b(1 : 2 : n_filas, 1 : 2 : n_columnas) = a(1 . . . n, 1. . . m) 8: // Predicción de las columnas. 9: for j= 1 hasta n_columnas a incrementos de 2do 10: b(1 . . . n_filas, j)=enosr_zoom(b(1:2:n_filas, j), n, orden)0 11: fin for 12: // Predicción de las filas. 13: for j= 1 hasta n_filas do 14: b(i, 1. . . n_columnas)=enosr_zoom(b(i, 1:2:n_filas), n, orden) 15: fin for 16: //Actualización de las variables. 17: n=n_filas;m=m_columnas 18: a=b 19: fin for 50 CAPÍTULO 3. TÉCNICAS DE INTERPOLACIÓN NO LINEALES Algoritmo 8 ENOSR Zoom Parte 1/3 Entrada: v,n, ,orden,met_particion ~v - Vector de datos. orden - Orden del polinomio interpolatorio. nLongitud del vector ~v. Salida: f ~ fVector con los valores interpolados. 1: ~ f=zeros(1,2∗n−1) 2: ~ f(1 : 2 : 2 ∗n−1) = ~v(1 . . . n) 3: //Obtenemos las máscaras correspondientes a cada uno de los stencils posibles a elegir. 4: m1=getMask(orden/2, orden/2) 5: m2=getMask(orden/2−1, orden/2 + 1) 6: m3=getMask(orden/2+1, orden/2−1) 7: m4=−5 16 ,21 16,−35 16 ,35 16 8: m5=35 16,−35 16 ,21 16,−5 16  9: //Predicción del primer elemento con máscara lineal. 10: ~ v0=~v[v1. . . vorden] 11: f2=~m2×~ v0 12: //Predicción del segundo elemento. Véase algoritmo ?? 13: s=selJerarquicaEno(~v, orden) 14: si s∈S2then 15: ~ v0=~v[v1. . . vorden] 16: f4=~m1×~ v0 17: else 18: ~ v0=~v[v2. . . vorden+1] 19: f4=~m2×~ v0 20: fin si 21: //Predicción del tercer elemento. Véase algoritmo ?? 22: s=selJerarquicaEno(~v, orden) 23: si S∈S2then 24: ~ v0=~v[v1. . . vorden] 25: f4=~m1×~ v0 26: else 27: ~ v0=~v[v2. . . vorden+1] 28: f4=~m2×~ v0 29: fin si 30: si s∈S1then 31: ~ v0=~v[v1. . . vorden] 32: f6=~m1×~ v0 3.3. MÉTODOS DE INTERPOLACIÓN NO LINEALES 51 Algoritmo 10 ENO-SR Zoom Parte 2/3 33: else si s∈S2then 34: ~ v0=~v[v3. . . vorden+2] 35: f6=~m3×~ v0 36: else 37: ~ v0=~v[v2. . . vorden+1] 38: f6=~m2×~ v0 39: fin si 40: for j= 8 hasta 2∗n−8a incrementos de 2do 41: //Véase algoritmo ?? 42: si=stencil(~v[vj 2−3. . . vj 2+2) 43: sd=stencil(~v[vj 2−1. . . vj 2+4) 44: // Si no están descentrados. 45: si Ssi=0I0and sd=0D0then 46: q11 = vj 2 47: q12 = −vj 2−3+ 4vj 2−2−6vvj 2−1+ 4vj 2 48: q21=4vj 2+1 −6vvj 2+2vj 2+3vj 2+4 49: q22 = vj 2+1 50: g1=q21−q11; g2=q22 −q12 51: si g1∗g2<0then 52: //Existe una discontinuidad en el intervalo. 53: ~v1=~v[vj 2−3. . . vj 2] 54: p1=~v1×~m4 55: ~v2=~v[vj 2+1 . . . vj 2+ 4] 56: p1=~v2×~m5 57: // Y se evalúa g en el punto medio. 58: g=p2−p1; 59: si g1∗g < 0then 60: fj=g2 61: else 62: fj=g1 63: fin si 64: fin si 65: ~ v0=~v[vi 2−1. . . vi 2+2] 66: fj=~m1×~ v0 67: //Si están descentrados, se aplica ENO jerárquico. 68: else 69: //Véase algoritmo 2 70: s=selJerarquicaEno(~v, orden) 52 CAPÍTULO 3. TÉCNICAS DE INTERPOLACIÓN NO LINEALES Algoritmo 12 ENO-SR Zoom Parte 2/3 71: si s∈S1then 72: ~ v0=~v[vj 2−1. . . vj 2+2] 73: fj=~m1×~ v0 74: else si s∈S2then 75: ~ v0=~v[vj 2−2. . . vj 2+1] 76: fj=~m3×~ v0 77: else 78: 79: ~ v0=~v[vj 2. . . vj 2+3] 80: fj=~m2×~ v0 81: fin si 82: fin si 83: fin for 84: // Predicción del antepenúltimo elemento 85: s=selJerarquicaEno(~v, orden) 86: si s∈S1then 87: ~ v0=~v[vn−4. . . vn−1] 88: fj=~m1×~ v0 89: else si s∈S2then 90: ~ v0=~v[vn−5. . . vn−2] 91: fj=~m3×~ v0 92: else 93: 94: ~ v0=~v[vn−3. . . vn] 95: fj=~m2×~ v0 96: fin si 97: // Predicción del penúltimo elemento 98: s=selJerarquicaEno(~v, orden) 99: si s∈S1then 100: ~ v0=~v[vn−3. . . vn] 101: fj=~m1×~ v0 102: else si s∈S2then 103: ~ v0=~v[vn−4. . . vn−1] 104: fj=~m3×~ v0 105: else 106: 107: ~ v0=~v[vn−3. . . vn] 108: fj=~m2×~ v0 109: fin si 110: //Predicción del último elemento 111: ~ v0=~v[vn. . . vn−orden+1] 112: f2n−2=~m3×~ v0 3.3. MÉTODOS DE INTERPOLACIÓN NO LINEALES 53 Algoritmo 13 Funcion Stencil Entrada: ~v ~v - Vector de seis elementos. Salida: s val - Descentramiento del vector ENO Jerárquico. 1: si |(v4−v3)−(v3−v2)|≤|(v5−v4)−(v4−v3)|then 2: si |((v5−v4)−(v4−v3)) −((v4−v3)−(v3−v2))|≤|((v4−v3)−(v3− v2)) −((v3−v2)−(v2−v1))|then 3: val =0C0 4: else 5: val =0I0 6: fin si 7: else 8: si |((v5−v4)−(v4−v3)) −((v4−v3)−(v3−v2))|≤|((v5−v4)−(v3− v2)) −((v5−v4)−(v4−v3))|then 9: val =C 10: else 11: val =0I0 12: fin si 13: fin si 54 CAPÍTULO 3. TÉCNICAS DE INTERPOLACIÓN NO LINEALES 3.3.3. Interpolación WENO La interpolación WENO (Weighed ENO) es otra técnica de reconstrucción no lineal que se presenta como otra mejora respecto a la interpolación ENO. Como se ha tratado anteriormente, ésta última técnica selecciona el stencil más adecuado para intervalo, consiguiendo una aproximación del orden de r+ 1 en aquellos intervalos que se hallan lo suficientemente aisladas. A pesar de ello, existen aspectos de la técnica que se pueden mejorar: En primer lugar, el proceso de selección del stencil es demasiado sensible a las perturbaciones, y un error de redondeo dado entre dos diferencias divididas muy próximas, incurriría en un cambio de selección stencil. En segundo lugar, otro punto a tener en cuenta es que, en aquellas regiones en las que las función es suave, no es necesario hacer esta selección del stencil, ya que la selección utilizada por cualquier método lineal obtendría el mismo conjunto de puntos. Por último, sería posible aumentar el orden de exactitud de la aproximación, ya que tomando stencils de rintervalos, la interpolación ENO se lleva a cabo seleccionando uno de entre rposibles, obteniendo como se ha mencionado, una aproximación del orden de r. Pero dado que existen 2r−1subintervalos contenidos en los rstencils, se pierde información proporcionada por r−1de estos stencils. Si la función es lo suficientemente suave en estas regiones, se podría llegar a obtener un orden de aproximación igual a 2rcomo máximo en estas regiones, utilizando la información dada por los 2r−1stencils. Para solucionar los dos primeros problemas, se presentó una estrategia de poda o sesgo, la cual consiste en tomar como base un stencil centrado en el intervalo donde se realiza la interpolación, y utilizarlo para modificar el criterio de selección del stencil con un parámetro de sesgo. En contraposición con la técnica ENO, que construye el interpolante seleccionando un stencil para cada subintervalo, el método WENO, asigna a cada uno de éstos subintervalos todos los stencils posibles, y el polinomio se calcula como una combinación convexa de los polinomios correspondientes a dichos stencils. Con este tipo de construcción, se prioriza a aquellos polinomios construidos a partir de stencils en donde la función es suave, de forma que aquellos que poseen alguna discontinuidad contribuyen al calculo de forma prácticamente nula. Por tanto, se conserva el efecto ENO (la interpolación en regiones cercanas a singularidades se obtiene mediante información solo de regiones donde la función es suave), y los errores cometidos por unos se pueden cancelar por otros (por la construcción mediante combinación convexa de polinomios), obteniéndose un orden de aproximación mayor. De manera formal, sean Si+k,k= 0,. ..,r−1los rstencils conteniendo al intervalo Ii, y pi+kel polinomio construido a partir de dicho stencil Si+k; el polinomio 3.3. MÉTODOS DE INTERPOLACIÓN NO LINEALES 55 k= 0 k= 1 k= 2 r= 2 1/2 1/2 1/2 r= 3 3/16 10/16 3/16 Tabla 3.1: Pesos óptimos interpolador queda como sigue: Si+k={xi+k−r, . . . , xi+k}, k = 0, . . . , r −1 pWENO i(x) = Pr−1 k=0 wi kpi+k(x) donde wi k≥0, k = 0, . . . , r −1, r−1 X k=0 wi k= 1 (3.11) Como se ha mencionado, el interpolante WENO toma información de 2rnodos con la pretensión de alcanzar una aproximación de este orden en aquellos intervalos donde la función sea suave. Sea ˜ρi2r−1la aproximación empleando 2s nodos {xi−r, . . . , xi+r−1}y sean ˜ρirlas aproximaciones obtenidas con los stencils Si+k. Dicha aproximación se puede expresar como una combinación lineal de las raproximaciones de orden r+ 1, esto es: ρ2r−1 i= r−1 X k=0 Cr kpr i+k(3.12) con Ckr≥0,∀kyPk= 0r−1Cr k= 1. Estas constantes se denominan pesos óptimos que los podemos encontrar en [5]. Pasamos a definir los pesos de modo que verifiquen, por un lado: wi k=Cr k+O(hr−1), k = 0, . . . , r −1(3.13) para imponer que la aproximación obtenida sea de orden 2r; y por otro lado para que la contribución de los polinomios que crucen alguna discontinuidad sea la menor posible. Para poder satisfacer esta última condición, definimos: wi k=αi k Pr−1 s=0 αs , k = 0, . . . , r −1,(3.14) 56 CAPÍTULO 3. TÉCNICAS DE INTERPOLACIÓN NO LINEALES con αi k=Cr k (+ISi+k)2.(3.15) Destacar que ISi+kes un indicador de suavidad de f(x)en el stencil Si+k, y  una constante positiva introducida para evitar la anulación del denominador. Llegados a este punto vemos que los pesos verifican la condición [3.11] independientemente del indicador de suavidad utilizado. Para que la contribución de los polinomios que cruzan alguna singularidad sea prácticamente nula, es suficiente con que ISi+k=O(1) en aquellos stencils que la función presenta alguna discontinuidad. El algoritmo propuesto queda como sigue:    Hacer k= 1 ...,L fk 2j=fk−1 j−1 fk 2j−1=wk−1 j−1sk−1 j−1+wk−1 jsk−1 j+wk−1 j+1 sk−1 j+1 donde      sk−1 j−1= 5fk−1 j−3+ 15fk−1 j−2+ 5fk−1 j−1+fk−1 j sk−1 j=−fk−1 j−2+9fk−1 j−1+9fk−1 j−fk−1 j+1 16 sk−1 j+1 = 5fk−1 j−3+ 15fk−1 j−2+ 5fk−1 j−1+fk−1 j Y los pesos wk−1 j−1,wk−1 jywk−1 j+1 se calculan mediante las expresiones wk−1 j−1=αk−1 j−1 αk−1 j−1+αk−1 j+αk−1 j+1 wk−1 j=αk−1 j αk−1 j−1+αk−1 j+αk−1 j+1 wk−1 j+1 =αk−1 j+1 αk−1 j−1+αk−1 j+αk−1 j+1 Los valores de αk−1 j−1αk−1 jαk−1 j−1son, respectivamente αk−1 j−1= 3 16 +ISk−1 j−1 αk−1 j= 10 16 +ISk−1 j αk−1 j+1 = 3 16 +ISk−1 j+1 siendo una constante positiva introducida para evitar la anulación del denominador, y que típicamente adquiere valores tales como = 10−5o= 10−6. 3.3. MÉTODOS DE INTERPOLACIÓN NO LINEALES 57 Finalmente solo queda definir cómo son calculados los indicadores de suavidad. En [9] se propone calcularlos como sigue: ISk−1 j−1=1 2hfxk−1 j−2, xk−1 j−1−fxk−1 j−3, xk−1 j−22+fxk−1 j−1, xk−1 j−fxk−1 j−2, xk−1 j−12i +fxk−1 j−1, xk−1 j−2fxk−1 j−2, xk−1 j−1+fxk−1 j−3, xk−1 j−22 ISk−1 j=1 2hfxk−1 j−1, xk−1 j−fxk−1 j−2, xk−1 j−12+fxk−1 j, xk−1 j+1 −fxk−1 j−1, xk−1 j2i +fxk−1 j, xk−1 j+1 −2fxk−1 j−1, xk−1 j+fxk−1 j−2, xk−1 j−12 ISk−1 j+1 =1 2hfxk−1 j, xk−1 j+1 −fxk−1 j−1, xk−1 j2+fxk−1 j+1 , xk−1 j+2 −fxk−1 j, xk−1 j+1 2i +fxk−1 j+11, xk−1 j+2 −2fxk−1 j, xk−1 j+1 +fxk−1 j−1, xk−1 j2 Algoritmo A continuación se propone el algoritmo de reconstrucción siguiendo la técnica WENO. 64 CAPÍTULO 3. TÉCNICAS DE INTERPOLACIÓN NO LINEALES Algoritmo 18 Racional Zoom Entrada: v,n, ,orden ~v - Vector de datos. nLongitud del vector ~v. orden - Orden del polinomio interpolatorio. Salida: f ~ fVector con los valores interpolados. 1: ~ f=zeros(1,2∗n−1) 2: ~ f(1 : 2 : 2 ∗n−1) = ~v(1 . . . n) 3: //Obtenemos las máscaras correspondientes a cada uno de los stencils posibles a elegir. 4: m1=getMask(orden/2, orden/2) 5: m2=getMask(orden/2−1, orden/2 + 1) 6: m3=getMask(orden/2+1, orden/2−1) 7: //Predicción del primer elemento con máscara lineal. 8: ~ v0=~v[v1. . . vorden] 9: f2=~m2×~ v0 10: si |f[v1. . . v4]|≤|f[v2. . . v5]|then 11: ~ v0=~v[v1. . . v4] 12: f2=~m1×~ v0 13: else 14: ~ v0=~v[v2. . . v5] 15: f2=~m2×~ v0 16: fin si 17: for i= 6 hasta 2∗n−6a incrementos de 2do 18: //Predicción del primer elemento con máscara lineal. 19: e=vj 2−1−vj 2−2 20: e2=vj 2−vj 2−1 21: ~ v0=~v[vn, vn−1, vn−2, vn−3] 22: f2=~m2×~ v0 23: w0=1+αe2 1 2+αe2 2+e2 1 24: w1=1+αe2 2 2+αe2 2+e2 1 25: fj=w0vj 2+w1vj 2+1 26: fin for 27: //Predicción del último elemento. 3.3. MÉTODOS DE INTERPOLACIÓN NO LINEALES 65 3.3.5. Interpolación PPH En esta sección se describe un esquema de interpolación a trozos polinomial denominada PPH (Piecewise Polynomial Harmonic). Esta técnica de reconstrucción posee varias características deseables como son: Por un lado, cada polinomio esta constituido por un stencil centrado de cuatro puntos. En segundo lugar, en las regiones suaves el la técnica produce el mismo resultado que en al utilizar algún método lineal, como Lagrange. Por último, la exactitud se ve reducida en regiones cerca de las singularidades pero sigue siendo mejor que en el caso lineal. A continuación se describe el operador de recontrstucción PPH de modo análogo al explicado para otras técnicas detalladas anteriormente. Sea IP k(x, fk)el operador de reconstrucción PPH y x∈R, tomemos jtal que x∈[xk j−1, xk j]. Entonces IP k(x, fk) = ˜ Pj(x, fk) = fk j, donde ˜ Pj(x, fk)es un polinomio formado a partir de los datos centrados, {fk j−2, fk j−1, fk j, fk j+1}, y tal que ˜ Pj(xk j−1, fk) = fk j−1y˜ Pj(xk j, fk) = fk j. Se dispone del conjunto de puntos{fk j−2, fk j−1, fk j, fk j+1}y se quiere realizar la predicción del punto medio, fk j−1/2. Como deducimos de lo anteriormente comentado, si la función no contiene singularidades en el intervalo [xk j−2, xk j+1], bastaría con una interpolación centrada para proporcionar una buena aproximación. Pero como hemos visto, cuando la señal muestra singularidades, dicha aproximación pierde exactitud. A continuación se discutirá la modificación propuesta cuando se detecta una singularidad en [xk j, xk j+1]. Supongamos que la diferencia dividida f[xj−1, xj, xj+1]es mayor o igual que f[xk j−2, xk j−1, xk j]en valor absoluto. Esto indica la posible presencia de una singularidad en un punto xd∈[xk j, xk j+1]. Se considera el trozo polinomial para [xk j−1, xk j] escrito como, Pj(x) = a0 + a1(x−xj−1 2) + a2(x−xj−1 2)2+a3(x−xj−1 2)3.(3.28) Para un esquema lineal centrado las cuatro condiciones de interpolación en los puntos xj−2,xj−1,xjyxj+1 son        a0−a13 2h+a2(3 2h)2−a3(3 2h)3=fj−2, a0−a11 2h+a2(1 2h)2−a3(1 2h)3=fj−1, a0−a11 2h+a2(1 2h)2−a3(1 2h)3=fj, a0−a13 2h+a2(3 2h)2−a3(3 2h)3=fj+1. Despejando, obtenemos que a1=fj−2−27fj−1+27fj−fj−1 24h. De este modo, el sistema 66 CAPÍTULO 3. TÉCNICAS DE INTERPOLACIÓN NO LINEALES anterior es equivalente a        a0−a13 2h+a2(3 2h)2−a3(3 2h)3=fj−2, a0−a11 2h+a2(1 2h)2−a3(1 2h)3=fj−1, a0−a11 2h+a2(1 2h)2−a3(1 2h)3=fj, a1=fj−2−27fj−1+27fj−fj−1 24h. Si introducimos los siguientes cambios de variable ej−3 2=f[xj−2, xj−1],ej−1 2= f[xj−1, xj],ej+1 2=f[xj, xj+1],Dj−1=f[xj−2, xj−1, xj]yDj=f[xj−1, xj, xj+1]; después de una serie de manipulaciones llegamos a la expresión: a1=ej−1 2+ 13ej−1 2 12 −1 12 Dj−1+Dj 2h. Se puede observar, que en presencia de una discontinuidad de salto en [xj, xj−1], a1=O(1 h), ya que Dj=O1 h2. Este comportamiento es debido a la mala aproximación de la reconstrucción en presencia de discontinuidades. Por ello, se sustituye la media aritmética por la armónica, y se obtiene la nueva expresión para a1modificada: ˜a1=ej−1 2+ 13ej−1 2 12 −1 12 2Dj−1+Dj 2Dj−1+Dj h. (3.29) La media armónica consigue adaptarse mejor a la presencia de singularidades porque cuando |Dj−1|es O(1) |Dj|es O(1 h2), la media armónica permanece siendo O(1) y, en consecuencia, ˜a1=O(1). Cabe destacar también que en las regiones suaves a1−˜a1=O(h3), ya que la diferencia entre la media armónica y la aritmética original es O(h2). Como resultado la interpolación es de cuarto orden, yfj−1 2−ˆ Pj(xj−1 2) = O(h4)La reconstrucción empleada entonces es de cuarto orden. a1=ej−1 2+ 13ej−1 2 12 . Se tiene entonces a1−˜ ˜a1=O(h)Aunque adaptada a las singularidades, el grado de exactitud se ha reducido a dos. De esta manera, el operador reconstrucción PPH estará constituido por la ecuación (3.29) así como por los nuevos operadores ˜a0,˜a1,˜a2˜a3si El algorimo queda como sigue:            Hacer k= 1 ...,L fk 2j=fk−1 j fk 2j−1=   fk j−1+fk j 2−1 4 Dfk j−1+Dfk j Dfk j−1+Dfk j ,si Dfk j−1Dfk j>0 fk j−1+fk j 2, en otro caso. 3.3. MÉTODOS DE INTERPOLACIÓN NO LINEALES 67 Algoritmo Algoritmo propuesto para la reconstrucción mediante la técnica PPH descrita en la sección anterior. 68 CAPÍTULO 3. TÉCNICAS DE INTERPOLACIÓN NO LINEALES Algoritmo 19 Interpolación PPH Entrada: niv,orden,im,met_particion niv - Niveles de zoom. orden - Orden del polinomio interpolatorio. im - Imagen original sobre la cual aplicar el zoom. Salida: b bImagen interpolada. 1: [n m] = size(a) 2: //Bucle para los niveles de zoom. 3: for k= 1 hasta niv do 4: n_filas = 2 ∗n−1 5: n_columnas = 2 ∗m−1 6: b=zeros(n_filas, n_columnas) 7: b(1 : 2 : n_filas, 1 : 2 : n_columnas) = a(1 . . . n, 1. . . m) 8: // Predicción de las columnas. 9: for j= 1 hasta n_columnas a incrementos de 2do 10: b(1 . . . n_filas, j)=pph_zoom(b(1:2:n_filas, j), n, orden, met_particion)0 11: fin for 12: // Predicción de las filas. 13: for j= 1 hasta n_filas do 14: b(i, 1. . . n_columnas)=pph_zoom(b(i, 1:2:n_filas), n, orden, met_particion) 15: fin for 16: //Actualización de las variables. 17: n=n_filas;m=m_columnas 18: a=b 19: fin for 3.3. MÉTODOS DE INTERPOLACIÓN NO LINEALES 69 Algoritmo 20 PPH Zoom Entrada: v,n, ,orden ~v - Vector de datos. nLongitud del vector ~v. orden - Orden del polinomio interpolatorio. Salida: f ~ fVector con los valores interpolados. 1: ~ f=zeros(1,2∗n−1) 2: ~ f(1 : 2 : 2 ∗n−1) = ~v(1 . . . n) 3: //Obtenemos la máscara correspondientes a cada uno de los stencils posibles a elegir. 4: m2=getMask(orden/2−1, orden/2 + 1) 5: //Predicción del primer elemento con máscara lineal. 6: ~ v0=~v[v1. . . vorden] 7: f2=~m2×~ v0 8: for i= 4 hasta 2∗n−4a incrementos de 2do 9: p1=vj 2−1−2vj 2+vj 2+ 1 10: p2=vj 2−2vj 2+1 +vj 2+ 2 11: aux =vj 2 +vj 2+1 2 12: si p1∗p2>0then 13: fj=aux −(p1∗p2) 4(p1+p2) 14: else 15: fj=aux 16: fin si 17: fin for 18: //Predicción del último elemento. 19: ~ v0=~v[vn, vn−1, vn−2, vn−3] 20: f2=~m2×~ v0 70 CAPÍTULO 3. TÉCNICAS DE INTERPOLACIÓN NO LINEALES Capítulo 4 Experimentos y resultados En este capítulo se pretende evaluar el grado de calidad que es posible obtener mediante la aplicación de los algoritmos anteriormente vistos. Para ello se comenzará con una breve descripción de aquellos indicadores útiles y medibles que son frecuentemente utilizados en el estudio de imágenes digitales. Acto seguido se detallará el procedimiento por el cual éstos serán medidos y evaluados. Por último se presentará las imágenes correspondientes a dos niveles de resolución con el fin de poder apreciar visualmente la calidad obtenida, y las tablas con los resultados del experimento para los restantes niveles. Para finalizar, se muestran una serie de gráficas con las que se pretende obtener conclusiones acerca del experimento. 4.1. Indicadores medibles para la evaluación de la calidad 4.1.1. Introducción Los siguientes indicadores de calidad en imágenes digitales que a continuación se presentan, son frecuentes en el ámbito de la manipulación y el procesamiento las mimas. A continuación se detalla, para cada uno de ellos, una breve descripción y la motivación de su uso para alguno de ellos. Para facilitar la descripción formal se ha considerado que se tienen dos imágenes AyA0definidas de la siguiente manera: A:= (ai,j)M×NyA0:= (a0 i,j)M×N, siendo A0una aproximación a A. 4.1.2. Error Cuadrático Medio El Error Cuadrático Medio (ECM) o Mean Squared Error, (MSE), es una medida del cuadrado del error entre dos imágenes; en el caso que nos ocupa entre 71 72 CAPÍTULO 4. EXPERIMENTOS Y RESULTADOS las imágenes AyA0. Este error mide el grado en que una imagen difiere con respecto a otra, y viene dado por la expresión siguiente: ECM =1 MN M−1 X i=1 N−1 X i=1 |(ai,j)−(a0 i,j)|2(4.1) En el registro de imágenes, es habitual hacer uso de esta medida para cuantificar el error que se produce entre el vector de movimiento real de la imagen objetivo y de la de referencia, y el vector de movimiento estimado en el procedimiento de registro. La presencia de este error se debe a que aveces la estimación no se calcula con una precisión suficiente. Otra medida directamente derivada del ECM es la raíz cuadrática media, o Root Mean Square, (RMS), calculada como la raíz cuadrada del ECM, de la siguiente manera: RMS =√ECM =v u u t M−1 X i=1 N−1 X j=1 |(ai,j)−(a0 i,j)|2(4.2) 4.1.3. Relación señal a ruido de pico (PSNR) La relación PSNR (Peak to Signal Noise Ratio) es utilizada para definir la relación entre la máxima energía posible de una señal y el ruido que afecta a su representación fidedigna. Debido a que muchas señales tienen un gran rango dinámico, el PSNR se expresa generalmente en escala logarítmica, utilizando como unidad el decibelio. Cabe recordar que un aumento de 20 dB corresponde a un decrecimiento de una décima parte en la diferencia RMS entre las dos imágenes. El uso más habitual del PSNR es como medida cuantitativa de la calidad de la reconstrucción en el ámbito de la compresión de imágenes. Esta medida se define como: PSNR = 20log b √ECM .(4.3) donde bes el mayor valor posible de la señal, y RMS es la raíz cuadrática media. Para una imagen en formato RGB, la definición del PSNR es la misma, pero el ECM se calcula como la media aritmética de los ECM de los tres colores (R, G y B). Los valores típicos que adopta este parámetro están entre 30 y 50 dB, siendo mayor cuanto mejor es la codificación. El comité MPEG emplea un valor umbral informal de 0,5 dB en el incremento del PSNR para decidir si se incluye una determinada mejora en un algoritmo de codificación, ya que se considera que este aumento del PSNR es apreciable visualmente. 4.1. INDICADORES MEDIBLES PARA LA EVALUACIÓN DE LA CALIDAD73 4.1.4. Correlación cruzada normalizada El Coeficiente de Correlación Cruzada o Normalized Cross Correlation, (NK), es una de las medidas de similitud más frecuentemente utilizada. Ésta se calcula entre parejas de bloques pertenecientes a la imagen de referencia y a la imagen objetivo, con el propósito de encontrar el máximo entre dicha medida. Aquél bloque con el que se consigue el máximo es el que determina la correspondencia finalmente establecida. Este coeficiente permite el alineamiento con precisión de imágenes que han sido trasladadas entre sí, aunque también es posible su aplicación entre imágenes que han sufrido rotaciones leves o escalados. Como desventajas, se suelen citar entre otras, el elevado coste computacional que requiere, aunque es inferior al de otras medidas frecuentes; y la excesiva planicidad de los máximos de similitud detectados, debido a la autosimilitud de las imágenes. NK =PM i=1 PN j=1(ai,j )·(a0 i,j ) PM i=1 PN j=1(ai,j )2. 4.1.5. Diferencia media La Diferencia Media o Average Difference, (AD) no es más que la media de las diferencias obtenidas entre píxeles de ambas imágenes. Formalmente AD =PM i=1 PN j=1 (ai,j )−(a0 i,j ) MN . 4.1.6. Diferencia Máxima La Diferencia Máxima o Maximum Difference, (MD) es la máxima diferencia existente entre píxeles de ambas imágenes. MD =Max |(ai,j)−(a0 i,j)|. 4.1.7. Error Absoluto Medio El Error Normalizado Medio (ENM), o Normalized Absolute Error, (NAE) de define como la media aritmética de los errores absolutos cometidos. NAE =PM i=1 PN j=1|(ai,j )−(a0 i,j )| PM i=1 PN j=1|(ai,j )|. 80 CAPÍTULO 4. EXPERIMENTOS Y RESULTADOS En la columna de la izquierda de las siguientes tablas se encuentra el nombre de las técnicas estudiadas en el tercer capítulo de esta memoria. A modo de recordatorio. ENO jerárquico óENOh Interpolación ENO mediante selección jerárquica del stencil utilizando cuatro puntos (3.3.1). ENO no jerárquico óENO Interpolación ENO mediante selección jerárquica del stencil utilizando cuatro puntos (3.3.1). WENO Interpolación WENO stencil utilizando cuatro puntos (3.3.3). ENO-SR Interpolación ENO-SR basada en la técnica ENO con selección jerárquica del stencil (3.3.2). RACIONAL Interpolación utilizando cuatro puntos en la selección del stencil y un α= 0,5(3.3.4). PPH Interpolación PPH (3.3.5). Aquellas siglas o abreviaturas utilizadas en las columnas de cada tabla corresponden a los indicadores de calidad de imagen mencionados con anterioridad, y son los que a continuación se enumeran. MSE Error Cuadrático Medio o Mean Squared Error (4.1.2). PSNR Relación Señal a Ruido de Pico, o Peak Signal to Noise Ratio (4.1.3). NCC Coeficiente de Correlación Normalizado, o Normalized Cross Correlation (4.1.4). AD Diferencia Media, oAverage Difference (4.1.5). MD Diferencia Máxima, o Maximum Difference (4.1.6). NAE Error Absoluto Normalizado, o Normalized Absolute Error (4.1.7). En la sección ?? se exponen las conclusiones acerca de los resultados observados. 4.3. RESULTADOS 81 4.3.1. Resultados con la imagen Dicom1 Técnica N MSE PSNR NCC AD MD NAE ENO no jerárquico 1 120,6907 27,3141 0,9145 0,2300 102 0,3874 2 141,1701 26,6334 0,8866 0,4619 128 0,4161 3 204,7422 25,0186 0,8501 0,7106 155 0,4746 4 275,9669 23,7222 0,7376 1,8971 182 0,5529 ENO jerárquico 1 134,8196 26,8310 0,9115 0,1981 102 0,4132 2 154,7419 26,2347 0,8911 0,3765 128 0,4337 3 217,0438 24,7653 0,8533 0,7455 149 0,4887 4 295,4249 23,4249 0,7386 1,7789 180 0,5662 WENO 1 136,6975 26,7732 0,8979 0,4118 121 0,4050 2 153,0548 26,2823 0,8476 0,8541 167 0,459 3 208,0133 24,9499 0,7485 2,0822 191 0,4722 4 353,2384 22,6624 0,5326 4,4272 229 0,5935 ENO-SR 1 136,9803 26,7642 0,9057 0,3680 105 0,4164 2 157,5326 26,1571 0,8824 0,5962 128 0,4372 3 219,9318 24,7079 0,8463 0,8533 151 0,4922 4 307,1795 23,2569 0,7319 1,661 181 0,5840 RACIONAL 1 115,6464 27,3316 0,9020 0,5802 104 0,3818 2 139,1557 26,6859 0,8586 1,0636 143 0,4034 3 193,8547 25,2565 0,8034 1,5401 154 0,491 4 318,4464 23,1004 0,6780 2,6364 189 0,5843 PPH 1 115,6464 27,4976 0,910 0,0660 101 0,3811 2 138,1151 26,7284 0,8937 0,2510 123 0,4116 3 196,3408 25,2007 0,8611 0,5099 151 0,4675 4 291,5029 23,4844 0,7389 1,7228 178 0,5635 Tabla 4.1: Resultados numéricos para la imagen Dicom1 mediante el uso de las diferentes técnicas con diferentes factores de zoom. 82 CAPÍTULO 4. EXPERIMENTOS Y RESULTADOS (a) (b) (c) (d) (e) (f) Figura 4.3: Imágenes Dicom1 con un nivel de zoom (x2) acotadas por una ventana de 175 ×151 píxeles. Reconstrucción con ENO no jerárquico (a), con ENO jerárquico (b), con ENO-SR (c), con WENO (d), con Racional (e), y con la técnica PPH (f). 4.3. RESULTADOS 83 (a) (b) (c) (d) (e) (f) Figura 4.4: Imágenes Dicom1 con un nivel de 2(x4) acotadas por una ventana de 175×151 píxeles. Reconstrucción con ENO no jerárquico (a), con ENO jérárquico (b), con ENO-SR (c), con WENO (d), con Racional (e), y con la técnica PPH (f). 84 CAPÍTULO 4. EXPERIMENTOS Y RESULTADOS 4.3.2. Resultados con la imagen Dicom2 Técnica N MSE PSNR NCC AD MD NAE ENO no jerárquico 1 113,1933 27,5926 0,9304 0,0602 117 0,3693 2 148,6390 26,4095 0,9043 0,0103 140 0,4133 3 217,7340 24,7515 0,8450 0,1768 171 0,4803 4 335,6572 22,8718 0,7044 1,1989 179 0,6042 ENO jerárquico 1 126,0167 27,1265 0,9283 0,00217 117 0,3942 2 158,3730 26,1340 0,9108 −0,0824 140 0,4286 3 231,6568 24,4824 0,8547 0,0439 187 0,4926 4 350,8093 22,6801 0,7155 0,9933 177 0,6179 WENO 1 130,4950 26,9749 0,9106 0,2102 141 0,3897 2 167,1330 25,9002 0,8567 0,3394 172 0,4181 3 240,1387 24,3262 0,7464 1,3369 176 0,4849 4 382,2118 22,3078 0,5660 2,6737 195 0,6401 ENO-SR 1 128,007 27,058 0,9227 0,1791 117 0,3977 2 161,3301 26,0536 0,9048 0,1283 140 0,4317 3 235,4351 24,4121 0,8472 0,2164 184 0,4973 4 361,3813 22,551 0,7048 1,0844 177 0,6321 RACIONAL 1 112,9673 27,6013 0,9187 0,3563 117 0,3651 2 147,5731 26,4407 0,8840 0,5129 155 0,4019 3 207,6204 24,9581 0,7930 1,0537 171 0,430 4 335,5185 22,8736 0,6203 2,4094 185 0,6066 PPH 1 108,5177 27,7758 −0,1370 −0,1370 117 0,3638 2 144,9833 26,5176 0,9135 −0,2287 144 0,4091 3 220,6592 24,6936 0,8538 −0,0245 174 0,4784 4 352,1596 22,6634 0,7211 0,6919 179 0,6186 Tabla 4.2: Resultados numéricos para la imagen Dicom2 mediante el uso de las diferentes técnicas con diferentes factores de zoom. 4.3. RESULTADOS 85 (a) (b) (c) (d) (e) (f) Figura 4.5: Imágenes Dicom2 con un nivel de zoom (x2) acotadas por una ventana de 174 ×151 píxeles. Reconstrucción con ENO no jerárquico (a), con ENO jérárquico (b), con ENO-SR (c), con WENO (d), con Racional (e), y con la técnica PPH (f). 86 CAPÍTULO 4. EXPERIMENTOS Y RESULTADOS (a) (b) (c) (d) (e) (f) Figura 4.6: Imágenes Dicom2 con un nivel 2 (x4) de zoom acotadas por una ventana de 174 ×151 píxeles. Reconstrucción con ENO no jerárquico (a), con ENO jérárquico (b), con ENO-SR (c), con WENO (d), con Racional (e), y con la técnica PPH (f). 4.3. RESULTADOS 87 4.3.3. Resultados con la imagen Lena Técnica N MSE PSNR NCC AD MD NAE ENO no jerárquico 1 61,7359 30,2254 0,9954 0,3050 189 0,0342 2 283,7605 23,6013 0,9836 0,5736 189 0,0349 3 712,4067 19,6035 0,9602 1,1324 216,076 0,1493 41410 16,6035 0,9228 1,7535 214,9427 0,2349 ENO jerárquico 1 60,9001 30,2846 0,9956 6,2969 189 0,0391 2287,0612 23,5511 0,9842 0,5570 189 0,0780 3744,4507 19,4124 0,9615 1,0174 216,0762 0,1535 41501 16,3656 0,9252 1,4454 216,5927 0,2422 WENO 1 157,5510 26,1565 1,0077 −0,2058 189 0,0636 2 455,4448 21,5464 1,0132 2,9024 202 0,1189 3 955,0030 18,3308 1,0143 5,3136 193 0,1850 4 1770 15,7553 0,9974 −7,5572 191 0,2714 ENO-SR 1 61,1501 30,2668 1,0077 −1,2058 189 0,0636 2 289,1822 23,5191 0,9839 0,5939 189 0,0882 3 744,3638 19,4130 0,9608 1,1265 196 0,1544 4 1513,7 16,3305 0,9246 1,4993 198 0,2456 RACIONAL 1 71,0149 29,6173 0,9968 0,2633 189 0,0353 2 330,3484 22,9411 0,9848 0,7057 194 0,0872 3 788,6687 19,1619 0,9634 1,1424 195 0,1484 4 1410 16,3871 0,9308 1,5220 191 0,273 PPH 1 61,2942 30,2566 0,9961 0,2147 189 0,0341 2 280,2441 23,6013 0,9840 0,4817 189 0,0863 3 704,4169 19,6525 0,9607 0,9674 195 0,1491 4 1397,2 16,678 0,9231 1,5848 191 0,2338 Tabla 4.3: Resultados numéricos para la imagen Lena.jpg mediante el uso de las diferentes técnicas con diferentes factores de zoom. 88 CAPÍTULO 4. EXPERIMENTOS Y RESULTADOS (a) (b) (c) (d) (e) (f) Figura 4.7: Imágenes Lena.jpg con un nivel de zoom (x2) acotadas por una ventana de 256 ×256 píxeles. Reconstrucción con ENO no jerárquico (a), con ENO jérárquico (b), con ENO-SR (c), con WENO (d), con Racional (e), y con la técnica PPH (f). 4.3. RESULTADOS 89 (a) (b) (c) (d) (e) (f) Figura 4.8: Imágenes Lena.jpg con un factor de zoom 2 (x4) acotadas por una ventana de 255 ×255 píxeles. Reconstrucción con ENO no jerárquico (a), con ENO jérárquico (b), con ENO-SR (c), con WENO (d), con Racional (e), y con la técnica PPH (f). 96 CAPÍTULO 4. EXPERIMENTOS Y RESULTADOS 4.3.6. Resultados con la imagen Tac2 Técnica N MSE PSNR NCC AD MD NAE ENO no jerárquico 1 39,3236 32,1843 0,9977 00,88 183 0,0455 2 255,7191 24,0532 0,9854 0,0401 195 0,1146 3 844,9424 24,0532 0,9490 0,2619 243 0,2150 4 1717 15,7816 0,8931 0,5474 249 0,3436 ENO jerárquico 1 38,3645 32,2915 0,9979 0,0047 183 0,0,456 2 251,9034 24,1185 0,4864 0,0336 195 0,1140 3 852,4802 18,8240 0,9504 0,3015 243 0,2190 4 1830 15,4913 0,8940 0,4101 255 0,354 WENO 1 154,4015 26,2443 1,602 −0,3554 221 0,0831 2 528,6831 20,8998 0,9914 −0,3181 246 0,1593 3 1250 17,1587 0,9648 −0,3181 246 0,1593 4 2870 14,9356 0,9302 −2,1281 225 0,3710 ENO-SR 1 38,5038 32,2758 0,9978 0,0108 183 0,0457 2 252,2411 26,1122 0,9864 0,0355 195 0,1142 3 654,9497 20,8114 0,95 0,3769 210 0,2198 4 1840 16,4653 0,8914 6,6097 250 0,3563 RACIONAL 1 47,8495 31,3324 0,9993 0,0384 185 0,0473 2331,3535 22,9279 0,9892 0,2011 238 0,1191 3 1043 17,9757 1,3236 0,9494 243 0,2223 4 1844 15,3568 0,9080 1,7185 247 0,3277 PPH 1 38,7738 32,2954 0,9981 −0,0609 183 0,0453 2 253,8404 26,0852 0,9854 −0,0182 195 0,1144 3 838,55 20,8955 0,9493 0,0624 243 0,2193 4 1748,6 16,7164 0,8406 0,528 249 0,3493 Tabla 4.6: Resultados numéricos para la imagen Tac2.jpg mediante el uso de las diferentes técnicas con diferentes factores de zoom. 4.3. RESULTADOS 97 (a) (b) (c) (d) (e) (f) Figura 4.13: Imágenes Tac2.jpg con un 1 nivel de zoom (x2) acotadas por una ventana de 512 ×512 píxeles. Reconstrucción con ENO no jerárquico (a), con ENO jérárquico (b), con ENO-SR (c), con WENO (d), con Racional (e), y con la técnica PPH (f). 98 CAPÍTULO 4. EXPERIMENTOS Y RESULTADOS (a) (b) (c) (d) (e) (f) Figura 4.14: Imágenes Tac2.jpg con un factor de zoom 2 (x4) acotadas por una ventana de 512×512 píxeles. Reconstrucción con ENO no jerárquico (a), con ENO jérárquico (b), con ENO-SR (c), con WENO (d), con Racional (e), y con la técnica PPH (f). 4.3. RESULTADOS 99 4.3.7. Resultados con la imagen Pet-Tc1 Técnica N MSE PSNR NCC AD MD NAE ENO no jerárquico 1 44,9343 31,6050 0,9939 0,0013 142 0,0783 2 224,1021 24,6263 0,9683 0,01184 238 0,1810 3 4740,2337 21,3709 0,9330 −0,0259 238 0,2713 4 781,2310 19,2030 0,87779 0,3312 245 0,3544 ENO jerárquico 1 44,3473 31,6221 0,9943 −0,000125 ∗10−5150 0,0780 2 244,3708 24,2503 0,9681 −0,0080 238 0,1882 3 518,0333 20,9872 0,9311 0,03377 248 0,845 4 837,5324 18,9008 0,881 0,2912 246 0,3726 WENO 1 174,4752 25,7135 0,9927 −0,5268 237 0,1475 2 368,7084 22,4640 0,9804 −0,09805 243 0,2286 3 587,9350 20,4375 0,9486 −,9213 249 0,2977 4 894,1853 18,6165 0,8838 −0,0439 252 0,3805 ENO-SR 1 44,5924 31,6382 0,9945 −0,80 150 0,0782 2 247,1329 24,2019 0,9673 0,0394 239 0,1894 3 524,3784 20,9344 0,9282 0,1813 247 0,3758 4 848,7539 18,8430 0,8746 0,5349 247 0,3758 RACIONAL 1 56,4816 30,6117 0,9961 0,0441 203 0,0836 2 257,9580 24,0153 0,9706 0,1312 239 0,1872 3 514,5613 21,0164 0,9343 0,2211 244 0,0705 4 831,2260 18,9336 0,8525 1,9148 255 0,355 PPH 1 44,1232 31,6841 0,9940 −0,0247 141 0,0779 2 221,1880 24,6832 0,9679 −0,0142 238 0,1803 3 416,7993 21,4209 0,9303 −0,0086 238 0,2699 4 772,9352 19,2499 0,8804 0,0891 244 0,3585 Tabla 4.7: Resultados numéricos para la imagen Pet-Tc1.jpg mediante el uso de las diferentes técnicas con diferentes factores de zoom. 100 CAPÍTULO 4. EXPERIMENTOS Y RESULTADOS (a) (b) (c) (d) (e) (f) Figura 4.15: Imágenes Pet1.jpg con un nivel de zoom (x2) acotadas por una ventana de 961 ×512 píxeles. Reconstrucción con ENO no jerárquico (a), con ENO jérárquico (b), con ENO-SR (c), con WENO (d), con Racional (e), y con la técnica PPH (f). 4.3. RESULTADOS 101 (a) (b) (c) (d) (e) (f) Figura 4.16: Imágenes Pet-Tc1.jpg con un factor de zoom 2 (x4) acotadas por una ventana de 961 ×512 píxeles. Reconstrucción con ENO no jerárquico (a), con ENO jérárquico (b), con ENO-SR (c), con WENO (d), con Racional (e), y con la técnica PPH (f). 102 CAPÍTULO 4. EXPERIMENTOS Y RESULTADOS 4.3.8. Resultados con la imagen Pet-Tc2 Técnica N MSE PSNR NCC AD MD NAE ENO no jerárquico 1 32,1185 38,0633 0,9960 0,0059 109 0,0393 2 170,4518 25,8141 0,9770 0,2311 214 0,0988 3 398,0473 25,8148 0,9508 0,3511 245 0,1837 4 983,8310 18,2016 0,8582 1,5086 250 0,3245 ENO jerárquico 1 27,8125 33,6884 0,9970 0,0285 101 0,0384 2 159,330 26,1077 6,9801 0,1955 214 0,0973 3 400,5415 22,1043 0,9531 0,3151 245 0,1856 4 962,0801 18,2987 0,8732 1,1633 250 0,3242 WENO 183,3407 28,6677 1,0072 −0,4396 180 0,0677 2 206,8424 27,9744 0,9837 0,2481 214 0,0998 3 618,6745 20,2162 0,9444 0,7119 250 0,1934 4 1345 16,8409 0,8633 1,5738 250 0,354 ENO-SR 1 27,9109 33,6731 ,09971 0,0255 101 0,0386 2 159,9576 26,0908 0,9799 0,2067 214 0,0978 3 405,7442 22,0483 0,9519 0,3626 245 0,1870 4 1006,2 18,1040 0,8666 1,2478 250 0,3304 RACIONAL 188,3704 28,6677 1,0072 −0,4396 180 0,0677 2 206,8428 24,9744 0,9837 0,2481 214 0,0998 3 618,6745 20,2162 0,9444 0,7119 250 0,1934 4 1345,8 16,8409 0,8633 1,57,38 250 0,3254 PPH 130,0271 33,3557 0,9971 −0,0182 85 0,0390 2 169,2704 25,84,50 ,09789 ,1069 214 0,0988 3 419,3260 21,9053 0,9509 0,0496 245 0,1877 4 1001,2 18,1255 0,8647 0,8665 250 0,3296 Tabla 4.8: Resultados numéricos para la imagen Pet-Tc2.jpg mediante el uso de las diferentes técnicas con diferentes factores de zoom. 4.3. RESULTADOS 103 (a) (b) (c) (d) (e) (f) Figura 4.17: Imágenes Pet-Ct2.jpg con un nivel de zoom (x2) acotadas por una ventana de 516×542 píxeles. Reconstrucción con ENO no jerárquico (a), con ENO jérárquico (b), con ENO-SR (c), con WENO (d), con Racional (e), y con la técnica PPH (f). 104 CAPÍTULO 4. EXPERIMENTOS Y RESULTADOS (a) (b) (c) (d) (e) (f) Figura 4.18: Imágenes Pet-Ct2.jpg con un factor de zoom 2 (x4). Reconstrucción con ENO no jerárquico (a), con ENO jérárquico (b), con ENO-SR (c), con WENO (d), con Racional (e), y con la técnica PPH (f). 4.3. RESULTADOS 105 4.3.9. Resultados con la imagen Pet1 Técnica N MSE PSNR NCC AD MD NAE ENO no jerárquico 1 251,8877 24,1187 0,9306 0,0585 255 0,2444 2 740,6606 19,4346 0,7822 0,4978 255 0,5137 3 120,4 17,3305 0,5941 0,5341 255 1,5484 4 2072,6 14,9636 0,3961 0,0265 255 0,1320 ENO jerárquico 1 239,6264 24,3355 0,9318 0,1152 255 0,2460 2 734,4960 19,7409 0,7874 0,4650 255 0,5162 3 1223,9 17,2533 0,5916 1,5453 255 0,7568 4 2124,7 14,8578 0,4020 −0,4020 255 1,1416 WENO 1 451,3236 21,5859 0,8571 1,3639 255 0,3543 2 937,8203 18,4096 0,6013 4,4571 255 0,5726 3 1390,6 15,2723 0,3690 8,5724 255 1,0293 4 2100 05,2723 0,2100 8,1692 255 1,0293 ENO-SR 1 240,1396 24,3262 0,9313 0,416 255 0,2469 2 735,6175 19,4643 0,7859 0,5404 255 0,7669 3 1225,6 17,1422 0,5877 1,5142 255 0,7669 4 2106,7 24,8149 0,3991 0,3621 255 1,1311 RACIONAL 1 248,0244 23,5972 0,9220 0,6031 255 0,2482 2 769,0261 19,2714 0,7611 2,1585 255 0,4986 3 1236,7 17,2082 0,4864 6,1446 255 06910 4 1924,4 15,2879 0,2692 7,7486 255 0,9739 PPH 1 237,5153 24,3262 0,9313 0,1416 255 0,2469 2 735,6175 19,4643 0,7822 0,1978 255 0,2137 3 1255,6 17,1422 0,5941 1,5484 255 0,7472 4 2106,7 14,8949 0,3991 0,3621 255 1,1311 Tabla 4.9: Resultados numéricos para la imagen Pet1.jpg mediante el uso de las diferentes técnicas con diferentes factores de zoom. 112 CAPÍTULO 4. EXPERIMENTOS Y RESULTADOS Figura 4.25: Gráfico que muestra los valores del PSNR obtenidos para cada técnica y para cada nivel, aplicados sobre la imagen Lena.jpg. Figura 4.26: Gráfico que muestra los valores del PSNR obtenidos para cada técnica y para cada nivel, aplicados sobre la imagen Geo.jpg. 4.3. RESULTADOS 113 Figura 4.27: Gráfico que muestra los valores del PSNR obtenidos para cada técnica y para cada nivel, aplicados sobre la imagen Tac1.jpg. Figura 4.28: Gráfico que muestra los valores del PSNR obtenidos para cada técnica y para cada nivel, aplicados sobre la imagen Tac2.jpg. 114 CAPÍTULO 4. EXPERIMENTOS Y RESULTADOS Figura 4.29: Gráfico que muestra los valores del PSNR obtenidos para cada técnica y para cada nivel, aplicados sobre la imagen Pet-Tc1.jpg. Figura 4.30: Gráfico que muestra los valores del PSNR obtenidos para cada técnica y para cada nivel, aplicados sobre la imagen Pet-Tc2.jpg. 4.3. RESULTADOS 115 Figura 4.31: Gráfico que muestra los valores del PSNR obtenidos para cada técnica y para cada nivel, aplicados sobre la imagen Pet1.jpg. Figura 4.32: Gráfico que muestra los valores del PSNR obtenidos para cada técnica y para cada nivel, aplicados sobre la imagen Pet2.jpg. 116 CAPÍTULO 4. EXPERIMENTOS Y RESULTADOS Figura 4.33: Gráficos de medias obtenidas para cada técnica en cada nivel. Figura 4.34: Gráficos de medias obtenidas para cada imagen en cada nivel. 4.3. RESULTADOS 117 Niv. ENO ENOh WENO ENOSR RAC PPH Dicom1 1 27,3141 26,8310 26,7732 26,7642 27,3316 27,4976 2 26,6334 26,2647 26,2823 26,1571 26,6959 26,7284 3 25,0186 24,7653 24,9499 24,7079 25,2565 25,2007 423,7222 23,4349 22,6624 23,2569 23,1004 23,4844 Dicom2 1 27,5926 27,1265 26,9749 27,0586 27,6013 27,7758 2 26,4095 26,1340 25,9002 26,0536 26,7407 25,5176 3 24,7515 24,4824 24,3262 24,4121 24,9581 24,6936 4 22,8718 22,6801 22,6634 22,5511 22,8736 22,6634 Lena 1 30,2254 30,2846 26,1565 30,2668 29,6117 30,2566 223,6013 23,5511 21,5464 23,5191 22,9411 23,6013 3 19,6035 19,4124 18,3308 19,4130 19,1619 19,6525 4 16,6035 16,3656 15,7553 16,3305 16,3871 16,6780 Geo 1 22,9247 22,2297 21,0658 22,1488 22,5190 23,1453 2 19,9252 19,4448 19,0043 19,3654 19,7008 20,0000 3 18,0696 17,6552 17,4473 17,5499 17,8960 18,0955 4 16,3071 15,9710 15,9357 15,8952 16,1294 16,4225 Tac1 1 28,125 28,6793 23,972 24,852 26,7670 28,1518 2 20,8502 21,1671 18,9375 19,5344 19,7666 20,9321 317,1625 17,1519 16,099 17,006 16,5243 17,0532 4 14,3698 14,4804 14,1893 14,0192 14,3547 14,2668 Tac2 1 32,1843 32,2915 26,2553 32,2758 31,3324 32,2454 2 24,0532 24,1185 20,8988 26,1122 22,9274 26,0852 3 18,8625 18,824 17,1587 20,8114 17,9754 20,8955 4 15,7816 15,4973 14,9356 16,4653 15,3568 16,7164 Pet-Tc1 131,6050 31,3221 25,7135 31,6382 30,6117 31,6841 2 24,6263 24,2503 22,4640 24,2019 24,0153 24,6832 3 21,3709 20,9872 20,4175 21,4209 21,0164 21,4209 4 19,2030 18,9008 18,6165 19,2494 18,9336 19,2494 Pet-Tc2 1 33,0633 33,6884 28,6677 33,9731 28,6677 33,3557 2 25,8148 26,1077 23,1062 26,0908 24,9744 25,845 322,1315 22,1043 20,0143 22,0483 20,2162 21,9053 4 18,2016 18,2987 18,1040 18,1040 16,8409 18,1255 Pet1 1 24,1187 24,3355 21,5859 24,3262 23,5972 24,9739 219,4946 19,4709 18,4096 19,4643 19,4709 19,3097 317,3303 17,2533 16,6807 17,1422 17,2533 17,244 4 14,9656 14,8575 15,2723 14,8149 14,8578 14,8117 Pet2 1 25,6257 25,6318 24,175 25,6304 25,5191 25,1684 2 21,6845 21,7042 20,4208 21,697 21,6110 21,6991 3 17,4864 17,5117 17,2709 17,4925 17,4562 17,437 4 13,7884 13,7519 14,8692 13,7267 13,8920 13,7129 Tabla 4.11: Tabla que resume los valores de PSNR obtenido para las distintas imágenes agrupados por tipo de técnica. 118 CAPÍTULO 4. EXPERIMENTOS Y RESULTADOS Tabla de medias y desviaciones por técnica Nivel ¯xi, σi2ENO ENOh WENO ENOSR RACIONAL PPH 1¯x125,7361 25,6925 23,2079 25,4120 25,0108 25,8169 σ2 15,9184 6,0309 4,5562 5,8133 5,4680 5,8018 2¯x222,3857 22,3043 21,0075 22,1795 22,0134 22,4327 σ2 23,8146 3,8441 3,4851 3,7390 3,8045 3,7456 3¯x320,3125 20,1801 19,4280 20,4305 19,5424 20,4060 σ2 33,4152 3,4402 3,2180 3,5103 3,6781 3,3880 4¯x418,5086 18,3912 18,0670 18,4025 18,6184 18,7667 σ2 418,5086 18,3912 18,0670 18,4025 18,6184 18,7667 Tabla 4.12: Media y desviación estándar de los valores del PSNR de todas las imágenes obtenidos para cada técnica y nivel de zoom. Tabla de medias y desviaciones por imagen (I) Nivel ¯xi, σi2Dicom1 Dicom2 Lena Geo Tac1 Tac2 1¯x127,0853 27,3550 29,4669 22,3389 26,7579 31,0975 σ2 10,3311 0,3402 1,6424 0,7341 1,9441 2,4006 2¯x226,2188 26,0759 23,1267 19,5734 20,1980 24,0326 σ2 20,7679 0,3437 0,8144 0,3757 0,9076 1,9793 3¯x325,2247 24,6040 19,2624 17,7856 16,8329 19,0879 σ2 30,4899 0,2383 0,4882 0,2742 0,4296 1,5045 4¯x423,2769 22,7172 16,3533 16,1102 14,2800 15,7922 σ2 40,3676 0,1289 0,3248 0,2157 0,1614 0,6808 Tabla de medias y desviaciones por imagen (II) Nivel ¯xi, σi2Pet-Tc1 Pet-Tc2 Pet1 Pet2 1¯x130,4291 31,9027 23,8229 25,2917 σ2 12,3447 2,5244 1,1818 0,5755 2¯x224,0402 25,3232 19,2700 21,4694 σ2 20,8139 1,1624 0,4268 0,5149 3¯x321,1056 21,4033 17,1506 17,4425 σ2 30,3911 1,0028 0,2379 0,0882 4¯x419,0255 17,9458 14,9300 13,9569 σ2 40,2542 0,5465 0,1767 0,4515 Tabla 4.13: Media y desviación estándar de los valores del PSNR obtenidos para todas las técnicas agrupadas por tipo de imagen y nivel de zoom. 4.3. RESULTADOS 119 Figura 4.35: Histograma que resume el porcentaje de veces que una técnica consigue mejores resultados que el resto para una configuración (<Imagen><Técnica><Nivel>). Cada barra aparece de un color diferente y atiende a la división en grupos establecida en base a la calidad obtenida. Aparecen en verde la técnica incluida en el grupo de alta calidad, en ámbar las clasificadas en el grupo de calidad media y en aquéllas del denominado grupo de baja calidad. Figura 4.36: Gráfico de dispersión obtenido en función del PSNR obtenidos para cada tipo de imagen. 120 CAPÍTULO 4. EXPERIMENTOS Y RESULTADOS Figura 4.37: Gráfico de dispersión obtenido en función del PSNR obtenidos para cada nivel de zoom. Figura 4.38: Gráfico de dispersión obtenido en función del PSNR obtenidos para cada técnica interpolatoria estudiada. 4.3. RESULTADOS 121 4.3.11. Análisis de los resultados obtenidos En esta sección se propone emitir una valoración acerca de los resultados obtenidos en la sección anterior. El análisis parte de las tablas de resultados de PSNR obtenidas por cada técnica y para cada imagen (tablas 4.1, 4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8. 4.9 y la 4.10)), obtenidas en los anteriores experimentos. En todas ellas han sido medidos los indicadores anteriormente comentados en función de la técnica y nivel utilizado. De ahora en adelante denominaremos configuración a la terna <Imagen-Técnica-Nivel>, para referirnos a la técnica empleada para ampliar una imagen con un nivel determinado de zoom. En las tablas de datos de cada imagen, se verifica el incremento del valor del PSNR, del Error Cuadrático Medio y el Error Absoluto Normalizado conforme aumenta el nivel de ampliación, aunque este aumento no es proporcional en todos los indicadores. Con el objetivo de poder observar la evolución del valor de PSNR obtenido en función de configuraciones diferentes, se han construido los gráficos 4.23, 4.24, 4.25, 4.26, 4.27, 4.28, , 4.30, 4.31 y 4.32. Todos ellos tratan de resumir la misma información obtenida a partir de imágenes de distinta naturaleza, y presentan en el eje de abcisas las diversas técnicas de reconstrucción empleadas, y en el de ordenadas un rango de PSNR dado. Para cada valor en el eje de abcisas, son mostradas cuatro barras representando los cuatro niveles de ampliación con los que se ha llevado a cabo el experimento. En primer lugar se muestra el gráfico correspondiente a la imagen Dicom1. Las reconstrucciones que mejores resultados obtienen para esta imagen de gammagrafía son, respectivamente la técnica PPH y la técnica ENO no jerárquico. A pesar de esto, los resultados entre distintas técnicas no son notables, ya que las técnicas disminuyen a razón de aproximadamente 1,2decibelios cuando se incrementa nivel de zoom empleado. Destacar la proximidad en decibelios que existe entre el nivel 1 de resolución con alguna técnica y el nivel 2 correspondiente a la misma técnica. En segundo lugar el gráfico 4.24 tampoco muestra diferencias significativas entre los resultados de las diferentes técnicas aplicadas sobre la imagen Dicom2. Para el primer nivel de ampliación (las diferencias se encuentran por debajo de un decibelio), los óptimos para ésta imagen se consiguen mediante la técnica PPH para el primer nivel y la técnica Racional para los restantes. En general, con esta imagen, la alternancia entre las distintas técnicas no produce un valor del PSNR medido sustancialmente diferente aunque visualmente unas imágenes parecen ser mucho mejor que otras. Esto puede ser producido por la propia naturaleza de la imagen: por un lado, influirá la resolución inicial de la imagen (debido al procedimiento inicial de reducción de la imagen original); y por otro, los contornos que describen la figura humana en la imagen original no son completamente nítidos, y la frontera entre la imagen humana y el resto no está 128 CAPÍTULO 4. EXPERIMENTOS Y RESULTADOS en este contexto. En segundo lugar, la técnica ENOSR solo ha resultado mejor que las restantes en dos de los casos estudiados. Se trata del cuarto nivel de ampliación en la imagen Pet-Tc1 y del valor obtenido en el primer nivel de la imagen Pet-Tc2. Dado que, por una parte, el valor obtenido en la primera de las imágenes es el mismo que el obtenido por la técnica PPH; y que por otra, el valor para el primer nivel que ha obtenido ésta técnica al ser aplicada sobre la imagen Pet-Tc2 dista en solo 0,3 del denominado óptimo para ese nivel, puede procederse a descartar ésta técnica como idónea para llevar a cabo una reconstrucción con resultados aceptables. Por todo lo expuesto y según los resultados obtenidos en este experimento, la técnica PPH se comporta de manera aceptable y obtiene mejores resultados que cualquier otra de las técnicas expuestas en este estudio. Esto ocurre en todos los tipos de imágenes empleadas, por lo que en esta memoria se considera que es la técnica que se debería utilizar para implementar un zoom digital, ya que se consiguen resultados en regiones suaves como lo haría un polinomio de Lagrange, y las regiones que presentan singularidades mediante la combinación de polinomios detallada anteriormente. Dependiendo del contexto en el cual se trabaje y de la calidad requerida, podría ser substituida por la técnica ENO con partición no jerárquica y en algún otro caso por la misma técnica con partición jerárquica del stencil; siempre y cuando las discontinuidades presentes en las imágenes con las cuales se trabaje, estén bien alejadas unas de otras. También se ha de descartar cualquier suposición relativa a la existencia de alguna correlación entre el tipo de imagen empleada y el la técnica idónea de interpolación. Después de observar los resultados obtenidos, a grandes rasgos todas las técnicas provén resultados que numéricamente fluctúan de manera similar entre distintos tipos de imágenes, salvo en las aspectos y casos comentados. Capítulo 5 Conclusiones El presente proyecto ha tratado de profundizar en las técnicas actuales para proporcionar un zoom digital aplicable a todo tipo de imágenes digitales, especialmente a imágenes médicas como los son las gammagrafías. Se ha detallado para todas ellas su motivación, idoneidad para el problema propuesto y sus característica, ademas de un pseudocódigo que muestra una posible implementación. En la sección de experimentos, han sido realizadas diversas pruebas con cada una de las técnicas y han sido evaluadas evaluadas mediante el uso de ciertos indicadores para imágenes digitales, y ello ha puesto de manifiesto las ventajas e inconvenientes que cada técnica produce, a la vista de los valores obtenidos y las imágenes de las reconstrucciones. Si bien, es verdad que las reconstrucciones no lineales tienden a disminuir tanto el fenómeno el fenómeno de Gibbs que aparece en las discontinuidades como la difusión de contornos, características típicas de los métodos lineales; la técnica que mejores resultados ha tenido presenta más similitudes en cuanto a resultados que los algoritmos lienales, que el resto de algoritmos para la reconstrucción de imágenes. Después de haber analizado los resultados, no se detecta ningún patrón o relación que establezca qué tipo de técnica utilizar para una determinada imagen. Dicho esto, la técnica que mejores resultados ha obtenido es la PPH tal y como se ha detallado, y que hacen de ésta una técnica idónea, fácil de implementar y eficiente de ser incluida en cualquier software que requiera un zoom digital. Por último, como posible ampliación de este trabajo, se propone el diseño de un experimento similar pero aplicado sobre imágenes en color. Ello solo conllevaría aplicar la técnica interpolatoria deseada a los diferentes canales de color de la imagen y quizás pudiera poner de manifiesto aspectos que en este estudio, o agudizar otros que hayan que los resultados de este estudio no hayan puesto de manifiesto. 129 130 CAPÍTULO 5. CONCLUSIONES Anexo i A continuación se muestra el código de la aplicación. Las archivos que han sido codificados y que aparecen a continuación son int_eno.m, int_enoSRm,int_weno.m,int_racional.m,int_pph.m,getMask.m, e indicadores.m. int_eno.m Se trata de una reconstrucción siguiendo la técnica ENO mediante una reconstrucción de cuatro puntos, tanto con seleción jerárquica, como no jerárquica. (Ver sección 3.3.1). Contiene las funciones: function [b] = int_eno(niv, a, met_particion) La función int_eno recibe como parámetros de entrada un nivel de zoom niv, la imagen de partida im, y un booleano met_particion que, como su propio nombre indica, establece el método de partición empleado. El valor 0 especifica un método no jerárquico de selección del stencil, y el 1 un método jerárquico de selección. El vector de salida contiene los valores de vademas de valores intercalados interpolados en forma de vector fila. function [f] = eno_zoom(v, n, met_particion) La función eno_zoom recibe como parámetros de entrada un vector, su longitud y un booleano denominado met_particion. El vector de salida es f, que contiene los valores de vademas de valores intercalados interpolados en forma de vector fila. int_enoSR.m 131 132 ANEXO I. Este fichero alberga las funciones int_enoSR,Stencil yenoSR_zoom. Implementa una reconstrucción ENO-SR de cuatro puntos basada en la técnica ENO jerárquica. (Ver sección 3.3.2) Implementa las funciones: function [b] = int_enoSR(niv, a) Esta función lleva a cabo la interpolación ENO Sucell Resolution. Recibe como parámetros los niveles de zoom y la imagen ay obtiene una ampliación de dicha imagen, b, de nivel niv. En ella, se llama a la función descrita a continuación, para cada fila y columna. function [f] = enoSR_zoom(v, n) La función enoSR_zoom recibe como parámetros de entrada un vector vy su longitud, y develve en fel vector de entrada con los valores intercalados interpolados. int_weno.m Este archivo contiene la implementación de las funciones int_weno yweno_zoom que llevan a cabo una interpolación WENO también de cuatro puntos. (Ver sección 3.3.3) En él se encuentran implementadas las funciones siguientes: function [b] = int_weno(niv, a) Esta función lleva a cabo la interpolación WENO mediante una reconstrucción de cuatro puntos. Recibe como parámetros los niveles de zoom y la imagen ay obtiene una ampliación de dicha imagen, b, de nivel niv. Para cada fila y columna, se realiza una llamada a la función que sigue. function [f] = weno_zoom(v, n) La función weno_zoom recibe como parámetros de entrada un vector vy su longitud, y develve en fel vector de entrada con los valores intercalados interpolados. int_racional.m Contiene las funciones int_racional yracional_zoom, necesarias para llevar a cabo una interpolación racional. (Ver sección 3.3.4) Contiene las siguientes funciones: 133 function [b] = int_racional(niv, a) Ésta funcion lleva a cabo la interpolación de una imagen a, basándose en la técnica Racional, y con los niveles de zoom niv pasados como parámtero. Para cada fila y columna de atiene lugar una llamada a la función racional_zoom. function [f] = racional_zoom(v, n) Ésta es la que lleva a cabo propiamente la reconstrucción Racional. Para cada vector fila v, de longitud n, se devuelve un vector fcon los valores de v, además de los respectivos valores intercalados obtenidos mediante ésta técnica. int_pph.m Fichero que implementa las funciones int_pph ypph_zoom Se trata también de una reconstrucción de cuatro puntos, pero llevada a cabo mediante una reconstrucción PPH. (Ver sección 3.3.5) En él se encuentran implementadas: function [b] = int_int_pph(niv, a) Análogamente a las anteriores, ésta función recibe los niveles de zoom deseados y la imagen de partida, a; y obtiene en bla imagen resultante del proceso. function [f] = pph_zoom(v, n) Ésta es la encargada de proceder a la reconstrucción de cada vector fila recibido como parámetro, v, a partir de éste y de su longitud n. Como salida obtiene f, que consta de los valores del vector vrecibido, ademas de aquellos que han sido obtenidos como resultado de la aplicación de la técnica PPH. getMask.m Esta función calcula las máscaras necesarias para el esquema de subdivisión basado en Lagrange. Su cabecera es function [mas,numer,deno]=getMask(l, r) 134 ANEXO I. Los argumentos de entrada son lyr, que son respectivamente la canitdad de puntos existentes a derecha e izquierda del punto de interés en el stencil seleccionado. Como salida se obtienen mas, indicadores.m En este archivo se encuentra el código de los diferentes indicadores implementados. Las funciones que los obtienen son, respectivamente, MeanSquareError, PeakSignaltoNoiseRatio,NormalizedCrossCorrelation,AverageDifference, MaximumDifference yNormalizedAbsoluteError. Todas ellas reciben como entrada dos imágenes: una de referencia, origImg, y otra imagen distImage que en este caso es una aproximación a la primera. function ECM = MeanSquareError(origImg, distImg) La salida que produce la función es el Error Cuadrático Medio, definido anteriormente en el punto4.1.2. function PSNR = PeakSignaltoNoiseRatio(origImg, distImg) Como resultado obtiene el PSNR según se ha definido en la sección 4.1.3. function NK = NormalizedCrossCorrelation(origImg, distImg) Produce el Coeficiente de Correlación Normalizada, tal y como se ha visto anteriormente en 4.1.4. function AD = AverageDifference(origImg, distImg) Calcula la diferencia media entre las dos imágenes de entrada. Ver sección 4.1.5. function MD = MaximumDifference(origImg, distImg) Obtiene la diferencia máxima entre dos píxeles como se ha visto en 4.1.6. function EAN = NormalizedAbsoluteError(origImg, distImg) Da como resultado el Error Absoluto Normalizado, definido en la sección ??. 15/09/11 10:29 C:\Users\Juan\Desktop\Indicadores.m 1 of 4 %****** % Main %****** % Programa para la obtención de indicadores de suavidad para imágenes % digitales. %ENTRADA: % origImg - Imagen de referencia. % origImg - Imagen interpolada. %SALIDA: % Se muestran por salida estándar los siguientes indicadores: % ECM - Error Cuadrático Medio % PSNR - Señal a ruido de pico % NK - Coeficiente de Correlación Cruzada Normalizada % AM - Diferencia Media % DM - Diferencia Máxima % EAN - Error Absoluto Normalizado function main(origImg, distImg); %Error Cuadrático Medio. ECM = MeanSquareError(origImg, distImg); disp('Error cuadratico medio = '); disp(ECM); %Peak Signal to Noise Ratio. PSNR = PeakSignaltoNoiseRatio(origImg, distImg); disp('Peak Signal to Noise Ratio = '); disp(PSNR); %Coeficiente de Correlación Normalizada. NK = NormalizedCrossCorrelation(origImg, distImg); disp('Coeficiente de Correlacion Cruzada Normalizado = '); disp(NK); %Diferencia Media. AD = AverageDifference(origImg, distImg); disp('Diferencia Media = '); disp(AD); %Máxima Diferencia. DM = MaximumDifference(origImg, distImg); disp('Diferencia Maxima = '); disp(DM); %Error Absoluto Normalizado EAN = NormalizedAbsoluteError(origImg, distImg); disp('Error Absoluto Normalizado = '); disp(EAN); %*********************** % Error Cuadrático Medio 15/09/11 10:29 C:\Users\Juan\Desktop\Indicadores.m 2 of 4 %*********************** %ENTRADA: % origImg - Imagen de referencia. % distImg - Imagen interpolada. %SALIDA: % ECM - Error Cuadrático Medio. %*********************** % Error Cuadrático Medio %*********************** %ENTRADA: % origImg - Imagen de referencia. % distImg - Imagen interpolada. %SALIDA: % ECM - Error Cuadrático Medio. function ECM = MeanSquareError(origImg, distImg) Image1 = double(origImg); Image2 = double(distImg); ECM =(sum((Image1(:)-Image2(:)).^2)/numel(Image1)); %******** % PSNR %******** %ENTRADA: % origImg - Imagen de referencia. % origImg - Imagen interpolada. %SALIDA: % PSNR - Peak Signal to Noise Ratio function PSNR = PeakSignaltoNoiseRatio(origImg, distImg) origImg = double(origImg); distImg = double(distImg); PSNR = 20*log10(255/sqrt(MeanSquareError(origImg, distImg))); end %*********************************************** % Coeficiente de Correlación Cruzada Normalizado %*********************************************** %ENTRADA: % origImg - Imagen de referencia. % origImg - Imagen interpolada. %SALIDA: 15/09/11 10:29 C:\Users\Juan\Desktop\Indicadores.m 3 of 4 % NK - Coeficiente de Correlación Cruzada Normalizado function NK = NormalizedCrossCorrelation(origImg, distImg) origImg = double(origImg); distImg = double(distImg); NK = sum(sum(origImg .* distImg)) / sum(sum(origImg .* origImg)); end %******************* % Diferencia Media %******************* %ENTRADA: % origImg - Imagen de referencia. % distImg - Imagen interpolada. %SALIDA: % DM - Diferencia Media. function DM = AverageDifference(origImg, distImg) origImg = double(origImg); distImg = double(distImg); [M N] = size(origImg); error = origImg - distImg; DM = sum(sum(error)) / (M * N); end %******************* % Diferencia Máxima %******************* %ENTRADA: % origImg - Imagen de referencia. % distImg - Imagen interpolada. %SALIDA: % MD - Diferencia Máxima. function DM = MaximumDifference(origImg, distImg) origImg = double(origImg); distImg = double(distImg); error = origImg - distImg; DM = max(max(error)); end 18/09/11 10:50 F:\UPV\Proyecto\ArchivosMa...\int_enoSR.m 2 of 6 % SALIDA: % f - Vector fila con valores interpolados. function [f] = enoSR_zoom(v, n) % Cramos el vector f, vector de salida de la funcion. f = zeros(1, 2*n); f(1:2:2*n-1) = v(1:n); % Obtencion de la mascaras. [m1] = getMask(orden/2, orden/2); [m3] = getMask(orden/2+1, orden/2-1); [m2] = getMask(orden/2-1, orden/2+1); % Prediccion del primer elemento con mascara lineal; f(2) = m2(1)*v(1) + m2(2)*v(2) + m2(3)*v(3) + m2(4)*v(4); % Prediccion del segundo elemento. % Cáculo de las diferencias dividas. elem1 = (v(2) - v(1)); elem2 = (v(3) - v(2)); elem3 = (v(4) - v(3)); elem4 = (v(5) - v(4)); dif1 = (elem2 - elem1); a_dif1 = abs(dif1); dif2 = (elem3 - elem2); a_dif2 = abs(dif2); if(a_dif1 <= a_dif2) f(4) = m1(1)*v(1) + m1(2)*v(2) + m1(3)*v(3) + m1(4)*v(4); else dif3 = (elem4 - elem3); t1 = dif2 - dif1; a_t1 = abs(t1); t2 = dif3 - dif2; a_t2 = abs(t2); if(a_t1 <= a_t2) f(4) = m1(1)*v(1) + m1(2)*v(2) + m1(3)*v(3) + m1(4)*v(4); else f(4) = m2(1)*v(2) + m2(2)*v(3) + m2(3)*v(4) + m2(4)*v(5); end end % Prediccion del tercer elemento con ENO jerarquico. elem1 = (v(2) - v(1)); elem2 = (v(3) - v(2)); elem3 = (v(4) - v(3)); elem4 = (v(5) - v(4)); elem5 = (v(6) - v(5)); dif1 = (elem2 - elem1); a_dif1 = abs(dif1); dif2 = (elem3 - elem2); a_dif2 = abs(dif2); dif3 = (elem4 - elem3); a_dif3 = abs(dif3); dif4 = (elem5 - elem4); a_dif4 = abs(dif4); if(a_dif2<=a_dif3) t1 = a_dif2 - a_dif1; a_t1 = abs(t1); 18/09/11 10:50 F:\UPV\Proyecto\ArchivosMa...\int_enoSR.m 3 of 6 t2 = a_dif3 - a_dif2; a_t2 = abs(t2); if(a_t2 <= a_t1) f(6) = m1(1)*v(2) + m1(2)*v(3) + m1(3)*v(4)+ m1(1)*v(5); else f(6) = m3(1)*v(1) + m3(2)*v(2) + m3(3)*v(3) + m3(4)*v(4); end else t2 = dif2 - dif2; a_t2= abs(t2); t3 = dif4 - dif3; a_t3 = abs(t3); if(a_t2 <= a_t3) f(6) = m1(1)*v(3) + m1(2)*v(4) + m1(3)*v(5)+ m1(1)*v(6); else f(6) = m3(1)*v(3) + m3(2)*v(4) + m3(3)*v(5) + m(4)*v(6); end end % Prediccion de valores intermedios for j=8:2:2*n-8 %Paso 1 : para cada intervalo determinamos los stencil a la izquierda %y a la derecha.. Si no esta descentrado aplicamos eno-h. s_izq = stencil(v(j/2-3:j/2+2)); s_der = stencil(v(j/2-1:j/2+4)); if(s_izq == 'I' & s_der == 'D') %Continuamos con el metodode resolucion subcelda comprobamos si %g(x_i*g(x_(i+1)) <0. q11 = v(j/2); q12 = -v(j/2-1) + 4*v(j/2-2) - 6*v(j/2-1) + 4*v(j/2); q21 = 3*v(j/2+1) - 6*v(j/2+2)+ 4*v(j/2+3) - v(j/2+4); q22 = (j/2+1); g1 = q21 - q11; g2 = q22 - q12; if(g1*g2 < 0) %Hay una discontinuidad en el intervalo. q1m = -5/16*v(j/2-3) + 21/16*v(j/2-2) - 35/16*v(j/2-1) + 35/16*v(j/2); q2m = 35/16*v(j/2+1) - 35/16*v(j/2+2) - 21/16*v(j/2+3) - 5/16*v(j/2+4); %funcion g evaluada en punto medio. gm = q2m - q1m; if(g1*gm < 0) f(j) = q2m; else f(j)= q1m; end else %ENO jerarquico. p1 = v(j/2-1) - v(j/2-2); p2 = v(j/2) - v(j/2-1); p3 = v(j/2+1) - v(j/2); p4 = v(j/2+2) - v(j/2+1); p5 = v(j/2+3) - v(j/2+2); dif1 = p2 - p1; dif2 = p3 - p2; dif3 = p4 - p3; 18/09/11 10:50 F:\UPV\Proyecto\ArchivosMa...\int_enoSR.m 4 of 6 dif4 = p5 - p4; if(as2 <= as3) aux1 = dif2 - dif1; a_aux1 = abs(aux1); aux2 = dif3 - dif2; a_aux2 = abs(aux2); if(a_aux2 <= a_aux1) f(j) = m1(1)*v(1) + m1(2)*v(2) + m1(3)*v(3) + m1(4)*v(4); else f(j) = m3(1)*v(j/2-2) + m3(2)*v(j/2-1) + m3(3)*v(j/2) + m3(4)*v (j/2+1); end else aux2 = dif3 - dif2; a_aux2 = abs(aux2); aux3 = dif3 - dif3; a_aux3 = abs(aux3); if(a_aux2 <= a_aux3) f(j) = m1(1)*v(1) + m1(2)*v(2) + m1(3)*v(3) + m1(4)*v(4); else f(j) = m2(1)*v(j/2) + m2(2)*v(j/2+1) + m2(3)*v(j/2+2) + m2(4)*v (j/2+3); end end end else % Aplicamos ENO jerárquico. elem1 = (v(j/2-1) - v(j/2-2)); elem2 = (v(j/2) - v(j/2-1)); elem3 = (v(j/2+1) - v(j/2)); elem4 = (v(j/2+1) - v(j/2+1)); elem5 = (v(j/2+3) - v(j/2+2)); dif1 = (elem2 - elem1); dif2 = (elem3 - elem2); a_dif2 = abs(dif2); dif3 = (elem4 - elem3); a_dif3 = abs(dif3); dif4 = (elem5 - elem4); if(a_dif2 <= a_dif3) aux1 = (elem2 - elem1); a_aux1 = abs(aux1); aux2 = (elem3 - elem2); a_aux2 = abs(aux2); if(a_aux2 <= a_aux1) f(j) = m1(1)*v(j/2-1) + m1(2)*v(j/2) + m1(3)*v(j/2+1)+ m1(4)*v(j/2+2); else f(j) = m3(1)*v(j/2-2) + m3(2)*v(j/2-1) + m3(3)*v(j/2) + m3(4)*v(j/2+1); end else aux2 = (dif3 - dif2); a_aux2 = abs(aux2); aux3= (dif4 - dif3); a_aux3 = abs(aux3); if(a_aux2 <= aux3) f(j) = m1(1)*v(j/2-1) + m1(2)*v(j/2) + m1(3)*v(j/2+1)+ m1(4)*v(j/2+2); else f(j) = m2(1)*v(j/2) + m2(2)*v(j/2+1) + m2(3)*v(j/2+2) + m2(4)*v(j/2+3); end end end end 18/09/11 10:50 F:\UPV\Proyecto\ArchivosMa...\int_enoSR.m 5 of 6 % Prediccion penultimo elemento. elem1 = (v(n-4) - v(n-5)); elem2 = (v(n-3) - v(n-4)); elem3 = (v(n-2) - v(n-3)); elem4 = (v(n-1) - v(n-2)); elem5 = (v(n) - v(n-1)); dif1 = (elem2 - elem1); a_dif1 = abs(dif1); dif2 = (elem3 - elem2); a_dif2 = abs(dif2); dif3 = (elem4 - elem3); a_dif3 = abs(dif3); dif4 = (elem5 - elem4); if(a_dif2 <= a_dif3) aux1 = (elem2 - elem1); a_aux1 = abs(aux1); aux2 = (elem3 - elem2); a_aux2 = abs(aux2); if(a_aux2 <= a_aux1) f(2*n-6) = m1(1)*v(n-4) + m1(2)*v(n-3) + m1(3)*v(n-2) + m1(4)*v(n-1); else f(2*n-6) = m3(1)*v(n-5) + m3(2)*v(n-4) + m3(3)*v(n-3) + m3(4)*v(n-2); end else aux2 = (elem3 - elem2); a_aux2 = abs(aux2); aux3 = (elem4 - elem3); a_aux3 = abs(aux2); if(a_aux2 <= a_aux3) f(2*n-6) = m1(1)*v(n-4) + m1(2)*v(n-3) + m1(3)*v(n-2) + m1(4)*v(n-1); else f(2*n-6) = m2(1)*v(n-3) + m2(2)*v(n-2) + m2(3)*v(n-1) + m2(4)*v(n); end end % Prediccion penultimo elemento. elem1 = (v(n-3) - v(n-4)); elem2 = (v(n-2) - v(n-3)); elem3 = (v(n-1) - v(n-2)); elem4 = (v(n) - v(n-1)); dif1 = (elem2 - elem1); a_dif1 = abs(dif1); dif2 = (elem3 - elem2); a_dif2 = abs(dif2); dif3 = (elem4 - elem3); a_dif3 = abs(dif3); if(a_dif3 <= a_dif2) f(2*n-4) = m1(1)*v(n-3) + m1(2)*v(n-2) + m1(3)*v(n-1) + m1(4)*v(n); else aux1 = (dif2 - dif1); a_aux1 = abs(aux1); aux2 = (dif3 - dif2); a_aux2 = abs(aux2); if(a_aux2 <= a_aux1) f(2*n-4) = m1(1)*v(n-3) + m1(2)*v(n-2) + m1(3)*v(n-1) + m1(4)*v(n); else f(2*n-4) = m3(1)*v(n-4) + m3(2)*v(n-3) + m3(3)*v(n-2) + m3(4)*v(n-1); end end % Prediccion del ultimo elemento. f(2*n-2) = m2(1)*v(n) + m2(2)*v(n-1) + m2(3)*v(n-2) + m2(4)*v(n-3); end % Función: stencil % ****************** % ENTRADA: 18/09/11 10:50 F:\UPV\Proyecto\ArchivosMa...\int_enoSR.m 6 of 6 % v - . % SALIDA: % r - . function [r] = stencil(v) elem1 = (v(2) - v(1)); elem2 = (v(3) - v(2)); elem3 = (v(4) - v(3)); elem4 = (v(5) - v(4)); elem5 = (v(6) - v(5)); dif1 = (elem2 - elem1); a_dif1 = abs(dif1); dif2 = (elem3 - elem2); a_dif2 = abs(dif2); dif3 = (elem4 - elem3); a_dif3 = abs(dif3); dif4 = (elem5 - elem4); if(a_dif2 <= a_dif3) t1 = dif2 - dif1; a_t1 = abs(t1); t2 = dif3 - dif2; a_t2 = abs(t2); if(a_t2 <= a_t1) r = 'C'; else r = 'I'; end else t2 = dif3 - dif2; a_t2 = abs(t2); t3 = dif4 - dif3; a_t3 = abs(t3); if(a_t2 <= a_t3) r = 'C'; else r = 'I'; end end end 18/09/11 10:52 F:\UPV\Proyecto\ArchivosMatlab\int_weno.m 1 of 3 %********************* % Interpolacion WENO %********************* %Función: int_weno %****************** %ENTRADA: % niv - Niveles de zoom. % a - Imagen a interpolar. %SALIDA: % b - Imagen interpolada. function [ b ] = int_weno(niv, a) [n m] = size(a); %Bucle para los niveles de zoom. for k=1:niv % Filas y columnas para la imagen interpolada. n_filas = 2*n; n_columnas = 2*m; % Creamos la nueva matriz. b=zeros(n_filas,n_columnas); b(1:2:n_filas,1:2:n_columnas) = a(1:n,1:m); %Algoritmo de prediccion para las columnas. for j=1:2:n_columnas b(1:n_filas,j) = weno_zoom(b(1:2:n_filas, j), n)'; end %Algorimo de preciccion para las filas. for i =1:n_filas b(i, 1:n_columnas) = weno_zoom(b(i,1:2:n_columnas),m); end %Actualización de variables. n = n_filas; m = n_columnas; % Inicializacion matriz de zoom. a = b; end b = uint8(b); end % Función: weno_zoom % ****************** 18/09/11 10:52 F:\UPV\Proyecto\ArchivosMatlab\int_weno.m 2 of 3 % ENTRADA: % v - Vector fila al que se le va a aplicar el algoritmo. % n - Longitud del vector v. % SALIDA: % f - Vector fila con valores interpolados. function [f] = weno_zoom(v, n) f = zeros(1, 2*n); f(1:2:2*n-1) = v(1:n); %Obtencion de la mascara ENO. [m1] = getMask(orden/2, orden/2); [m3] = getMask(orden/2-1, orden/2+1); [m2] = getMask(orden/2+1, orden/2); %Prediccion del primer elemento. f(2) = m1(1)*v(1) + m1(2)*v(2) + m1(3)*v(3) + m1(4)*v(4); %Prediccion del segundo elemento. elem2 = abs(-v(1) + 3*v(2) - 3*v(3) + v(4)); elem3 = abs(-v(2) + 3*v(3) - 3*v(4) + v(5)); if(elem2 < elem3) f(4) = m1(1)*v(1) + m1(2)*v(2) + m1(3)*v(3) + m1(4)*v(4); else f(4) = m2(1)*v(2) + m2(2)*v(3) + m2(3)*v(4) + m2(4)*v(5); end %Prediccion valores intermedios epsilon = 10^(-6); for j=6:2:2*n-6 elem1 = (v(j/2-1) - v(j/2-2)); elem2 = (v(j/2) - v(j/2-1)); elem3 = (v(j/2+1) - v(j/2)); elem4 = (v(j/2+2) -v(j/2+1)); elem5 = (v(j/2+3) - v(j/2+2)); %Calculo de los indicadores de suavidad: i_izq = 1/2 * ((elem2 - elem1)^2 + (elem3 - elem2)^2)+(elem3 - 2*elem2 + elem1)^2; i_central = 1/2 * ((elem3 - elem2)^2 + (elem4 - elem3)^2) + (elem4 - 2*elem3 + elem2)^2; i_der = 1/2 * ((elem4-elem3)^2 + (elem5-elem4)^2) + (elem5 - 2*elem4 + elem3)^2; %Calulo de los numeradores delos pesos. alfa_izq = 3./16/(epsilon + i_izq)^3; alfa_der = 10./16/(epsilon + i_der)^3; alfa_central = 3./16/(epsilon + i_central)^3; peso_izq = alfa_izq/(alfa_izq + alfa_der + alfa_central); peso_der = alfa_der/(alfa_izq + alfa_der + alfa_central); peso_central = alfa_central/(alfa_izq + alfa_der + alfa_central); %Suma de los 3 polinomios por sus pesos. q_central = m1(1)*v(j/2-1) + m1(2)*v(j/2) + m1(3)*v(j/2+1)+ m1(4)*v(j/2+2); q_izquierdo = m3(1)*v(j/2-2) + m3(2)*v(j/2-1) + m3(3)*v(j/2) + m3(4)*v(j/2+1); q_derecho = m2(1)*v(j/2) + m2(2)*v(j/2+1) + m2(3)*v(j/2+2) + m2(4)*v(j/2+3); f(j) = peso_izq*q_izquierdo + peso_der*q_derecho + peso_central*q_central; 18/09/11 10:52 F:\UPV\Proyecto\ArchivosMatlab\int_weno.m 3 of 3 end %Predicción del pnultimo elemento. elem1 = abs(-v(n-4) + 3*v(n-3) - 3*v(n-2) + v(n-1)); elem2 = abs(-v(n-3) + 3*v(n-2) - 3*v(n-1) + v(n)); if(elem2 <= elem1) %f(2*n-4) = (-v(n-3)+9*v(n-2)+9*v(n-1)-v(n))/16; f(2*n-4) = m1(1)*v(n-3) + m1(2)*v(n-2) + m1(3)*v(n-1) + m1(4)*v(n); else f(2*n-4) = m3(1)*v(n-4) + m3(2)*v(n-3) + m3(3)*v(n-2) + m3(4)*v(n-1); end %prediccion del ultimo elemento f(2*n-2) = m2(1)*v(n) + m2(2)*v(n-1) + m2(3)*v(n-2) + m2(4)*v(n-3); end 18/09/11 10:52 F:\UPV\Proyecto\Archivo...\int_racional.m 1 of 2 %****************** % Interpolacion Racional %****************** %Función: int_racional %****************** %ENTRADA: % niv - Niveles de zoom. % a - Imagen a interpolar. %SALIDA: % b - Imagen interpolada. function [ b ] = int_racional(niv, a) [n m] = size(a); %Bucle para los niveles de zoom. for k=1:niv % Filas y columnas para la imagen interpolada. n_filas = 2*n; n_columnas = 2*m; % Creamos la nueva matriz. b=zeros(n_filas,n_columnas); b(1:2:n_filas,1:2:n_columnas) = a(1:n,1:m); % Algoritmo de prediccion para las columnas. for j=1:2:n_columnas b(1:n_filas,j) = racional_zoom(b(1:2:n_filas, j), n)'; end % Algorimo de preciccion para las filas. for i =1:n_filas b(i, 1:n_columnas) = racional_zoom(b(i,1:2:n_columnas),m); end %Actualización de variables. n = n_filas; m = n_columnas; % Inicializacion matriz de zoom. a = b; end b = uint8(b); end % Función: racional_zoom % ****************** % ENTRADA: % v - Vector fila al que se le va a aplicar el algoritmo. 18/09/11 10:52 F:\UPV\Proyecto\Archivo...\int_racional.m 2 of 2 % n - Longitud del vector v. % SALIDA: % f - Vector fila con valores interpolados. function [f] = racional_zoom(v, n) f = zeros(1, 2*n); f(1:2:2*n-1) = v(1:n); % Obtencion de la mascara ENO. [m1] = getMask(orden/2, orden/2); [m3] = getMask(orden/2-1, orden/2+1); [m2] = getMask(orden/2+1, orden/2); % Prediccion del primer elemento. f(2) = m1(1)*v(1) + m1(2)*v(2) + m1(3)*v(3) + m1(4)*v(4); % Prediccion del segundo elemento. elem2 = abs(-v(1) + 3*v(2) - 3*v(3) + v(4)); elem3 = abs(-v(2) + 3*v(3) - 3*v(4) + v(5)); if(elem2 < elem3) f(4) = m1(1)*v(1) + m1(2)*v(2) + m1(3)*v(3) + m1(4)*v(4); else f(4) = m2(1)*v(2) + m2(2)*v(3) + m2(3)*v(4) + m2(4)*v(5); end % Prediccion valores intermedios epsilon = 10^(-6); for j=6:2:2*n-4 elem1 = (v(j/2+1)-v(j/2+2)); elem2 = (v(j/2-1) - v(j/2)); %pesos a = 1/2; w0 = (1 + a*(elem1^2))/(2+a*(elem2^2+elem1^2)); w1 = (1 + a*(elem2^2))/(2+a*(elem2^2+elem1^2)); f(j) = w0*v(j/2) + w1*v(j/2+1); end % Prediccion del ultimo elemento. f(2*n-2) = m2(1)*v(n) + m2(2)*v(n-1) + m2(3)*v(n-2) + m2(4)*v(n-3); end