Repositorio Institucional de Documentos

Proyecto Fin de Carrera Desarrollo de una aplicación software para la fusión robusta de mapas de profundidad y reconstrucción 3D densa Autor Sergio Ibáñez Sanahuja Director Pedro Piniés Rodríguez Escuela de Ingeniería y Arquitectura 2013 Dedicado a toda mi familia. Para los que est´an, y para los que est´an por venir. La investigaci´on es para ver lo que todo el mundo ha visto, y pensar en lo que nadie ha pensado. Albert Szent-Gy¨orgi Agradecimientos Quisiera expresar mi agradecimiento a todo el mundo que me ha apoyado a lo largo de la realizaci´on del proyecto, a mi familia, especialmente a mi madre y mi pareja ya que su apoyo ha sido incondicional, a mis compa˜neros de clase durante estos maravillosos cinco a˜nos, a mis compa˜neros de piso por agradarme este largo verano especialmente a mi compa˜nero en el d´ıa a d´ıa Francisco Javier, con el que he compartido un verano de duro trabajo. Por ´ultimo al director del proyecto, Pedro, por su paciencia y apoyo durante todo el proyecto. Gracias a todos porque sin vuestro apoyo esto no hubiera sido posible. ´ Indice 1. Introducci´on 1 1.1. Contexto en el que se enmarca el proyecto . . . . . . . . . . . . . . 1 1.1.1. Localizaci´on de la c´amara . . . . . . . . . . . . . . . . . . . 1 1.1.2. Mapas de profundidad . . . . . . . . . . . . . . . . . . . . . 2 1.2. Estadodelarte ............................. 3 1.3. Procesopropuesto............................ 3 1.3.1. Fusi´on de mapas . . . . . . . . . . . . . . . . . . . . . . . . 3 1.4. Algoritmo Primal Dual . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.5. Organizaci´on .............................. 5 2. Image Denoising 7 2.1. Introducci´on............................... 7 2.2. TV-ROF................................. 10 2.2.1. Proceso ............................. 10 2.3. Huber-ROF ............................... 11 2.3.1. Proceso ............................. 12 2.4. TGV-ROF................................ 13 2.4.1. Proceso ............................. 14 2.5. Elecci´on de par´ametros y resultados . . . . . . . . . . . . . . . . . . 16 2.5.1. Par´ametros para la convergencia . . . . . . . . . . . . . . . . 16 2.5.2. Par´ametros de ponderaci´on . . . . . . . . . . . . . . . . . . 18 3. Fusi´on de mapas 25 3.1. Mapas de profundidad . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.2. Mapasvirtuales............................. 25 3.3. Fusi´on .................................. 27 3.3.1. Proceso ............................. 27 3.4. Resultados de la fusi´on . . . . . . . . . . . . . . . . . . . . . . . . . 29 3.4.1. Generaci´on de mapas virtuales . . . . . . . . . . . . . . . . . 29 3.4.2. Fusi´on.............................. 29 3.4.3. Medida de tiempos . . . . . . . . . . . . . . . . . . . . . . . 32 i ´ INDICE ´ INDICE 4. Conclusiones 35 4.1. LineasFuturas ............................. 35 4.2. Valoraci´on personal . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 A. Optimizaci´on convexa: Algoritmo de Primal Dual 41 A.1. Algoritmo Primal Dual . . . . . . . . . . . . . . . . . . . . . . . . . 41 A.2. Transformada de Legendre-Fenchel .................. 43 A.2.1.Valorabsoluto.......................... 44 A.2.2.NormadeHuber ........................ 47 A.2.3.TGV............................... 48 A.3. Resoluci´on del problema . . . . . . . . . . . . . . . . . . . . . . . . 50 A.3.1. Optimizaci´on convexa de primer orden . . . . . . . . . . . . 50 A.3.2. M´etodo del Gradiente . . . . . . . . . . . . . . . . . . . . . 51 A.3.3. Otros m´etodos de optimizaci´on . . . . . . . . . . . . . . . . 51 A.4.Proximalmap.............................. 52 A.4.1. Elecci´on de par´ametros . . . . . . . . . . . . . . . . . . . . . 53 B. Gesti´on del proyecto 55 B.1.Metodolog´ıa............................... 55 B.2. Tecnolog´ıas empleadas . . . . . . . . . . . . . . . . . . . . . . . . . 56 B.3. Herramientas utilizadas . . . . . . . . . . . . . . . . . . . . . . . . . 60 C. Programaci´on de Primal Dual 63 D. Interfaces 67 D.1. Interfaz de prueba de par´ametros . . . . . . . . . . . . . . . . . . . 67 D.1.1.Requisitos............................ 67 D.1.2. Resultado de la interfaz . . . . . . . . . . . . . . . . . . . . 68 D.1.3. Detalles de utilidad . . . . . . . . . . . . . . . . . . . . . . . 70 D.1.4. Detalles de implementaci´on . . . . . . . . . . . . . . . . . . 72 D.2. Interfaz de visualizaci´on . . . . . . . . . . . . . . . . . . . . . . . . 77 D.2.1.Requisitos............................ 77 D.2.2. Resultado de la interfaz . . . . . . . . . . . . . . . . . . . . 78 D.2.3. Detalles de implementaci´on . . . . . . . . . . . . . . . . . . 79 D.3.Interfazdefusi´on ............................ 83 D.3.1.Requisitos............................ 83 D.3.2. Resultado de la interfaz . . . . . . . . . . . . . . . . . . . . 84 D.3.3. Detalles de implementaci´on . . . . . . . . . . . . . . . . . . 87 ii ´ INDICE ´ INDICE E. Fusi´on 91 E.1.Mapasvirtuales............................. 91 E.1.1. Lectura del cubo . . . . . . . . . . . . . . . . . . . . . . . . 95 E.1.2. Programaci´on.......................... 96 iii Cap´ıtulo 1 Introducci´on La obtenci´on autom´atica de modelos 3D densos de un determinado entorno usando una c´amara como ´unico sensor tiene un gran inter´es debido a su bajo coste y al amplio abanico de aplicaciones: reconstrucci´on y modelado de obras art´ısticas, obtenci´on de modelos densos de ciudades (Google Maps), asistencia en cirug´ıas no invasivas como laparoscopias u otras aplicaciones que no dependen estrictamente de la reconstrucci´on como en la industria del entretenimiento donde se puede usar para desarrollar juegos para consolas Xbox o Wii que permiten interactuar al jugador con el mundo real usando realidad aumentada, detecci´on de obst´aculos en navegaci´on de veh´ıculos, segmentaci´on, etc. El proceso para obtener una reconstrucci´on densa a partir de un conjunto de im´agenes es laborioso y complejo. Desde un punto de vista de alto nivel podemos identificar distintos bloques, a continuaci´on se explica todo este proceso. 1.1. Contexto en el que se enmarca el proyecto En la figura 1.1 se muestran los distintos bloques en los se divide el proceso desde el c´alculo de los mapas de profundidad hasta las distintas aplicaciones que pudieran hacer uso de ellos. A continuaci´on se detallan de forma resumida estos bloques. 1.1.1. Localizaci´on de la c´amara La clave para poder conseguir extraer un depth map1de cualquier escena es conseguir informaci´on de la profundidad, para ello se requiere m´as informaci´on que una sola imagen. 1Depth map: Mapa de profundidad 1 Secci´on 1.1 1. Introducci´on Figura 1.1: Esquema de todo el proceso desde la consecuci´on de los mapas de profundidad hasta hasta las aplicaciones de este Los mapas de profundidad usados como entrada del algoritmo implementado en este proyecto, se han construido a partir de informaci´on capturada por una c´amara monocular. Puesto que de una sola imagen no es posible extraer informaci´on de la profundidad, como si que ser´ıa posible en una c´amara est´ereo, se requiere informaci´on redundante para conocerla. Esta se consigue al tomar varias im´agenes y conocer con precisi´on el desplazamiento de la c´amara al tomar estas, de este modo se podr´a deducir la profundidad a la que se encuentra el punto en cuesti´on. Este proceso tambi´en se puede realizar a la inversa, lograr a partir un mapa de profundidad la localizaci´on espacial de la c´amara en el momento de realizar las distintas capturas. 1.1.2. Mapas de profundidad En todos los aspectos, un mapa de profundidad es equivalente a una imagen, s´olo que en vez de almacenar valores de color, almacena distancias. A partir de un conjunto de im´agenes en niveles de gris y las correspondientes posiciones de la c´amara desde donde se tomaron se puede crear un mapa de profundidad mediante el algoritmo Dense Tracking and Mapping[RAND11]. La idea intuitiva es que si la profundidad de un punto en el espacio 3D est´a bien calculada, los niveles de gris de los p´ıxeles proyectados en todas las im´agenes en las que se vio deber´ıa coincidir. Desafortunadamente esta t´ecnica produce im´agenes de profun2 1. Introducci´on Secci´on 1.2 didad ruidosas y con espurios que pueden afectar a la calidad de las aplicaciones que hace uso de ellas. El objetivo del este proyecto es eliminar o al menos reducir estos efectos en los depth maps construidos como se explica en el apartado 1.3. 1.2. Estado del arte En cuanto al estado del arte las investigaciones en las que se basa este trabajo son relativamente recientes. Por un lado art´ıculos sobre sistemas que se encargan de todos los procesos necesarios para obtener una reconstrucci´on 3D usando una c´amara monocular, es decir abarcan todo el proceso mostrado en la figura 1.1, por ejemplo Dense Tracking and Mapping in Real-Time [RAND11] o Dense Reconstruction On–the–Fly [AWB]. Adem´as se encargan de realizarlo en tiempo real. Existe tambi´en otros art´ıculos que tratan este problema pero usando c´amaras Kinect que dan mucha informaci´on ya que consiguen de manera sencilla cada uno de las mapas de profundidad, por ejemplo Kinect-Fusion [SI11]. Por ´ultimo otro art´ıculo muy relacionado es el TGV-Fusion [TPB11], pero con la diferencia de que este se centra en la reconstrucci´on a partir de im´agenes a´ereas, por lo que parte del proceso realizado en este proyecto no es necesaria. En este caso las im´agenes de entrada van a ser im´agenes de corta distancia donde la escena cambia mucho con un leve movimiento del punto vista, por lo que ser´a necesaria la programaci´on de un mecanismo de virtualizaci´on como se explicar´a m´as adelante. Todos estos trabajos son la fuente de multitud de aplicaciones tales como las comentadas al inicio del cap´ıtulo. 1.3. Proceso propuesto La aportaci´on propuesta en este proyecto con respecto al esquema anterior aparece en la figura 1.2. En concreto se ha incluido un m´odulo que es le que se va a encargar que realizar la fusi´on. Este m´odulo se basa en el art´ıculo de TGV-fusion [TPB11] para im´agenes a´ereas pero con la particularidad de que se realizar´a un proceso intermedio para la transformaci´on de los mapas de entrada por lo dicho en el apartado anterior. 1.3.1. Fusi´on de mapas Como se observa en la figura 1.2, para realizar la fusi´on ser´a necesario m´as de un depth map, el nuevo bloque se encarga de fusionar estos mapas de profundidad que presenten cierto solapamiento (es decir que todos los mapas tengan una misma regi´on de la escena en com´un) para generar uno nuevo de mayor calidad. 3 Secci´on 1.4 1. Introducci´on Figura 1.2: Nuevo esquema del proceso desde la consecuci´on de los mapas de profundidad hasta hasta las aplicaciones de este pasando por la fusi´on Los mapas de profundidad resultantes de la etapa anterior a este nuevo m´odulo tienen informaci´on redundante, ya que se han generado varios con distintas im´agenes, este hecho da la capacidad de fusionar estos mapas para lograr uno de mayor calidad que los anteriores. Esta fusi´on se realiza aplicando t´ecnicas variacionales de optimizaci´on convexa. En este punto se centra el proyecto, concretamente en la aplicaci´on de este algoritmo basado en Primal Dual para la fusi´on de estos mapas de profundidad de forma eficiente y robusta. Como se indic´o en el resumen de la memoria, en una primera etapa del proyecto se estudiar´a el problema del denosising, es decir la eliminaci´on de ruido de una imagen, ya que es un problema bastante similar al de la fusi´on, ambos tienen la misma base, por lo que una vez programadas todas las normas2del denoising la programaci´on de la fusi´on ser´a m´as r´apida ya que se utilizar´an exactamente los mismos conceptos. 1.4. Algoritmo Primal Dual La optimizaci´on de funciones es una rama muy importante de estudio en las matem´aticas, en general, no es posible encontrar el m´ınimo o el m´aximo absoluto de un problema de optimizaci´on aunque s´ı existen, para determinados tipos de funciones, soluciones ´optimas o buenas heur´ısticas [NY04]. 2Normas: Cada una de las variantes del algoritmo de optimizaci´on utilizado en el denoisisg que se estudiar´an m´as adelante. 4 1. Introducci´on Secci´on 1.5 Afortunadamente en este proyecto nos encontramos en uno de esos casos, en los que las funciones a minimizar que nos vamos a encontrar, presentan una estructura especial que permite obtener la soluci´on global (´optima) del problema. La estructura gen´erica de los problemas que vamos a resolver, tanto en denoising como en TGV-fusion, es la siguiente: m´ın uF(Ku)+G(u) (1.1) Donde F y G son funciones convexas aunque no tienen por que ser diferenciables y K es un operador lineal. La soluci´on a esta optimizaci´on se basa en el algoritmo Primal Dual que como veremos m´as adelante, cuando se aplica a problemas de visi´on por computador, se puede acelerar considerablemnte usando tarjetas gr´aficas (GPUs) y programaci´on paralela en CUDA [CUD10]. 1.5. Organizaci´on A lo largo de memoria se van a ir abarcando conceptos de manera gradual para ayudar a la comprensi´on. Todo el proceso que se sigue a continuaci´on se basa en el algoritmo Primal dual cuyo an´alisis est´a realizado en el anexo A.Una vez comprendido este algoritmo se abarcar´an los distintos problemas de denoising que se trata, como se ha dicho, de una aplicaci´on muy directa e intuitiva del algoritmo Primal Dual (Cap´ıtulo 2), a continuaci´on se abarcar´a el problema de la fusi´on as´ı como todos los pasos previos para realizarla (Cap´ıtulo 3), por ´ultimo habr´a un peque˜no apartado de conclusiones de todo el proyecto donde se expondr´an las impresiones personales, y se explicar´an posibles l´ıneas futuras de trabajo en relaci´on al proyecto (Cap´ıtulo 4). 5 Secci´on 1.5 1. Introducci´on 6 Cap´ıtulo 2 Image Denoising Como se ha comentado, las funciones de la forma 1.1 se pueden aplicar a muchos problemas, relacionados o no con el tema de este proyecto. Una aplicaci´on casi directa, y muy intuitiva del problema es el uso de esta para eliminaci´on de ruido en im´agenes. En este cap´ıtulo, se van a detallar tres algoritmos con la forma 1.1 que nos dar´an pie a, en primer lugar conocer en mas detalle el algoritmo de optimizaci´on Primal Dual ya que se trata de una aplicaci´on mas sencilla que la fusi´on de mapas, y adem´as se aprender´an conceptos muy ´utiles para apartados posteriores. 2.1. Introducci´on Todos los problemas que se van a tratar en este cap´ıtulo tienen pr´acticamente la siguiente estructura: m´ın uZ|ru|+ 2kugk2 2d(⌦) (2.1) Esta esta integral afecta al espacio continuo de la imagen, y como tal se trata de un problema variacional [LIRF92]. Esta compuesta de dos funciones, para lograr la minimizaci´on habr´a que encontrar un resultado ´optimo para ambas partes. Cada una de estas partes tiene un significado valora unas caracter´ısticas concretas que deber´a tener la imagen resultado. |ru|: Se trata del regularizador, en este caso es el valor absoluto del gradiente de la imagen propio del m´etodo TV-ROF que se ver´a mas adelante. Este t´ermino tratar´a de valorar soluciones cuya variaci´on entre p´ıxeles sea peque˜na, para poder eliminar el ruido de alta frecuencia. 7 Secci´on 2.1 2. Image Denoising Figura 2.1: Ejemplo de una mala elecci´on de ,alaizquierdasehaescogidoun demasiado peque˜no, a la derecha uno demasiado grande, y en el centro uno que podr´ıa ser ´optimo  2kugk2 2: Esta parte se denomina data term. Ser´a com´un a todos los m´etodos estudiados en el denoising. Valora el hecho de que la soluci´on se parezca lo m´aximo a la imagen ruidosa. En este caso en particular se trata una norma cuadr´atica. Existen otras normas que se podr´ıan colocar en este punto tales como la norma L1consistente en una funci´on valor absoluto o la norma Huber que ser´a utilizada para la fusi´on y se explicar´a m´as adelante. En resumen, una colaboraci´on entre ambas eliminar´a detalles no deseados, mientras que preservar´a los detalles importantes tales como bordes El par´ametro ser´a el que pondere la relaci´on entre ambos t´erminos, por lo que un demasiado grande no eliminar´a el ruido por completo, mientras que uno demasiado peque˜no lo har´a en exceso eliminando excesiva informaci´on. En la figura 2.1 se observa este efecto. Como la expresi´on 2.1 es continua, ser´a necesario discretizarla, resultando la siguiente expresi´on. m´ın uF(Ku)+G(u) (2.2) Donde Kes un operador discreto equivalente al gradiente, en este caso se trata del c´alculo hacia delante del mismo. Es decir, para cada uno de los p´ıxeles calcula la diferencia entre este y el siguiente en filas y en columnas. Kui,j =(ui+1,j ui,j ,u i,j+1 ui,j) La expresi´on 2.2 se va a sustituir ahora por una equivalente con que la poder trabajar, ya que las funciones F(Ku)yG(u)notienenporqueserdiferenciables, por lo que no es posible la aplicaci´on directa de m´etodos cl´asicos. Para ello se realiza una transformaci´on que logra el prop´osito de convertir la expresi´on en diferenciable, pero no a cualquier precio ya que la minimizaci´on inicial se ha convertido en una maximizaci´on y minimizaci´on simult´aneas. 8 2. Image Denoising Secci´on 2.1 En resumen, este proceso realiza una dualizaci´on de la funci´on no diferenciable, para transformarla a diferenciable. Todo el proceso detallado para alcanzar la siguiente expresi´on a partir de la 2.2 se detalla en el anexo A. m´ın usup p <Ku,p>+G(u)F⇤(p)(2.3) m´ın usup ps.t. |p|1 <Ku,p>+G(u)(2.4) m´ın usup p C(u, p)(2.5) En primer lugar la expresi´on 2.3 es directamente la que resulta en el anexo A. En la 2.4, se encuentra la misma pero indicando la funci´on dual como restricci´on de la maximizaci´on. Por ´ultimo la expresi´on 2.5 es una simplificaci´on de las dos anteriores que facilitar´a las operaciones que vienen a continuaci´on. El algoritmo que se usar´a para la resoluci´on del problema tiene la siguiente estructura b´asica. 8 < : yn+1 =(I+@F⇤)1(yn+K¯xn)(1) xn+1 =(I+⌧@G)1(xn⌧K⇤yn+1)(2) ¯xn+1 =xn+1 +✓(xn+1 xn)(3) (2.6) Se trata de un algoritmo llamado proximal map, muy similar al gradiente ascendente/descendente pero con unas peculiaridades que lo hace m´as eficientes. Se resumen a continuaci´on las operaciones implicadas en cada uno de los pasos. 1. Se trata de la aplicaci´on directa del m´etodo del gradiente ascendente a la expresi´on 2.4, con la particularidad que se realiza la evaluaci´on de y(pen este caso) en el paso siguiente, y que se aplica el operador proximal map, para tener en cuenta la restricci´on que genera la funci´on dual. 2. Esta parte es incluso m´as parecida a un gradiente descendente cl´asico ya que para ella no existe restricci´on. 3. Esta parte es un factor de suavizado que se aplica a la soluci´on para lograr una convergencia mas r´apida. Para conocer m´as detalles de este m´etodo de optimizaci´on as´ı como saber en detalle las operaciones del proceso, dir´ıjase al anexo A. 9 Secci´on 2.5 2. Image Denoising 2.5. Elecci´on de par´ametros y resultados La soluci´on que se obtiene para cada uno de los m´etodos depende de dos tipos de par´ametros. Par´ametros de convergencia ⌧y(pasos del proximal map), que controlan la velocidad de convergencia de la soluci´on. Par´ametros de ponderaci´on (,↵,↵1y↵2) que controlan el comportamiento entre ajuste a los datos originales y la regularizaci´on. En la literatura sobre Primal-Dual no hay una soluci´on clara en la elecci´on de estos par´ametros, por lo que en este proyecto se han desarrollado una serie de aplicaciones que permiten tener una idea de su valor de forma razonable. (Estos valores dependen de la imagen sobre la que act´ua el m´etodo pero se va a suponer que los par´ametros ´optimos obtenidos de estas pruebas son extrapolables a cualquier imagen) En cuanto a las pruebas realizadas a lo largo de todo el cap´ıtulo se han desarrollado sobre una tarjeta gr´afica nVidia Tesla M2090. 2.5.1. Par´ametros para la convergencia Para asegurar la convergencia los pasos ⌧ydel proximal map deben cumplir la siguiente expresi´on. ⌧L21 (2.15) Donde L=kKksiendo kKk=p8, este valor es independiente de la imagen. Se puede encontrar una explicaci´on m´as detallada de esta elecci´on en el anexo A secci´on 4.1. Para buscar los valores ⌧yque mejoran la velocidad de convergencia, es decir para reducir el n´umero de iteraciones para lograr la soluci´on ´optima, se ha desarrollado el siguiente proceso. Se eligen ⌧ytal que se cumple la ecuaci´on 2.15 y como se trata de un problema convexo simplemente se sobreitera el m´etodo para llegar a´un resultado ´optimo u⇤. Esta soluci´on ser´a utilizada como criterio de parada. El criterio de parada se basa en la siguiente ecuaci´on. ek=v u u t nP ixels X 0 (uk ij u⇤ ij)2 Donde kes la iteraci´on actual. Se ha establecido un umbral para cuando e1e3 el algoritmo se detenga. 16 2. Image Denoising Secci´on 2.5 (a) Imagen de la prueba 0 0.5 1 1.5 2 2.5 3 3.5 4 20 40 60 80 100 120 140 160 180 σ Iteraciones TV Huber TGV (b) Resultados Figura 2.3: Resultados generados por la aplicaci´on de prueba de par´ametros para  El resultado de esta prueba para cada uno de los m´etodos estudiados se observa en la figura 2.3, as´ı como con la imagen que se ha utilizado para conseguirlo. El valor ´optimo de que hacen m´ınimo el numero de iteraciones es cualquiera entre 0,5y1,25 para todos los m´etodos, se puede obtener el valor de ⌧simplemente igualando a 1 la ecuaci´on 2.15. Tiempo de ejecuci´on Para los valores ´optimos de ⌧yelegidos en el apartado anterior, se quiere comprobar el tiempo de ejecuci´on para ver si es posible la ejecuci´on en tiempo real, ya que este depende de las iteraciones necesarias para converger. TV Iteraciones 50 100 150 200 250 Tiempo (ms) 11 17 25 31 38 Huber Iteraciones 25 50 100 150 200 250 Tiempo (ms) 710 16 23 30 35 TGV Iteraciones 150 200 250 Tiempo (ms) 59 76 94 Se puede observar que para TV y Huber es posible aumentar el numero de iteraciones hasta 200, y seguir´ıa ejecut´andose en tiempo real, si es que se trabaja a 30 17 Secci´on 2.5 2. Image Denoising frames por segundo. Para TGV el minimo numero de iteraciones para que converja es 150, es decir 59 milisegundos. Si se fuese m´as flexible en el tiempo real suponiendo que cada 100 ms fuera posible actualizar la imagen, se podr´ıa aumentar hasta 250 el n´umero de iteraciones y conseguir la ejecuci´on en tiempo real. Se ha implementado una una aplicaci´on que realiza este procedimiento, en el anexo D hay m´as informaci´on. 2.5.2. Par´ametros de ponderaci´on Para elegir los mejores par´ametros de ponderaci´on (TV, Huber, TGV), ↵1y↵2 (TGV) que regulan el compromiso entre la proximidad de la soluci´on a la imagen de entrada, y regularizaci´on del resultado, se ha dise˜nado el siguiente procedimiento. 1. Para cada algoritmo se realiza el n´umero de iteraciones calculado en la secci´on 2.5.1, es decir TV 50, Huber 25 y TGV 150. 2. Se obtiene para un rango de los par´ametros a estudio, por ejemplo k2 [0...n], la soluci´on uktras correr el algoritmo el n´umero de iteraciones comentado en el punto anterior. Puesto que se dispone de la imagen original sin ruido uGT (Ground truth) que aparece en la figura 2.3 a, se calcula la relaci´on se˜nal ruido SNR entre la soluci´on uky la imagen sin ruido uGT mediante la siguiente f´ormula. SNRk=10log 10 Pi,j(uGT i,j )2 Pi,j(uk i,j uGT i,j )2 Se elige el par´ametro que maximiza la SNR. Elecci´on de  Tanto TV como Huber como TGV dependen de este par´ametro. Se recuerda que pondera la importancia de la se˜nal de entrada (altos) frente a la regularizaci´on (bajos). Para el caso del TGV que depende de dos valores adicionales (↵1y↵2) se han fijado a ↵1=0,5y↵2=1,5. Para realizar este estudio se toma la imagen uGT y se a˜nade ruido gausiano de media 0 y =0,1 (Se recuerda que la imagen uGT esta normalizada entre 0 y 1). Aplicando el proceso explicado al inicio de la secci´on se obtiene la gr´afica de la figura 2.4 a. Para el TV el mejor es 12, para Huber 8 y para TGV 5. La SNR para la imagen con ruido (=0,1) y la SNR para cada una de las normas utilizando el valor ´optimo de se puede observar en la figura 2.5. 18 2. Image Denoising Secci´on 2.5 0 2 4 6 8 10 12 14 16 18 20 12 14 16 18 20 22 24 26 λ Decibelios TV Huber TGV (a) Ruido gausiano con =0,1 0 2 4 6 8 10 12 14 16 18 20 5 10 15 20 λ Decibelios TV Huber TGV (b) Ruido gausiano con =0,3 0 2 4 6 8 10 12 14 16 18 20 5 10 15 20 λ Decibelios TV Huber TGV (c) Ruido sal y pimienta en el 15 % de la imagen Figura 2.4: Valores ´optimos de para distintos tipos y valores de ruido. 19 Secci´on 2.5 2. Image Denoising (a) Imagen ruidosa snr = 14,4(b) TV snr = 24,4 (c) Huber snr = 23,7(d) TGV snr = 24,3 Figura 2.5: Prueba de comportamiento para los tres m´etodos estudiados con un ruido gausiano con =0,1 20 2. Image Denoising Secci´on 2.5 (a) Imagen ruidosa snr =6 (b) TV snr = 18,7 (c) Huber snr = 16,6(d) TGV snr = 19 Figura 2.6: Prueba de comportamiento para los tres m´etodos estudiados con un ruido gausiano con =0,3 Se observa que todas las normas realizan un gran trabajo. A pesar de que la norma TV logra una mejor SNR, el resultado de TGV es visualmente mejor ya que no se observa starcasing problem. Se estudia tambi´en el efecto de subir la desviaci´on est´andar del ruido a =0,3. Los valores ´optimos de para este caso se muestran en la figura 2.4 b, y son para TV es 6, para Huber 4, y para TGV 2. La SNR de la imagen con ruido y para cada norma con los valores ´optimos de en este caso se observan en la figura 2.6. Para este caso se observa que a pesar de que exista mejora entre la imagen con ruido y los resultados, estos no son ni de lejos tan buenos como en el caso anterior. Por ´ultimo se estudia el efecto del ruido de sal y pimienta. En la figura 2.4 c se muestran los resultados obtenidos para para cada una de las normas. Se observa que para TV es 6, para Huber 5 y para TGV 2. Se muestran en la figura 2.7 los resultados obtenidos de esta prueba. En este caso, como en los anteriores la mejor opci´on es TGV, pero aun as´ı, nunca 21 Secci´on 2.5 2. Image Denoising (a) Imagen ruidosa snr =8 (b) TV snr = 19,1 (c) Huber snr = 16,1(d) TGV snr = 19,3 Figura 2.7: Prueba de comportamiento para los tres m´etodos estudiados con un ruido sal y pimienta que afecta al 15 % de los p´ıxeles 22 2. Image Denoising Secci´on 2.5 Figura 2.8: Espectrograma resultante del an´alisis se˜nal ruido de ↵1y↵2,ambas entre 0 y 5, para =5 llega a suprimir por completo el ruido sal y pimienta, adem´as elimina muchos de los detalles de la imagen al intentarlo. La norma id´onea para este tipo de ruido hubiera sido la TV-L1, aunque para este proyecto no se ha estudiado. Elecci´on de ↵1y↵2 Para el caso del TGV, para un ruido gausiano de =0,1 y para el valor ´optimo de calculado en el apartado anterior (= 5), se estudia el efecto de modificar ↵1y↵2.Elresultadoeselqueapareceenlafigura2.8.Seobservaqueexisteun amplio rango de valores de ↵1y↵2para los que la se˜nal ruido es alta (franja roja), incluidos los valores ↵1=0,5y↵2=1,5 utilizados para las pruebas de la elecci´on de . Para el estudio de todos estos par´ametros se ha implementado una aplicaci´on cuyos detalles se encuentran en el anexo D. 23 Secci´on 2.5 2. Image Denoising 24 Cap´ıtulo 3 Fusi´on de mapas La parte central del proyecto es la aplicaci´on de los conceptos explicados anteriormente basados en Primal Dual para la fusi´on de mapas de profundidad. La base de todo el proceso que se va a explicar a continuaci´on son los mapas de profundidad. M´as concretamente un conjunto de ellos, sobre los que se realizar´an distintas operaciones para lograr una mejora en la calidad del mapa resultante. 3.1. Mapas de profundidad Como se ha dicho en la introducci´on del cap´ıtulo, los mapas de profundidad son el principal componente de todo este proceso. Un mapa de profundidad es simplemente una imagen que contiene informaci´on relacionada con la distancia de las superficies de los objetos de la escena desde un solo punto de vista. Estos mapas se pueden construir de mediante distintas t´ecnicas, una de las m´as actuales y que m´as se est´an utilizando es el sistema Kinect, que captura las distancias a la escena que tiene delante gracias a un sensor infrarrojo que hay junto a la c´amara. Por otro lado se pueden utilizar otras t´ecnicas para la construcci´on de mapas de profundidad a partir de im´agenes tomadas mediante una c´amara monocular como puede ser DTAM (Dense Tracking and Mapping in Real-Time) [RAND11]. En este caso se ha utilizado una t´ecnica similar a esta para el c´alculo de los mapas. 3.2. Mapas virtuales Para realizar la fusi´on de distintos mapas de profundidad deben de estar generados desde el mismo punto de vista. Para ello se ha de realizar un proceso que se detalla en la figura 3.1. En resumen, lo que hace el sistema es proyectar los datos contenidos en el mapa de profundidad que se quiere trasladar en un cubo 3D, y leer esos datos de vuelta 25 Secci´on 3.4 3. Fusi´on de mapas (a) Mapa estimado. snr= 17.5 (b) Resultado. snr=19.1 Figura 3.6: Resultado de la fusi´on de 12 mapas estimados. uno m´as parecido al ground truth, por lo que la SNR deber´ıa aumentar en este caso. Los resultados de esta prueba se muestran en la figura 3.6.Se confirma el aumento de la SNR tras la fusi´on de los mapas. Tanto en los bordes de la mesa como en los de la impresora se aprecia una mejora a simple vista. En la parte del monitor no corrige este efecto debido a que todos los mapas estimados tienen una mala estimaci´on para esta parte, por lo que la soluci´on ser´ıa mejorar el m´etodo de creaci´on de depth maps. Para comprobar que la fusi´on puede resolver este problema se introducen mapas de profundidad de una buena calidad dentro de la fusi´on el resultado deber´ıa mejorar en gran medida.En la figura 3.7(a) se muestra el mapa obtenido tras fusionar un 85 % de mapas estimados y un 15 % de mapas de mas calidad y la SNR del resultado, y en la figura 3.7(b) el resultado y el SNR para un 50 % de mapas estimados y un 50 % de mapas de mas calidad. La figura 3.8 muestra el efecto que tiene introducir mapas de mas calidad en la SNR. Como se esperaba se puede decir que la fusi´on mejora conforme la cantidad de mapas buenos introducidos aumenta. Otro efecto curioso es que la mejora depende del orden en el que se introducen debido al nivel de solapamiento respecto al depth map de referencia. Se observa con claridad que cuando se introduce el depth map numero 7 se produce una gran mejora como se observa en la pendiente de la gr´afica. 3.4.3. Medida de tiempos En cuanto a la medida de tiempos el inter´es se centra en dos procesos. Cuanto cuesta crear un solo mapa virtual, ya que que para una aplicaci´on real este tiempo esta espaciado y depende del n´umero de im´agenes en niveles de 32 3. Fusi´on de mapas Secci´on 3.4 (a) Resultado. snr=19.53 (b) Mapa estimado. snr= 21.8 Figura 3.7: Resultado de la fusi´on habiendo a˜nadido distinto porcentaje de mapas de mas calidad. 0 2 4 6 8 10 12 19 20 21 22 23 24 25 26 27 Mapas perfectos SNR Figura 3.8: gr´afica de mapas perfectos 33 Secci´on 3.4 3. Fusi´on de mapas gris necesarias para construir el depth map. Este tiempo es de una media de 200 milisegundos, ya que depende en cierto modo del mapa que se introduce en el cubo. Cabe destacar que esta prueba se ha realizado para una resoluci´on de 1024x1024x1024. Una bajada de esta resoluci´on reducir´ıa significativamente el tiempo de virtualizaci´on, por ejemplo para 512x512x512, el tiempo desciende hasta los 32 milisegundos (Tiempo real). Por otro lado es interesante saber el tiempo dedicado a la fusi´on en funci´on de el n´umero de mapas que se fusionan. En la figura 3.9 se observa este dato. 0 2 4 6 8 10 12 50 60 70 80 90 100 110 120 Número de mapas de profundidad Tiempo (ms) Figura 3.9: Tiempos de realizaci´on de la fusi´on en funci´on del n´umero de mapas implicados. Esta prueba se ha realizado para 100 iteraciones un valor que asegura un buen resultado. Estas pruebas, as´ı como todas las de la memoria se han realizado sobre una tarjeta gr´afica nVidia Tesla M2090. 34 Cap´ıtulo 4 Conclusiones A lo largo de proyecto se han ido cumpliendo cada uno de los objetivos planteados al inicio del mismo. A continuaci´on, se listan las aportaciones realizadas por parte de este PFC. Explicaci´on sencilla de los c´alculos matem´aticos involucrados en todo el proceso, desde el algoritmo utilizado para la resoluci´on del Primal Dual, hasta su aplicaci´on en concreto para m´etodos de denoising y fusi´on. Estudio detallado de todos los par´ametros implicados en el denioising, tanto de convergencia como de ponderaci´on. De este modo se tiene una idea m´as clara para futuros trabajos sobre las elecciones a tomar en este aspecto, ya que la documentaci´on existente no hace ´enfasis en este tema. Estudio de la posible aplicaci´on en tiempo real de los algoritmos de denoising. ´ Util para aplicaciones posteriores. Aplicaci´on del proceso de TGV-Fusion para im´agenes proyectivas y no solamente ortogr´aficas como se hab´ıa realizado hasta la fecha. Estudio de los resultados de la fusi´on sobre la influencia del punto de vista y de distintos valores y tipos de ruido. 4.1. Lineas Futuras En cuanto a posibles trabajos posteriores con relaci´on de este proyecto: Estudio de la posible mejora del proceso de fusi´on mediante aplicaci´on de distintos m´etodos de denosing despu´es de la virtualizaci´on de los mapas, justo antes de la fusi´on. 35 Secci´on 4.2 4. Conclusiones Realizaci´on de un promedio de los mapas directamente en el cubo que realiza la virtualizaci´on. En vez de transformar cada uno por separado. Comparar este resultado con el obtenido mediante la t´ecnica utilizada en este proyecto. Realizaci´on de un m´etodo variacional con los datos directamente cargados en el cubo de virtualizaci´on, y no a posteriori como se realiza en este caso. Estudio de la distribuci´on del trabajo enviado a la tarjeta gr´afica, para optimizarlo seg´un la arquitectura en la que se ejecuta. Por ´ultimo, hubiera sido interesante la integraci´on de la fusi´on realizada dentro de todo el proceso, desde la captura de im´agenes hasta la reconstrucci´on de la escena. Y estudiar la posible ejecuci´on en tiempo real de todo el proceso. 4.2. Valoraci´on personal Respecto a la valoraci´on personal, este proyecto ha supuesto todo un reto, ya que se ha tratado del proyecto m´as grande que he realizado, gracias a esto he aprendido mucho, muchas veces de los propios errores. Tambi´en a supuesto un reto en cuanto a las herramientas empleadas, ya que antes de realizarlo carec´ıa de experiencia en programaci´on en Cuda. Tampoco nunca hab´ıa realizado un documento en L A T EXlocualhahechomuchom´aslenta la etapa de escritura de la memoria. La realizaci´on de este proyecto me ha permitido conocer el mundo de la visi´on por computador, un ´area que a lo largo de la carrera no hab´ıa tocado demasiado, pero que siempre me hab´ıa llamado la atenci´on. Adem´as este proyecto se ha basado en la investigaci´on, un campo en mi opini´on muy interesante para cualquier ingeniero ya que es donde puedes aplicar todo tu ingenio, trabajar con tecnolog´ıas punteras, as´ı como estar al d´ıa de todo, gracias a las distintas publicaciones de centros de investigaci´on y universidades, en definitiva, el l´ımite te lo pones t´u mismo. En cuanto al desarrollo del proyecto han sido aplicados conceptos de muy diversas materias aprendidas durante los ´ultimos 5 a˜nos de carrera. Dentro de estos conceptos, por ejemplo los relacionados con la arquitectura de computadores, fueron muy ´utiles para comprender la arquitectura de la GPU, y de este modo programar de forma m´as eficiente. Por otro lado, se han aplicado multitud de conceptos relacionados con el c´alculo y el c´alculo num´erico. En definitiva me he dado cuenta de que toda la formaci´on recibida a lo largo de la carrera tiene una gran aplicaci´on. En cuanto a mi valoraci´on personal propiamente dicha, ha sido una experiencia muy buena, ya que el tema del proyecto me gustaba, aunque en ocasiones lo vi 36 4. Conclusiones Secci´on 4.2 demasiado te´orico, m´as tarde me di cuenta de que la teor´ıa matem´atica es la base de todos los problemas relacionados con este tema, y que es necesaria una buena comprensi´on de la misma para poder abarcar cualquiera de los problemas resueltos a lo largo del proyecto. 37 Secci´on 4.2 4. Conclusiones 38 Bibliograf´ıa [AWB] Gottfried Graber Thomas Pock Andreas Wendel, Michael Maurer and Horst Bischof. Dense reconstruction on-the-fly. [Bra00] G. Bradski. The OpenCV Library. Dr. Dobb’s Journal of Software Tools, 2000. [BV04] Elsevier Stephen Boyd and Lieven Vandenberghe. Convex Optimization. Cambridge University Press, Marzo 2004. [CP11] A. Chambolle and T. Pock. A first-order primal-dual algorithm for convex problems with applications to imaging. Journal of Mathematical Imaging and Vision,40(1),MayoMayo2011. [CUD10] CUDA by Example: An Introduction to General-Purpose GPU Programming. 2010. [LIRF92] S. Osher L. I. Rudin and E. Fatemi. Nonlinear total variation based noise removal algorithms. Proc. of the 11th annual Int. Conf. of the Center for Nonlinear Studies on Experimental mathematics : computational issues in nonlinear sciences, 1992. [NY04] Nesterov and Yurii. Introductory lectures on convex optimization: A basic course, volume 87. Springer, 2004. [RAND11] S. J. Lovegrove R. A. Newcombe and A. J. Davison. Dtam: Dense tracking and mapping in real-time. Proceedings of the 2011 International Conference on Computer Vision ICCV, 2011. [Roc97] R. T. Rockafellar. Convex Analysis. Princeton Math. Series, Princeton Univ. Press, 1997. [SI11] Otmar Hilliges David Molyneaux Richard Newcombe Pushmeet Kohli1 Jamie Shotton1 Steve Hodges1 Dustin Freeman1 5 Andrew Davison2 Andrew Fitzgibbon Shahram Izadi, David Kim. Kinectfusion: Realtime 3d reconstruction and interaction using a moving depth camera. 39 BIBLIOGRAF´ IA BIBLIOGRAF´ IA Proceedings of the 24th annual ACM symposium on User interface software and technology., 2011. [Sum10] Mark Summerfield. Advanced Qt Programming. Prentice Hall, 2010. [TPB11] Lukas Zebedin Thomas Pock and Horst Bischof. Tgv-fusion. 2011. 40 Ap´endice A Optimizaci´on convexa: Algoritmo de Primal Dual La optimizaci´on de funciones es una rama muy importante de estudio en las matem´aticas, en general, no es posible encontrar el m´ınimo o el m´aximo absoluto de un problema de optimizaci´on aunque s´ı existen, para determinados tipos de funciones, soluciones ´optimas o buenas heur´ısticas. [NY04] Afortunadamente en este proyecto nos encontramos en uno de esos casos, en los que las funciones a minimizar que nos vamos a encontrar, presentan una estructura especial que permite obtener la soluci´on global (´optima) del problema. La estructura gen´erica de los problemas que vamos a resolver, tanto en denoising como en TGV-fusion, es la siguiente: m´ın uF(Ku)+G(u) (A.1) Donde F y G son funciones convexas aunque no tienen por que ser diferenciables y K es un operador lineal. La soluci´on a esta optimizaci´on se basa en el algoritmo Primal Dual que como veremos m´as adelante, cuando se aplica a problemas de visi´on por computador, se puede acelerar considerablemnte usando tarjetas gr´aficas (GPUs) y programaci´on paralela en CUDA. A.1. Algoritmo Primal Dual Citando a [Roc97] ”La gran l´ınea divisoria en optimizaci´on no es entre problemas lineales y no lineales sino problemas convexos y no convexos”. Para un n´umero considerable de problemas convexos con cierta estructura es posible encontrar una soluci´on ´optima con buena precisi´on y en un tiempo razonable, independientemente de la inicializaci´on. En nuestro caso la expresi´on 2.1 est´a compuesta por dos funciones convexas que admite una soluci´on ´optima [CP11]. 41 Secci´on A.2 A. Optimizaci´on convexa: Algoritmo de Primal Dual −2−1.5 −1−0.5 0 0.5 1 1.5 2 −0.5 0 0.5 1 1.5 2 p>1 p <= 1 Figura A.6: Gr´afica ´util para comprender el dual de la funci´on regularizadora de Huber F⇤(p)=⇢↵p2 2|p|1 1|p|>1(A.8) A.2.3. TGV En el caso del TGV como en los m´etodos anteriores existe un cambio en la parte del regularizador, pero no se trata de un cambio de funci´on, en este caso introduce una funci´on nueva vque lo que pretende es dar cabida a las variaciones lineales del gradiente. Esto es, en superficies planas de la imagen en las que la variaci´on de color sea lineal no habr´a penalizaciones que le obliguen a generar peque˜nas zonas de igual color como en el TV, si no que permitir´a esta variaci´on gradual del color, en la figura A.7, se observa este hecho en funciones sencillas. Esta nueva variable vrepresenta la variaci´on del gradiente a lo largo de la imagen, por lo que como se ha dicho antes, minimizar esta funci´on tiende a generar ´areas de gradiente constate. En la figura se observan cuatro gr´aficas de funciones significativas para entender el problema, la primera gr´afica es la supuesta funci´on de entrada con ruido, la segunda la soluci´on que generar´ıa el m´etodo TV, como se ha dicho antes se observa que tiende a generar superficies planas, en las que no varia el gradiente. En la 48 A. Optimizaci´on convexa: Algoritmo de Primal Dual Secci´on A.2 0 5 10 15 20 25 30 0 2 4 6 U 0 5 10 15 20 25 30 0 5 Rof 0 5 10 15 20 25 30 0 5 TGV 0 5 10 15 20 25 30 −2 0 2 v Figura A.7: Comparaci´on de los distintos m´etodos en una funci´on sencilla tercera gr´afica se observa el resultado que resultar´ıa del m´etodo TGV, donde se conservan las pendientes como deber´ıa. Por ´ultimo se representa la funci´on v, que es la cual permite que el TGV se comporte de esa forma, porque el m´etodo minimiza la variaci´on total de esta funci´on como lo hacia TV para la propia imagen, por lo que generar´a en ella sectores planos que se corresponder´an variaciones constantes del gradiente en la imagen resultado. La parte del regularizador queda pues de la siguiente forma. |ru|TGV =↵1|ruv|+↵2|rv|(A.9) Al haber ahora dos t´erminos existir´an dos funciones a dualizar, la ventaja es que ambas son la funci´on valor absoluto, por lo que esas funciones son sencillas como se ha visto en el caso del TV. A continuaci´on se muestran ambas, como en casos anteriores pes el gradiente de u, y en este caso en particular qva a ser el gradiente de v. F⇤(p)=⇢0|p|↵1 1|p|>↵ 1 F⇤(q)=⇢0|q|↵2 1|q|>↵ 2 49 Secci´on A.3 A. Optimizaci´on convexa: Algoritmo de Primal Dual A.3. Resoluci´on del problema A partir del apartado anterior, podemos concluir que la ecuaci´on A.1 de la que se part´ıa, tras aplicar la transformaci´on de Lengendre-Fenchel para el t´ermino F(Ku) resulta la siguiente expresi´on. m´ın nsup p <Ku,p>+G(u)F⇤(p) (A.10) + m´ın nsup p <Ku,p>+G(u)sujetoa|p|1 Ahora la expresi´on ya es diferenciable, ya que la parte F(Ku) ha quedado transformada, y la parte G(u)eradiferenciabledeiniciocomosehasupuestoen el apartado de caracter´ısticas. Por lo que resulta un problema de minimizaci´on, maximizaci´on de funciones diferenciables. En este punto es posible aplicar m´etodos de primer orden (basados en que la funci´on es diferenciable al menos una vez) para la optimizaci´on. A.3.1. Optimizaci´on convexa de primer orden En primer lugar se va a definir la optimizaci´on mediante m´etodos de primer orden ya que es el algoritmo en que se basar´a la soluci´on del problema Primal-Dual para los problemas que se desean resolver. Interesa resolver el siguiente problema: Dada f:Rn! R, campo escalar una vez diferenciable con continuidad, encontrar ¯x2Rnsoluci´on de: m´ın x2Rnf(x) (A.11) Encontrar una soluci´on de este problema corresponde a encontrar un punto ¯xque satisfaga: f(¯x⇤)f(x)8x2Rn(A.12) El punto ¯x⇤se denomina m´ınimo absoluto de fsobre Rn Sabiendo que la funci´on a estudiar es convexa y diferenciable, si pensamos en funciones de una variable, habremos encontrado un m´ınimo cuando f0(x)=0, pero al extrapolar este problema a Rn, podemos demostrar de forma te´orica que nos encontramos en un m´ınimo cuando se cumple la siguiente condici´on necesaria y suficiente. [BV04] 50 A. Optimizaci´on convexa: Algoritmo de Primal Dual Secci´on A.3 ¯xes un punto estacionario de f rf(¯x)=~ 0() 2 6 6 6 6 4 @f(¯x) @x1 @f(¯x) @x2 . . . @f(¯x) @xn 3 7 7 7 7 5 =2 6 6 6 4 0 0 . . . 0 3 7 7 7 5 Acontinuaci´onsevaadescribirelm´etododelgradientedescendentequeser´ael que se utilizar´a con ciertos matices para resolver el problema. A.3.2. M´etodo del Gradiente El m´etodo del gradiente descendente es uno de los m´etodos sencillo de optimizaci´on de funciones diferenciables. Se trata de un m´etodo iterativo que dado un punto inicial busca la direcci´on de m´aximo descenso y elije un paso para acercarse a la soluci´on. Matem´aticamente lo se puede ver como que en cada iteraci´on se busca una direcci´on d2Rndada por la siguiente expresi´on. m´ın d2Rnf(x)dsujeto a kdk2=1 Con lo que resulta la siguiente expresi´on. d=rf(x) Es decir, la direcci´on local de m´aximo descenso es la direcci´on del gradiente negativo. En resumen el m´etodo del gradiente descendente obedece a la siguiente expresi´on donde en cada paso hay que calcular tau usando un algoritmo de line search. xn+1 =xn⌧rf(xn) (A.13) A.3.3. Otros m´etodos de optimizaci´on Como se ha dicho anteriormente el m´etodo del gradiente es el m´etodo m´as sencillo, y como tal en condiciones normales es de los m´as lentos en alcanzar la soluci´on, existen muchos otros m´etodos a priori m´as r´apidos y eficientes. Por ejemplo uno de los m´as conocidos es el m´etodo de Newton,quesebasaenelmismo concepto, pero realiza una aproximaci´on cuadr´atica a la hora de elegir la direcci´on de descenso. 51 Secci´on A.4 A. Optimizaci´on convexa: Algoritmo de Primal Dual La clave para elegir el de gradiente y no el de newton es que en la pr´actica dado que cada componente del gradiente se puede calcular de forma independiente, se puede paralelizar totalmente, dada la independencia de los datos. Si se piensa en problemas relacionados con la imagen, para hilo de ejecuci´on har´a las operaciones para cada pixel, mientras que en el caso de m´etodos de aproximaci´on cuadr´atica los datos no son independientes y resultan operaciones m´as complejas que involucran a m´as datos, por lo que no aprovechar´ıamos la potencia de la programaci´on gr´afica. A.4. Proximal map A pesar de que es posible la utilizaci´on del m´etodo del gradiente para la optimizaci´on de la expresi´on A.10, se va a utilizar un m´etodo llamado proximal map.Se trata de un m´etodo muy similar al del gradiente descendente pero incluye varias caracter´ısticas que lo hace perfecto para la resoluci´on del Primal Dual. Como se ha dicho en anteriores ocasiones se parte de la siguiente expresi´on general. m´ın xsup y <Kx,y>+G(x)sujetoa|y|1 + m´ın nsup p C(x, y)sujetoa|y|1 Se muestra en funci´on de una funci´on de energ´ıa standard para facilitar la comprensi´on. Acontinuaci´onsemuestralaestructuradeestem´etodo. 8 < : yn+1 =(I+@F⇤)1(yn+K¯xn)(1) xn+1 =(I+⌧@G)1(xn⌧K⇤yn+1)(2) ¯xn+1 =xn+1 +✓(xn+1 xn)(3) (A.14) Este es el esquema de la resoluci´on del problema Primal Dual, en el se realiza una maximizaci´on y una minizaci´on para cada una de las iteraciones. El la parte 1 del proceso se realiza la maximizaci´on. Como se hac´ıa para el m´etodo del gradiente descendente se deriva en funci´on de y. ˜yn+1 =yn+ryC(xn,yn+1) En este punto se observa la primera de las ventajas del proximal map,eval´uaay en el punto siguiente del que se encuentra en la iteraci´on actual, lo que acelerar´a el proceso. 52 A. Optimizaci´on convexa: Algoritmo de Primal Dual Secci´on A.4 La siguiente ventaja viene dada por el operador (I+@F⇤)1,sulaborconsiste en realizar una proyecci´on de modo que cumpla la restricci´on que impone la funci´on dual, de lo que resulta lo siguiente. yn+1 =˜yn m´ax(1,˜yn+1) La parte 2 del proceso realiza la minimizaci´on en x,estecasoenalgomas sencillo que el anterior ya que al realizar la primera derivada desaparece la funci´on dual as´ı como su restricci´on por lo que el operador (I+⌧@G)1es trivial. En este caso sigue apareciendo la ventaja de que la fuci´on es evaluada en el punto siguiente al realizar los c´alculos por lo que la derivada tendr´a la siguiente forma. xn+1 =xn⌧rxC(xn+1,yn+1) En este caso ytambi´en se eval´ua en el punto siguiente ya que puesto que este resultado es conocido de antemano. En cuanto al tercer elemento no tiene una importancia vital para el problema, se trata de una relajaci´on para lograr que el m´etodo converja con mayor rapidez. En cuanto al paso, para esta forma de funci´on de energ´ıa en concreto se puede demostrar que tomando ⌧yde forma adecuada se asegura la convergencia. Por lo que no ser´a necesario realizar el c´alculo del paso en cada iteraci´on como si que lo ser´ıa con el m´etodo del gradiente. En resumen el m´etodo proximal map utilizado se basa en el m´etodo del gradiente descendente y ascendente evaluados en el punto siguiente y con un paso conocido de antemano haciendo que en realidad la convergencia sea mas r´apida y los c´alculos m´as sencillos que en un m´etodo de gradiente tradicional. A.4.1. Elecci´on de par´ametros Una particularidad de la funci´on de energ´ıa de este proyecto es que se puede asegurar que con una correcta elecci´on de los par´ametros el m´etodo de gradiente para este problema siempre converge. Esta correcta elecci´on se refiere a que los distintos pasos ⌧ydeben cumplir la siguiente expresi´on. ⌧L21 (A.15) Siendo L=kKkes decir, la norma de K.Enelart´ıculoenelquesebasaesta parte del proyecto [CP11], escoge ambas con valor aproximado de 0,3536 ya que se puede demostrar que kKk=p8, por lo que haciendo en la expresi´on anterior ⌧=el resultado es el comentado anteriormente. 53 Secci´on A.4 A. Optimizaci´on convexa: Algoritmo de Primal Dual 54 Ap´endice B Gesti´on del proyecto En este cap´ıtulo del anexo se va a detallar todo el proceso que se ha ido realizando a lo largo del proyecto. B.1. Metodolog´ıa Dadas las caracter´ısticas del proyecto no se aplic´o ninguna metodolog´ıa concreta para su realizaci´on. A pesar de esto el desarrollo del proyecto se puede dividir en distintas fases. Aprendizaje En principio se realiz´o una etapa de aprendizaje, por un lado del algoritmo Primal Dual, y por otro del primero de los problemas de denoising, ya que se trata de el m´as sencillo, pero en el que se aplican todos los conceptos que se necesitar´ıan para los siguientes problemas. Alavezqueseaprend´ıanlosconceptoste´oricosymatem´aticos,secomenz´oa programar los distintos problemas en distintas tecnolog´ıas gradualmente hasta llegar a la el ´ultimo escal´on, una vez all´ı, se fueron implementando los dem´as m´etodos de denoising as´ı como la fusi´on. Se detallan a continuaci´on todos los pasos llevados a cabo. Matlab Como se ha dicho antes se tom´o el m´as sencillo de los problemas de denoising, el TV-Rof, para ello se estudi´o el problema propuesto en el Pock [CP11], y se program´o de forma literal en primer lugar en lenguaje matlab,yaquemediantesu uso las operaciones complejas con matrices se tornan muy sencillas. Por lo que no se tuvo en este primer acercamiento una visi´on a nivel de pixel, y por tanto era 55 Secci´on B.2 B. Gesti´on del proyecto imposible una interpretaci´on paralela del m´etodo. Aun as´ı esta etapa como se ha dicho ayud´o a la compresi´on del problema y del algoritmo. C++ Una vez que el problema estaba programado y funcionando en matlab,secomenz´o a implementar en C++. La principal novedad, y punto importante de C++frenteamatlab,huboque cambiar la concepci´on del problema, hab´ıa que pensar a nivel de pixel. El paso de pensar en matrices a pensar en p´ıxeles, fue lo mas complejo de esta etapa pero a la vez lo m´as interesante ya que se ve´ıa claramente que se pod´ıa paralelizar de forma sencilla. CUDA El siguiente paso, y ´ultimo, fue la implementaci´on en CUDA. La principal ventaja del uso de GPU’s es, como se comenta en la memoria, la gran paralelizaci´on que se logra, ya que las funciones primal y dual se ejecutan una vez por cada pixel, yencadaunodeestoskernels 1se realizan operaciones sencillas en las que las GPU’s marcan la gran diferencia. Una vez se llegados a este punto, se fueron realizando los siguientes algoritmos, despu´es de tenerlos todos implementados, se inici´o la realizaci´on de las aplicaciones que les daban soporte a estos m´etodos. M´as adelante en los anexos correspondientes se detalla el proceso seguido para la implementaci´on de cada una de las interfaces. B.2. Tecnolog´ıas empleadas La elecci´on de las tecnolog´ıas empleadas se bas´o principalmente en las directrices dadas por el director del proyecto, ya que todos los trabajos en ese grupo de investigaci´on se basaban en las mismas tecnolog´ıas por lo que ten´ıan experiencia en posibles problemas que pudieran surgir. La principal tecnolog´ıa empleada, y que marc´o la diferencia es Cuda/C ++ utilizada para el backend de la aplicaci´on, es decir, la programaci´on de todos los problemas propuestos en los distintos Papers. Como se ha comentado en otras ocasiones a lo largo de la memoria, gracias a la independencia de las operaciones para cada pixel en el c´alculo del Primal Dual, se logra una aceleraci´on, y una gran mejora en el tiempo de ejecuci´on respecto a versiones en C++omatlab,porejemploparamatlab, solo dos iteraciones, le 1kernel: Funci´on que se ejecuta en la gpu de forma paralela 56 B. Gesti´on del proyecto Secci´on B.2 costaba casi treinta minutos aunque la mayor parte de este tiempo estaba creando una matriz para realizar el c´alculo del gradiente. En el caso de C++ sin paralelizar para unas 200 iteraciones tardaba sobre 4 segundos. Por ´ultimo en la versi´on m´as r´apida de cuda llega a hacer las 200 iteraciones en treinta milisegundos, es decir, en tiempo real. A continuaci´on se explica de forma resumida que es Cuda,yquecaracter´ısticas tienen la arquitectura Nvidia. Cuda es una plataforma dise˜nada conjuntamente a nivel software y hardware para aprovechar la potencia de una GPU en aplicaciones de proposito general, permite una comunicaci´on sencilla entre CPU yGPU, y su hardware esta orientado al c´alculo paralelo. El paradigma de programaci´on de cuda consiste en programar ciertas funciones a paralelizar tal y como se har´ıa si se quisiese ejecutar en un solo hilo, y el sistema se encarga autom´aticamente de ejecutar esta funci´on paralelamente en tantos procesadores de la GPU como se requiera. En la figura B.1 se puede observar una visi´on global de la arquitectura Fermi, la tercera generaci´on del hardware paralelo de Nvidia,enellasedistinguena primera vista una serie de multiprocesadores, y los distintos niveles de memoria. Cada uno de los multiprocesadores esta compuesto en este caso de 32 n´ucleos que podr´an ejecutarse en paralelo, as´ı como memoria cach´e local y distintas herramientas para planificar la ejecuci´on de los hilos. En cuanto a la programaci´on en Cuda como se ha comentado antes, es muy parecida a la programaci´on secuencial en CPU,lanovedadest´aenelcompilador nvcc que es el que selecciona las partes que han sido programadas para la gpu (Funciones kernel) y las compila para su ejecuci´on ella. Las interacciones de datos entre la CPU ylaGPU se realizan mediante un bus PCI. Otro concepto importante para la programaci´on en cuda es el uso de distintos tipos de memoria, que ser´an ´utiles para acelerar el proceso. En la figura B.2 se observa la jerarqu´ıa entre los distintos niveles de memoria. Para comprender esto hay que tener claro la forma que cuda divide el trabajo para enviarlo a los distintos n´ucleos de la GPU.Cuda estructura el trabajo en varios niveles. El hilo o thread es el nivel m´as bajo, cada uno de ellos se ejecutar´a en un n´ucleo. Estos threads se agrupan en bloques generalmente de 32. Cada uno de estos bloques se ejecutar´a exclusivamente en un multiprocesador. Por ´ultimo el nivel m´as alto se conoce como grid, ser´a quien agrupe todos los bloques de una ejecuci´on. De todas las memorias de la figura se usar´an a lo largo del proyecto las que son visibles para todos los multiprocesadores que trabajan con el mismo grid, es decir: 57 C. Programaci´on de Primal Dual 5//Declaracion de los parametros del metodo 6float gamma = 0 . 7 f ⇤lambda ; 7float theta ; 8 9//Reserva espacio en gpu para los datos necesarios 10 11 cudaMalloc (...) ; 12 ... 13 14 // Inicializacion de los datos 15 16 cudaMemcpy( . . . , cudaMemcpyHostToDevice ) ; 17 cudaMemset ( . . . ) ; 18 19 //Carga de las texturas 20 cudaBindTexture2D ( . . . ) ; 21 ... 22 23 //Parametros de las dimensiones de trabajo 24 dim3 nThreads (NTHREADS BLOCK, NTHREADS BLOCK) ; 25 dim3 nBlocks (divUp( Inoisy . cols , nThreads . x) , divUp( Inoisy . rows , nThreads .y)) ; 26 27 //Bucle de ejecucion 28 for (int iter = 0; iter <iteraciones ; ++iter){ 29 30 updateDualGPU <<< nBlocks , nThreads >>> (...) ; 31 32 theta = 1.0 f/sqrtf (1.0 f+2⇤gamma⇤tau) ; 33 34 updatePrimalGPU <<< nBlocks , nThreads >>>(...) ; 35 36 tau = theta⇤tau ; 37 sigma = sigma/theta ; 38 } 39 40 cudaMemcpy( I r e s u l t . data , gpu u, nbytes, cudaMemcpyDeviceToHost); 41 42 //Liberacion de las memorias de textura 43 cudaUnbindTexture ( . . . ) ; 44 45 //Liberacion la memoria en la gpu 46 cudaFree (..) ; 47 48 return Iresult ; 49 } 64 C. Programaci´on de Primal Dual El c´odigo anterior es bastante autoexplicativo, en ´el principalmente hay tres zonas. La primera de declaraci´on de variables ya carga de datos, lo ´unico a destacar en esta parte son las instrucciones finales donde inicializa los par´ametros para la llamada al kernel. En resumen lo que hace en esta parte es asignar en primer lugar el n´umero de threads que ser´an generados para cada bloque, para ello declara una malla de 16*16 threads, despu´es calcula el n´umero de bloques que ser´an necesarios para recorrer toda la imagen. M´as tarde en la llamada a la funci´on kernel se le indica estos dos valores, para que se realice la paralelizaci´on deseada. En la parte central esta el bucle de minimizaci´on. En cuanto a este bucle solo ejecuta la actualizaci´on del primal y del dual una vez por cada iteraci´on, indicando los bloques y threads necesarios para una m´axima paralelizaci´on como se acaba de explicar. La estructura solo cambia ligeramente en el caso de m´etodo Huber,yaque como se ha dicho con anterioridad este m´etodo admite una variaci´on que hace que sea mas r´apida la convergencia, lo ´unico que cambia para este caso es que no se actualizan los par´ametros en cada iteraci´on si no que se inicializan al principio de forma diferente a la anterior, y mantienen el mismo valor durante todo el proceso. Por ´ultimo se encuentra la zona de liberar memoria y devolver resultado. La estructura de las funciones kernel, para la actualizaci´on Primal Dual es la siguiente 1global void update primalOdual GPU ( . . . ) { 2 3int i=blockDim.x⇤blockIdx .x + threadIdx .x; 4int j=blockDim.y⇤blockIdx .y + threadIdx .y; 5 6if(i<width && j<height ){ 7//Coordenadas en memoria de textura 8float tx=i+0.5f ; 9float ty=j+0.5f ; 10 /⇤ 11 Realizacion de los calculos 12 ⇤/ 13 } 14 } 65 C. Programaci´on de Primal Dual En primer lugar lo que hace es averiguar para que pixel se esta ejecutando el kernel, una vez hecho esto se comprueba que ese supuesto pixel esta dentro de la imagen, ya que cabe la posibilidad que se hayan ejecutado mas hilos que p´ıxeles tiene la imagen. Tras esto calcula las coordenadas tx yty,queser´anlas coordenadas donde leer´a en la memoria de textura, esto es porque la memoria textura lee la imagen como algo constante, y no discreto como la memoria global, por lo que realiza una interpolaci´on, y al sumarle ese 0.5f se desplaza al centro del pixel en cuesti´on pudiendo as´ı comprobar realmente el valor deseado. En cuanto a los c´alculos para cada uno de los kernel de cada uno de los m´etodos es literalmente programar los resultados de las deducciones matem´aticas que se encuentran en el cap´ıtulo de denoising o de fusi´on en la memoria. 66 Ap´endice D Interfaces D.1. Interfaz de prueba de par´ametros El origen de esta interfaz fue improvisado, ya que no se pensaba hacer de antemano, pero surgi´o de la necesidad, ya que en los art´ıculos en los que se basaba el proyecto eran algo ambiguos con el tema de los par´ametros que hab´ıa que introducir. D.1.1. Requisitos En cuanto a requisitos, se definieron sencillamente los siguientes: Se introducir´a una imagen desde archivo para las pruebas, que ser´a mostrada en la interfaz. Existir´a un apartado para introducir ruido en la imagen, una vez introducido, se deber´a mostrar en la interfaz. Se tendr´a la opci´on de elegir entre cada uno de los algoritmos de denoising del proyecto. Seg´un el algoritmo elegido se mostrar´a la posibilidad de hacer unas pruebas u otras. Se podr´an seleccionar los rangos en los que hacer las pruebas y el valor de otros par´ametros que afecten al algoritmo aunque la prueba no se centre en ellos. Se mostrar´an los resultados en una gr´afica embebida en la interfaz. Se mostrar´a tambi´en al final de la prueba la imagen resultado. 67 Secci´on D.1 D. Interfaces Figura D.1: Aspecto de la interfaz de prueba de par´ametros al inicio de su ejecuci´on. Las pruebas se realizaran en un thread aparte para que la interfaz no quede bloqueada durante la prueba. Se mostrar´a con una barra de carga o similar el proceso de la prueba para tener noci´on del tiempo que va a necesitar. D.1.2. Resultado de la interfaz Se muestra a continuaci´on el dise˜no de la interfaz, por un lado, el estado inicial en la figura D.1, y por otro un esquema de la interfaz completa, con todos los elementos visibles para facilitar la explicaci´on en la figura D.2. 1. Se trata de el widget que mostrar´a en cada momento la imagen con la que se este trabajando, al principio mostrar´a la imagen que se cargue, y una vez introducido el ruido, mostrar´a la imagen ruidosa. 2. Es la zona de introducir ruido en la imagen, en ella se puede elegir el tipo de ruido que se desea mediante el combobox,as´ıcomosuvalor. 3. En esta zona simplemente est´a el bot´on para cargar la imagen desde archivo, yelcomboBox para la elecci´on del algoritmo sobre el que se quiere realizar la prueba. 4. Esta es la zona que m´as espacio ocupa en la interfaz, pero tambi´en es la m´as importante. En ella aparecen dos zonas diferenciadas. 68 D. Interfaces Secci´on D.1 Figura D.2: Esquema de todos los elementos de la interfaz de par´ametros. a) Se trata de la zona de introducci´on de los intervalos a muestrear de las distintas variables as´ı como, por un lado par´ametros relacionados con las pruebas, y por otro par´ametros relacionados con el m´etodo. M´as adelante se detallar´an todos ellos. b)En esta zona representan los resultados que se han obtenido en la prueba, cabe destacar que para la mayor´ıa de las pruebas esta representaci´on es en tiempo real gracias a que la interfaz y las pruebas se ejecutan threads independientes. Tambi´en existe en esta zona la opci´on de hacer HoldOn,esdecir,la opci´on de que se conserven los datos en la gr´afica de una prueba, al realizar la siguiente, ´util para poder comparar los resultados con otros m´etodos, o elegir otro rango de acci´on de la prueba sin que tengas que repetir parte de ella. 5. Se trata simplemente del bot´on de ejecutar que desata todas las pruebas. 6. Esta zona solo se muestra despu´es de que se haya realizado una prueba, en ella se muestra el antes y el despu´es de la imagen ruidosa, siendo el despu´es el mejor resultado conseguido por todas las muestras ejecutadas en la prueba. El aspecto de esta sencilla ventana emergente se muestra en la figura D.3. 7. Se trata de la barra de carga, calcula el n´umero de muestras procesadas para dar informaci´on sobre el estado de toda la prueba. 69 Secci´on D.1 D. Interfaces Figura D.3: Aspecto de la ventana emergente que se genera tras pulsar el bot´on mostrar resultado despu´es de ejecutar una prueba. D.1.3. Detalles de utilidad Se van a detallar a continuaci´on que pruebas se pueden realizar para cada algoritmo, no de forma cualitativa como se ha hecho en la memoria, sino que se van a describir los par´ametros requeridos para cada prueba as´ı como los resultados que genera. Prueba ⌧y Los par´ametros que se requieren para realizar esta prueba, son, por un lado el rango de valores de ⌧oen los que se va a realizar la prueba y por otro lado los par´ametros relacionados con las prueba, como el n´umero de iteraciones a realizar al principio para generar una soluci´on que haya convergido y las iteraciones m´aximas a realizar a lo largo de la prueba, ya que si este n´umero no est´a acotado, para una mala elecci´on de los par´ametros puede tardar mucho a converger y la prueba se eternizar´ıa. Por ´ultimo hay que introducir tambi´en el valor de , lo que ofrece un mayor rango de posibilidades en las pruebas. Cabe destacar que la interfaz es ´unica para pruebas de ⌧yde, por lo que existe un radioButton que da la opci´on de elegir una u otra, tampoco ser´ıa necesario que hubiera dos pruebas diferenciadas, ya que una depende de la otra. En este caso solo se mide la velocidad de convergencia por lo que no existe una imagen resultado mejor que las dem´as por lo que en la ventana de mostrar el resultado se mostrar´a la soluci´on sobreiterada. 70 D. Interfaces Secci´on D.1 Figura D.4: Pseudodiagrama de clases de parte de la aplicaci´on de prueba de par´ametros. Prueba de  Esta prueba est´a disponible para todos los algoritmos, es bastante simple, s´olo es necesario introducir el intervalo de datos que se desean probar, y las iteraciones que se han de realizar para cada muestra. Prueba de ↵1y↵2 Esta prueba es exclusiva del algoritmo TGV,eslamascomplejadelasrealizadas ya que estas dos variables no dependen de s´ı mismas, por lo que el n´umero de muestras se multiplica. Para esta prueba hay que introducir los intervalos en los que realizar la prueba, as´ı como las iteraciones a realizar por muestra, y el valor de . El resultado de esta prueba es distinto a todos los anteriores ya que ahora el valor del an´alisis se˜nal ruido a realizar depende de dos variables, por lo que el resultado se muestra en un espectrograma. 71 Secci´on D.1 D. Interfaces D.1.4. Detalles de implementaci´on En este apartado se va a detallar lo m´as representativo de la parte de la implementaci´on tanto lo relacionado con la interfaz en QT como la parte de las pruebas en C++. En toda la aplicaci´on se pueden diferenciar tres clases principales, una es la que soporta la interfaz (MainWindow), otra se utilizar´a para ejecutar todas las pruebas (Lanzador), y otra en la que est´an programados los distintos algoritmos de denoising (Normas). En la figura D.4 se muestra un pseudodiagrama de clases en el que se muestra la jerarqu´ıa de las dos primeras clases. Se observa que todo esta heredado de la clase QObject que es la principal clase de QT,graciasaellaporejemplosepuedenhacerusodelossignalsyslots. En este punto el diagrama se divide en dos ramas, una llega ya a la clase lanzador, y otra a QWidget,esta´ultimaclaseser´alaquedoteasushijosde distintas propiedades para su representaci´on en la interfaz. Tras ella se encuentra la clase QMainWindow que simplemente es una distribuci´on predefinida en la interfaz, con posibilidad de crear de forma sencilla barras de herramientas, men´us, barras de estado, etc. Por ´ultimo al final de esta rama se encuentra la clase MainWindow,estaesla propia interfaz rellena ya con todos los elementos. A continuaci´on se detalla las distintas clases principales nombradas anteriormente. MainWindow Es la clase principal de la aplicaci´on maneja todo lo relacionado con la interfaz, esto queda claro al observar los distintos slots que la componen. Solo con leer el nombre se observa que cada uno de ellos atiende a una de las posibles acciones que se pueden desatar en la interfaz, cargar una imagen (loadImage()), introducirle el ruido (gereraRuido(), ejecutar las pruebas (Ejecutar()), etc. Uno de los requisitos m´as laboriosos fue el hecho de que las pruebas se ejecutaran en otro thread para que la interfaz no quedara bloqueada durante ese tiempo, para ello se construy´o la clase Lanzador, que esta asociada a la variable lanzadera de clase QThread.El funcionamiento de este sistema es el siguiente: 1lanzadera=new QThread() ; 2lanzador .moveToThread(lanzadera) ; Se declara el Thread y se indica que la instancia de la clase Lanzador se ejecute en ´el. 72 D. Interfaces Secci´on D.1 1connect(&lanzador , SIGNAL(done()) , this ,SLOT(finEjecucion())); 2connect(this ,SIGNAL(TVtauSigma(float ,float ,float ,int ,int ,int , float)), &lanzador , SLOT(TVtauSigma(float ,float ,float ,int ,int , int ,float))); En las l´ıneas anteriores se indica que cuando la clase lanzador emita la se˜nal done,esdecir,afinalizadolaspruebas,seejecuteelslotfinEjecucion(). La segunda l´ınea le indica que cuando la interfaz emita la se˜nal TV tausigma con todos los par´ametros, en el lanzador se ejecutar´a el slot con el mismo nombre, quien realizar´a la prueba. Esta segunda l´ınea se repite para todas las pruebas que se puede realizar. Todas estas interaciones entre las clases Mainwindow y Lanzador,est´anrepresentadaseneldiagramaporlasflechasazules. 1lanzadera>start (); 2emit TVtauSigma(/⇤parametros de la prueba ⇤/); Ahora en el momento que quiere ejecutar una prueba en concreto se inicializa el thread y se emite la se˜nal que capturar´a la el lanzador, quien ejecutara la prueba, adem´as lo har´a en un thread paralelo, por lo que la interfaz no quedar´a bloqueada esperando el c´alculo. Una vez realizada la prueba y le´ıdos los resultados se ejecutar´a lo siguiente. 1lanzadera>quit () ; Con lo que se detendr´a el thread a la espera de la siguiente ejecuci´on. Al igual la interfaz emite se˜nales con los valores de una prueba para que la realice el lanzador, este devuelve los valores de la prueba, estas se˜nales ser´an detalladas en el siguiente apartado donde se define la clase Lanzador. Dentro de la interfaz tambi´en cabe destacar la inclusi´on de widgets1de la librer´ıa qwt, que se trata de una librer´ıa que incluye nuevos componentes y utilidades para aplicaciones t´ecnicas, en este caso han sido usados dos de estos widgets,ambos para la representaci´on de datos, uno m´as sencillo el qwtPlot que se encarga simplemente de representar en una gr´afica 2D los datos generados en las pruebas, y otro m´as complejo que realiza una representaci´on 3D, mediante un espectrograma, ´util para la prueba de ↵1y↵2. 1Widget: D´ıcese de cada uno de los elementos de la interfaz 73 Secci´on D.2 D. Interfaces Figura D.6: Esquema de las distintas zonas de la interfaz de visualizaci´on. Aunque pudieran parecer iguales el funcionamiento de estos dos threads es muy distinto, para entenderlo hay que observar la figura D.7. Como se puede observa la estructura es muy parecida a la de la aplicaci´on de prueba de par´ametros con el a˜nadido del thread de lectura de la c´amara. Antes de la implementaci´on de este thread, se barajaron dos posibilidades. Por un lado, se baraj´o el hecho de hacerlo de igual modo que la clase Lazador, es decir, que la interfaz cada cierto tiempo ejecutar´ıa el thread para conseguir un frame de v´ıdeo y trabajar con ´el. El principal defecto de esta t´ecnica es que se ha de inicializar el thread a cada frame que pide la interfaz, y esto sucede con mucha frecuencia unas 30 veces por segundo. La otra opci´on fue el programar una clase que heredara directamente de QThread,yqueseejecutaracadavezquelaaplicaci´onsepusieraenmodo v´ıdeo, y durante este tiempo lea frames de la c´amara. Esta ´ultima fue la opci´on elegida como se puede observar en la figura D.7, su funcionamiento en resumen consiste en que cada vez que el usuario indica el modo v´ıdeo se activa el thread de la siguiente forma. 80 D. Interfaces Secci´on D.2 Figura D.7: Pseudodiagrama de clases de la aplicaci´on de visualizaci´on 1camera . start () ; 2 3timer = new QTimer( this); 4connect(timer ,SIGNAL(timeout()) , this ,SLOT(leerCamara())); 5timer>start(15); // se puede poner 0 En primer lugar activa el thread simplemente llamando la funci´on start(), tras esto se activa un timeout que leer´a una imagen del thread cada 15 milisegundos en este caso. Esta lectura se realiza de la siguiente forma. 1void MainWindow : : leerCamara () { 2Candado . lock () ; 3imgOriginal = camera.getImagen() ; 4Candado . unlock () ; 5imgRuido=addNoise ( imgOriginal ,0.0 , valorRuido ) ; 6ui>mostrar original>setImage(imgRuido) ; 7ui>mostrar original>updateGL() ; 8} Se llama a la funci´on getImagen(),quedevuelvelaimagenqueenesemomento 81 Secci´on D.2 D. Interfaces captura la c´amara. Una vez capturada cuando el usuario pulsa ejecutar se realiza el denosing de esa imagen, si al finalizar la ejecuci´on no se ha salido del modo v´ıdeo, se vuelve a realizar el denoising de la ´ultima imagen que la interfaz haya le´ıdo del thread de captura. Una vez se haya salido del modo de v´ıdeo, por ejemplo para cargar una imagen desde archivo, o directamente para cerrar la aplicaci´on, se detiene la ejecuci´on del thread de captura de la siguiente forma. 1timerDenoising>stop() ; 2timer>stop() ; 3camera . stop () ; Adem´as de detener el thread se detiene el timeout que se hab´ıa activado para la captura de las im´agenes, por lo que la interfaz no pedir´a al thread m´as im´agenes, lo que producir´ıa un error. A continuaci´on se resume la implementaci´on de las distintas clases implicadas. MainWindow Poco mas se puede de decir de esta clase, la ´unica miga que pudiera tener es el tratamiento de los distintos threads, de la que ya se ha hablado. M´as all´a de eso su estructura y funcionamiento de esta clase es similar que en la aplicaci´on de prueba de par´ametros. Lanzador El funcionamiento de esta clase es similar que en la aplicaci´on de prueba de par´ametros. Capture En cuanto a la implementaci´on de la clase Capture,esunaclasequehereda de QThread,porloquetieneunaspropiedadesespeciales,supartecentralse encuentra en la funci´on run(), que es la que se ejecuta autom´aticamente cuando el thread se inicia. Se muestra a continuaci´on esta funci´on run() de forma simplificada. 1cv : : VideoCapture capture (0) ; 2while (/⇤Fin del thread⇤/){ 3capture >> frame ; 82 D. Interfaces Secci´on D.3 4cv : : cvtColor(frame ,imGrey ,CVRGB2GRAY) ; 5Candado . lock () ; 6imagenParaInterfaz=imGrey; 7Candado . unlock () ; 8} Al principio inicia el canal de captura, tras ello comienza el bucle de captura, dentro de ´el captura un frame, convierte a escala de grises, y almacena el resultado para que la interfaz lo lea cuando lo desee. Para evitar que la interfaz lea la imagen justo cuando esta clase la esta escribiendo, se ha colocado un mutex en el momento en el escribe la imagen, y cuando la lee en la interfaz. Normas Esta clase es sencillamente una reuni´on de todos los m´etodos de denosing del proyecto. Para mas informaci´on acerca de la implementaci´on de esta clase dir´ıjase al anexo correspondiente. D.3. Interfaz de fusi´on Esta interfaz es el soporte de todo el proceso de fusi´on de mapas explicado en la memoria, desde la elecci´on de los distintos mapas de profundidad y elecci´on del ruido pasando por la virtualizaci´on hasta la realizaci´on de la fusi´on. D.3.1. Requisitos A continuaci´on se muestran una serie de requisitos generados antes de realizar la aplicaci´on, con directrices de su funcionamiento. Se deber´an introducir los mapas profundidad en la aplicaci´on y simult´aneamente las posiciones y par´ametros de las c´amaras de referencia. Ser´a posible la visualizaci´on de cada uno de los mapas una vez introducidos. Se deber´a de poder introducir ruido gausiano as´ı como datos espurios a cada mapa de profundidad por separado. En todo momento se tendr´a la posibilidad de utilizar para la virtualizaci´on y posterior fusi´on los mapas estimados, o los mapas perfectos (ground truth) indistintamente. 83 Secci´on D.3 D. Interfaces Figura D.8: Aspecto inicial de la interfaz de fusi´on habiendo cargado ya un lote de mapas de profundidad Una vez haya mapas cargados se habilitar´a la opci´on de realizar la virtualizaci´on. Antes de realizar la virtualizaci´on se elegir´an las resoluciones de los distintos ejes del cubo. Una vez virtualizados se podr´an consultar los mapas resultantes. Una vez realizada la fusi´on se mostrar´a el resultado as´ı como el resultado del an´alisis se˜nal ruido respecto del mapa en ground truth visto desde la posici´on sobre la que se encuentra el resultado. La fusi´on de los mapas se realizar´a en un thread independiente de modo que la interfaz no quede bloqueada durante este proceso. D.3.2. Resultado de la interfaz En la figura D.8 se observa el aspecto inicial de la aplicaci´on una vez cargados un lote de mapas de profundidad, y en la figura D.9 la situaci´on final despu´es de haber realizado la virtualizaci´on y posterior fusi´on de los mapas de profundidad. Como en las anteriores ocasiones la figura D.10 muestra numeradas cada una de las partes en las que se divide la interfaz para facilitar la explicaci´on. 84 D. Interfaces Secci´on D.3 Figura D.9: Aspecto final de la interfaz de fusi´on habiendo realizado ya la fusi´on de los mapas de profundidad virtualizados 1. Se trata de los botones para cargar los mapas de profundidad al comienzo de la ejecuci´on y para generar los mapas virtuales en el momento en el que es posible. 2. Esta zona es la de elecci´on del tipo de los mapas as´ı como de introducci´on de ruido a estos, que ser´an la entrada del siguiente paso. a) Aqu´ı se muestran todos lo mapas introducidos al comienzo, es posible navegar entre ellos para poder visualizarlos todos mediante los botones anterior y siguiente. Adem´as en esta parte existe la opci´on de usar mapas en ground truth o mapas estimados, esta elecci´on se aplicar´a solamente al mapa que se este mostrando en ese momento. b) Se trata de la zona de generaci´on de ruido en los mapas de profundidad introducidos al inicio, las opciones disponibles son la introducci´on, por un lado de ruido gausiano, y por otro de datos espurios, para lo que se indica el porcentaje de pixeles afectados y la variaci´on m´axima en estos. Como en el caso anterior estas modificaciones se ejecutan solamente en el mapa que se muestra en ese momento. 3. Es la pesta˜na de visualizaci´on de los mapas de profundidad ya colocados 85 Secci´on D.3 D. Interfaces Figura D.10: Esquema de las distintas partes que conforman la interfaz de la fusi´on de mapas. Figura D.11: Di´alogo para la introducci´on de la ruta de los mapas de profundidad todos respecto de la misma referencia, esta parte se muestra en la figura D.9. Como en el caso del visualizador de mapas inicial da la opci´on de navegar por todos ellos mediante los botones anterior y siguiente. 4. Se trata de la zona donde se mostrar´a el resultado final del experimento. 5. Por ´ultimo, en esta zona se encuentra, por un lado el bot´on que ordenar´a la ejecuci´on de la fusi´on siempre y cuando sea posible, y por otro la caja de texto donde se devolver´a el valor del an´alisis se˜nal ruido del resultado obtenido respecto de el ground truth correspondiente a la c´amara de referencia del resultado. Existen partes de la interfaz que no se muestran en la figura D.10, ya que se tratan de las distintas ventadas emergentes que se muestran para que el usuario introduzca los datos oportunos. 86 D. Interfaces Secci´on D.3 Figura D.12: Di´alogo para la introducci´on de la resoluci´on en las distintas dimensiones del cubo Figura D.13: Di´alogo de introducci´on de par´ametros para la fusi´on de mapas de profundidad. En la figura D.11, se muestra la primera de ellas, se acciona cuando el usuario pulsa el bot´on de cargar los mapas, en ella se debe introducir la ruta donde se almacenan los mapas de profundidad de entrada, as´ı como los par´ametros y posici´on de las c´amaras desde las que fueron construidos. El di´alogo de la figura D.12, se muestra antes de realizar la virtualizaci´on de los mapas de entrada, una vez se a pulsado el bot´on correspondiente, en ´el se introduce la resoluci´on del cubo en las distintas dimensiones. Por ´ultimo, el di´alogo mostrado en la figura A.5 entra en acci´on justo antes de realizar la fusi´on, en ´el el usuario tiene a posibilidad de introducir nuevos par´ametros para el m´etodo que generar´a el resultado final. D.3.3. Detalles de implementaci´on En cuanto a la implementaci´on la estructura es similar a las interfaces anteriores en cuanto al uso de threads, en este caso, como en la interfaz de prueba de par´ametros solo existen dos threads, el de la propia interfaz, y el que se encarga de realizar la fusi´on, y su estructura es la misma que en caso anterior, en la figura D.14 se observa esta distribuci´on. 87 Secci´on D.3 D. Interfaces Figura D.14: Pseudodiagrama de clases de la interfaz de Fusi´on MainWindow En cuanto a la clase MainWindow es similar a las anteriores ocasiones. Se basa en una serie de slots que se act´uan seg´un las se˜nales emitidas por los elementos de la interfaz. En cuanto a los par´ametros almacenados en la clase, en otras ocasiones se almacenaban las distintas im´agenes del proceso, la original, la ruidosa, etc. En este caso esas im´agenes se han convertido en vectores de im´agenes, adem´as no hay una entrada como en los problemas de denoising, hay dos, por un lado est´an los mapas de profundidad estimados, y por otro est´an los ground truth correspondientes, por lo que se duplican las variables, adem´as de eso, hay que que tener en cuenta que para representar los mapas han de estar normalizados, lo que significa otro vector m´as. Por ´ultimo tambi´en almacena el resultado de la virtualizaci´on normalizado y sin normalizar. Lanzador El lanzador es incluso m´as sencillo que en los casos anteriores ya que solo existe una funci´on a la que llamar, se trata la que hace la fusi´on. Esta funci´on tiene como entrada una serie de par´ametros que se le env´ıan al lanzador mediante la se˜nalfusion generada por la interfaz, adem´as de eso necesita, como es l´ogico los 88 D. Interfaces Secci´on D.3 mapas ya virtualizados, para ello la interfaz antes de emitir la se˜nal, ejecutar´a la funci´on cargaMapasVirtuales(),deestemodoyasepodr´arealizarlafusi´onyen este caso se realizar´a en un thread paralelo, con lo que la interfaz no se quedar´a inutilizada mientras tanto. Di´alogos En cuanto a las distintas clases de di´alogos su funcionamiento en bien sencillo. Se ejecutan desde la clase MainWindow, una vez aceptados la interfaz consulta los datos que ha introducido el usuario mediante las distintas funciones get que existen en las distintas clases di´alogo. Programaci´on en GPU Hasta este punto se ha detallado el funcionamiento de la interfaz, pero lo realmente importante es lo que queda detr´as de ella, se trata de las distintas clases que realizan la virtualizaci´on de los mapas y la fusi´on de los mismos. La programaci´on de la creaci´on de los mapas virtuales esta detallada en el anexo correspondiente a la fusi´on. En cuanto a la fusi´on propiamente dicha, se ha creado una clase llamada tgvFusion que simplemente tiene implementado el m´etodo de forma similar a la clase Normas utilizada para el denoising. 89 Secci´on E.1 E. Fusi´on Una vez hecho esto se realiza el muestreo, comenzando por la primera intersecci´on del cubo con el rayo, y usando como paso la m´axima distancia que asegura que se van a recorrer todos los v´oxeles intersectados, esta es. m´ın( dx resx ,dy resy ,dz resz )⇤0,5 Siendo dla distancia de la arista del cubo en cada una de las dimensiones. Para cada una de las muestras se tiene un valor en el espacio, el siguiente paso ser´a encontrar el voxel al que corresponde ese punto, se trata de una operaci´on sencilla ya que la coordenada a consultar se encuentra respecto de sistema de referencia del cubo, por lo que s´olo ser´a necesario dividir estas coordenadas reales entre la resoluci´on de cada eje del cubo. Una vez consultado el valor del voxel en cuesti´on se compara con el consultado en la muestra anterior, de modo que si se detecta que ha habido un cambio de signo, significar´a que se a atravesado la superficie. Tras esto a´un se realiza una interpolaci´on consultando el intervalo en el que se ha detectado la superficie para estimar mejor la distancia real. Una vez conseguida esta distancia se ha de trasladar de nuevo a la coordenada origen. Hay que recordar que en el mapa de profundidad se almacenan las distancias en el eje z, y no en l´ınea recta desde el origen que es el dato que se conoce en este momento, por lo que habr´a que calcular esta distancia. El valor obtenido de este c´alculo ser´a el introducido en el mapa de profundidad trasladado ya a la referencia deseada. Una vez realizados todos los pasos anteriores para cada uno de los mapas de entrada ya est´an dispuesto para realizar la fusi´on. E.1.2. Programaci´on Acontinuaci´onseenumeranlasclasesutilizadaspararealizarelprocesoanterior, en cuanto a su programaci´on interna es seguir los pasos descritos anteriormente por lo que no se va a hacer hincapi´e en ella. ReadDepthInfo : Se encarga de la lectura de mapas de profundidad de un directorio determinado as´ı como de la informaci´on de las c´amaras desde las que est´an construidos estos detph maps. Para lograr esto se requiere que los mapas y los par´ametros de las c´amaras tengan un formato concreto. 96 E. Fusi´on Secci´on E.1 RangeImageOp Se trata de una clase con operaciones muy ´utiles con mapas de profundidad como las siguientes: •Normalizar y desnormalizar mapas de profundidad, ya que los mapas para realizar operaciones con ellos como a˜nadir ruido o realizar la fusi´on necesitan no estar normalizados, sin embargo para representarlos si que hay que normalizarlos, de ah´ı la gran utilidad de este tipo de funciones. •A˜nadir ruido, esta es una utilidad tambi´en muy ´util de esta clase, con ella se puede introducir en la imagen de forma sencilla ruido gausiano, o datos espurios en un porcentaje de la imagen que se requiera. •Por ´ultimo permite el calculo del ´ındice de se˜nal ruido de un mapa de profundidad. En resumen es una clase muy ´util para este proceso, tanto como lo ser´a para la fusi´on. MatrixTransf Se trata de una clase que define una matriz de 3x4, as´ı como ciertas operaciones sobre ella ´utiles para c´alculos con matrices de transformaci´on. CudaCube Es la m´as importante de las anteriores es la que se encarga de realizar el proceso anterior propiamente dicho. Este se es totalmente paralelizable por lo que el coraz´on del mismo esta programado en CUDA Por ello esta clase en realidad est´a compuesta de dos una en C++ yotraen CUDA,enlapartedeC++ simplemente realiza las llamadas oportunas a la clase en CUDA del mismo nombre que ser´a la que realizar´a los c´alculos en paralelo. Las principales funciones de esta clase son las siguientes: •setDimFromDepthMap:Estaesla´unicadelastresfuncionesquesevan a describir que se ejecuta enteramente en la cpu, se encarga simplemente de calcular los par´ametros del cubo, este c´alculo lo realiza usando la c´amara del primer depth map que se ha le´ıdo. •fillCube: Como su propio nombre indica se encarga de rellenar el cubo que se ha generado mediante la funci´on anterior, con los datos de un depth map cualquiera. Esta funci´on esta paralelizada para que se ejecute una vez por cada pixel de la cara XY, y dentro de cada kernel se recorra toda la profundidad del cubo. 97 Secci´on E.1 E. Fusi´on •getVirtualDepthMap:Realizael´ultimopasodelproceso,esdecirleelos datos del cubo desde la posici´on deseada para conseguir el depth map virtual. Esta funci´on realiza el mismo proceso una vez para cada pixel del nuevo depth map, por lo que la paralelizaci´on es evidente. 98