Modelado general de ruido en imágenes. Aplicaciones médicas

ESCUELA DE INGENIERÍA DE TELECOMUNICACIÓN Y ELECTRÓNICA PROYECTO FIN DE CARRERA MODELADO GENERAL DE RUIDO EN IMÁGENES. APLICACIONES MÉDICAS Titulación: Ingeniero de Telecomunicación Autor: Guillermo Valentín Socorro Marrero Tutor: Dr. Eduardo Rovaris Romero Fecha: Junio de 2016 ESCUELA DE INGENIERÍA DE TELECOMUNICACIÓN Y ELECTRÓNICA PROYECTO FIN DE CARRERA MODELADO GENERAL DE RUIDO EN IMÁGENES. APLICACIONES MÉDICAS HOJA DE FIRMAS Alumno Fdo.: Guillermo Valentín Socorro Marrero Tutor Fdo.: Dr. Eduardo Rovaris Romero Fecha: Junio de 2016 ESCUELA DE INGENIERÍA DE TELECOMUNICACIÓN Y ELECTRÓNICA PROYECTO FIN DE CARRERA MODELADO GENERAL DE RUIDO EN IMÁGENES. APLICACIONES MÉDICAS HOJA DE EVALUACIÓN Calificación: ______________________________ Presidente Fdo.: Vocal Secretario/a Fdo.: Fdo.: Fecha: Junio de 2016 iii Agradecimientos La redacción de esta memoria y lo que supone el hecho de que usted la esté leyendo no hubiese sido posible sin la ayuda y la colaboración de muchas personas. Vaya desde estas líneas mi agradecimiento especial a las siguientes: A mis tutores, Karl y Eduardo, sin cuya orientación y ayuda nunca hubiese terminado este proyecto. A Javier Miranda, Suni, Samuel, Mari Carmen, Tomás y Ángel Capote, compañeros y profesores de Teleco a los que admiro y con los que me resisto a que dejen de serlo. A mis compañeros en el SIANI y sobre todo a Rafael Montenegro que en su momento cometió la insensatez de incorporarme a su grupo, que ahora es también el mío. Allí estoy creciendo a una velocidad que creía reservada a la adolescencia. Gracias a todos por el apoyo y la infinita paciencia. A mi familia: mi madre, a Pili, Domingo, Carolina y Carlos, por hacerme sentir querido y arropado. A Yurena, Jesús, Romina, y Laura que saben lo importante que es para mí la amistad y que han “sufrido” con resignación mis frecuentes ausencias durante la carrera. Mención cum laude a Yanis y Omar quienes, aparte de ser buenos amigos, fueron también compañeros de fatigas. Y por último, por qué no decirlo, a mi cerebro, por haber estado siempre ahí a pesar de nuestras diferencias. A todos, muchas gracias. iv v Índice de contenidos 1. INTRODUCCIÓN ......................................................................................................... 3 1.1 Motivación .............................................................................................................. 4 1.2 Antecedentes ........................................................................................................... 4 1.3 Líneas de investigación actuales ............................................................................. 7 1.4 Descripción del proyecto ........................................................................................ 8 1.4.1 Objetivos .......................................................................................................... 8 1.4.2 Fases del proyecto ........................................................................................... 9 1.5 Estructura de la memoria ...................................................................................... 10 2. BASE DE DATOS Y PREPROCESAMIENTO DE LAS IMÁGENES .................... 13 2.1 Base de datos ........................................................................................................ 13 2.1.1 Descripción de la base de datos ..................................................................... 14 2.1.2 Imágenes de prueba ....................................................................................... 15 2.1.3 Imágenes de carácter general ......................................................................... 16 2.1.4 Imágenes médicas .......................................................................................... 18 2.2 Preprocesamiento de las imágenes ....................................................................... 21 2.2.1 Formato de almacenamiento de las imágenes ............................................... 23 3. MODELO DE RUIDO ................................................................................................ 25 3.1 Introducción .......................................................................................................... 25 3.2 Clasificación de ruido en imágenes ...................................................................... 26 3.2.1 Ruido aditivo ................................................................................................. 27 3.2.2 Ruido multiplicativo ...................................................................................... 29 3.2.3 Modelo generalizado de ruido dependiente de la señal de Selva-Alparone .. 31 3.2.4 Otros tipos de ruido ....................................................................................... 32 3.3 Ruido en aplicaciones médicas ............................................................................. 34 3.3.1 Ruido en imágenes de rayos X ...................................................................... 34 3.3.2 Ruido en imágenes de ultrasonidos ............................................................... 35 3.3.3 Ruido en imágenes de resonancia magnética ................................................ 36 3.4 Modelo de ruido propuesto ................................................................................... 39 4. FILTROS DE RESTAURACIÓN DE IMÁGENES ................................................... 42 4.1 Introducción .......................................................................................................... 42 4.2 Estrategias de filtrado ........................................................................................... 44 4.2.1 Filtros en el dominio espacial ........................................................................ 44 xii Índice de tablas Tabla 2.1 Volúmenes de la base de datos del USC-SIPI. ................................................... 17 Tabla 2.2 Características de las colecciones del TCIA incluidas en la base de datos. ........ 19 Tabla 2.3 Descripción de los campos del formato de archivo IMD. ................................... 24 Tabla 3.1 Valores del exponente  para distintos sistemas de imágenes. ........................... 29 Tabla 3.2 Valores del modelo propuesto para ruido en imágenes de rayos X. ................... 40 Tabla 3.3 Valores del modelo propuesto para imágenes de resonancia magnética. ........... 41 Tabla 4.1 Expresiones matemáticas de las medias más empleadas en filtros MF. ............. 45 Tabla 7.1 Medida del error en la estimación basada en filtros (prefiltrado) ....................... 88 Tabla 7.2 Valores de error para la estimación sin suavizar. ................................................ 91 Tabla 7.3 Valores de error para la estimación suavizada. ................................................... 91 Tabla 7.4. Valores de error en la estimación (𝜎𝑟𝑒𝑓 = 1). ................................................... 93 Tabla 7.5 Valores de error en la estimación (𝜎𝑟𝑒𝑓 = 2). .................................................... 93 Tabla 7.6 Valores de error en la estimación (𝜎𝑟𝑒𝑓 = 5). .................................................... 95 Tabla 7.7 Valores de error en la estimación (𝜎𝑟𝑒𝑓 = 20). .................................................. 97 Tabla 7.8 Valores de error en la estimación de 𝜎 en imágenes de rayos X. ....................... 99 Tabla 7.9 Valores de error en la estimación de 𝜎 en imágenes de resonancia magnética. 101 Tabla 12.1 Costes de amortización de los recursos software. ........................................... 117 Tabla 12.2 Costes de amortización de los recursos hardware. .......................................... 117 Tabla 12.3 Factor de corrección en función del tiempo de realización del proyecto. ....... 119 Tabla 12.4 Cálculo del presupuesto de ejecución material. .............................................. 119 Tabla 12.5 Coeficiente de ponderación en función del presupuesto. ................................ 120 Tabla 12.6 Cálculo del coste del material fungible. .......................................................... 121 Tabla 12.7 Presupuesto total del proyecto. ...................................................................... 122 xiii Siglas y Acrónimos CCD, Charge Coupled Device. CDF, Cummulative Density Function. CR, Computed Radiolography. CT, Computed Tomography. COIT, Colegio Oficial de Ingenieros de Telecomunicación. DICOM, Digital Imaging and Communications in Medicine. DCT, Discrete Cosine Transform. DFT, Discrete Fourier Transform. DWT, Discrete Wavelet Transform. FDP, Función de Densidad de Probabilidad. FBTD, Filering-based Transform Domain. GIMET, Grupo de Imagen Médica. ICA, Independent Component Analysis. IMD, Image Data. LMMSE, Linear Minimum Mean Square Error. MF, Mean Filter. MR, Magnetic Resonance. NM, Nuclear Medicine. NLM, Non-local Means. PCA, Principal Component Analysis. PET, Positron Emission Tomography. PDE, Partial Differential Equation. SAR, Synthetic Aperture Radar. SDN, Signal-Dependent Noise. SIPI, Signal and Image Processing Institute. SNF, SUSAN Neighbourhood Filter. SNR, Signal-to-Noise Ratio. SURE, Stein’s Unbiased Risk Estimate SUSAN, Smallest Univalue Segment Assimilating Nucleus. TCIA, The Cancer Imaging Archive. US, Ultrasound. YNF, Yaroslavsky Neighbourhood Filter. xiv 1 PARTE I MEMORIA 2 3 1. INTRODUCCIÓN El procesamiento digital de imágenes ha experimentado un enorme desarrollo en las últimas décadas, tanto desde un punto de vista teórico, con la publicación de una gran cantidad de artículos en distintas ramas relacionadas con esta área del conocimiento, como desde un enfoque más práctico, con el desarrollo de nuevos sistemas de adquisición y modalidades de imagen. Este crecimiento ha tenido su reflejo en ámbitos que van desde las telecomunicaciones, la biología, la ciencia de los materiales o la robótica hasta la medicina, en los que la imagen digital se ha consolidado como una herramienta imprescindible. Independientemente del sistema que se considere, uno de los principales problemas asociados al procesamiento digital de imágenes es el tratamiento del ruido. No en vano, las imágenes siempre están contaminadas en mayor o menor medida por señales indeseadas o presentan algún tipo de distorsión. Con la intención de caracterizar las fuentes de ruido y mejorar las prestaciones de los algoritmos de restauración, se está potenciando el desarrollo de modelos más precisos del ruido en imágenes digitales. En esta línea se enfoca este Proyecto Fin de Carrera en el que se plantea la definición de un modelo general de ruido con especial interés en el comportamiento del mismo en las distintas modalidades de imágenes médicas. 4 1.1 Motivación El estudio del ruido se ha considerado siempre un asunto relevante en el procesamiento digital de imágenes y el desarrollo de estrategias para su cancelación o atenuación sigue siendo un desafío. Desde un punto de vista práctico, la presencia del ruido representa un serio problema en el tratamiento de las imágenes y dificulta no sólo la visualización de las mismas sino también procesos posteriores que se desee llevar a cabo con ellas como pueden ser la segmentación, el realce de contraste, el reconocimiento de patrones o el diagnóstico por imagen. Existe una gran variedad de filtros disponibles para la restauración, pero muchos de ellos carecen de adaptatividad ya que no tienen en cuenta los niveles de ruido en la imagen, o que han sido diseñados suponiendo de antemano unos valores concretos. Esto penaliza de manera notable las prestaciones de tales sistemas cuando las condiciones no son las previstas y fuerza a que su uso quede restringido a aplicaciones concretas. Por otro lado, se dispone de filtros adaptativos preparados para funcionar de forma óptima si se les proporciona una estimación del nivel de ruido. Debe tenerse en cuenta que el nivel de ruido no sólo puede cambiar de una modalidad de imagen a otra, ni entre imágenes de la misma naturaleza, sino que puede presentar variaciones suaves dentro de la propia imagen. Surge por tanto la necesidad de desarrollar estrategias de estimación del nivel de ruido presente en imágenes, que proporcionen a estos filtros adaptativos la información necesaria para conseguir una mayor atenuación del efecto del ruido. 1.2 Antecedentes Un buen punto de partida para analizar cómo ha evolucionado el desarrollo de sistemas para la estimación del nivel de ruido presente en imágenes es el estudio publicado por S. I. Olsen [41] en 1993. En él se recogen los seis métodos que a criterio del autor conformaban la punta de lanza en esta rama de procesamiento de imágenes. Pasados treinta años, las ideas que inspiraron aquellos algoritmos siguen teniendo vigencia y continuamente aparecen publicados nuevos trabajos que recuperan algunas de las soluciones propuestas. 5 Se mencionan en el artículo algoritmos elementales basados en filtros, como el uso de medias aritméticas o filtros de mediana. Los resultados obtenidos en aplicaciones prácticas distaban mucho de ser satisfactorias por lo que comenzó a explorarse la conveniencia de no considerar la señal en su conjunto sino identificar y analizar regiones concretas donde la estimación devolvía mejores resultados, desechando el resto de píxeles de la imagen, sentando las bases de lo que hoy constituye la categoría de estrategias basadas en el procesamiento por bloques. Este cambio supuso una gran mejora, especialmente en señales con alto contraste o con abundancia de líneas delgadas que marcaban los contornos de los objetos observados, a la vez que ponía de manifiesto la necesidad de los criterios de selección de los bloques. Paralelamente a estos primeros ensayos, Bracho y Sanderson [15] se centraron en el estudio del gradiente y concretamente en su magnitud, publicando sus conclusiones en 1983. Se percataron de que dicha magnitud seguía una distribución de Rayleigh con parámetros dependientes del nivel de ruido presente en la imagen, siempre que la distribución del ruido fuese gaussiana y la imagen fuese relativamente homogénea. Se podía entonces inferir la varianza del ruido analizando la función de densidad acumulada de la magnitud del gradiente, un análisis que estaba basado en el cómputo de histogramas locales. Para mejorar la robustez del sistema frente a muestras alejadas de los valores normales se aplicaba un suavizado previo al histograma obtenido. En una tercera línea de investigación J.S. Lee [31] y G.A. Mastin [38], desarrollaban técnicas de filtrado adaptativo que requerían, para mejorar sus prestaciones, de un conocimiento aproximado del nivel de la relación señal a ruido en la imagen. La conclusión definitiva fue que una buena manera de estimar la varianza era calcularla localmente para todos los píxeles y considerar como valor estimado el mínimo de todos los valores de varianza del entorno. La idea subyacente es que los promedios locales bastaban para protegerse ante valores muy alejados de la media, mientras que la elección de la varianza mínima en el entorno descontaba el efecto de la estructura de la imagen, que siempre provocaba una sobreestimación del estadístico. Seis años más tarde, en 1989, el propio Lee y K. Hoppel [33] propusieron un método alternativo basado en la transformada de Hough que permitía discernir las componentes de ruido dependientes de la imagen de aquéllas que no lo eran. 6 La última estrategia recogida por Olsen en su compendio era la propuesta por P. Meer y J. Jolion [39] en 1990. Consistía en la estimación de la varianza a distintas escalas, considerando progresivamente bloques de tamaño mayor, en una estructura piramidal. En un análisis posterior, se contrastaban las curvas de evolución de las estimaciones en función de las escala con curvas teóricas características para niveles concretos de ruido. Un algoritmo iterativo permitía obtener el valor que proporcionaba el mejor ajuste. En 1993, J. Immerkaer [28] publicó un método de estimación rápida de la varianza, utilizando simplemente un filtro de convolución con una máscara de pequeñas dimensiones (3x3). La mejora en la exactitud de las estimaciones fue notable. La máscara implementa el operador laplaciano, un operador que se mostró muy eficaz para discriminar la estructura de la imagen y proporcionar una señal de características estadísticas semejantes a las del ruido. El método es rápido, fácil de implementar y robusto incluso para valores de ruido elevado. Sin embargo, cuando la varianza es pequeña, el filtro de Immerkaer siempre tiende a sobreestimar su valor, considerando como ruido los pequeños detalles presentes en la escena observada en la imagen. En 1999, K. Rank, M. Lendl y R. Unbehauen [48] propusieron otra estrategia para estimar la varianza del ruido. El algoritmo consta de tres pasos: 1. La descomposición del operador laplaciano en dos operadores diferenciales, vertical y horizontal, para procesar la imagen primero en una dirección y luego en la otra. 2. Cálculo del histograma de estimaciones locales de la varianza y 3. Análisis de los histogramas para inferir las funciones de distribución acumulada subyacentes y, a partir de ellas, estimar la varianza. Posteriormente, R. C. Bilcu y M. Vehvilainen [14] introdujeron como paso intermedio la detección de bordes en la imagen. Comprobaron que identificando los píxeles correspondientes a los contornos y no considerando los valores del laplaciano en eso puntos, el valor estimado para la varianza se acerca más a su valor verdadero. Con esta mejora, unida a un procedimiento más sofisticado en el análisis para el cálculo de los histogramas, se solventó parcialmente el problema de la sobreestimación provocada al considerar los detalles como ruido. La obtención de esos histogramas conlleva un incremento notable en las necesidades de cómputo, lo que constituye la principal desventaja de esta estrategia. 7 En el año 2002, A. Amer, A. Mitiche y E. Dubois [8] publicaron un artículo en el que se detalla la manera de identificar con mayor fiabilidad las regiones homogéneas en la imagen. Para ello definieron un conjunto de ocho máscaras específicas que se corresponden con patrones de regiones de alta variabilidad en la intensidad (esquinas, fronteras –alineadas con los ejes u oblicuas–, etc.), obteniendo mejores resultados que aplicando exclusivamente los detectores de bordes desarrollados hasta el momento. El trabajo posterior, publicado por D.-H. Shin, R.-H. Park, S. Yang y J.-H. Jung [19] en 2005, fue innovador en el sentido de que integraba en una misma solución las ideas de las dos familias de estrategias (las basadas en filtros y las basadas en bloques). En concreto, consta de una etapa selección inicial de bloques homogéneos, descartando los de alta variabilidad. En un segundo paso filtra la imagen en los bloques seleccionados obteniendo una aproximación inicial de la imagen ideal, que al sustraerla de la imagen ruidosa proporciona una estimación del ruido. El último trabajo de alto impacto en este período fue el publicado en 2008 por S.-C. Tai y S.-M.Yang [49]. En él se retoma el método de Immerkaer –con el filtrado laplaciano seguido de la estimación de la varianza– pero sugiriendo una detección de bordes al comienzo del proceso y no tras el filtrado como proponían Bilcu y Vehvilainen. 1.3 Líneas de investigación actuales En la actualidad se siguen empleando los métodos comentados en el apartado anterior, añadiendo mejoras principalmente en dos aspectos: diseñado estimadores más robustos y perfeccionado las medidas de homogeneidad de los bloques. Ejemplos de métodos englobados en el primer grupo son los que vienen desarrollando Yang y Tai [56], desde 2010, en los que se aprovechan las relaciones existentes entre estadísticos para estimar el que resulte conveniente en cada situación y deducir a partir de él los restantes. En la misma línea, se ha incorporado el análisis de componentes principales (PCA, Principal Component Analysis) a los estimadores clásicos. También son relevantes los trabajos de Santiago Aja-Fernández y Karl Krissian, aplicando estimadores lineales de mínimo error cuadrático medio (LMMSE, Linear 14 2.1.1 Descripción de la base de datos Debido al desarrollo y abaratamiento de los sistemas de captura de imágenes y la facilidad que las tecnologías de la comunicación ofrecen para su publicación, existe en la actualidad una enorme disponibilidad de imágenes. Son numerosas las aplicaciones que se apoyan o se sustentan en el uso de imágenes, abarcando ámbitos muy dispares que van desde la fotografía artística, los sistemas de vigilancia o autenticación, o la observación por satélite, hasta el diagnóstico por imagen. Para abordar el desarrollo de este Proyecto Fin de Carrera, en su doble vertiente –general, pero sesgada hacia las aplicaciones médicas– se hace necesario seleccionar, de entre todas las disponibles, un conjunto reducido de imágenes que sea suficientemente representativo de las situaciones que se pretende estudiar y que, por otro lado, contenga un número considerable de imágenes médicas. En concreto, se han definido tres clases bien diferenciadas de imágenes, a saber, imágenes de prueba, imágenes de carácter general e imágenes médicas. En la Figura 2.1 se presenta la estructura de la base de datos propuesta mientras que las características de cada clase de imagen se detallan en los siguientes apartados. Figura 2.1 Estructura de la base de datos del proyecto. 15 2.1.2 Imágenes de prueba La primera clase definida en la base de datos es la de imágenes de prueba. Estas imágenes no son el resultado de la observación de ningún objeto o fenómeno, sino que han sido construidas ex profeso con el fin de recrear escenarios de test. Facilitan el desarrollo y la evaluación de los algoritmos en un entorno controlado y, por otro lado, permiten diseñar experimentos que pongan de manifiesto las principales virtudes y debilidades de las distintas estrategias. (a) (b) (c) (d) Figura 2.2 Ejemplos de imágenes de prueba incluidas en la base de datos. (a) Cilindro y (b) su representación como superficie; (c) composición y (d) su imagen inversa. La base de datos de imágenes de prueba tiene un tamaño de 10.4 MB y contiene 11 imágenes en escala de grises, con distintas resoluciones y niveles de intensidad. En la Figura 2.2 se representan ejemplos de este tipo de imágenes, que incluyen figuras regulares u ovaladas, útiles para el estudio de comportamiento ante contornos rectos o curvos, respectivamente, y las composiciones de figuras giradas para evaluar, por ejemplo, la sensibilidad de los resultados de un algoritmo frente a la orientación respecto de las direcciones principales. 16 2.1.3 Imágenes de carácter general El segundo grupo de imágenes lo conforman las de carácter general, es decir, aquéllas que sí son el resultado de una observación, pero sin especializarse en un ámbito o aplicación concretos (como podría ser la imagen de rostros, de huellas digitales o radioastronómicas). Se trata por tanto de una miscelánea de imágenes de origen dispar. Este es el tipo de conjuntos de imágenes que se emplea habitualmente en la literatura para comparar las características de algoritmos de procesamiento digital. De entre las numerosas bases de datos disponibles para este propósito, se ha elegido la proporcionada por el Departamento de Ingeniería Eléctrica Ming Hsieh de la Universidad del Sur de California, bajo la denominación de Signal and Image Processing Institute (SIPI) ImageDatabase [61], por ser una de las más ampliamente utilizadas. En la Figura 2.3 se muestra el aspecto del portal del USC-SIPI, a través del cual se puede acceder a la base de datos y descargar libremente cualquier paquete temático de imágenes. Figura 2.3. Portal de acceso a la base de datos del USC-SIPI. La base de datos USC-SIPI se divide en cuatro volúmenes, atendiendo a los escenarios representados. La temática de cada volumen se detalla en la Tabla 2.1. Las resoluciones disponibles son 256x256, 512x512 y 1024x1024 muestras, mientras que la codificación de niveles de intensidad es de 8 bits por muestra para imágenes en escala de grises y 24 bits por muestra para las imágenes en color. De entre las disponibles, se han seleccionado para incluir en la base de datos de este proyecto un total de 12 imágenes, con un tamaño resultante de 40.3 MB. Cuatro de estas imágenes se reproducen en la Figura 2.4, en la que queda constancia de la variedad en los grados de contraste y en la textura de las diferentes regiones. 17 Tabla 2.1 Volúmenes de la base de datos del USC-SIPI. Nombre del volumen Temática de las imágenes Textures Texturas variadas y mosaicos Aerials Imágenes aéreas tomadas desde gran altura Miscellaneous Miscelánea Sequences Fotogramas de filmaciones en movimiento (a) (b) (c) (d) Figura 2.4 Ejemplos de imágenes de la base de datos del USC-SIPI. (a) Cameraman, (b) Lena, (c) Bárbara y (d) Mandrill. Este lote de imágenes se ha completado con una serie de fotografías, de características similares a las de la base de datos del SIPI, captadas expresamente para este trabajo con la finalidad de probar los filtros con imágenes naturales. 18 2.1.4 Imágenes médicas La tercera y última clase la constituyen las imágenes médicas. Tener acceso a este tipo de imágenes no es fácil, debido principalmente a la información sensible que contienen y al celo con el que debe aplicarse la legislación en lo referente a la protección de datos del paciente. Por este motivo lo frecuente es encontrar imágenes aisladas o colecciones con un reducido número de pruebas, pobremente documentadas, esto es, sin información acerca del tipo de paciente, la modalidad de captura concreta ni el dispositivo empleado y su configuración. Sin embargo, existen proyectos de investigación que publican y comparten los datos objeto de su estudio en forma de colecciones bien documentadas. El portal The Cancer Imaging Archive (TCIA) [63], centraliza y facilita el acceso a las bases de datos de varios de estos proyectos, proporcionando no sólo una ingente cantidad de imágenes médicas sino también otros datos complementarios accesibles desde la propia web del TCIA o a través de los repositorios de cada proyecto particular. En la Figura 2.5 se muestra el aspecto de la página de entrada al portal y el menú de acceso las diferentes colecciones. Figura 2.5 Portal de The Cancer Imaging Archive (TCIA). 19 Tabla 2.2 Características de las colecciones del TCIA incluidas en la base de datos. Colección Zona anatómica Modalidades Nº de pacientes Breast Diagnosis Mama MR, CT, MG 88 TCGA-GBM Cerebro MR, CT, DX 262 QIN-Head Neck Cabeza y cuello PT, CT, SR, SEG 156 RIDER Lung PET-CT Pulmón PT, CT 244 TCGA-PRAD Próstata CT, PT, MR 7 TCGA-LIHC Hígado MR, CT, PT 65 REMBRANDT Cerebro MR 130 De los 60 proyectos actualmente accesibles desde el portal, se han escogido para su inclusión en la base de datos de este trabajo imágenes de las colecciones que aparecen en la Tabla 2.2. En la selección de las colecciones ha primado que estuviera representado el mayor número posible de modalidades de imagen (rayos X, tomografía computerizada, resonancia magnética, etc.) y, a la vez, un número suficiente de pacientes. En la Figura 2.6 se muestran ejemplos de imágenes de las diferentes modalidades consideradas. Toda la información se comparte siguiendo el protocolo DICOM (Digital Imaging and Communication in Medicine) [64], un estándar ampliamente extendido que establece un formato de fichero y un protocolo de comunicación de red basado en TCP/IP, normalizando de este modo la gestión, visualización y el intercambio de pruebas médicas. El protocolo permite incluir información adicional (metadatos) en la imagen y de entre ellos y algunos de ellos como las resoluciones en cada dimensión son imprescindibles para la interpretación geométrica de las imágenes médicas. MATLAB incorpora rutinas que permiten gestionar este tipo de archivo y en especial cargar sus muestras como matrices y acceder a los metadatos mencionados. Sus interfaces se muestran a continuación. dicomImage = dicomread(inputDicomFilename) imageInfo = dicominfo(inputDicomFilename) 20 (a) (b) (c) (d) (e) (f) Figura 2.6 Ejemplos de imágenes de la base de datos TCIA de distintas modalidades médicas: (a) MR, resonancia magnética, (b) CR, radiografía computerizada, (c) CT, tomografía computerizada, (d) US, ultrasonidos (e) PET, tomografía por emisión de positrones y (f) NM, medicina nuclear. 21 2.2 Preprocesamiento de las imágenes La base de datos contiene imágenes de una dimensión arbitraria D, normalmente 2 ó 3, que representan exclusivamente el nivel de una determinada magnitud 𝑓, también denominada intensidad. Se trata, por tanto, de imágenes en escala de grises, matemáticamente equivalentes a un campo escalar definido sobre una rejilla regular del espacio de dimensión D, de coordenadas 𝒙=(𝑥1,𝑥2, …, 𝑥𝐷). Admiten pues las dos representaciones de la Figura 2.7 y esta equivalencia justifica que, en adelante, nos refiramos a ellas indistintamente como imágenes o señales. (a) (b) Figura 2.7 Dos representaciones admisibles para una imagen 2D: (a) convencional, en escala de grises y (b) como superficie sobre un dominio bidimensional. De este modo, las imágenes en color, que contienen varias componentes y son, por tanto, matemáticamente equivalentes a campos vectoriales, deben ser transformadas antes de su inclusión en la base de datos del proyecto. La función rgb2gray de MATLAB reduce las tres componentes (roja, verde y azul) de la imagen en color RGB, a la única componente de la señal en escala de grises. El resultado de su aplicación a la imagen Mandrill se muestra en la Figura 2.8. imageGray = rgb2gray(imageRGB) 22 (a) (b) Figura 2.8 Conversión de la imagen en color a escala de grises. (a) Imagen en color con las tres componentes RGB y (b) su conversión a escala de grises. Es asimismo frecuente que las imágenes no presenten valores continuos de magnitud, sino que éstos aparezcan discretizados por un valor concreto, el cuanto, determinado en general por el número de bits empleados para la codificación de cada muestra. También suele fijarse el rango de variación permitido para la intensidad, siendo dos los intervalos más ampliamente utilizados como referencia, [0, 1] y [0, 256]. Las imágenes que hayan sido definidas atendiendo a otro rango de variación, pueden normalizarse al intervalo de referencia escalando linealmente, tal como se representa en la Figura 2.9. En cualquier caso, los algoritmos propuestos en este trabajo no imponen estas limitaciones de carácter discreto y valor acotado, por lo que la cuantificación o normalización de las señales queda restringida a aquellos experimentos en los que sea estrictamente necesaria con fines de comparación con otros algoritmos que sí requieran estas propiedades. (a) (b) Figura 2.9 FDP de la distribución de niveles para la imagen Mandrill normalizada a dos intervalos de referencia distintos. (a) Intervalo [0, 256] y (b) intervalo [0, 1]. 23 2.2.1 Formato de almacenamiento de las imágenes Por todo lo mencionado anteriormente, y con el fin de homogeneizar la información procedente de fuentes tan dispares y codificadas de distinta manera, se ha decidido establecer un formato de archivo común para todas las imágenes incorporadas en la base de datos, al que nos referiremos como formato IMD (Image Data) para el que se ha reservado la extensión .imd. Se ha optado por un formato en texto plano, sacrificando la velocidad de procesamiento y la capacidad de compresión que proporciona la codificación binaria, en favor de una mayor versatilidad y simplicidad en la gestión de los datos de los archivos de texto, directamente editables por parte del usuario. Esta característica es deseable en las tareas de investigación para facilitar la inspección, creación y modificación de las imágenes y la interoperabilidad de las distintas herramientas que potencialmente pudieran emplearse en este trabajo o en futuras extensiones. En la Figura 2.10, se esquematiza la estructura general del formato propuesto, describiéndose en la Tabla 2.3 el significado de cada uno de sus campos. Para la gestión de este formato de archivos se han implementado el conversor imdConverter y las rutinas saveImage y loadImage, para el almacenamiento en disco y su posterior recuperación, respectivamente. Las interfaces de estas rutinas se muestran a continuación. Para una descripción más detallada se remite al lector a la Parte II, Planos y Programas, de este documento. imageImd = imdConverter(imageNonImd) saveImage(imageImd, outputFilename) imageImd = loadImage(inputFilename) La base de datos resultante contiene un total de 495 imágenes, con un tamaño de 416 MB. 30 recibe este nombre porque la señal ruidosa adopta el aspecto de los patrones característicos de color en la película fotográfica de haluro de plata. Por su parte, el ruido speckle es característico de sistemas coherentes, como los de ultrasonidos o los radares de apertura sintética, y se manifiesta como un patrón de interferencias que viene provocado por las sumas constructivas y destructivas de reflexiones difusas que se producen en el objeto observado [52] [59], alternándose, respectivamente, puntos brillantes con píxeles oscuros. (a) (b) (c) (d) Figura 3.3 Ruido multiplicativo de media 𝝁𝝃 nula y varianza 𝝈𝝃𝟐= 𝟎.𝟏𝟑. Imágenes ruidosas (a) Cameraman y (b) Barbara, y las respectivas señales de ruido (c) y (d). 31 3.2.3 Modelo generalizado de ruido dependiente de la señal de SelvaAlparone Con el fin de obtener un modelo más completo del ruido dependiente de la imagen, capaz de describir el comportamiento de la mayoría de sistemas de adquisición, M. Selva y L. Alparone [7], generalizaron el caso multiplicativo proponiendo en 2009 el siguiente modelo paramétrico: 𝑔(𝒙)=𝑓(𝒙)+ 𝑓(𝒙)𝛾 ∙ 𝜉(𝒙)+𝜂(𝒙)=𝑓(𝒙)+ 𝑣(𝒙)+𝑤(𝒙) (3.5) (a) (b) (c) (d) Figura 3.4 Ruido según el modelo generalizado con parámetros 𝝁𝝃 = 0, 𝝈𝝃 = 0.12, 𝝁𝜼 = 0 y 𝝈𝜼 = 3.6. Imágenes ruidosas (a) Cameraman y (b) Barbara, y las respectivas señales de ruido (c) y (d). 32 en el que 𝑓 es la imagen libre de ruido, modelada como un proceso correlado y no estacionario, 𝜉 es un proceso aleatorio incorrelado de media nula y varianza 𝜎𝜉2, independiente de 𝑓, que modela la componente multiplicativa del ruido, y 𝜂 modela el ruido electrónico aditivo considerado como un proceso aleatorio gaussiano de media nula y varianza 𝜎𝜂2, independiente de los dos anteriores. Anulando 𝑣 o 𝑤 se obtienen, respectivamente, los modelos de ruido aditivo y multiplicativo. 3.2.4 Otros tipos de ruido Al margen de los descritos anteriormente, existe una gran cantidad de tipos de ruido en imágenes. Algunos ejemplos notables son los siguientes:  Ruido sal y pimienta: la señal ruidosa toma valores extremos en píxeles distribuidos aleatoriamente en la imagen. Su presencia desvirtúa totalmente el valor de la imagen en ese píxel y lo hace irrecuperable y sólo queda decidir si es ruido o no y sustituir su valor por otro acorde a píxeles adyacentes, potencialmente también ruidosos. (a) (b) Figura 3.5 Ruido de sal y pimienta. (a) Imagen ruidosa y (b) imagen filtrada.  Artefactos por movimiento: provocan el difuminado de la imagen si el movimiento es aleatorio, o las llamadas réplicas fantasma (ghost) si es periódico, presentando el ruido en este último caso una fuerte correlación con la propia imagen y tomando la forma de una versión desplazada y atenuada de una región de la imagen original. Puede reducirse con estrategias de decorrelación de señales. 33 (a) (b) Figura 3.6 Ruido por réplicas fantasma en una imagen de resonancia magnética (a) antes y (b) después del filtrado.  Ruido periódico: la componente de ruido puede seguir siendo independiente de la intensidad de la imagen, pero sí presenta una fuerte correlación consigo misma. Suele tratarse con algoritmos específicos que lo atenúan operando directamente en el dominio de la frecuencia. (a) (b) (c) (d) Figura 3.7 Ruido periódico (a) señal ruidosa, con espectro representado en (c); (d) supresión de las componentes espectrales dominantes del ruido y (b) imagen filtrada en el dominio espacial, con una reducción apreciable del nivel de ruido. 34 En cualquier caso, para todos ellos existen algoritmos específicos capaces de reducir sus efectos con mayor eficacia que integrándolos como parte de los procesos aleatorios en los modelos de ruido anteriores. Por este motivo, no se considerarán como objeto de estudio de este Proyecto Fin de Carrera, recomendándose en cada caso un preprocesamiento con el filtro específico apropiado. 3.3 Ruido en aplicaciones médicas 3.3.1 Ruido en imágenes de rayos X El ruido en la obtención de imágenes de rayos X presenta dos fuentes principales: ruido cuántico y ruido térmico. La parte dominante es el ruido cuántico, que proviene de fluctuaciones en el número de fotones en la trayectoria de los rayos (los que salen de la fuente, los que atraviesan el cuerpo observado, los que alcanzan el detector y los que excitan la película). Matemáticamente se modela como una componente de ruido multiplicativo directamente proporcional a la intensidad de la imagen, en la que el proceso aleatorio sigue la distribución de Poisson característica de estos fenómenos de naturaleza cuántica. Por otro lado, también debe considerarse el ruido térmico de los circuitos electrónicos del sistema de adquisición, que admite una modelización como ruido incorrelado gaussiano de media nula independiente de la imagen. La expresión del modelo completo de ruido para imágenes de rayos X queda de la siguiente manera: 𝑀(𝒙)=𝐴(𝒙)+ 𝐴(𝒙) ∙ 𝑢(𝒙)+𝑤(𝒙) (3.6) donde 𝑀 es la magnitud en la imagen observada, 𝐴 es la imagen ideal, libre de ruido, 𝑢 es un proceso aleatorio incorrelado e independiente de 𝐴, cuyos niveles siguen una distribución de Poisson de parámetro 𝜆𝑢, y 𝑤 es un proceso incorrelado e independiente de 𝐴 y 𝑢, que sigue una distribución gaussiana de media nula y varianza 𝜎𝑤2. Algunos estudios [35] defienden un modelo más elaborado del término multiplicativo, sustituyendo la distribución de Poisson por una Poisson compuesta (superposición de un número de distribuciones de Poisson, siendo este número un proceso 35 (a) (b) Figura 3.8 Ruido en imagen de rayos X. (a) Imagen ideal y (b) imagen ruidosa. aleatorio). El modelo resultante no constituye una particularización del general propuesto por Selva-Alparone, pero admite la siguiente aproximación definiendo los valores equivalentes de 𝛾 y 𝑢: 𝑀(𝒙)=𝐴(𝒙)+𝐴(𝒙)𝛾𝑒𝑞 ∙ 𝑢𝑒𝑞(𝒙) (3.7) 3.3.2 Ruido en imágenes de ultrasonidos En imágenes obtenidas mediante ultrasonidos, el ruido dominante es el speckle, producido por el patrón de interferencias de las ondas reflejadas en los tejidos y otros elementos de dimensiones comparables a la longitud de onda de trabajo. En la literatura se diferencia entre dos fuentes de ruido. La primera de ellas es la presencia de un número considerable de elementos dispersos que provocan reflexiones difusas de la señal transmitida por el sistema de adquisición, como ocurre con los glóbulos de la sangre. Por otro lado también se genera ruido speckle cuando la señal es reflejada de forma difusa por los propios tejidos [23]. En ambos casos, el nivel de ruido multiplicativo es muy superior al de ruido aditivo, por lo que este último puede despreciarse. En esa situación, estableciendo los valores equivalentes de 𝛾 y 𝑢, el modelo de ruido para imágenes de ultrasonidos queda definido por la siguiente expresión: 36 (a) (b) Figura 3.9 Ruido en imágenes de ultrasonidos. (a) Imagen original y (b) imagen ruidosa. 𝑀(𝒙)=𝐴(𝒙)+𝐴(𝒙)𝛾𝑒𝑞 ∙ 𝑢𝑒𝑞(𝒙) (3.8) que coincide formalmente con la obtenida para el caso de ruido en imágenes de rayos X. 3.3.3 Ruido en imágenes de resonancia magnética Los sistemas modernos de adquisición de imágenes por resonancia magnética cuentan con una o varias bobinas. En ellas se detecta el débil campo generado por la rotación de los momentos magnéticos angulares de los núcleos atómicos cuando el paciente, o el objeto de análisis, es sometido a un determinado campo magnético externo, tal como se muestra en la Figura 3.10. Existe cierto consenso en que es precisamente el paciente la fuente principal de ruido en esta modalidad de imagen [6] [10] [24], debido a movimientos propios de funciones fisiológicas como la respiración o el proceso digestivo. En lo meramente tecnológico, otras fuentes contaminantes significativas son el ruido térmico en la entrada de radiofrecuencia de las bobinas, el acoplamiento en ellas de las corrientes parásitas de Foucault y el ruido electrónico producido por los circuitos de la cadena de recepción. 37 (a) (b) Figura 3.10 Adquisición de imágenes de resonancia magnética. (a) Ubicación de las bobinas en el escáner y (b) alineación de los momentos magnéticos nucleares con el campo magnético exterior. El proceso de adquisición y obtención de la imagen se describe en la Figura 3.11. Los datos obtenidos en los detectores de cada bobina se corresponden matemáticamente con la transformada de Fourier de la imagen, definidas en el dominio espectral, que en este ámbito se conoce como espacio k. Aplicando la transformada inversa de Fourier se obtienen las correspondientes señales en el dominio espacial, que finalmente se combinan para componer la imagen definitiva. Figura 3.11 Etapas del procesamiento de la señal en imágenes MR. De izquierda a derecha, señal detectada en cada bobina en el espacio k, correspondiente señal el espacio físico y combinación de los resultados de las bobinas para obtener la imagen. 38 (a) (b) Figura 3.12 Ruido en imágenes de resonancia magnética. (a) Imagen original y (b) imagen ruidosa. En el modelo de ruido simplificado se considera que cada componente de la transformada (parte real y parte imaginaria) está contaminada por ruido normal gaussiano de media nula y una determinada varianza, 𝜎𝐾2, que se supone igual para todas las bobinas. Tras aplicar la transformada inversa, la imagen obtenida presenta también contaminación con ruido normal gaussiano de media nula en cada componente, con el valor de varianza equivalente, 𝜎2, para el dominio espacial. Los valores de intensidad en la imagen final se corresponden con la magnitud de estas señales complejas, siendo su expresión matemática: 𝑀(𝒙)=√ [𝐴(𝒙)+𝑛𝑟(𝒙)]2+ [𝑛𝑖(𝒙)]2 (3.9) La existencia de estas componentes gaussianas en cuadratura implica que el ruido presente en el valor de la magnitud (𝑀) sigue una distribución riciana, dependiente además del valor de intensidad ideal libre de ruido (𝐴), es decir, del sujeto observado. Este comportamiento ha sido ampliamente descrito en la literatura [6] [47] haciendo hincapié en la dificultad que representa trabajar con este tipo de distribuciones y la complejidad de recuperar con ciertas garantías la señal original cuando está sujeta a ruido de esta naturaleza. Por este motivo, a la hora de abordar de forma práctica la estimación del parámetro 𝜎 suele recurrirse a modelos aún más simplificados como el ruido aditivo gaussiano, considerando semejante el comportamiento de ambos modelos cuando el nivel de ruido es bajo. Sin embargo, esta hipótesis deja de ser cierta en imágenes con alto 39 contraste en las que la disparidad en el valor de intensidad de los píxeles es significativa y se refleja en la parte multiplicativa de la señal de ruido que deja de ser despreciable frente a la aditiva, poniendo de manifiesto la necesidad de modelos de ruido más precisos como el que se presenta en el siguiente apartado. 3.4 Modelo de ruido propuesto El modelo propuesto para el ruido en imágenes, aplicable tanto en un contexto general como particularizado al ámbito de las modalidades médicas, es el siguiente 𝑔(𝒙)=𝑓(𝒙)+ ∑[𝑓(𝒙)𝛾𝑟· 𝜉𝑟(𝒙)] 𝑅 𝑟=1 (3.10) en el que 𝑔 es la imagen ruidosa, 𝑓 es la imagen libre de ruido, modelada como un proceso correlado y no estacionario y R es el número de términos de ruido. Para cada uno de esos términos 𝛾𝑟 es el exponente que determina la relación entre las magnitudes de imagen y ruido, y 𝜉𝑟 es un proceso aleatorio incorrelado e independiente de 𝑓 y del resto de procesos aleatorios. Conceptualmente, no se diferencia mucho del modelo para ruido dependiente de la señal propuesto por Selva-Alparone, en el que está basado. Sin embargo, desde el punto de vista matemático, su expresión es más compacta, debido a que todos los términos de ruido responden a la misma expresión. Asimismo es más general en el sentido de que no fija ni impone restricciones en las funciones de distribución que gobiernan dichos términos. Estas dos propiedades son deseables desde el punto de vista de la automatización de los algoritmos para la estimación de los parámetros del modelo, simplificando su implementación y reduciendo la necesidad de intervención por parte del usuario. Como muestra de flexibilidad que proporciona este modelo, se aplica a continuación a dos de los ejemplos presentados anteriormente que no tenían encaje en el modelo generalizado de ruido dependiente de la señal. 46 (a) (b) (c) (d) Figura 4.2 Compromiso entre reducción de ruido y conservación de detalles para filtros de promedios. (a) Imagen original, (b) versión ruidosa, (c) MF con entorno de radio 5 y (d) MF con entorno de radio 12. 4.2.1.2 Filtros lineales de convolución Como generalización de los filtros de media aritmética puede considerarse un subconjunto de filtros lineales (LF, Linear Filters) de convolución, también llamados filtros de máscara. En este caso, se sustituye el promedio por la suma ponderada de valores de intensidad de la ecuación (4.2), en la que los coeficientes 𝑤𝑖 se corresponden con el peso relativo de la intensidad del píxel vecino 𝒙𝑖 en el cálculo de la salida del filtro para la posición 𝒙. 47 𝑦(𝒙) =𝐿𝐹𝜌 𝑓(𝒙)= ∑ [𝑤𝑖·𝑓(𝒙𝑖)] 𝒙𝑖 ∈ 𝐵(𝒙,𝜌) (4.2) Con el fin de que la suma resultante sea representativa del verdadero valor de intensidad en la imagen libre de ruido, los coeficientes 𝑤𝑖 deben cumplir dos propiedades: 1.) 0≤ 𝑤𝑖≤1 (4.3) 2.) ∑ 𝑤𝑖 = 1 𝒙𝑖 ∈ 𝐵(𝒙,𝜌) (4.4) Estos filtros quedan totalmente determinados por su máscara de convolución, en la que los elementos son precisamente los coeficientes de los respectivos píxeles del entorno. Los coeficientes, así como la máscara, son constantes e independientes del píxel procesado por lo que no se modifican en función de los distintos valores de intensidad presentes en la imagen. En la Figura 4.3, se representan tres máscaras distintas. En la primera de ellas los elementos de la máscara tienen el mismo valor, por lo que la expresión de la suma ponderada se reduce al caso particular de media aritmética. Las dos restantes se corresponden con máscaras de filtros en los que se han considerado entornos cuadrados (b) o circulares (c) y con un valor de los coeficientes que decae según el cuadrado de la distancia al píxel central. (a) (b) (c) Figura 4.3 Máscaras de filtros de convolución. (a) Filtro promedio, (b) ponderado con máscara cuadrada y (c) ponderado con máscara circular. 48 Es habitual establecer una ponderación decreciente con la distancia radial al píxel central, empleándose pesos inversamente proporcionales al cuadrado de la distancia o en forma de campana gaussiana de anchura prefijada. En la ecuación (4.5) se particulariza la expresión general de filtros lineales de la ecuación (4.2) para una ponderación gaussiana en la que el parámetro 𝑏 permite ajustar dicha anchura. Disminuyendo el valor de este parámetro se reducen los pesos relativos de los píxeles más alejados del central, obteniendo un promedio más local. 𝑦(𝒙) = ∑ [(𝑒− ||𝒙𝑖 − 𝒙||2 𝑏2 𝐶)·𝑓(𝒙𝑖)] 𝒙𝑖 ∈ 𝐵(𝒙,𝜌) (4.5) En esta expresión 𝐶 es el factor de normalización que asegura que los pesos cumplan las propiedades (4.3) y (4.4). Dada la regularidad de la rejilla sobre la que se define la imagen, este coeficiente es independiente del píxel central 𝒙 considerado, y su valor es: 𝐶= ∑ 𝑒− ||𝒙𝑖 − 𝒙||2 𝑏2 𝒙𝑖 ∈ 𝐵(𝒙,𝜌) (4.6) (a) (b) (c) Figura 4.4 Filtro lineal con ponderación gaussiana. (a) Imagen original, (b) imagen ruidosa y (c) salida del filtro lineal de parámetros ρ = 5 y b = 10. 49 En la Figura 4.4 se representa la salida del filtro lineal de radio 𝜌 = 5 con ponderación gaussiana de parámetro b = 10. La reducción del ruido es comparable a la obtenida para filtros de promedio, pero con el filtro lineal el efecto de difuminado de los contornos y la pérdida de detalles es mucho menor. 4.2.1.3 Filtro de mediana Como alternativa a los filtros lineales y con el fin de reducir el efecto negativo que éstos tienen sobre los bordes de la imagen, surgen los filtros de mediana. El algoritmo es similar a los presentados en los apartados anteriores, considerando el mismo conjunto de píxeles del entorno. Sin embargo, en esta nueva estrategia no se ponderan ni promedian los valores de intensidad, sino que se escoge directamente la mediana del conjunto, lo que equivale a ordenar los valores de intensidades y tomar el que ocupa la posición central. Si, tal como se comentó en el apartado de filtros de promedios, se consideran entornos simétricos centrados en el píxel, el número de valores es impar y por tanto no hay ambigüedad en la elección de la mediana. Nótese que si éste no fuera el caso y el número de píxeles fuese par, la mediana podría definirse bien como uno de los dos valores centrales o como la media aritmética de ambos, en cuyo caso se recomienda la primera de las soluciones. El motivo de esta elección es la conveniencia de que el valor de intensidad resultante sea uno de los ya presentes en la imagen, evitando así el suavizado de los bordes. (a) (b) (c) Figura 4.5 Comparativa de los filtros de convolución (b) y de mediana (c) en respuesta a la señal ruidosa (a). 50 La Figura 4.5 permite comparar las respuestas del filtro lineal de convolución y el filtro de mediana. Con los tamaños de entorno utilizados en aplicaciones reales, la capacidad del filtro de mediana para reducir el ruido es comparable a la de los filtros de máscara, a la vez que evita el difuminado de la imagen, respetando los bordes y conservando en la medida de lo posible los detalles. Sin embargo, si bien la respuesta a saltos abruptos en el valor de intensidad es apropiada, el filtro presenta muy mal comportamiento en presencia de regiones angulosas, filamentosas o con contornos no redondeados. En la Figura 4.6 se observan las similitudes en la respuesta del filtro ante un elemento de contorno circular y otro de contorno cuadrado de semejantes proporciones e intensidades. La supresión de las esquinas resulta evidente para este segundo elemento y constituye el mayor inconveniente de los filtros de mediana, lo que hace desaconsejable su utilización en aplicaciones en las que se requiere respetar fielmente los contornos presentes en la imagen, como sucede en el ámbito del diagnóstico oncológico por imagen [47]. (a) (b) Figura 4.6 Comportamiento del filtro de mediana en contornos angulosos: (a) imagen ruidosa y (b) respuesta del filtro de mediana. 4.2.1.4 Filtro de Vecindad de Yaroslavsky (YNF) Esta nueva clase de filtros ha sido desarrollada paralelamente por L. Yaroslavsky, con el nombre de filtros de vecindad de Yaroslavsky (YNF, Yaroslavsky Neighbourhood Filter) [57], y por J. S. Lee, bajo la denominación de Filtros Sigma [32]. En adelante nos referiremos a ellos como YNF, ya que es el término más extendido. A semejanza de los filtros lineales, en los YNF la intensidad resultante se obtiene como una suma ponderada 51 de intensidades de los píxeles cercanos pero, a diferencia de aquéllos, los valores de los coeficientes sí dependen del píxel considerado y de los valores de intensidades implicados en la suma. Su respuesta viene dada por: 𝑦(𝒙) =𝑌𝑁𝐹ℎ,𝜌𝑓(𝒙)= 1 𝐶(𝒙)∑ [𝑒− |𝑓(𝒙𝑖)− 𝑓(𝒙)|2 ℎ2 𝑓(𝒙𝑖)] 𝒙𝑖 ∈ 𝐵(𝒙,𝜌) (4.7) en la que 𝜌 y h, permiten controlar respectivamente el tamaño del entorno y la sensibilidad de los pesos con la diferencia en los valores de intensidad. El factor de normalización 𝐶(𝒙) depende en este caso del píxel central 𝒙 y las intensidades de los píxeles del entorno, siendo su valor: 𝐶(𝒙)= ∑ 𝑒− |𝑓(𝒙𝑖)− 𝑓(𝒙)|2 ℎ2 𝒙𝑖 ∈ 𝐵(𝒙,𝜌) (4.8) El filtro YNF introduce un cambio importante en el concepto de cercanía entre píxeles ya que a la hora de calcular el peso relativo de los píxeles adyacentes se consideran dos criterios: proximidad geométrica y proximidad en el valor de intensidad. La proximidad geométrica se impone, ya que se considera nulo el coeficiente para todos los píxeles no contenidos en el entorno 𝐵(𝒙,𝜌). Ajustando el parámetro 𝜌 se regula el tamaño de este entorno. Por otro lado, la proximidad en valores de intensidad se fuerza indirectamente, debido a que la expresión de los coeficientes penaliza a los píxeles vecinos cuya intensidad se aleja de la del píxel central. El parámetro h del filtro permite controlar la sensibilidad de los coeficientes a diferencias de intensidades, penalizando tanto más cuanto menor sea el valor del parámetro. Empleando esta estrategia se reduce la probabilidad de que el valor de intensidad resultante para un píxel diste mucho de su valor inicial, disminuyendo el efecto de suavizado de bordes y la pérdida de detalles respecto de los filtros lineales. Por el mismo motivo, se evita la supresión de contornos angulosos de las regiones presentes en la imagen, propia de los filtros de mediana, ya que los coeficientes para píxeles que no comparten la misma región homogénea con el píxel central son prácticamente despreciables frente a los que sí la comparten. En la Figura 4.7 se ilustra esta propiedad, resaltando en oscuro las posiciones con coeficientes dominantes para cada escenario. 52 (a) (b) Figura 4.7 Distribución de los valores relativos de los pesos en el YNF en dos escenarios: (a) contorno curvo y (b) promediando en regiones no conexas. La región oscura es la dominante. En la Figura 4.8 se muestran de manera conjunta los resultados del filtrado con YNF y los obtenidos empleando el filtro de mediana. En lo referente a la reducción de ruido, el filtro de mediana tiene un comportamiento ligeramente mejor que el YNF en regiones de intensidad prácticamente constante ya que este último replica en ocasiones la propia variación del ruido. Esto se debe a la naturaleza aleatoria de la señal ruidosa, que provoca que en determinados entornos la intensidad del píxel central tome valores extremos de esa región, de forma que las ponderaciones de los píxeles vecinos sean bajas y, por tanto, estos tengan poco efecto en el valor devuelto por el filtro para ese píxel central. Sin embargo, la notable capacidad de los filtros YNF para conservar los detalles de la imagen y en especial para respetar los contornos de las regiones, los hace recomendables frente a todos los estudiados anteriormente. (a) (b) (c) Figura 4.8 Comparativa de la respuesta de los filtros de mediana y YNF: (a) imagen ruidosa, y la respuesta de los filtros (b) de mediana y (c) YNF. 53 4.2.1.5 Filtros SUSAN y Bilaterales Los filtros YNF son menos conocidos que versiones que evolucionaron a partir de ellos como son los filtros SUSAN [51] (también referidos como filtros de vecindad SUSAN (SNF, SUSAN Neighbourhood Filters) y los filtros bilaterales [54]. En ambos algoritmos, se sustituye la restricción de proximidad geométrica impuesta al definir el entorno 𝐵(𝒙,𝜌), por una ponderación en función de la distancia radial al píxel central, que se añade a la ya existente debida a la diferencia de intensidades. Para la ponderación en distancia suele escogerse una función gaussiana, tal como se recoge en las ecuaciones (4.9) y (4.10), en las que el parámetro 𝜌 determina el ancho de la campana y 𝐶(𝒙) es el nuevo factor de normalización. 𝑦(𝒙) =𝑆𝑁𝐹ℎ,𝜌𝑓(𝒙)= 1 𝐶(𝒙)∑[𝑒− ||𝒙𝑖 − 𝒙||2 𝜌2𝑒− |𝑓(𝒙𝑖)− 𝑓(𝒙)| ℎ2 𝑓(𝒙𝑖)] 𝒙𝑖 (4.9) 𝐶(𝒙)=∑[𝑒− ||𝒙𝑖 − 𝒙||2 𝜌2𝑒− |𝑓(𝒙𝑖)− 𝑓(𝒙)|2 ℎ2] 𝒙𝑖 (4.10) De esta manera, para obtener el valor de intensidad de cada píxel debe calcularse una suma ponderada que se extiende, al menos de manera teórica, a todos los píxeles de la imagen, lo que lo convierte en un método muy costoso desde el punto de vista computacional. En cuanto a resultados, el comportamiento es en esencia el mismo que el de los filtros YNF, evitando de igual modo el difuminado de los bordes. Cuantitativamente, la sustitución del promediado original por el gaussiano proporciona una ligera mejora en la reducción del ruido. 4.2.1.6 Filtros de promedios no locales Los filtros de promedios no locales (NLM, Non-local Means) pueden considerarse como una versión de los filtros de promedios SUSAN que redefinen la proximidad entre píxeles incorporando el concepto de autosimilitud. La autosimilitud hace referencia al parecido existente entre pequeñas regiones de la imagen, potencialmente distantes entre sí, que presentan valores parecidos en su intensidad y disposición espacial. Están especialmente pensados para situaciones en las que existe un alto grado de redundancia en 54 Figura 4.9 Identificación de regiones autosimilares. la imagen, como ocurre en la Figura 4.9, con las regiones similares enmarcadas en recuadros del mismo color. La respuesta del filtro NLM viene dada por las expresiones (4.11) y (4.12) 𝑦(𝒙) =𝑁𝐿𝑀ℎ 𝑓(𝒙)= 1 𝐶(𝒙)∑[𝑒− 𝑑𝑓2(𝑓𝒙𝑖, 𝑓𝒙) ℎ2 ·𝑓(𝒙𝑖)] 𝒙𝑖 (4.11) 𝑑𝑓2(𝑓𝒙𝑖, 𝑓𝒙)= ‖𝑓𝒙𝑖− 𝑓𝒙‖𝐺𝑏 2 (4.12) en las que 𝑑𝑓2(𝑓𝒙𝑖, 𝑓𝒙) es una medida de la similitud entre los entornos 𝑓𝒙 𝑦 𝑓𝒙𝑖, centrados respectivamente en el punto 𝒙, donde se evalúa la respuesta, y otro punto de la imagen 𝒙𝑖, arbitrariamente distante del primero. Matemáticamente se calcula como la norma L2 de la diferencia entre valores de intensidad de los elementos de ambos entornos, aplicando previamente una ponderación gaussiana de parámetro b, en la manera descrita en el apartado dedicado a los filtros lineales de convolución. El coeficiente de normalización 𝐶(𝒙) para este tipo de filtros toma el valor: 𝐶(𝒙)=∑𝑒− 𝑑𝑓2(𝑓𝒙𝑖, 𝑓𝒙) ℎ2 𝒙𝑖 55 Experimentalmente se demuestra que, salvo para imágenes con una marcada autosimilitud, la respuesta de los filtros de promedios no locales no difiere mucho de la obtenida con los filtros SUSAN o YNF, mientras que la complejidad computacional es sensiblemente mayor, incrementándose de forma cuadrática con el número de píxeles de la imagen. 4.2.1.7 Filtros basados en ecuaciones en derivadas parciales Otra familia de métodos aplicados en el ámbito del filtrado de imágenes es la que surge del análisis de comportamiento del sistema visual humano. Los estudios sobre el sistema de captación visual [46] y la capacidad de reconocimiento de los detalles en presencia de ruido, concluyeron que un modelo simplificado para reproducir este procesamiento es la aplicación directa de la ecuación de difusión. 𝜕Ψ 𝜕t=𝑑𝑖𝑣 (𝑔(|∇Ψ|)· ∇Ψ) expresión que, si se discretiza en el tiempo y se aplica un esquema en diferencias finitas para las derivadas espaciales, admite su interpretación como un filtro iterativo de restauración en el que t es el índice de iteración, Ψ es el campo escalar de los valores de intensidad de la imagen y 𝑔 una función monótona decreciente. Su funcionamiento se basa en suponer que el error se difunde por la imagen pero sólo entre regiones con niveles semejantes –es decir, gradiente bajo–, mientras que no se difunde en regiones donde el gradiente es alto como contornos y regiones con texturas muy marcadas. Ejemplos clásicos de este tipo de filtros son el de Perona-Malik [42] y el Histace-Rousseau [27]. Desde el punto de vista algorítmico podría clasificarse como un filtro de realce de bordes aunque la señal filtrada es también útil como aproximación de la imagen ideal no degradada. Sin embargo, esta reducción del nivel de ruido es básicamente cualitativa, es decir, la señal filtrada se percibe visualmente como una versión mejorada de la ruidosa en la que los contornos de las regiones y los detalles se aprecian con mayor nitidez. Sin embargo, el análisis cuantitativo de la señal resultante pone de manifiesto que la SNR decrece en amplias regiones de la imagen y, lo que es peor, el ruido presente en la imagen filtrada no conserva las propiedades estadísticas del ruido original. Por este motivo, los 62 Figura 4.13 Etapas del filtro propuesto. Finalmente cabe destacar el comportamiento de los filtros basados en ecuaciones en derivadas parciales, como el propuesto por Histance-Rousseau [27], en los que la restauración se modela como un fenómeno de difusión del ruido. En ellos, se establece un umbral en el gradiente de intensidades que limita la difusión en determinadas direcciones, lo que permite conservar y hasta realzar los contornos de la imagen incluso en presencia de potencias de ruido elevadas. Presentan la ventaja de homogeneizar los valores de SNR en la imagen pero, desafortunadamente, alteran completamente las características estadísticas del ruido. Esto los convierte en ideales para una mejora desde el punto de vista cualitativo, es decir, cuando las imágenes van a ser directamente observadas, pero en poco recomendables para el análisis automático del ruido con fines de caracterizarlo. Por todo ello se propone un filtrado basado en promedios en los que la selección de píxeles considerados atiende a dos criterios: la proximidad geométrica y una limitación en la variación de la intensidad. En lo relativo a la proximidad geométrica, se impondrá que las regiones seleccionadas como entorno de un píxel sean conexas, tal como se recomienda en estrategias semejantes como la de entornos crecientes ideada por Jiři Jan [29] o Li [34], o la de crecimiento adaptativo propuesta por Balafar [10]. Por otro lado, la limitación en la en la variación admisible de los valores de intensidad en píxel vecinos se lleva a cabo estimado las derivadas direccionales vistas desde píxel central y estableciendo un valor umbral para dichas derivadas. Los pasos del algoritmo, descritos en la Figura 4.13, se detallan en los siguientes apartados. 4.3.1 Definición de los entornos Para la obtención de la magnitud asociada a un píxel se consideran inicialmente los píxeles en su vecindad, de forma análoga a la empleada en los filtros de promedio. Debe tenerse en cuenta que la selección de la forma del entorno y su tamaño condicionan las 63 propiedades del filtro (isótropo o anisótropo, sesgado, lineal, etc.) y su respuesta. En principio, nada hace favorecer a priori ninguna dirección, motivo por el cual el entrono debe ser simétrico y centrado en el píxel analizado en cada momento. Se busca además que el conjunto de valores escogidos sea representativo de la región a la que pertenece el píxel y, por tanto, debe imperar un criterio de proximidad, lo que apunta hacia entornos que sigan criterios de distancia geométrica en detrimento de estrategias de análisis por bloques rectangulares. En estas últimas, pueden quedar excluidos píxeles más cercanos al central que otros píxeles que sí han sido seleccionados. Por todo lo anterior, el entorno escogido es una bola de radio 𝜌, considerando distancia euclídea. Un ejemplo de este tipo de entornos se representa en la Figura 4.14. En ella aparecen marcados en azul los píxeles vecinos al central, destacado en rojo, considerando un entorno de radio 𝜌 = 5. Nótese que en este caso la rejilla, aunque uniforme, tiene resoluciones diferentes para los distintos ejes. Figura 4.14 Píxeles vecinos (azul) en el entorno del píxel central (rojo). 4.3.2 Aproximación de las derivadas direccionales Una vez escogidos los píxeles del entorno, se analizan todas las direcciones definidas tomando como origen el píxel central y pasando por alguno de los píxeles vecinos. Cada una de las semirrectas asociadas a estas direcciones puede contener uno o varios píxeles vecinos pero, de existir varios, todos ellos –incluido el píxel central– aparecen equiespaciados en la semirrecta, debido a la regularidad de la rejilla sobre la que se define la imagen. En la Figura 4.15 se muestra el análisis de estas direcciones empleando un entorno circular de radio 5 sobre una rejilla con resoluciones distintas en los 64 Figura 4.15 Semirrectas definidas entre el píxel central y sus vecinos. dos ejes. El número de píxeles presentes en cada dirección varía desde 2 (semirrecta verde) hasta 5 (semirrecta roja). Sobre cada una de estas orientaciones, o visto de forma equivalente, sobre cada semirrecta definida por el procedimiento explicado anteriormente, se puede estimar la derivada direccional del valor de intensidad evaluada en el píxel central del entorno. Este valor de la derivada proporciona información sobre la variación local de la intensidad, que servirá de discriminante para diferenciar distintas regiones en la imagen. Para la estimación se emplea el esquema en diferencias finitas de la ecuación (4.14), en el que 𝑓󰆹′(𝒙𝟎) es la estimación considerando 𝒙𝟎 como píxel central, 𝒓 es el vector radial en la semirrecta estudiada, cuyo módulo 𝑟 coincide con la distancia entre píxeles de la misma, siendo 𝑖 𝑟 la distancia radial al píxel vecino 𝒙𝒊. Los coeficientes 𝑐𝑘 se obtienen mediante la resolución del sistema de ecuaciones correspondiente. 𝑓󰆹′(𝒙𝟎)= ∑𝑐𝑖 𝑓(𝒙𝒐+ 𝑖 𝒓) 𝑀 𝑖=1 = ∑𝑐𝑖 𝑓(𝒙𝒊) 𝑀 𝑖=1 (4.14) La aplicación de este esquema permite obtener tantas estimaciones como píxeles vecinos estén incluidos en la semirrecta. La primera estimación (M = 1) se obtiene considerando dos muestras: el píxel central y el vecino más cercano y es a éste último al que se asocia el valor resultante. Añadiendo progresivamente nuevas muestras, en orden creciente de distancia radial respecto del píxel central, se obtienen las aproximaciones 65 restantes, de forma que a cada píxel le corresponda un valor de la variación en intensidad en la orientación determinada por la semirrecta a la que pertenecen. 4.3.3 Selección del umbral en la variación de la intensidad Con el fin de determinar qué píxeles serán analizados conjuntamente para obtener el valor de intensidad filtrado correspondiente al píxel central, se establece una selección por umbral. Este umbral representa la máxima variación admisible entre píxeles cercanos para considerar que comparten una región homogénea con las mismas propiedades estadísticas. La determinación de su valor constituye el principal problema en la implementación del filtro, ya que el valor óptimo se considera en general dependiente de la aplicación y las características particulares de las imágenes consideradas. Con el fin de evitar la necesidad de intervención por parte del usuario, se fija este umbral como aquel valor en la estimación de derivadas direccionales bajo el cual se encuentra un porcentaje prefijado del total de estimaciones computadas para la imagen. Esta estrategia ha sido aplicada con éxito en algoritmos con necesidades semejantes como el propuesto por Olsen [41] o Tai [49] que requieren el establecimiento de un umbral para la varianza en el valor de intensidades. Una vez establecido el valor umbral, debe seleccionarse el criterio de cribado de los píxeles en función de su localización geométrica respecto del píxel central. En principio, podrían descartarse exclusivamente los píxeles para los cuales se obtuvo una variación estimada superior al umbral. Esta estrategia presenta el inconveniente de considerar regiones homogéneas que no son conexas. Para evitar estas situaciones se ha optado por un tratamiento recurriendo nuevamente al procesamiento por orientaciones. En concreto, para la semirrecta asociada a cada orientación se van comprobando progresivamente los píxeles en orden creciente de distancia radial respecto del píxel central. Si la estimación de la derivada correspondiente a un píxel no supera la máxima variación admisible, entonces se considera como píxel vecino del central. En caso contrario, ese píxel y todos los restantes de la misma semirrecta son descartados. Procediendo de igual forma con todas las semirrectas del entorno, se restringe el concepto de regiones homogéneas a dominios conexos. 66 4.3.4 Cálculo de los valores de intensidad El último paso consiste en asignar un valor de intensidad a la imagen filtrada integrada en la posición del píxel central, considerando que la región es relativamente homogénea y que un promedio de los valores de intensidades de todos los píxel del entorno proporciona más representativo de la región en su conjunto. En concreto, por su sencillez y buenas prestaciones, se ha optado por emplear la media aritmética. Procediendo de igual forma en cada uno de los píxeles de la imagen, considerándolos como píxeles centrales en sus respectivos entornos, y repitiendo el proceso de selección de vecinos potenciales, cribado y promediado, se completa la imagen filtrada. 4.4 Estructura de la imagen y prefiltrado Como se adelantó en la introducción de este capítulo, los filtros de restauración de imágenes pueden emplearse tanto para obtener directamente la imagen recuperada, como para realizar un prefiltrado que proporcione una primera aproximación de la imagen buscada –denominada estructura– y luego seguir el esquema de la Figura 4.1 (b) con el fin de obtener valores aproximados de los parámetros de ruido que servirán como entradas de posteriores etapas de filtrado. El error cometido al considerar como estimación del ruido la diferencia entre la señal ruidosa y la estructura de la imagen puede llegar a ser muy alta, especialmente en las regiones en las que la imagen ideal presenta cambios abruptos de intensidad. La estimación local de la varianza de este ruido estimado proporciona en dichas regiones valores que superan ampliamente el valor real. Se corre además el riesgo de extender esta sobreestimación a toda la imagen si se plantea una etapa de posprocesado para suavizar el campo de varianza obtenido. También debe tenerse en cuenta que la estructura extraída es en esencia una versión suavizada de la imagen ruidosa en la que se entremezclan el verdadero valor de la imagen 67 y la alteración provocada por el valor medio del ruido, cuando su media no sea nula, como ocurre en el caso de la modalidades médicas estudiadas (rayos X y resonancia magnética). En el Capítulo 6, dedicado a la caracterización del ruido, se analiza el procedimiento para obtener estimaciones de los parámetros del modelo de ruido, en el que éstos y otros problemas son tenidos en cuenta, detallando las soluciones adoptadas para solventarlos. 68 5. PROCESAMIENTO POR BLOQUES Todos los filtros analizados en el Capítulo 4, a pesar de operar de forma local, tienen en cuenta la imagen en su conjunto y proporcionan valores en la salida para todos los píxeles de la imagen. En el procesamiento por bloques el enfoque es diferente ya que, una vez dividida la imagen en bloques, y establecido un criterio de selección, se descartan todos los bloques que no lo cumplan y no se devuelve ningún valor para los píxeles que lo componen. En el caso concreto de caracterización del ruido, interesa realizar el análisis en aquellas regiones que se consideren suficientemente homogéneas y desechar los bloques en los que la variabilidad sea acusada, evitando de esta manera las regiones en las que se sabe de antemano que el error cometido por los estimadores será mayor. 69 5.1 Descomposición en bloques La definición de los bloques juega un papel crucial en esta estrategia, ya que las decisiones que se tomen en tamaño, forma y distribución de los bloques condicionarán las prestaciones de los filtros. Estas decisiones implican normalmente compromisos entre dos características deseables y habrá que primar las que se consideren más deseables en función de la aplicación concreta y los objetivos perseguidos. En las siguientes secciones se analizan estas cuestiones sobre la descomposición de la imagen en bloques. 5.1.1 Tamaño Una decisión crítica en el procesamiento de imágenes por bloques es la del tamaño de dichos bloques, especialmente cuando se trata de caracterizar estadísticamente fenómenos de naturaleza aleatoria como el ruido. Desde el punto de vista de la localidad, interesa definir bloques del menor tamaño posible, analizando en esa situación las relaciones entre valores de píxeles muy próximos entre sí. De esta manera se pueden detectar con mayor resolución los cambios en la imagen. Por otro lado, la exactitud de los estimadores, especialmente si se trata de estimadores consistentes –como suele ser habitual– está directamente relacionada con el número de muestras consideradas, de forma que interesa tomar bloques de tamaño considerable. El número de valores debe ser suficiente para representar de forma fidedigna el fenómeno aleatorio, evitando en la medida de lo posible que la existencia de valores aberrantes (outliers) alteren sensiblemente los valores estimados. En la Figura 5.1 se representa la esperanza del error relativo cometido al estimar la desviación estándar de un ruido blanco gaussiano de media nula, definido según la ecuación (5.1), en la que 𝜎 representa el valor estimado y 𝜎 el valor verdadero. Se comprueba que este valor es independiente de la desviación considerada y que su comportamiento es monótono decreciente respecto del número de muestras tenidas en cuenta. En concreto, se observa que son necesarias al menos 34 muestras para que el error quede por debajo del 10% y 131 si se fija el umbral en el 5%. 𝜀𝑟=𝜎−𝜎 𝜎 (5.1) 70 Figura 5.1 Esperanza del error relativo en la estimación de la desviación estándar en función del número de muestras disponibles. Este compromiso, conjuntamente con el coste computacional añadido que representa trabajar con bloques de mayor tamaño, determina las dimensiones finales de los bloques. Recurriendo a los estudios realizados al respecto [49], se comprueba que en la práctica está muy extendido el uso de bloques con un número relativamente elevado de píxeles (7x7) cuando se estiman estadísticos, mientras que se opta por otros mucho más ajustados al píxel central (3x3) para el cálculo de operaciones diferenciales. En este punto se recomienda un tamaño de bloque suficiente para mantenerse por debajo del primer umbral mencionado anteriormente, 𝜀𝑟 = 10%, es decir, considerar al menos 35 vecinos. 5.1.2 Forma En lo referente a la forma de los bloques, apenas existen ejemplos de geometrías distintas de la dos más habituales: el bloque rectangular y el disco. Suele imponerse además la simetría en torno al píxel central, según los ejes coordenados para el caso de bloques rectangulares que pasan así a convertirse en cuadrados, y simetría radial en el caso de bloques circulares. 71 En el apartado anterior se destacó la importancia de considerar un número mínimo de píxeles para no incurrir en errores de estimación altos. Si se analiza cuál debe ser el tamaño de los bloques para contener, al menos, un número prescrito de muestras, se obtienen los resultados mostrados en la Figura 5.2, en la que se observa que el bloque en forma de disco permite un mejor ajuste al número de vecinos requeridos, mientras que con bloques cuadrados a menudo se excede notablemente ese valor prescrito. Figura 5.2 Adaptación de los bloques al número de vecinos requeridos en función de su forma. Sin embargo, y tal como demuestran los resultados reflejados en la Figura 5.3, el valor estimado para la desviación estándar de un proceso gaussiano de media nula se aproxima más al valor real utilizando bloques cuadrados, motivo por el que es preferible esta forma frente a la de bloque circular. Adicionalmente, desde el punto de vista de la implementación, resulta más cómodo trabajar con este tipo de estructuras. Figura 5.3 Exactitud en la estimación de la desviación estándar en función de la forma del bloque. 78 Figura 6.1 Diagrama de bloques para estimación empleando el operador laplaciano. (a) (b) (c) (d) Figura 6.2 Caracterización de ruido aplicando prefiltrado. (a) Imagen ruidosa, (b) estructura de la imagen obtenida con el filtrado basado en derivadas direccionales, (c) estimación del ruido y (d) estimación de la desviación estándar del ruido. 79 (a) (b) Figura 6.3 Caracterización de ruido con filtro laplaciano. (a) Imagen ruidosa y (b) estimación de la desviación estándar a partir de la salida del filtro. 6.3 Supresión de bordes A pesar de la alta insensibilidad del operador laplaciano a la imagen libre de ruido, la señal que devuelve siempre presenta parte de la estructura de dicha imagen. Este remanente de la imagen ideal se presenta con frecuencia como líneas delgadas, principalmente en las regiones en las que el contraste de intensidades entre píxeles cercanos es alta, es decir, allí donde se encuentran los contornos de los objetos observados. Esta distorsión provoca una sobreestimación de los valores locales de la varianza que, tal y como se ha comentado en capítulos anteriores, se extiende al resto de la imagen si se practica un suavizado del campo de variaciones resultantes de la estimación. Para reducir en la medida de lo posible este efecto, se recurre a la detección de bordes y al posterior descarte de píxeles próximos a los contornos detectados antes de abordar la estimación de los parámetros, tal como aparece en el diagrama de la Figura 6.4. Los valores puntuales de varianza para esos píxeles se obtienen como promedios de los obtenidos para píxeles cercanos no descartados o en caso de ser necesario, con estrategias de reparación como la propuesta en el apartado 5.3. En lo referente al método de detección de bordes, destacan dos estrategias. En la primera se estima la magnitud del gradiente de la imagen a partir de máscaras sencillas (Sobel o Prewitt) y se considera perteneciente al contorno a todo píxel en el que la magnitud de este gradiente sea elevada. La segunda está 80 Figura 6.4 Diagrama de bloques con la incorporación de la supresión de bordes. basada en el algoritmo Canny [55] de detección de bordes en la que también se estima el gradiente pero que incorpora unas etapa previa de filtrado y una posterior de seguimiento de las líneas de contornos atendiendo a la conectividad entre píxeles. Requiere de una mayor carga computacional pero los resultados se ajustan mucho mejor a la localización real de los contornos. (a) (b) Figura 6.5 Resultados de la detección de bordes utilizando (a) la magnitud del gradiente y (b) el algoritmo Canny. 81 En la Figura 6.5 se muestran conjuntamente los resultados de la detección de bordes para las dos estrategias: magnitud del gradiente y algoritmo Canny. En ambos casos se ha aplicado una dilatación artificial posterior para aumentar el grosor del borde detectado, añadiendo también los píxeles vecinos. El comportamiento del algoritmo Canny es manifiestamente mejor. 6.4 Estimación de estadísticos En lo referente a la estimación de estadísticos, tres son las posibles entradas a esta etapa: la salida del operador laplaciano, la estimación del ruido proveniente de la etapa de prefiltrado o la señal ruidosa. La señal considerada dependerá del estimador que se utilice. El primer método es el más sencillo pero válido sólo se desea estimar el valor de la varianza, ya que toma como entrada la señal devuelta por el operador laplaciano en la que se ha filtrado la información de la intensidad de la imagen. Si L es la salida del filtro laplaciano para un píxel concreto, la estimación de la varianza en él viene dada por la expresión: 𝜎=16 |𝐿| (6.2) Una segunda estrategia consiste en considerar la estimación de la señal de ruido devuelta por la etapa de prefiltrado, sin pasar por el filtro laplaciano, lo que permite calcular no sólo la varianza sino cualquier estadístico. Sin embargo, debe tenerse en cuenta que al no utilizar el operador laplaciano, la estimación del ruido está contaminada por los detalles de la imagen libre de ruido que la etapa de prefiltrado eliminó sólo parcialmente. Sobre esta señal se puede aplicar directamente el estimador del estadístico que interese (media, varianza, kurtosis, etc.), o alternativamente, calcular el histograma y utilizar esa información para realizar un ajuste con la función de densidad acumulada propia del modelo de ruido, obteniendo los parámetros que definen esa distribución y pudiendo calcular a partir de ellos el estadístico de interés. La última estrategia consiste en tomar directamente la señal ruidosa, sin etapa de prefiltrado, procesamiento por bloques ni operador laplaciano, y aplicar un estimador de 82 máxima verosimilitud de los parámetros del modelo de ruido. Debe tenerse en cuenta que cada estimación representa en este caso un problema de optimización y, por tanto, el coste computacional es muy alto. En un trabajo reciente, Liu [36] propone un estimador que proporciona directamente la terna de parámetros del modelo de ruido de Selva-Alparone pero de forma global, es decir, una terna para toda la imagen. Podría plantearse procesar la imagen por bloques, estimar la terna en cada uno de ellos e interpolar los valores restantes empleando la técnica descrita en el apartado 5.3, pero el coste computacional se dispara convirtiéndolo en impracticable. Sea cual sea la señal de entrada considerada y el estimador escogido, la señal devuelta por esta etapa es una imagen completa para el estadístico deseado, donde los valores han sido estimados de forma local y por lo tanto puede existir una diferencia acusada entre valores adyacentes, contraviniendo la hipótesis de que los estadísticos del ruido varían de forma suave. En el siguiente apartado se aborda este problema. 6.5 Suavizado espacial de estadísticos Por la propia naturaleza aleatoria del ruido, el campo de valores estimados para cualquier estadístico, y en concreto para la varianza, presenta siempre variaciones de alta frecuencia, lo que contradice la hipótesis de trabajo de fenómenos aleatorios en los que sus estadísticos varían de forma suave en la imagen. Estas variaciones de alta frecuencia pueden justificarse si se tiene en cuenta que, aunque se usen estimadores consistentes y por grande que sea el entorno o bloque escogido, el número de muestras disponibles para la estimación local es siempre reducido, lo que mantiene un grado de incertidumbre en los valores obtenidos. Con el fin de alcanzar un resultado consistente con las hipótesis de suavidad de los estadísticos del ruido, puede incluirse una etapa final de suavizado, siendo propuesto el suavizado local iterativo con un número prefijado de iteraciones. Alternativamente, si se dispone de información adicional sobre la forma que adoptan en la realidad estos parámetros del ruido (parábolas centradas en el centro de la imagen, nulo en las esquinas, o cualquier otro modelo geométrico) puede sustituirse el suavizado por un ajuste de curvas en el sentido de mínima energía del error cometido. El resultado será entonces tan suave como lo sea el modelo geométrico de referencia, y la aproximación 83 será mejor, pero desde el punto de vista práctico carece de interés puesto que la forma en la que los estadísticos evolucionan en el dominio de la imagen no es conocida a priori. En la Figura 6.6 la imagen de la izquierda representa la estimación directa con el estimador insesgado de la varianza, antes de la etapa de suavizado. En la derecha se muestra el mismo ejemplo tras aplicar la estrategia de suavizado descrita. La imagen suavizada no está exenta de variaciones considerables de sus valores, pero estas variaciones son menos frecuentes y acusadas que en la versión sin suavizar. (a) (b) Figura 6.6 Efecto del suavizado en los valores estimados: imagen (a) antes y (b) después del suavizado. 84 7. EXPERIMENTOS Y RESULTADOS En capítulos anteriores se ha descrito el modelo de ruido propuesto, el conjunto de filtros disponibles para restauración de imágenes ruidosas y prefiltrado, así como las distintas estrategias de estimación caracterización del ruido, englobadas en dos categorías: basadas en filtros y basadas en el procesamiento por bloques. El objetivo de este capítulo es evaluar el comportamiento de los algoritmos y obtener experimentalmente el error cometido en la estimación de parámetros del ruido. 85 7.1 Descripción del experimento En las siguientes secciones de este capítulo se evaluará el error cometido en la estimación local de la desviación estándar de ruido para varias imágenes, estrategias de estimación y niveles de ruido. El que se supondrá como verdadero valor de la desviación estándar (𝜎) se genera aleatoriamente en cada ensayo, procediendo de la siguiente manera: 1. Se establece un valor de referencia 𝜎𝑟𝑒𝑓 y un rango de variación Δ𝜎. 2. Se genera para cada vértice de la imagen un valor aleatorio de 𝜎 en el intervalo [𝜎𝑟𝑒𝑓−Δ𝜎 2, 𝜎𝑟𝑒𝑓+Δ𝜎/2]. 3. Se calcula el valor de 𝜎 mediante interpolación bilineal, imponiendo así las hipótesis referentes a la variación suave del parámetro en la imagen. (a) (b) Figura 7.1 Aspecto de los campos de 𝝈 generados aleatoriamente. Los campos de 𝜎 obtenidos de esta manera tienen el aspecto de los mostrados en la Figura 7.2. Con estos valores se genera la señal de ruido y se suma a la señal original. Una vez obtenido el campo de 𝜎 estimado, esto es, 𝜎, se evalúa el error cometido aplicando cinco medidas:  Raíz del error cuadrático medio (rmse).  Error absoluto medio (mae). 86  Error absoluto máximo (maxae).  Error relativo medio (mre).  Error relativo máximo (maxre). En las tablas se muestran los valores promedio para 100 realizaciones. 7.2 Resultados para estimación basada en filtros (prefiltrado) En este apartado se muestran los resultados obtenidos para la estimación de la desviación estándar del ruido aditivo añadido, utilizando un sistema con prefiltrado. El filtro empleado es el basado en derivadas direccionales, propuesto en el Capítulo 4. El valor verdadero de esta desviación se ha generado siguiendo el procedimiento descrito en la sección anterior, con un valor de referencia 𝜎𝑟𝑒𝑓 = 10 y con un margen de variación del 20%. El campo de 𝜎 obtenido es el mostrado en la Figura 7.4, mientras que la señal ruidosa y la estructura de la imagen obtenido con el prefiltrado se muestran, respectivamente en las Figuras 7.2 y 7.3. Figura 7.2 Imagen con ruido aditivo (𝝈 = 10) 87 Figura 7.3 Estructura de la imagen (prefiltrado) Figura 7.4 Verdadero valor puntual de la desviación estándar (𝝈𝒓𝒆𝒇 = 10) 94 Figura 7.12 Estimación final de la desviación estándar (𝝈𝒓𝒆𝒇 = 2).  Estimación con 𝜎𝑟𝑒𝑓=𝟓. Figura 7.13 Verdadero valor puntual de la desviación estándar (𝝈𝒓𝒆𝒇 = 5). 95 Figura 7.14 Estimación final de la desviación estándar (𝝈𝒓𝒆𝒇 = 5). Tabla 7.6 Valores de error en la estimación (𝝈𝒓𝒆𝒇 = 5). Medida del error en la estimación Valor Raíz del error cuadrático medio (rmse) 0.9630 Error absoluto medio (mae) 0.7467 Error absoluto máximo (maxae) 4.2143 Error relativo medio (mre) 0.1566 Error relativo máximo (maxre) 0.8739 96  Estimación con 𝜎𝑟𝑒𝑓=𝟐𝟎. Figura 7.15 Verdadero valor puntual de la desviación estándar (𝝈𝒓𝒆𝒇 = 20). Figura 7.16 Estimación final de la desviación estándar (𝝈𝒓𝒆𝒇 = 20). 97 Tabla 7.7 Valores de error en la estimación (𝝈𝒓𝒆𝒇 = 20). Medida del error en la estimación Valor Raíz del error cuadrático medio (rmse) 4.8781 Error absoluto medio (mae) 4.0700 Error absoluto máximo (maxae) 9.0572 Error relativo medio (mre) 0.2241 Error relativo máximo (maxre) 0.4882 En general, se observa un ligero incremento en el error de estimación conforme aumenta el nivel de ruido, manteniéndose en cualquier caso por debajo del 25% de error relativo. Sí es significativo el aumento en el error cometido para valores muy bajos de la desviación estándar como, por ejemplo, 𝜎𝑟𝑒𝑓 = 1, situación en la cual el error relativo se aproxima al 40%. Este valor tan elevado puede explicarse, entre otros factores, por la sobreestimación del ruido que se produce al utilizar el filtro laplaciano cuando el nivel de ruido es bajo [28]. 7.5 Resultados para la estimación de ruido en imágenes de rayos X Procediendo de forma análoga a como se hizo con ruido aditivo en las secciones anteriores, se evalúa a continuación el comportamiento del método propuesto para la estimación de ruido en imágenes de rayos X. Para esta modalidad de imagen médica, el modelo de ruido apropiado es el puramente multiplicativo, que queda totalmente definido por el parámetro 𝜎𝑢𝑒𝑞. El valor de este parámetro es significativamente menor que el de la desviación estándar del ruido, puesto que esta última depende del valor de intensidad de la imagen libre de ruido. En este escenario, el hecho de que en la estimación de estadísticos del ruido se pueda identificar la escena representada en la imagen es fruto de la dependencia multiplicativa y no un efecto indeseado provocado por la decorrelación insuficiente entre ruido y señal. En la Figura 7.17 se muestran las imágenes de rayos X original y ruidosa, el valor verdadero de la desviación estándar del ruido y las estimaciones del estadístico antes y después de la etapa de suavizado. 98 (a) (b) (c) (d) (e) Figura 7.17 Estimación de la desviación estándar del error en imágenes de rayos X. (a) Imagen original, (b) imagen ruidosa, con 𝝈𝒖𝒆𝒒= 0.05, (c) valor verdadero de la desviación estándar, (d) valor estimado sin suavizar y (e) tras la etapa de suavizado. 99 Tabla 7.8 Valores de error en la estimación de 𝝈 en imágenes de rayos X. Medida del error en la estimación Valor (sin suavizar) Valor (suavizando) Raíz del error cuadrático medio (rmse) 6.3014 6.2918 Error absoluto medio (mae) 3.8474 3.8297 Error absoluto máximo (maxae) 63.7250 63.6675 Error relativo medio (mre) 0.3289 0.3246 Error relativo máximo (maxre) 0.9002 0.8981 Las medidas de error en la estimación de la desviación estándar del ruido, presentadas en la Tabla 7.8, muestran un incremento en el error relativo medio, que crece respecto del obtenido para ruido aditivo. La mayor complejidad del modelo de ruido en imágenes de rayos X y, en concreto, la dependencia de éste con el valor de intensidad de la imagen, dificulta la extracción del ruido y la estimación de sus estadísticos. En cualquier caso, y al margen de los valores numéricos, se observa que cualitativamente el método de estimación es capaz de seguir la evolución de la desviación estándar en la imagen. 7.6 Resultados para la estimación de ruido en imágenes de resonancia magnética El último conjunto de experimentos lo constituyen las estimaciones de ruido en imágenes de resonancia magnética, en las que rige el modelo descrito en el Capítulo 3. A diferencia de en imágenes de rayos X, el ruido no puede considerarse como puramente multiplicativo, sino que incluye una componente aditiva independiente de la intensidad de la imagen original. En escáneres de resonancia magnética, la varianza del error introducido en ambas componentes, real e imaginaria, es la misma, motivo por el cual el modelo de ruido queda totalmente determinado por un único parámetro, 𝜎𝑛. De nuevo, la dependencia del ruido con la intensidad de la imagen original provoca que en la representación de la desviación estándar del ruido puedan identificarse los contornos de la escena observada. En la Figura 7.18 se muestra el resultado de un experimento sobre imágenes de resonancia magnética, incluyendo las versiones original y ruidosa de la imagen, el verdadero valor de 100 (a) (b) (c) (d) (e) Figura 7.18 Estimación de la desviación estándar del error en imágenes de resonancia magnética. (a) Imagen original, (b) imagen ruidosa, con 𝝈𝒏= 0.12, (c) valor verdadero de la desviación estándar, (d) valor estimado sin suavizar y (e) tras la etapa de suavizado. 101 Tabla 7.9 Valores de error en la estimación de 𝝈 en imágenes de resonancia magnética. Medida del error en la estimación Valor (sin suavizar) Valor (suavizando) Raíz del error cuadrático medio (rmse) 6.1592 6.2000 Error absoluto medio (mae) 3.5821 3.6004 Error absoluto máximo (maxae) 24.2970 24.1543 Error relativo medio (mre) 0.6509 0.6642 Error relativo máximo (maxre) 10.4356 20.6900 la desviación estándar y los valores estimados de ese parámetros antes y después de la etapa de suavizado. Los valores de error obtenidos en los experimentos para ruido en imágenes de resonancia magnética, detallados en la Tabla 7.9, muestran un comportamiento discreto del método, con un error relativo en la estimación de la desviación estándar superior al 65% y que no mejora con la etapa de suavizado. La existencia de dos componentes de ruido: la multiplicativa, similar a la presente en imágenes de rayos X, y la aditiva, compromete aún más la estimación de los estadísticos del ruido, lo que provoca el aumento en el error cometido. En cuanto a la capacidad del estimador de seguir la evolución del parámetro del ruido a través de la imagen, ésta es apreciable, aunque con menor acierto que en el caso de los dos modelos anteriores (aditivo y multiplicativo). 102 8. CONCLUSIONES A tenor de los resultados obtenidos en los distintos ensayos y experimentos, pueden extraerse las conclusiones que se detallan a continuación:  Se ha propuesto un modelo de ruido en imágenes satisfactorio en el sentido de que es general y sencillo en su formulación ya que todos los términos responden a la misma expresión matemática. Además, modelos específicos de ruido en modalidades de imagen médicas, como por ejemplo resonancia magnética o rayos X, pueden ser tratados ahora como casos particulares del modelo general.  El filtro de imágenes propuesto, basado en derivadas direccionales, aunque fue pensado originalmente como etapa de prefiltrado, tiene por sí solo un buen comportamiento como filtro restaurador, suavizando la imagen para reducir el efecto del ruido a la vez que se respetan los contornos de los objetos representados. Este comportamiento es especialmente atractivo en ámbitos como la imagen médica en los que a menudo resulta crucial conservar estructuras muy delgadas.  El prefiltrado se ha mostrado como una estrategia poco conveniente para la caracterización del ruido ya que proporciona una estimación del mismo con una alta correlación con la señal ruidosa.  El procesamiento por bloques ofrece buenos resultados en la estimación de los parámetros del ruido, principalmente para la varianza ya que en ese caso puede emplearse el filtro laplaciano y cancelar así la correlación con la señal ruidosa, que provoca una sobreestimación del parámetro.  El algoritmo diseñado para reparar los mapas incompletos de estimaciones de parámetros (aquéllos con valores no definidos para algunos píxeles), constituye una potente herramienta puesto que es capaz de reconstruir señales a partir de un reducido número de valores puntuales uniformemente distribuidos. Piénsese en su utilidad para acelerar el procesamiento por bloques de la imagen cuando el tiempo sea un factor crítico. 103  Para ruidos multiplicativo y general, presentes en las imágenes de modalidades médicas como rayos X, ultrasonidos y resonancia magnética, el error de estimación crece pero no llega a los niveles del obtenido con estrategias basadas en filtros.  Para el modelo general, la mejor alternativa parece ser el uso de estimadores de máxima verosimilitud que, de momento, son válidos para estimar los parámetros del modelo en toda la imagen en su conjunto, sin poder seguir las variaciones de los mismos en el dominio de la imagen. Se plantean, además, los siguientes trabajos futuros:  En principio, el modelo de ruido se considera satisfactorio. Investigar sobre modelos existentes que no tengan cabida en él podría ser interesante.  Para el filtro basado en derivadas direccionales, resultaría conveniente revisar el criterio de cribado para que no se degraden las prestaciones en presencia de niveles altos de ruido, lo que en este momento representa su principal debilidad.  En lo concerniente a la estimación de los parámetros del modelo de ruido, podría plantearse el uso combinado de procesamiento por bloques, estimadores de máxima verosimilitud y el algoritmo propuesto de reparación de estimaciones incompletas de parámetros con el fin de probar una estrategia que proporcionara la terna de parámetros del modelo de Selva-Alparone y que fuera capaz de capturar las variaciones de los mismos en la imagen. 110 10. PLANOS Y PROGRAMAS Junto con esta memoria, se facilita un CD con la base de imágenes utilizada y las implementaciones en código MATLAB de todos los algoritmos descritos en ella. Contiene además otros algoritmos y experimentos varios, no mencionados en el documento pero que fueron objeto de estudio durante el desarrollo de este Proyecto Fin de Carrera. Para la ejecución de los mismos, se recomienda al lector ubicarse en la carpeta donde se encuentran el código fuente y la base de datos y ejecutar la siguiente instrucción: addpath(genpath('.'), '-end'); que incluye el árbol completo de directorios del proyecto en la ruta de búsqueda de MATLAB, facilitando la ejecución de las rutinas. 111 PARTE III PLIEGO DE CONDICIONES 112 11. PLIEGO DE CONDICIONES Para la realización de este proyecto se ha hecho uso de un conjunto de herramientas software y equipos hardware cuyas características principales se detallan en los siguientes apartados. 11.1 Recursos Software Las siguientes herramientas software se  Windows® 7 Professional: sistema operativo empleado en el proyecto, especialmente en los apartados de desarrollo de los algoritmos y la redacción de la documentación.  Linux (Ubuntu 12.04 LTS 64-bits): sistema operativo empleado principalmente para la realización de ensayos con código compilable, así como la recopilación y gestión de las imágenes de las bases de datos.  MATLAB® versión 7.10.0.499 (R2010a): o Image Processing Toolbox versión 7.0 o Statistics Toolbox versión 7.3  Microsoft Office® 2013: paquete de herramientas que incluye Microsoft Word, Micrososft Excel, y Microsoft Power Point, utilizadas para elaboración de la memoria y la presentación del proyecto. 11.2 Recursos Hardware  Estación portátil Dell Precission M4800, dotada con procesador Intel® Core™i7-4800MQ CPU @2.70 GHz, 8GB de memoria RAM y 500 GB de disco duro.  PC de sobremesa Acer Aspire M3641, con procesador Intel® Core™2 Quad CPU Q6600 @ 2.40 GHz, 3 GB de memoria RAM y 500 GB de disco duro. 113 PARTE IV PRESUPUESTO 114 12. PRESUPUESTO Don Guillermo Valentín Socorro Marrero, autor del presente Proyecto Fin de Carrera, declara que: El Proyecto Fin de Carrera con título “Modelado de Ruido en Imágenes. Aplicaciones Médicas”, desarrollado en el Grupo de Imagen, Tecnología Médica y Televisión (GIMET) de la Escuela de Ingeniería de Telecomunicación y Electrónica de la Universidad de Las Palmas de Gran Canaria, en un periodo de 10 meses, tiene un coste total de desarrollo de 51.707,60 euros, correspondiente a la suma de las cantidades consignadas a los apartados que se detallan en este presupuesto, calculadas siguiendo las recomendaciones del Colegio Oficial de Ingenieros de Telecomunicación (COIT). Fdo: Guillermo Valentín Socorro Marrero. Autor del proyecto. Junio de 2016. 115 12.1 Desglose del Presupuesto Para la realización del presupuesto se han seguido las recomendaciones del Colegio Oficial de Ingenieros de Telecomunicación (COIT) sobre los baremos orientativos para trabajos profesionales, publicadas en el documento “Tarifas de Derechos para Visado y Baremos de Honorarios Orientativos para 2012”. Los honorarios y umbrales monetarios que figuran en dicho documento se han actualizado aplicando los coeficientes de corrección correspondientes en función de la variación del IPC, publicada por el Instituto Nacional de Estadística (INE). El presupuesto se ha desglosado en varias secciones en las que se detallan los distintos costes asociados al desarrollo del proyecto. Estos costes se dividen en:  Recursos materiales.  Trabajo tarifado por tiempo empleado.  Costes de redacción del proyecto.  Material fungible.  Derechos de visado del COIT.  Gastos de tramitación y envío.  Aplicación de impuestos. 12.2 Recursos Materiales Esta partida incluye los gastos derivados de la amortización del hardware empleado en el desarrollo del proyecto, así como el relativo a las licencias de las distintas herramientas software utilizadas, tanto para el desarrollo de los algoritmos como para la redacción de la memoria. Se estipula el coste de amortización para un período de 3 años. Para ello, se utiliza un sistema de amortización lineal o constante, en el que se supone que el inmovilizado material se deprecia de forma constante a lo largo de su vida útil. La cuota de amortización anual se calcula usando la siguiente fórmula: 𝐶𝐴=𝑉𝐴−𝑉𝑅 𝑉𝑈 (12.1) 116 donde: 𝐶𝐴 es la cuota de amortización anual. 𝑉𝐴 es el valor de adquisición. 𝑉𝑅 es el valor residual que se supone que tendrá el elemento al final de su vida útil. 𝑉𝑈 es el número de años de vida útil. Los costes de amortización se calculan para el primer año, teniendo en cuenta la duración del proyecto. El valor final de dicha amortización se rige por la siguiente expresión: 𝐴𝑀=𝐶𝐴· 𝑇𝑈 12 (12.2) siendo: 𝐴𝑀 el valor de la amortización. 𝐶𝐴 la cuota de amortización anual. 𝑇𝑈 el número de meses correspondientes al tiempo de uso. 12.2.1 Recursos software Las herramientas software utilizadas en este proyecto, que incluyen el sistema operativo, los paquetes para el desarrollo de algoritmos y las herramientas de edición de la documentación, son:  Windows® 7 Professional.  Linux (Ubuntu 12.04 LTS 64-bits).  MATLAB® versión 7.10.0.499 (R2010a). o ImageProcessingToolbox versión 7.0. o StatisticsToolbox versión 7.3.  Microsoft Office® 2013. En la Tabla 12.1 se desglosan los costes de los recursos software, calculados mediantes las ecuaciones (12.1) y (12.2), teniendo en cuenta que el periodo de amortización lineal es de 3 años y la duración aproximada del proyecto se ha establecido en 10 meses. 117 Tabla 12.1 Costes de amortización de los recursos software. Concepto VA (€) VR (€) CA (€) AM (€) Windows® 7 Professional 185,00 0,00 61,67 51,39 Ubuntu 0,00 0,00 0,00 0,00 MATLAB 2000,00 0,00 666,67 555,56 Image Toolbox 1000,00 0,00 333,33 277,78 Statistics Toolbox 1000,00 0,00 333,33 277,78 Microsoft Office® 2013 269,00 0,00 89,67 74,72 Total (€) 1237,23 Los costes de amortización total del software ascienden a mil doscientos treinta y siete euros con veintitrés céntimos (1237,23 €). 12.2.2 Recursos hardware Los recursos hardware empleados para el desarrollo del proyecto se corresponden con el equipamiento informático, que consta de los siguientes ordenadores:  Estación portátil Dell Precission M4800, dotada con procesador Intel® Core™i7-4800MQ CPU @2.70 GHz, 8GB de memoria RAM y 500 GB de disco duro.  PC de sobremesa Acer Aspire M3641, con procesador Intel® Core™2 Quad CPU Q6600 @ 2.40 GHz, 3 GB de memoria RAM y 500 GB de disco duro. En la Tabla 12.2 se muestran los costes de amortización de los recursos hardware, considerando de nuevo amortización lineal a 3 años y tiempo de uso de 10 meses. Tabla 12.2 Costes de amortización de los recursos hardware. Concepto VA (€) VR (€) CA (€) AM (€) Dell Precission M4800 1585,00 0,00 528,33 440,28 Acer Aspire M3641 360,00 0,00 120,00 100,00 Total (€) 540,28 118 Los costes de amortización total del hardware ascienden a quinientos cuarenta euros con veintiocho céntimos (540,28 €). 12.3 Trabajo tarifado por tiempo empleado Atendiendo a los baremos orientativos del COIT, el importe correspondiente a las horas de trabajo del ingeniero empleadas en la realización del proyecto, se calcula mediante la siguiente fórmula: 𝐻=𝐶𝑡·(74,88·𝐻𝑛+ 96,72·𝐻𝑒)· 𝐶𝑝 € (12.3) expresión en la que: 𝐻 son los honorarios totales por el tiempo dedicado. 𝐻𝑛 son las horas normales trabajadas, es decir, aquéllas dentro de la jornada laboral. 𝐻𝑒 son las horas especiales, fuera de la jornada laboral. 𝐶𝑡 es el factor de corrección en función del número de horas trabajadas. 𝐶𝑝 es el factor de corrección por evolución del IPC. El tiempo empleado para la realización de este proyecto es de aproximadamente 10 meses, equivalente a 1.500 horas (10 meses x 4 semanas/mes x 37,5 horas/semana). Todas estas horas se consideran normales, siendo nulo el número de horas especiales. Por otro lado, consultando la Tabla 12.3 se comprueba que el factor de corrección en función del tiempo empleado es 𝐶𝑡=0,40, como corresponde a trabajos que superen las 1.080 horas. Además, la variación del IPC desde la finalización del año en el que fueron publicadas las recomendaciones del COIT hasta diciembre de 2015, facilitada por el INE, es de -0,5%, resultando un factor de corrección por evolución de los precios 𝐶𝑝= 0,995. La expresión de los honorarios totales queda entonces como: 𝐻=0,40·(74,88·1.500+ 96,72·0)· 0,995=44.703,36 € Los honorarios totales por tiempo dedicado antes de impuestos ascienden a cuarenta y cuatro mil setecientos tres euros con treinta y seis céntimos (44.703,36 €). 119 Tabla 12.3 Factor de corrección en función del tiempo de realización del proyecto. Tiempo empleado Factor de Corrección 𝑪𝒕 Hasta 36 horas 1,00 De 36 a 72 horas 0,90 De 72 a 108 horas 0,80 De 108 a 144 horas 0,70 De 144 a 180 horas 0,65 De 180 a 360 horas 0,60 De 360 a 540 horas 0,55 De 540 a 720 horas 0,50 De 720 a 1.080 horas 0,45 Más de 1080 horas 0,40 12.4 Costes de redacción del proyecto Los gastos asociados a la redacción del proyecto se obtienen aplicando la siguiente expresión: 𝑅=0,07·𝑃·𝐶ℎ (12.4) en la que: R son los gastos de redacción del proyecto. 𝑃 es el valor presupuestado de ejecución material, considerando los gatos de amortización y personal. 𝐶ℎ es el coeficiente de ponderación de costes de redacción en función del valor presupuestado. Tabla 12.4 Cálculo del presupuesto de ejecución material. Concepto Importe (€) Recursos software 1237,23 Recursos hardware 540,28 Tarifa por tiempo empleado 44.703,36 Total 46.480,87