Repositorio Institucional de Documentos

Proyecto Fin de Carrera de Ingenier´ıa Industrial Procesamiento de im´agenes a trav´es de m´etodos variacionales y de optimizaci´on convexa Francisco Javier Casanaba Bened´e Director: Lina Mar´ıa Paz Departamento de Inform´atica e Ingenier´ıa de Sistemas Centro Polit´ecnico Superior Septiembre de 2013 A mi familia, muy especialmente a mi hermano, quien seguro superar´a todas las piedras que la ingenier´ıa ponga en su camino. Hay una fuerza motriz m´as poderosa que el vapor, la electricidad y la energ´ıa at´omica: la voluntad. Albert Einstein Agradecimientos Me gustar´ıa que estas l´ıneas sirvieran para expresar mi enorme gratitud a todas aquellas personas que me han ayudado durante la realizaci´on de este proyecto, con menci´on especial a mi tutora Lina Mar´ıa Paz, quien con una generosidad fuera de lo com´un siempre tuvo un hueco para m´ı cuando lo necesit´e. Sergio, mi compa˜nero de fatigas, se merece un apartado exclusivo en este cap´ıtulo, pues su compa˜n´ıa y conocimientos de programaci´on durante estos meses han hecho mucho m´as agradable mi trabajo. Finalmente, agradecer a mi familia y amigos su comprensi´on y apoyo durante una etapa en la que no he podido atenderles todo lo se merecen. A todos ellos, de coraz´on, muchas gracias. ´ Indice 1. Introducci´on 1 1.1. Conceptosgenerales........................... 2 1.1.1. La transformaci´on de Legendre-Fenchel . . . . . . . . . . . . 3 1.1.2. Ladualidad........................... 3 1.1.3. Funciones no siempre diferenciables . . . . . . . . . . . . . . 4 1.1.4. Funciones convexas . . . . . . . . . . . . . . . . . . . . . . . 4 1.2. La funci´on de energ´ıa . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.3. ElGAP ................................. 6 1.4. Organizaci´on de la memoria . . . . . . . . . . . . . . . . . . . . . . 7 2. Denoising 11 2.1. Motivaci´on................................ 11 2.2. ElmodeloTV-ROF........................... 11 2.2.1. C´alculo del GAP . . . . . . . . . . . . . . . . . . . . . . . . 13 2.2.2. Los par´ametros del algoritmo . . . . . . . . . . . . . . . . . 14 2.3. El modelo Huber-ROF . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.3.1. El par´ametro α......................... 17 2.4. El modelo TV L1............................ 18 3. Zooming 21 3.1. Motivaci´on................................ 21 3.2. El modelo de energ´ıa . . . . . . . . . . . . . . . . . . . . . . . . . . 21 4. Image Deconvolution 25 4.1. Motivaci´on................................ 25 4.2. Modelodeenerg´ıa............................ 26 5. Image Inpainting 29 5.1. Motivaci´on................................ 29 5.2. Elalgoritmo............................... 29 i ´ INDICE ´ INDICE 6. Conclusiones 31 A. Gesti´on del proyecto 33 B. Deducciones matem´aticas 35 B.1.TV-ROF................................. 35 B.1.1. C´alculo de ∂pE(u, p) ...................... 35 B.1.2. C´alculo de ∂uE(u, p) ...................... 35 B.2.HUBER-ROF.............................. 36 B.2.1. C´alculo de ∂pE(u, p) ...................... 36 B.3.TV-L1.................................. 36 B.3.1. C´alculo de ∂uE(u, p) ...................... 36 B.4. Derivaci´on de la energ´ıa para el problema del zooming . . . . . . . 36 B.5. Derivaci´on de la funci´on de energ´ıa para el problema de la deconvoluci´on ................................. 38 C. Tipos de ruido 41 C.1.ElruidoGaussiano ........................... 41 C.2. El ruido de sal y pimienta . . . . . . . . . . . . . . . . . . . . . . . 42 D. Software 45 D.1.ElCMakelist.txt ............................ 45 D.2. Archivos de interfaz . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 D.3.CUDA.................................. 48 E. Programaci´on 51 F. Resultados 55 F.1. Resultados Denoising . . . . . . . . . . . . . . . . . . . . . . . . . . 55 F.1.1. Evaluaci´on del modelo TV-ROF . . . . . . . . . . . . . . . . 55 F.1.2. Evaluaci´on del modelo Huber-ROF . . . . . . . . . . . . . . 58 F.1.3. Evaluaci´on del modelo TV-L1 . . . . . . . . . . . . . . . . . 59 F.1.4. Comparaci´on de modelos . . . . . . . . . . . . . . . . . . . . 61 F.2. Resultados Zooming . . . . . . . . . . . . . . . . . . . . . . . . . . 64 F.3. Resultados Deconvolution . . . . . . . . . . . . . . . . . . . . . . . 68 F.4. Resultados Image Inpainting . . . . . . . . . . . . . . . . . . . . . . 73 G. Manual de usuario 81 G.1.VentanaPrincipal............................ 81 G.2.VentanaDenoising ........................... 82 G.3.VentanaZooming............................ 85 G.4. Ventana Deconvolution . . . . . . . . . . . . . . . . . . . . . . . . . 89 ii ´ INDICE ´ INDICE G.5. Ventana Inpainting . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 iii Secci´on 1.1 1. Introducci´on to se har´a referencia a los diferentes problemas de acuerdo a los nombres recibidos en ingl´es: denoising, zooming, deconvolution (motion deblurring) e inpainting. En este proyecto se aborda el modelado de cada problema desde el punto de vista de los modelos continuos variacionales. El porqu´e de esta selecci´on se debe al gran potencial que tienen dichos m´etodos para tratar la informaci´on densa de forma robusta. As´ı mismo, las soluciones proporcionadas para dichos m´etodos est´an basadas en los avances desarrollados dentro del campo de la optimizaci´on convexa. Desde este punto de vista, los algoritmos derivados permiten calcular soluciones globales de forma eficiente debido al tratamiento independiente de los p´ıxeles. Como resultado, los algoritmos usados en este proyecto permiten la explotaci´on del c´alculo masivo paralelo en GPGPU (General Purpose Graphich Processing Unit) siendo de gran importancia en las aplicaciones de tiempo real. En este cap´ıtulo se introducen las bases te´oricas de los algoritmos empleados en el proyecto. En la secci´on 1.1, se expondr´an los conceptos matem´aticos sobre los que se sustentan los algoritmos empleados. As´ı, se comienza definiendo la transformaci´on de Legendre-Fenchel, c´omo y por qu´e se aplica en la resoluci´on del algoritmo. Se contin´ua definiendo la dualidad, qu´e es el problema dual y c´omo y por qu´e resolverlo. Para finalizar, se tratan los problemas de funciones no siempre diferenciables y c´omo poder transformarlos en una funci´on convexa resoluble. En la secci´on 1.2, se incluye la definici´on de la funci´on de energ´ıa, eje principal de los modelos que se han utilizado en este proyecto. Al tratarse de un m´etodo num´erico, se debe definir un criterio de parada, as´ı surge el concepto del GAP descrito en la secci´on 1.3. Finalmente, la secci´on 1.4 muestra la organizaci´on de la memoria del proyecto, suponiendo el fin de este cap´ıtulo de introducci´on. 1.1. Conceptos generales Los modelos usados para cada problema concreto se basan en normas que son convexas pero no diferenciables, de forma que no se pueden utilizar algoritmos de optimizaci´on convencionales. Adem´as, el problema general consiste en la minimizaci´on de una funci´on de energ´ıa, funci´on que no es necesariamente convexa y por tanto no se puede alcanzar una soluci´on global. Con el fin de solucionar esta situaci´on, se aplicar´a la transformaci´on de Legendre Fenchel (Apartado 1.1.1). Esta transformaci´on requiere dualizar una funci´on, motivo por el que se explica el concepto de dualizaci´on y c´omo se debe aplicar. M´as tarde se expone un ejemplo de funciones no siempre diferenciables donde se dualiza con el fin de generar un problema diferenciable en todo su rango. Para finalizar, se comenta la necesidad de partir de una funci´on convexa para el algoritmo de soluci´on, pues de lo contrario no se podr´ıa obtener una soluci´on general. [2]. 2 1. Introducci´on Secci´on 1.1 1.1.1. La transformaci´on de Legendre-Fenchel La transformada de Legrendre-Fenchel (LF en adelante) [7] de una funci´on continua pero no necesariamente diferenciable se define como f∗(p) = sup x∈R{px −f(x)}(1.1) La figura 1.1 ilustra la transformaci´on LF de una funci´on cuadr´atica. Geom´etricamente se trata de la b´usqueda de un punto (x, f(x)) para el cual una recta de pendiente pproduzca el m´aximo corte con el eje vertical. El conjunto de rectas tangentes forma una envolvente que representa a s´ı misma la funci´on. p=f0(x) (1.2) Figura 1.1: Recta tangente a una funci´on cuadr´atica La representaci´on general vectorial de la transformaci´on LF para funciones multivariables se define como: f∗(p) = sup x∈R{xtp−f(x)}(1.3) 1.1.2. La dualidad La dualidad es el principio por el que observamos una misma funci´on desde dos diferentes perspectivas: primal y dual. Si se supone que la transformaci´on LF es reversible, se puede afirmar que (x, f(x)) ⇐⇒ (p, f∗(p)) (1.4) 3 Secci´on 1.1 1. Introducci´on donde pes la pendiente y f*(p) es el llamado conjugado convexo de la funci´on f(x). Un conjugado nos permite construir un problema dual que puede ser m´as f´acil de resolver que el problema primal. Observando la figura 1.1, se puede afirmar que cada punto de la funci´on primal se puede representar por la pendiente de la tangente a ese punto, obteniendo el dual. De esta manera, los puntos del dual son las pendientes del primal, y los puntos del primal son las pendientes del dual. El conjugado de LF es siempre convexo. 1.1.3. Funciones no siempre diferenciables Observemos la ecuaci´on [2] definida al principio, donde pes la pendiente tangente a f(x) en el punto xdefinida como la derivada f0(x) ¿Qu´e sucede si la pno es diferenciable en todo rango de x? La ecuaci´on com´un de una recta es la siguiente: y=px −c(1.5) Si se modifica esta ecuaci´on para que se refiera a un punto x∗no diferenciable, se puede reescribir como f(x∗) = px∗−c(1.6) El problema dual de este punto no diferenciable se convierte en una funci´on lineal en p. A pesar de la simplicidad de este ejemplo lineal, este mismo principio puede emplearse en el caso de funciones no diferenciables ya que el resultado es una nueva funci´on dual en pque puede ser diferenciable. Una de las grandes ventajas del problema dual es que pese a que el primal puede ser no diferenciable, el problema dual s´ı que lo es. 1.1.4. Funciones convexas Los problemas de procesamiento de im´agenes tratados en este proyecto se modelan como un problema de minimizaci´on de una funci´on de energ´ıa. De forma general, el punto ˆxes el m´ınimo global de una funci´on siempre y cuando ∇f(ˆx) = 0 y la funci´on sea convexa. Sin embargo, existen muchas funciones que no son siempre diferenciables o convexas, por lo que no se puede realizar el c´alculo de su gradiente. La figura 1.2 ilustra un conjunto de funciones convexas. 4 1. Introducci´on Secci´on 1.2 Figura 1.2: Ejemplos de funciones convexas La funci´on cuadr´atica es continua y diferenciable en todo su dominio. No as´ı la funci´on valor absoluto, que en el punto x= 0 no es diferenciable. En este ´ultimo caso, se puede demostrar que la transformada de LF es convexa. Por tanto, basta con realizar la transformaci´on de LF a una funci´on para obtener una funci´on convexa. 1.2. La funci´on de energ´ıa Los m´etodos de procesamiento de im´agenes (Image Processing en ingl´es) se basan en la minimizaci´on de una funci´on de energ´ıa. Esta funci´on est´a formada por dos t´erminos: El primero de ellos es el llamado Regularizador. Este t´ermino es el que cuantifica la variaci´on dentro de una propia imagen, tratando que las diferentes superficies tengan una textura uniforme y las esquinas y bordes queden lo m´as definidas posibles. El segundo t´ermino es el denominado ”Data Term”, que es el que cuantifica cu´an distinta es la imagen soluci´on (imagen arreglada) de la imagen original (imagen defectuosa). Al aplicar estos m´etodos se parte de la base de que la imagen defectuosa contiene a´un suficiente informaci´on para ser restaurada, aunque no en su totalidad, pero s´ı con un alto porcentaje de mejora sobre su estado inicial. En general, la ecuaci´on de la funci´on de energ´ıa a minimizar es la siguiente: m´ın x∈XE(x) (1.7) m´ın x∈XF(Kx) + G(x) (1.8) 5 Secci´on 1.3 1. Introducci´on Donde K:X→Yes un mapa lineal entre dos espacios vectoriales X e Y equipados con un operador de producto escalar <·,·>y una norma |·| =|·,·|. G :X→ [0,+´ınf) y F∗:X→[0,+´ınf) son funciones propias convexas, semi-continuas inferiormente [6]. El primer t´ermino de la ecuaci´on es el regularizador, el cual usa como ´unica informaci´on la variable primal. El segundo cuantifica la variaci´on de la variable primal respecto a su estado inicial. En cada caso particular de restauraci´on de im´agenes, la funci´on de energ´ıa ser´a redefinida para una variable primal que corresponder´a a la imagen soluci´on definida en el dominio Ω ∈R. A continuaci´on se dualiza esta funci´on del que resulta un problema de minimizaci´on y maximizaci´on sobre las variables primal y dual respectivamente. m´ın x∈Xm´ax y∈Y{hKx, yi−F∗(y) + G(x)}(1.9) Este resultado es fundamental para derivar los algoritmos de optimizaci´on empleados en este proyecto. 1.3. El GAP Una vez definida la funci´on, es necesario generar un par´ametro que nos muestre de una forma clara y concisa cu´anto ha sido maximizada y minimizada la funci´on de energ´ıa. Ese valor toma el nombre de GAP. Para generar una funci´on que defina el GAP se hace uso de la definici´on de transformada LF, ya mencionada con anterioridad. As´ı, si se observa la ecuaci´on anterior, se puede apreciar que F(Kx) = m´ax y∈Y{hKx, yi−F∗(y)}(1.10) En esta ecuaci´on se puede modificar el t´ermino hKx, yiy transformarlo a hx, K∗yi, obteniendo por tanto la ecuaci´on m´ın x∈Xm´ax y∈Y{hx, K∗yi−F∗(y) + G(x)}(1.11) Si se realiza el siguiente procedimiento de la misma forma que el anterior: m´ax x∈X{hx, K∗yi−G(x)}=G∗(K∗y) (1.12) m´ın x∈X{−hx, K∗yi+G(x)}=−G∗(K∗y) m´ın x∈X{−hx, −K∗yi+G(x)}=−G∗(−K∗y) m´ın x∈X{hx, K∗yi+G(x)}=−G∗(−K∗y) 6 1. Introducci´on Secci´on 1.4 Por lo que aparecen dos optimizaciones derivadas de la funci´on de energ´ıa: m´ın x∈X{F(Kx) + G(x)}(1.13) m´ax y∈Y{−G∗(−K∗y)−F∗(y)} El GAP se define como la diferencia entre ambas optimizaciones: GAP = m´ın x∈X{F(Kx) + G(x)}−m´ax y∈Y{−G∗(−K∗y)−F∗(y)}(1.14) Es el t´ermino m´as preciso de optimizaci´on, pues la funci´on de energ´ıa dualizada debe ser minimizada respecto a una variable y maximizada respecto a otra (ver figura 1.3). El GAP es decreciente (salvo en los instantes finales de convergencia, donde aparece algo de ruido) y aporta una idea muy sencilla de cu´anto se est´a optimizando la funci´on. No obstante, el c´alculo de GAP en ciertas aplicaciones puede ser muy costoso o incluso imposible, ya que requiere una dualizaci´on que no siempre se podr´a llevar a cabo. Cuando ´esto sucede, simplemente se comprueba que la funci´on que se minimiza efectivamente es decreciente y alcanza un m´ınimo. Figura 1.3: Evoluci´on de los procesos de minimizaci´on y maximizaci´on para las variables Primal y Dual respectivamente. As´ı, el objetivo es alcanzar el punto silla donde la funci´on Primal es m´ınima y la funci´on Dual es m´axima. 1.4. Organizaci´on de la memoria La memoria de este proyecto se estructurar´a de forma que, tras un primer cap´ıtulo de introducci´on, se expondr´an los diferentes problemas a tratar con sus diferentes secciones. De esta forma: Denoising: Se comienza explicando la motivaci´on por la que se desea resolver el problema, dando ejemplos pr´acticos donde podr´ıa ser necesaria. Denoising puede 7 Secci´on 1.4 1. Introducci´on ser resuelto de diferentes formas posibles, dependiendo de la funci´on de energ´ıa utilizada. As´ı, los m´etodos toman un nombre asociado a su resoluci´on. El regularizador se denomina Total Variation cuando se resuelve seg´un la norma 2. En caso de que se resuelva mediante la norma de Huber, se denominar´a ”Huber”. El t´ermino G(x) se denomina L1 si se resuelve seg´un la norma 1 o ROF si se resuelve seg´un la norma cuadr´atica. As´ı, los tres posibles m´etodos de resoluci´on que se manejar´an son TV-ROF, HUBER-ROF y TV-L1. En los casos en los que sea posible, se calcular´a el GAP. Una vez explicados los tres m´etodos, se finaliza haciendo referencia a los resultados situados en el Ap´endice F. Zooming: Se comienza con una motivaci´on explicativa ofreciendo los motivos por los que se plantea la resoluci´on del problema. Despu´es se da una breve introducci´on al m´etodo num´erico de Jacobi, necesario para la resoluci´on del algoritmo. Una vez explicado, se deduce, explica y resuelve el algoritmo, incidiendo en las diferencias respecto al caso anterior. Image Deconvolution: Nuevamente, se introduce el problema mediante un apartado de motivaci´on explicando por qu´e se produce este fen´omeno y por qu´e es importante solucionarlo. El siguiente apartado explica c´omo obtener la m´ascara que genera el blurring en la imagen. Es muy importante el c´alculo de esta m´ascara porque ser´a necesaria en la resoluci´on del algoritmo. A continuaci´on se resuelve dicho algoritmo. Inpainting: Es el ´ultimo de los problemas a tratar. La estructura seguida es similar a la de los casos anteriores. Se comienza con una explicaci´on en la que se comenta porqu´e puede surgir y por qu´e solucionarlo. A continuaci´on, se deduce y resuelve el algoritmo. Anexos: Donde se encuentran tanto explicaciones acerca de la programaci´on, de matem´aticas que no se introducen en la memoria por su complejidad pero que son necesarias para la resoluci´on del algoritmo, resultados de los diferentes algoritmos y un manual de usuario final, que es la conclusi´on definitiva del trabajo. En este proyecto se ha trabajado en la soluci´on de distintos problemas, obteniendo importantes datos en su resoluci´on tales como el gap, el tiempo de ejecuci´on, las iteraciones que tarda en converger o c´omo de fiable es el resultado final respecto a imagen original, la que no sufre modificaci´on. Estos resultados dependen del algoritmo obviamente, pero tambi´en de diferentes par´ametros que pueden estar acotados por la convergencia del m´etodo, pero que poseen cierta holgura y pueden variar lo ´optimo del m´etodo. As´ı, se ha generado una aplicaci´on donde se permite al usuario escoger entre los cuatro problemas citados y proceder a su resoluci´on. Cada uno de ellos tiene dos posibles modos de resoluci´on. El primero de ellos, llamado Programa Principal permite al usuario introducir la imagen, da˜narla en los casos que sea necesario y arreglarla con los par´ametros que se requieran. Como resultado se obtiene la imagen resultado y una gr´afica de evoluci´on de GAP o 8 1. Introducci´on Secci´on 1.4 funci´on de energ´ıa respecto a iteraciones seg´un el caso. La segunda pesta˜na es la de Evaluaci´on de Par´ametros que es capaz de calcular el par´ametro ´optimo para una imagen concreta. 9 Cap´ıtulo 2 Denoising 2.1. Motivaci´on Denoising se refiere al proceso de filtrado del ruido de una imagen. En la industria es muy probable encontrar im´agenes que deban someterse a un tratamiento de este tipo. Por ejemplo, en una cadena de manufactura de muebles hay gran cantidad de virutas esparcidas por el aire, que pueden perturbar la imagen. Tambi´en en el env´ıo de datos se puede a˜nadir ruido que la perturbe. Mediante un algoritmo primal-dual somos capaces de recuperar razonablemente la imagen original, pudiendo someterla as´ı a cualquier algoritmo de reconocimiento, de medici´on o de control de calidad. A continuaci´on se va a tratar el problema utilizando tres modelos diferentes que permiten analizar el impacto al cambiar la norma tanto en el regularizador como en el data term: modelo ROF, el huber ROF y el TVL1. 2.2. El modelo TV-ROF El modelo ROF se define como el problema variacional m´ın uZΩ|Du|+λ 2ku−gk2 2(2.1) Donde ues la imagen resultado, gla imagen original (la que tiene ruido) y λes un par´ametro utilizado para definir la acomodaci´on entre los dos t´erminos de la ecuaci´on. Esta ecuaci´on es continua, como una imagen se divide en p´ıxeles y ´estos son discretos, se reescribe de la siguiente forma m´ın u∈Uk ∇uk1+λ 2ku−gk2 2(2.2) 11 Secci´on 2.4 2. Denoising (a) α=0.1 (b) α=0.5 (b) α=1 Figura 2.2: Variaci´on del par´ametro αde la norma de Huber αes un par´ametro que optimiza el resultado, no la convergencia. Por este motivo, al igual que λ, se debe evaluar el SnR para determinar su valor ´optimo. 2.4. El modelo TV L1 Al igual que en el modelo ROF, se usa la variaci´on total para el t´ermino de regularizaci´on pero esta vez se introduce la norma L1para el t´ermino de datos. m´ın u∈Uk ∇uk1+kλ(u−g)k1(2.29) Se dualiza de nuevo el regulador, obteniendo las mismas expresiones que en el caso del modelo ROF. El problema que surge ahora es que la norma L1del t´ermino de comparaci´on con la imagen original implica una nueva dualizaci´on. kλ(u−g)k1= m´ax q∈Q(hq, λ(u−g)i−δQ(q)) (2.30) Se tendr´a que definir un nuevo convexo conjugado, realizar una ∂qde la funci´on y adem´as a˜nadir nuevos t´erminos. La derivaci´on es un poco m´as complicada y requiere del uso del concepto de subgradiente para derivar el t´ermino de datos respecto a la variable dual. La norma L1puede ser definida como una funci´on a trozos, donde: ∂u|f(u)|   f0(u)si u > 0 −f0(u)si u < 0 [-f’(u), f’(u)] si u = 0 (2.31) 18 2. Denoising Secci´on 2.4 (a) Funci´on valor absoluto (b) Funci´on derivada del valor absoluto Figura 2.3: Ilustraci´on del concepto de subgradiente para la funci´on valor absoluto Para el punto x= 0, puede haber un n´umero infinito de rectas tangentes a la funci´on en el intervalo p∈[−f0(x), f0(x)]. Por este motivo aunque la funci´on es convexa, no es diferenciable. En nuestro caso concreto, se puede escribir lo siguiente ∂uτλ|u−g|1   τλ si u −g > τλ −τλ si u −g < −τλ indef si |u−g| ≤ τλ (2.32) Este resultado se utiliza para resolver el algoritmo a partir del c´alculo de las derivadas sobre la siguiente funci´on de energ´ıa dualizada: m´ın u∈Um´ax p∈P(hp, ∇ui−δp(p) + λku−gk1) (2.33) Donde: ∂pE(u, p) = ∇u−αpn+1 =pn+1 −pn σ(2.34) ∂uE(u, p) = −(−divp +∂uλku−gk1=un+1 −un τ) (2.35) Tras este desarrollo matem´atico, se obtiene el valor de un+1. Gracias a la ecuaci´on 2.35, la actualizaci´on del primal se ejecuta de la siguiente manera: un+1    un+τdivpn+1 −τλ si u −g > τλ un+τdivpn+1 +τλ si u −g < −τλ g si |u−g| ≤ τλ (2.36) 19 Secci´on 2.4 2. Denoising El algoritmo empleado es el mismo que en el caso del ROF. Los resultados y conclusiones de los diferentes m´etodos del algoritmo se encuentran en el Ap´endice F.1 20 Cap´ıtulo 3 Zooming 3.1. Motivaci´on Cu´antas veces al reducir una fotograf´ıa y volver a ampliarla ha aparecido pixelada. Habitualmente, al expandir una fotograf´ıa cada p´ıxel se repite por un factor de ampliaci´on, lo que al final resulta una imagen muy poco realista, donde se aprecia con demasiada claridad que ha sido ampliada y no es una imagen original. La soluci´on m´as habitual consiste en la aplicaci´on de una simple interpolaci´on lineal. Si se ampl´ıa cuatro veces una imagen (el doble de alto y el doble de ancho), los p´ıxeles generados no tienen el valor constante de su predecesor, sino una interpolaci´on de sus vecinos. As´ı se consigue una sensaci´on m´as homog´enea en la imagen. Sin embargo, una interpolaci´on lineal podr´ıa producir un efecto de difuminaci´on de la imagen, eliminando los detalles. La aplicaci´on de un algoritmo primal dual para la minimizaci´on de una funci´on de energ´ıa puede ser una soluci´on al problema. Este apartado trata de su desarrollo matem´atico, implementaci´on y conclusiones. 3.2. El modelo de energ´ıa La ecuaci´on 3.1 define el problema zooming como la minimizaci´on de una energ´ıa. m´ın u∈Uk ∇uk1+λ 2k(Au −g)k2 2(3.1) Una diferencia b´asica respecto al problema de denoising es el uso del operador lineal representado por la matriz A. En el caso general del denoising esta matriz puede considerarse como la matriz identidad. Sin embargo, en el problema del zooming es de vital importancia, pues la matriz Ainfluye directamente en t´ermino de datos de la energ´ıa para transformar las dimensiones de la imagen ampliada u 21 Secci´on 3.2 3. Zooming en las dimensiones de la imagen de partida g. Consideremos el siguiente ejemplo, en el que la imagen de entrada es: g=1 3 2 4  Y, con un factor de ampliaci´on de 2 (2 de ancho y 2 de alto), la salida es de la forma: u=    1 5 9 13 2 6 10 14 3 7 11 15 4 8 12 16     Donde los valores interiores de la matriz hacen referencia a la posici´on de sus t´erminos en un vector de la matriz organizado por columnas. La matriz Aque transforma dimensiones de ues por tanto de la siguiente forma: A=    X X 0 0 X X 0000000000 0 0 X X 0 0 X X 00000000 00000000X X 0 0 X X 0 0 0000000000X X 0 0 X X     La primera fila hace referencia al primer t´ermino del vector g, donde las posiciones 1, 2, 5 y 6 apuntan a la matriz u. F´ısicamente se puede interpretar como que el p´ıxel 1 de la matriz g se expande hacia los p´ıxeles 1, 2, 5 y 6 de la matriz u. Lo siguiente a analizar es cu´anto vale la X. Hay tantas X en una fila como factor de ampliaci´on elevado al cuadrado. Si cada X valiese 1, el producto de Au ser´ıa de aproximadamente, factor al cuadrado veces el valor de g. El t´ermino de datos valora la variaci´on entre la imagen de salida y la de entrada, por lo que para que sean del mismo orden se requiere que X=1 s2, donde ses el factor de ampliaci´on. m´ın u∈Um´ax p∈Php, ∇ui+λ 2kAu −gk2 2−δp(p) (3.2) El problema dual ser´a similar al de denoising ROF, pues el ´unico cambio respecto a ese problema est´a en un t´ermino que no depende de p. El procedimiento para derivar el paso de actualizaci´on de la variable primal se describe en el Apendice B.4 un+1 = (I+τλATA)−1(un+τdivpn+1 +τλATg) (3.3) La expresi´on 3.3 requiere de la soluci´on de un sistema lineal de ecuaciones de la forma Mx =bque en muchos casos puede llevar al desarrollo de un algoritmo muy 22 3. Zooming Secci´on 3.2 lento con dificultad en su paralelizaci´on mediante GPU. Para evitar este problema se utiliza el m´etodo de soluci´on de ecuaciones de Jacobi. [1]. Si se sustituye µ=1 τ, se obtiene el siguiente resultado: un+1(τI +λATA) = µun+divpn+1 +λATg(3.4) Para solucionar el sistema de ecuaciones se toma M=µI +λATAyb=µun+ divpn+1 +λATg. La matriz Mde Jacobi se debe descomponer en DyR: D= (µ+λ s4)I, ya que los elementos de A est´an divididos por s2y los valores de la diagonal de AAestar´an divididos por s4. N´otese que Des un valor constante. R=λATA−λ s4I b=µun+divpn+1 +λATg Por tanto, el resultado de la actualizaci´on primal es el siguiente: un+1 ={λATA−λ s4I}un+µun+divpn+1 +λATg µ+λ s4 (3.5) Es importante decir que el m´etodo de Jacobi es iterativo, por lo que esta actualizaci´on deber´ıa realizarse varias veces en cada actualizaci´on del primal. Es decir, que el orden de las iteraciones se elevar´ıa al cuadrado. Dada la r´apida convergencia del algoritmo primal dual, tan s´olo se realizar´a una iteraci´on de Jacobi. Se puede observar que la matriz ATApromedia los p´ıxeles de la imagen grande que hacen referencia a la imagen peque˜na mediante el factor 1 s4. Ese valor medio es el mismo para todos los p´ıxeles de la imagen grande que tienen en com´un un mismo p´ıxel de la peque˜na. Este resultado permite la implementaci´on sencilla del algoritmo. Los resultados se encuentran en el Ap´endice F.2. 23 Cap´ıtulo 4 Image Deconvolution 4.1. Motivaci´on La calidad de las im´agenes capturadas mediante c´amaras digitales dependen del tiempo de exposici´on del sensor. Si aparecen objetos m´oviles en la escena o la c´amara se encuentra en movimiento, un tiempo de exposici´on excesivo puede afectar las im´agenes ya que en cada instante de tiempo el sensor recibe y promedia diferente informaci´on de intensidad. Este fen´omeno es tambi´en conocido como motion blurring o convoluci´on de movimiento. A partir de una imagen que ha sufrido un proceso de convoluci´on se puede recuperar la imagen original en unas condiciones muy aceptables. S´olo se necesitan otros dos par´ametros adem´as de la imagen: la longitud en p´ıxeles que se ha movido la imagen y la direcci´on en que ´esto ha sucedido. La ecuaci´on que relaciona la velocidad del objeto respecto a la c´amara y la velocidad expresada en p´ıxeles, se define como e=vT d(4.1) Donde v es la velocidad en el plano de la imagen, f´acil de obtener a partir de la velocidad real del objeto y la distancia al mismo, T es el tiempo de exposici´on y d el tama˜no de un p´ıxel. La velocidad en el plano de la imagen se puede obtener gracias a una sencilla relaci´on trigonom´etrica representada en la figura 4.1. 25 Secci´on 4.2 4. Image Deconvolution Figura 4.1: Transformaci´on entre el vector de velocidad del objeto y el vector de velocidad en el plano de imagen Suponiendo la velocidad como un vector, hay que calcular el vector proporcional en el plano imagen. Denotando vcomo la velocidad en el plano imagen, Vcomo la velocidad real, dcomo la distancia focal y Dcomo la distancia desde la lente hasta el objeto en movimiento, se deduce que: v=dV D+d(4.2) Colocando una c´amara de bajas prestaciones en una l´ınea de manufactura, sabiendo a la velocidad y direcci´on en que se mueve, y la distancia de la c´amara a la l´ınea, se puede predecir la convoluci´on que sufrir´a la imagen obtenida y se podr´a arreglar con gran exactitud. 4.2. Modelo de energ´ıa El caso de la deconvoluci´on es an´alogo al del zooming, con la salvedad de que la matriz Arepresenta el movimiento en p´ıxeles aplicado a la imagen. Esta matriz se puede modelar mediante un kernel local o matriz de convoluci´on que representa el movimiento lineal de los p´ıxeles. Dicha matriz es f´acil de calcular requiriendo s´olo como par´ametros la longitud en p´ıxeles y su ´angulo de inclinaci´on de la recta. La figura 4.2 muestra un ejemplo de una m´ascara calculada para una longitud de 7 p´ıxeles y un ´angulo de 45ode inclinaci´on. 26 4. Image Deconvolution Secci´on 4.2 A=           000000,0144942 0 00000,0375741 0,128286 0,0144944 0 0 0 0,0375741 0,128286 0,0375742 0 0 0 0,0375741 0,128286 0,0375741 0 0 0 0,0375742 0,128286 0,0375741 0 0 0 0,0144944 0,128286 0,0375741 0 0 0 0 0 0,0144942 0 0 0 0 0           Figura 4.2: C´alculo de una m´ascara de motion blurring. Los valores obtenidos est´an normalizados, de manera que la suma global de todos los elementos es 1. Figura 4.3: Representaci´on del proceso de convoluci´on sobre la imagen de entrada (derecha), dado un operador de motion blurring desplazado circularmente (izquierda). Para modelar el problema de la deconvoluci´on, se parte de la ecuaci´on 4.3 que describe el modelo de minimizaci´on de energ´ıa: m´ın u∈Uk ∇uk1+λ 2k(Au −g)k2 2(4.3) El paso m´as importante es la derivaci´on de la actualizaci´on de la variable pri27 A. Gesti´on del proyecto implementar el algoritmo en la tarjeta gr´afica. Como cada thread de la tarjeta es capaz de realizar una operaci´on sencilla, se puede sustituir ese bucle ”for”por una paralelizaci´on completa, y as´ı optimizar el tiempo. La soluci´on es programar en CUDA, un conjunto de herramientas desarrollado por nVidia para codificar algoritmos en GPU nVidia. Como no se dispon´ıa de un ordenador personal con GPU nVidia, se solicit´o una cuenta en el servidor ”Hermes”de unizar donde se ejecutar´ıan los algoritmos. As´ı, se procedi´o a la realizaci´on de los diferentes algoritmos en Hermes. Se implementaron TV-ROF, Huber-ROF, TV-L1 y Zooming. Una vez que se quer´ıa realizar el algoritmo de Image deconvolution, era necesaria una librer´ıa que calculase la FFT (Fast Fourier Transform). Al no encontrarse instalada en el servidor y, sabiendo que ´este iba a estar apagado durante el mes de Agosto, se gener´o una cuenta en el ordenador de la universidad de mi tutora, Lina Mar´ıa Paz, con el fin de que pudiese seguir trabajando ah´ı. Cuando se implantaron los diferentes algoritmos, surgi´o la duda de cu´ales eran los par´ametros ´optimos que se deb´ıan colocar en cada uno de dichos algoritmos. Una vez entendido el significado de cada par´ametro y qu´e optimiza, la manera de actuar era relativamente sencilla: recorrer un rango de valores con un paso relativamente peque˜no y mostrar por pantalla el resultado de cada caso (ratio de se˜nal/ruido o iteraciones hasta convergencia). Este hecho provocaba la necesidad de realizar un nuevo programa, de forma que para ordenar todos los algoritmos, se decidi´o realizar una interfaz de usuario donde se pod´ıa seleccionar el tratamiento al que se quer´ıa someter a la imagen y poder observar la imagen modificada, la arreglada, c´omo evoluciona el GAP o la funci´on de energ´ıa y los la evoluci´on de los diferentes par´ametros. 34 Ap´endice B Deducciones matem´aticas B.1. TV-ROF B.1.1. C´alculo de ∂pE(u, p) ∂pE(u, p) = ∂p(hp, ∇ui+λ 2ku−gk2 2−δp(p)) (B.1) Se simplifica t´ermino a t´ermino: ∂p(hp, 5ui) = 5u, ya que se trata de un producto escalar. ∂p(λ 2ku−gk2 2) = 0, ya que no depende de p. ∂p(δp(p)) = 0, Se coloca como 0 para permitir su c´alculo. Sin embargo, este t´ermino impone la restricci´on de que la variable dual en el resultado final no sea mayor que 1. Por tanto se concluye que: ∂pE(u, p) = ∇u(B.2) B.1.2. C´alculo de ∂uE(u, p) ∂uE(u, p) = ∂u(hp, ∇ui+λ 2ku−gk2 2−δp(p)) (B.3) Se simplifica t´ermino a t´ermino ∂u(hp, 5ui) = ∂u(−hu, divpi) = −divp. En el primer paso se convoluciona el producto de vector por gradiente en un producto negado de vector por divergencia. Se trata de una propiedad matem´atica aplicable al operador ∇. Una vez que es un producto derivable, se obtiene el resultado. ∂u(λ 2ku−gk2 2) = λ(u−g), al ser una norma cuadr´atica simple. 35 Secci´on B.2 B. Deducciones matem´aticas ∂u(δp(p)) = 0, ya que no depende de u. Por tanto se concluye que ∂uE(u, p) = −divp +λ(u−g) (B.4) B.2. HUBER-ROF B.2.1. C´alculo de ∂pE(u, p) ∂pE(u, p) = ∂p(hp, ∇ui−δp(p)−α 2kpk2+λ 2ku−gk2 2) (B.5) ∂p(hp, ∇ui) = ∇u, ya que se trata de un producto escalar. ∂p(δp(p)) = 0, Es el mismo caso que el anterior. Se coloca como 0 pero su efecto se a˜nadir´a al final del desarrollo. ∂p(α 2kpk2) = αp Una derivada com´un de un valor al cuadrado. ∂p(λ 2ku−gk2 2) = 0, ya que no depende de p. Por tanto, obtenemos como resultado ∂pE(u, p) = 5u−αp (B.6) B.3. TV-L1 B.3.1. C´alculo de ∂uE(u, p) ∂uE(u, p) = ∂u(hp, ∇ui−δp(p) + λku−gk1) (B.7) ∂u(hp, ∇ui) = ∂u(−hu, divpi) = −divp, ∂u(δp(p)) = 0 ∂u(λ 2ku−gk2 2) = λ(u−g) ∂uE(u, p) = −divp +∂uλku−gk1(B.8) B.4. Derivaci´on de la energ´ıa para el problema del zooming Partiendo del modelo de energ´ıa, se desea calcular su derivada respecto a la variable primal: 36 B. Deducciones matem´aticas Secci´on B.4 ∂uE(u, p) = ∂u(hp, ∇ui+λ 2kAu −gk2 2−δp(p)) (B.9) En primera instancia, se obtiene que: ∂u(hp, ∇ui) = ∂u(−hu, divpi) = −divp ∂u(δp(p)) = 0 ∂u(λ 2kAu −gk2 2)), requiere un desarrollo m´as amplio: ∂ukAu −gk2 2=∂u(Au −g)T(Au −g) (Au −g)T(Au −g) = ((Au)T−gT(Au −g)) ((Au)T−gT(Au −g)) = (uTAT−gT)(Au −g) (uTAT−gT)(Au −g) = uTATAu −uTATg−gTAu +gTg se realiza el cambio B=ATAy se calcula as´ı el primer t´ermino: ∂u(uTATAu) = ∂uuTBu ∂u(uTBu)=(B+BT)u Se invierte el cambio: ∂u(uTATAu) = (ATA+ (ATA)Tu ∂u(uTATAu)=(ATA+ATA)u ∂u(uTATAu) = 2ATAu Se calcula el resto de t´erminos: ∂u(−uTATg−gTAu) = ∂u(−VTg−gTV) = ∂u(−2gTV) ∂u(−2gTV) = ∂u(−2gTAu) = 2ATg gTg= 0 Al final por tanto: λ 2kAu −gk2 2) = λ(ATAu −ATg) ∂uE(u, p) = −divp +λ(ATAu −ATg) (B.10) Se sustituye la derivada por su definici´on: un−un+1 τ=−divp +λ(ATAun+1 −ATg) (B.11) (un−un+1) = −τdivpn+1 +τλ(ATAun+1 −ATg) (B.12) Finalmente, la actualizaci´on puede derivarse de la siguiente expresi´on: un+1(I+τλATA) = un+τdivpn+1 +τλATg(B.13) 37 Secci´on B.5 B. Deducciones matem´aticas B.5. Derivaci´on de la funci´on de energ´ıa para el problema de la deconvoluci´on Partiendo de la definici´on de funci´on de energ´ıa, y considerando que la ∂uE(u, p) = 0, se obtiene que: ∂uE(u, p) = ∂u(−hp, ∇ui+λ 2kAu −gk2 2−δp(p)) (B.14) Sustituyendo t´erminos por su valor: ∂uE(u, p) = divp +λAT(Au −g) (B.15) Si se coloca ˆu=u−τdivp en el algoritmo, se puede escribir que divp =u−ˆu τ ∂uE(u, p) = u−ˆu τ+λAT(Au −g) = 0 (B.16) Se realiza el cambio de la matriz Aque como ya se ha comentado es muy complicada de implementar y se sustituye por k∗, que es la convoluci´on. Tiene el significado f´ısico de pasar la m´ascara por la imagen. u−ˆu τ+λk∗∗(k∗u−g) = 0 (B.17) Donde k∗es la convoluci´on y k∗∗es el conjugado que convoluciona, es la equivalencia a la ATanterior. En este momento se aplica la Transformada de Fourier a la izquierda y derecha de la ecuaci´on. F(u)−F(ˆu) τ+λF(k)∗(F(k)F(u)−F(g)) = F(0) (B.18) Al aplicar Fourier, las convoluciones se convierten en productos. Se sigue operando hasta dejar a un lado de la ecuaci´on F(u): F(u)−F(ˆu) + τλF(k)2F(u)−τλF(k)∗F(g)) = 0 (B.19) F(u)(1 + τλF(k)2) = F(ˆu) + τλF(k)∗F(g) (B.20) F(u) = F(ˆu) + τλF(k)∗F(g) (1 + τλF(k)2)(B.21) Se realiza la transformada inversa y as´ı se obtiene la suci´on final: u=F−1F(ˆu) + τλF(k)∗F(g) (1 + τλF(k)2)(B.22) 38 B. Deducciones matem´aticas Secci´on B.5 Es necesaria una kdel tama˜no de la imagen para poder realizar la transformada de Fourier y que su producto tenga sentido. De esta manera, siguiendo el algoritmo de actualizaci´on primal dual habitual, y mediante el c´alculo de la FFT en CUDA, se puede realizar la deconvoluci´on de una forma muy r´apida. 39 Ap´endice C Tipos de ruido C.1. El ruido Gaussiano El ruido Gaussiano asigna a cada p´ıxel de la imagen una funci´on gaussiana centrada en el valor de dicho p´ıxel. Figura C.1: Comparaci´on de la funci´on gaussiana para diferentes valores de σ De esta forma, para valores de σpeque˜nos, el valor del p´ıxel cambiar´a muy poco el la mayor´ıa de los casos. Para valores de σgrandes, los p´ıxeles cambiar´an mucho su valor en media. El c´odigo que crea el ruido gaussiano es el siguiente: 1Mat addNoise (Mat &I , float mean , float std ) { 3Mat Ino isy ; 5Mat I n o i s e ( I . rows , I . c ol s , CV 32FC1) ; randn ( Inoise , mean , std ) ; 7 41 Secci´on C.2 C. Tipos de ruido I no i sy = I+I n o i s e ; 9cv : : max( Inoisy , 0.0 , Inoisy ) ; cv : : min ( Inoisy , 1.0 , Inoisy ) ; 11 13 return I noi sy ; } La funci´on randn genera un valor aleatorio de media mean y de σ std. Se impone mean igual a cero, de forma que la matriz Inoise creada es una funci´on centrada en cero y de la desviaci´on est´andar asignada. Luego se suma Inoise a la matriz original para crear el ruido, acotando los valores m´aximos y m´ınimos que pueda tener. C.2. El ruido de sal y pimienta Recibe este nombre por su similitud a la sal y la pimienta al colocar los p´ıxeles de color blanco y negro. Es una forma de generar espurios y comprobar lo robusto del algoritmo ante su aparici´on. El par´ametro que recibe como entrada es el porcentaje relativo de p´ıxeles espurios que se desea aparezcan en la imagen. Se supone misma cantidad de p´ıxeles generados blancos que negros. El c´odigo para generarlo es el siguiente: Mat Add s al y pim ient a (Mat o r i g i n a l , float porcentaje ) 2{int cols = o r i g i n a l . c ol s ; 4int rows = o r i g i n a l . rows ; bool rotar = f a l s e ; 6cv : : Mat I = o r i g i n a l ; f o r (int i = 0; i<c o l s ; i++) 8f o r (int j =0; j<rows ; j++) { 10 int valor = rand () %100; i f ( ( float) valor <= porcentaje ) 12 {i f ( rot a r == false) 14 {I . at<float >(i , j ) = 0 ; 16 rotar = true ; } 18 e l s e { 20 I . at<float >(i , j ) = 1 ; rotar = f a l s e ; 42 C. Tipos de ruido Secci´on C.2 22 } } 24 } return I ; 26 } Para cada p´ıxel se genera un valor aleatorio que al hacer %100 se calcula el resto de ese valor con 100. As´ı, el n´umero resultante es un valor entre 00 y 99. Si el n´umero aleatorio generado es inferior al porcentaje, significa que ese p´ıxel va a ser un espurio. La variable booleana rotar es quien se encarga de decidir si el espurio ser´a blanco o negro, cambiando su valor para el siguiente espurio. 43 Ap´endice E Programaci´on La estructura b´asica de el archivo donde se ejecutan los algoritmos matem´aticos deducidos en los apartados 2,3,4 y 5 sigue una secuencia como la expuesta en el apartado anterior. En este apartado se explica c´omo se realizan las actualizaciones del Primal y del Dual. Se comienza con un ejemplo: 1global void kern el upda te dua l ( f l o a t 2 ∗p dev , float ∗u hat dev , in t cols , in t rows , float sigma ) { 3in t x = blockDim . x∗blockIdx . x + threadIdx . x ; in t y = blockDim . y∗blockIdx . y + threadIdx . y ; 5 i f ( x <c o l s && y <rows ) 7{float tx = x + 0.5 f ; 9float ty = y + 0.5 f ; in t o f f s e t = y∗c o l s+x ; 11 float gx = tex2D ( tex u hat , tx +1, ty )−tex2D ( tex u hat , tx , ty ) ; 13 float gy = tex2D ( tex u hat , tx , ty+1)−tex2D ( tex u hat , tx , ty ) ; 15 i f (x == ( cols −1)) { 17 gx = 0.0 f ; } 19 i f ( y == ( rows −1)) 21 {gy = 0.0 f ; 23 } 25 f l o a t 2 p = tex2D ( tex p , tx , ty ) ; p . x = p . x + sigma∗gx ; 51 E. Programaci´on 27 p . y = p . y + sigma∗gy ; 29 float norma = s q r t f (p . x∗p . x+p . y∗p . y ) ; float den = fmaxf ( 1 . 0 f , norma) ; 31 p . x = p . x/den ; p . y = p . y/den ; 33 p dev [ o f f s e t ] = p ; 35 } } El primer paso es la definici´on de los valores xeyen funci´on del bloque y el Thread concreto de ese bloque en el que se vaya a ejecutar el c´odigo. Esta funci´on actualiza el dual, quien, seg´un el algoritmo expuesto, requiere del gradiente de ˆu, que ha sido guardado en textura. Como la textura hab´ıa sido definida previamente con dos dimensiones, es muy sencillo tomar los valores del gradiente ya que se hace como se har´ıa en una matriz de p´ıxeles: El p´ıxel de la derecha menos el p´ıxel en el que se ejecuta la funci´on es el gradiente en x y el p´ıxel de abajo menos el p´ıxel actual es el gradiente en y. ptambi´en se encuentra en textura, lo que hace m´as r´apido el algoritmo de actualizaci´on. Una vez calculada la norma, la p.x y la p.y, se debe enviar mediante la instrucci´on pdev[offset] = p; a la variable pdev, para que su nuevo valor se almacene en textura y en la siguiente iteraci´on se pueda proceder de la misma manera. global void kernel update Primal ( float ∗u dev , float ∗u hat dev , float ∗g dev , f l o a t 2 ∗p dev , float tau , float lambda , float theta , in t cols , in t rows ) 2 { 4 in t x = blockDim . x∗blockIdx . x + threadIdx . x ; 6in t y = blockDim . y∗blockIdx . y + threadIdx . y ; 8float uhat ; float u ; 10 i f ( x <c o l s && y <rows ) 12 {float tx = x + 0.5 f ; 14 float ty = y + 0.5 f ; in t o f f s e t = y∗c o l s+x ; 16 float k i j x = 0 , k i j y = 0 , k i 1 j x = 0 , k i j 1 y = 0 ; 18 f l o a t 2 p y i j 1 = tex2D ( tex p , tx , ty −1) ; 52 E. Programaci´on f l o a t 2 p i j = tex2D ( tex p , tx , ty ) ; 20 f l o a t 2 p x i 1 j = tex2D ( tex p , tx −1,ty ) ; 22 i f (y >0) { 24 k i j 1 y = p y i j 1 . y ; } 26 i f (y <rows −1) 28 {k i j y = p i j . y ; 30 } 32 i f (x >0) { 34 k i 1 j x = p x i 1 j . x ; } 36 i f (x <cols −1) 38 {k i j x = p i j . x ; 40 } 42 float u old = tex2D ( tex u , tx , ty ) ; float g = tex2D ( tex g , tx , ty ) ; 44 float dip = kijx −k i 1 j x+kijy −k i j 1 y ; uhat = u old +tau∗dip ; 46 u = ( uhat+tau∗lambda∗g ) /(1+ tau ∗lambda) ; uhat = u+theta ∗(u−u old ) ; 48 u dev [ o f f s e t ] = u ; u hat dev [ o f f s e t ] = uhat ; 50 } } La manera de actualizar el Primal cambia sustancialmente en funci´on del m´etodo. Este ejemplo corresponde al TV-ROF. El primer paso com´un a todos los m´etodos es la realizaci´on de la divergencia de p. Para ello se definen las variables kijx, kijy ,ki1jx ykij1y, que se refieren al valor de pen el p´ıxel actual en p.x y en p.y, al p´ıxel de la derecha en el caso de p.x y al p´ıxel de abajo en el caso de p.y. La divergencia se calcula como se muestra en el c´odigo. A continuaci´on se debe guardar el valor de uold para la futura actualizaci´on de ˆu. A partir de ah´ı, se debe seguir el algoritmo. Finalmente, tanto ˆucomo udeben enviarse a la variable que est´a en textura para que la actualizaci´on sea efectiva. 53 Ap´endice F Resultados Se exponen a continuaci´on los resultados y conclusiones obtenidas a partir de los diferentes algoritmos. F.1. Resultados Denoising F.1.1. Evaluaci´on del modelo TV-ROF El primer experimento se realiza sobre una imagen a la que se le ha adicionado ruido gaussiano con desviaci´on est´andar σr. El resultado de la figura F.1 se ha obtenido al ejecutar el algoritmo TV-ROF. (a) Imagen original (b) Imagen ruidosa σr= 0,1 (c) Imagen resultado λ= 8,0 TV-ROF Figura F.1: Primer experimento del modelo TV-ROF Como es de esperar, el uso de la norma del gradiente no penaliza las discontinuidades en la imagen pudi´endose mantener los bordes. As´ı mismo, las regiones con intensidades cercanas se suavizan. Durante la ejecuci´on del algoritmo se ha 55 Secci´on F.1 F. Resultados llevado a cabo un an´alisis de verificaci´on del GAP definido para cada modelo. La figura F.2 muestra el resultado para un numero fijo de iteraciones. Se puede notar que el GAP tiende a cero con apenas 30 iteraciones. Figura F.2: Evoluci´on del GAP respecto a las iteraciones Figura F.3: Imagen ruidosa σ= 0,3 Figura F.4: Imagen resultado λ= 8,0 TV-ROF Se puede observar que al aumentar dr´asticamente el porcentaje de ruido en la imagen el algoritmo no puede recuperar totalmente la imagen original. No obstante, ´esto ocurrir´a para cualquier algoritmo de filtrado. Los par´ametros ´optimos del m´etodo se pueden obtener mediante una simulaci´on. As´ı, se obtienen las gr´aficas: 56 F. Resultados Secci´on F.1 (a) Dependencia SnR respecto a λen TV-ROF (b) Dependencia iteraciones hasta convergencia respecto a τen TV-ROF Figura F.5: Evaluaci´on de par´ametros λyτ De acuerdo a los resultados, el valor ´optimo para λser´ıa de 12 tal como se observa en F.5 a) y un valor ´optimo de τser´ıa de 0.25, de acuerdo a la gr´afica F.5 b). Los valores ´optimos tambi´en depender´an de la cantidad de ruido introducido, siendo λespecialmente sensible a esta variaci´on. La figura F.4 se ha calculado para σr= 0.3. Si se realiza la evaluaci´on de SnR de los par´ametros para este ruido, el valor de λvar´ıa como se muestra en la figura . Figura F.6: dependencia SnR respecto a λen TV-ROF para un ruido de σ= 0.3. El ´optimo de λha variado a 4.2. Con ese nuevo valor se obtiene una resultado distinto. La figura F.7 muestra una comparaci´on entre el resultado previamente obtenido con λ= 8 y el valor ´optimo: 57 Secci´on F.1 F. Resultados (a) Imagen resultado λ= 8,0 TV-ROF (b) Imagen resultado λ= 4,2 TV-ROF Figura F.7: Comparaci´on de resultados cuando se incrementa el ruido en la imagen y se aplican valores diferentes de λ. Cuanto menor es el error, mayor es la fiabilidad de la imagen sobre la que el algoritmo act´ua, por lo que el resultado ser´a muy similar a dicha imagen de forma que λser´a muy elevada. Cuanto mayor es el error, m´as importante ser´a el t´ermino de gradiente, as´ı que el λ´optimo deber´a ser menor. F.1.2. Evaluaci´on del modelo Huber-ROF En el caso del modelo de Huber-ROF, observar que el algoritmo primal dual converge mucho m´as r´apido a la soluci´on en comparaci´on con el modelo TV-ROF. La figura F.8 ilustra el resultado obtenido. Las figuras F.9 y F.10 representan c´omo var´ıa el GAP y la relaci´on SnR respecto a un par´ametro. Se puede concluir por tanto que este par´ametro siempre converge en las mismas iteraciones. El valor ´optimo de λes de 7.5 y el de αes de 0.025. (a) Imagen ruidosa σr= 0,1 (b) Imagen resultado λ= 5,0 Huber-ROF Figura F.8: Resultado por el algoritmo primal dual para el modelo Huber-ROF. 58 F. Resultados Secci´on F.1 Figura F.9: Evoluci´on del GAP por iteraci´on para el modelo Huber-ROF. (a) Dependencia SnR respecto a λ (b) Dependencia SnR respecto a α Figura F.10: Evaluaci´on de par´ametros para el modelo Huber-ROF F.1.3. Evaluaci´on del modelo TV-L1 El algoritmo TV-L1 contiene una diferencia esencial respecto a los dos anteriores: La dificultad del c´alculo del GAP. A diferencia de los otros modelos, se ha calculado el valor de la funci´on de energ´ıa. Dicha funci´on de energ´ıa debe alcanzar un valor m´ınimo una vez el algoritmo haya convergido a la soluci´on ´optima. La figura F.11 muestra la imagen de resultado obtenida mientras la figura F.12 muestra la evoluci´on de la energ´ıa para el modelo TV-L1. 59 Secci´on F.2 F. Resultados (a) SnR, s = 4 (b) SnR, s = 8 Figura F.20: Evaluaci´on de precisi´on respecto a la variaci´on del parametro λ, para s = 4 y s = 8. En el experimento se ha tomado µ= 100. Aplicando la ecuaci´on anterior, se deduce que el valor m´aximo de λque asegura convergencia es de 1828.57. En el experimento se ha superado ese valor y el m´etodo ha convergido, pero m´as lentamente. Observado la gr´afica, el valor que se escoger´ıa es el ´optimo m´as lejano de la no convergencia, es decir, λ= 800. En el caso de factor igual a 8, el valor m´aximo de λes de 6606.45, por lo que el experimento se acota al valor 6000. Se encuentra el valor 1800 como ´optimo te´orico en este caso. En este algoritmo λ decrece muy lentamente una vez alcanzado el ´optimo, sin embargo conviene no tomar λmuy alto, pues el m´etodo podr´ıa no converger. A continuaci´on, se comprueba el valor ´optimo de τ. De nuevo se va realizar el mismo experimento para dos escalados diferentes y se interpretar´an los resultados obtenidos: 66 F. Resultados Secci´on F.2 (a) SnR, s = 4 (b) SnR, s = 8 Figura F.21: Evaluaci´on de precisi´on para el par´ametro τ, para s = 4 y s= 8 En el primer caso, el valor ´optimo de τest´a en torno a 0.025, donde el m´etodo converge en unas 800 iteraciones. ´ Esto se debe al error m´ınimo establecido para alcanzar la convergencia. Si es error admisible se aumenta, el m´etodo es muy r´apido pues disminuye el numero de iteraciones. En el segundo caso, el valor ´optimo de τes similar, de 0.02, sin embargo las iteraciones hasta la convergencia aumentan considerablemente. Al tener un factor mayor, las iteraciones se multiplican. θes el ´ultimo par´ametro a analizar de acuerdo a su influencia en la aceleraci´on de la convergencia. En la figura F.22 se comprueba c´omo evoluciona su valor para s = 4. 67 Secci´on F.3 F. Resultados Factor de ampliaci´on Escala Tiempo m´ınimo de convergencia 4 100x100 243ms 4 50x50 90ms 8 50x50 1589ms Tabla F.2: Tiempo de ejecuci´on para el algoritmo de zooming (a) Resultado ´optimo (b) SnR Figura F.22: Evaluaci´on de la precisi´on para el par´ametro θ, para s = 4. Se observa que θdebe ser lo m´as peque˜na posible para que disminuya las iteraciones. El m´ınimo para asegurar la convergencia es de 0.5, por lo que se mantiene ese valor como ´optimo. Adicionalmente se han realizado tres experimentos para obtener un valor aproximado del tiempo de c´omputo del algoritmo. La tabla F.2 resume los resultados para diferentes valores de escala. Se puede por tanto comprobar que el tiempo de c´omputo aumenta sustancialmente con el tama˜no de la imagen, pero es m´as sensible ante variaciones en el factor de ampliaci´on. F.3. Resultados Deconvolution Dada una imagen, se aplica una m´ascara de motion blurring de 10 p´ıxeles de longitud de movimiento y 45ode inclinaci´on. La figura F.23 expone la imagen original, la imagen degradada con blurring y la imagen de resultado obtenida despu´es de la deconvoluci´on. 68 F. Resultados Secci´on F.3 (a) Imagen original (b) Imagen con blurring, dp= 10, α= 45o (c) Imagen resultado λ= 1000 Figura F.23: Resultado obtenido tras la deconvoluci´on aplicando el algoritmo primal dual. Para evaluar la eficacia del algoritmo de deconvoluci´on, se ha incrementado el efecto de blurring. La figura muestra el resultado para una m´ascara lineal de dp=100 (a) Imagen original (b) dp= 100, α= 45o(c) λ= 1000 Figura F.24: Resultado de deconvoluci´on para una imagen degradada con dp=100 Una primera conclusi´on que se puede obtener tras observar la evoluci´on de la funci´on de energ´ıa (ver figura F.25) es que ´esta disminuye muy poco debido a la utilizaci´on de la transformada de Fourier que transforma el resultado del dominio espacial al de la frecuencia. 69 Secci´on F.3 F. Resultados (a) Evoluci´on de FdU, dp= 100 (b) Evoluci´on de λ(c) Evoluci´on de τ Figura F.25: Evaluaci´on de par´ametros Si se observan las gr´aficas del par´ametro λse deduce que lo m´as correcto es tomar una λinfinita, pues ´esta no deja de aumentar. Sin embargo, ¿qu´e sucede en el caso de que una imagen se vea afectada al mismo tiempo por blurring y ruido? Imaginemos un 5 % de la imagen con p´ıxeles blancos y negros. Al transformar la imagen al dominio frecuencial y optimizar ah´ı el algoritmo, se intentar´a mantener a toda costa la cantidad de p´ıxeles ruidosos. Adem´as, al aplicar la m´ascara de forma inversa en la resoluci´on del algoritmo, el ruido se expande. Como conclusi´on, se utilizar´a λ= 1000 como par´ametro ´optimo. En cuanto al par´ametro τ, se observa que para un valor de 0.012 es ´optimo. El par´ametro θno se incluye pues no provoca variaci´on en las iteraciones hasta convergencia en su rango de aplicaci´on. Para poder mostrar c´omo evoluciona el algoritmo en caso de que haya ruido en la imagen, se ha habilitado la posibilidad de generarlo en la interfaz. De esta forma, se realiza un blurring sobre la imagen y adem´as se a˜nade un ruido gaussiano de σ= 0.1. La imagen da˜nada es la siguiente: Figura F.26: Imagen con blurring y ruido de σ= 0.1 70 F. Resultados Secci´on F.3 El ´optimo te´orico de la deconvoluci´on es λ=∞, se asigna λ= 10000. El ´optimo te´orico para solucionar el ruido est´a en torno a 10. De esta forma, para comprender mejor el algoritmo, se hacen cuatro experimentos: λ= 10, λ= 100, λ= 1000 y λ = 10000. (a) λ= 10 (b) λ= 100 (c) λ= 1000 (d) λ= 10000 Figura F.27: Resultado de deconvoluci´on para diferentes valores de λal adicionar ruido gaussiano. Si se observa la imagen F.27 a), no ha reparado su blurring. Sin embargo, ya no aparece nada de ruido y las superficies han sido homogeneizadas. En la imagen F.27 b) se ha removido la mayor parte del ruido y la degradaci´on por movimiento ha sido casi reparada. La imagen F.27 c) muestra lo que sucede si λtoma un valor demasiado alto, el ruido no s´olo no desaparece, sino que es deconvolucionado con 71 Secci´on F.3 F. Resultados el m´etodo y se ha expandido. La imagen F.27 d) muestra el resultado para un λ = 10000 para el cual el blurring ha sido eliminado, pero el ruido persiste. El problema que se trata de resolver en este apartado es el del blurring, no el ruido, por lo que una σ= 0.1 es muy grande para este punto. Se a˜nade un ruido de σ= 0.01, inapreciable para el ser humano y se realiza el mismo experimento. Simulando una situaci´on real, se empieza con una λ= 1000 que hab´ıa sido sugerida en el an´alisis de par´ametros, ya que no se sabe que existe ruido. Figura F.28: Imagen con gran blurring y ruido de σ= 0.01 (a) λ= 1000 (b) λ= 10000 72 F. Resultados Secci´on F.4 (c) λ= 100000 (d) λ= 1000000 Figura F.29: Resultado de deconvoluci´on para differentes valores de λal adicionar ruido gaussiano. Como el resultado de λ= 1000 no repara el blurring completamente, se piensa en aumentar el valor de λ. Conforme se aumenta en un factor de 10, el resultado cada vez es peor. Este ´ultimo an´alisis demuestra que pese a que un an´alisis te´orico muestre que el λ´optimo es ∞, conviene no aumentarlo por encima de cierto valor, ya que puede haber algo de ruido inapreciable que puede empeorar dr´asticamente el resultado. En cuanto a la eficiencia del m´etodo, se ha evaluado tambi´en el tiempo de ejecuci´on. N´otese que los tiempos no var´ıan demasiado en funci´on de la longitud de movimiento en p´ıxeles (Ver tabla F.3). Longitud de movimiento Tiempo m´ınimo de convergencia 10 17 ms 100 55 ms Tabla F.3: Tiempo de ejecuci´on para el m´etodo de convoluci´on. F.4. Resultados Image Inpainting La imagen utilizada en este apartado vuelve a ser Lena.jpg, la misma que se utiliz´o en el apartado de denoising. Como primera simulaci´on, se da˜na el 50 % de la imagen y se recupera seg´un el algoritmo de inpainting. 73 Secci´on F.4 F. Resultados Figura F.30: Imagen da˜nada Figura F.31: Imagen arreglada λ= 640 Figura F.32: Evoluci´on de FdU El algoritmo es muy potente, ya que es capaz de recuperar la imagen pr´acticamente en su totalidad. Hay que decir que en la realidad no se encontrar´ıa un resultado tan fiel a la realidad, puesto que en la funci´on que repara la imagen se debe introducir qu´e p´ıxeles son espurios. En esta aplicaci´on se maneja la informaci´on perfecta de qu´e p´ıxeles son defectuosos, lo que en la realidad no es posible. Para que este algoritmo tenga una aplicaci´on m´as realista, se debe crear un algoritmo que sea capaz de calcular puntos espurios con precisi´on. Como ejemplo de la potencia del algoritmo si maneja informaci´on perfecta, se introduce una imagen da˜nada al 85 %. 74 F. Resultados Secci´on F.4 Figura F.33: Imagen da˜nada Figura F.34: Imagen arreglada λ= 640 Una vez se han mostrado los resultados del algoritmo, se expone la optimizaci´on de par´ametros: Figura F.35: Evoluci´on de λFigura F.36: Evoluci´on de τ Al igual que en el caso anterior, vuelve a aparecer una λinfinita como ´optima. De nuevo, hay que tomar este resultado con precauci´on, porque en esta aplicaci´on se maneja informaci´on perfecta. Si no se detectan todos los p´ıxeles espurios. Para demostrarlo, se presenta a continuaci´on una imagen da˜nada al 50 %, donde tan s´olo se ha detectado la mitad de da˜nados y se aplica una λmuy elevada: 75 Secci´on G.2 G. Manual de usuario G.2. Ventana Denoising Cuando se escoge la opci´on Denoising se abre la siguiente interfaz de usuario: Figura G.2: Interfaz denoising al inicio de ejecuci´on Punto 1. El bot´on lleva asignada la elecci´on de una imagen. Al pulsarlo, se abre una ventana que permite navegar por las carpetas del sistema y escoger la imagen deseada. Autom´aticamente esta imagen se guarda en la memoria del programa y se muestra en el hueco habilitado para Imagen original. Siempre que se introduce una imagen, se da˜na o se arregla, cuando se muestra en su lugar de la interfaz, se modifica su tama˜no para que se ajuste al hueco que le corresponde. Punto 2. En el desplegable inferior al bot´on Adicionar ruido se escoge el tipo de ruido que se desea: un ruido gaussiano de valor σel que se introduce en el espacio habilitado para ello, o el ruido de sal y pimienta, al que se le introduce como par´ametro el porcentaje de imagen da˜nada que se desea. Este porcentaje se introduce en el hueco donde ahora se ve la palabra sigma. Al cambiar el tipo de ruido, se cambia el nombre sigma por % para hacerlo m´as intuitivo para el usuario. La imagen da˜nada se guarda en memoria y se muestra en el widget Imagen ruidosa. Punto 3. La ejecuci´on del programa depende del algoritmo utilizado. El desplegable que se encuentra en esta zona permite elegir entre TV-ROF, Huber-ROF y TV-L1. Una vez que se escoge un m´etodo se habilitan espacios para introducir los par´ametros en funci´on del m´etodo (por ejemplo, en TV-ROF, θdepende de τyλ, por lo que s´olo se habilitan estos dos ´ultimos) con unos par´ametros que convergen cargados por defecto. Una vez se pulsa el bot´on, se ejecuta el algoritmo en la tarjeta gr´afica y muestra por pantalla la imagen resultado. 82 G. Manual de usuario Secci´on G.2 Punto 4. La gr´afica muestra en escala logar´ıtmica o decimal en funci´on del desplegable la evoluci´on del GAP o funci´on de energ´ıa seg´un el caso de los algoritmos. As´ı, se puede observar que en cada iteraci´on el GAP decrece hasta 0, valor que toma cuando el m´etodo ha convergido. Por su parte, la funci´on de energ´ıa nunca valdr´a 0, alcanzar´a un m´ınimo en el que se mantendr´a cuando haya convergido. Punto 5. Tras la ejecuci´on del algoritmo se muestra por pantalla la imagen resultado, la evoluci´on del GAP o funci´on de energ´ıa y el tiempo de c´omputo. Este tiempo de c´omputo es el que invierte el algoritmo en realizar los c´alculos propios del algoritmo, descontando el tiempo invertido en calcular el GAP, ya que ´este es m´as de cien veces mayor y es absurdo calcular el GAP en cada iteraci´on para saber si el m´etodo ha convergido. Si no se sabe con exactitud las iteraciones necesarias para converger, es computacionalmente m´as barato aumentar en uno el orden de las iteraciones previstas y no calcular el GAP, que es muy costoso. Punto 6. Las pesta˜nas de la intefaz. Aqu´ı se puede cambiar entre las diferentes opciones. A continuaci´on se muestra la interfaz una vez se ha reparado una imagen. El ruido generado es de sal y pimienta y se soluciona mediante el ´unico algoritmo que lo soluciona: El TV-L1. Se puede observar c´omo σse ha cambiado por % al cambiar el tipo de ruido y el resultado en logar´ıtmico. En TV-L1 se calcula la funci´on de energ´ıa y no el GAP as´ı que tambi´en ha cambiado en los ejes de la gr´afica. El TV-L1 es un m´etodo muy r´apido y con tan s´olo 100 iteraciones y ante tan poco ruido, tarda muy poco tiempo en finalizar la ejecuci´on. Figura G.3: Interfaz denoising tras simulaci´on Se muestra ahora la pesta˜na de Evaluaci´on de Par´ametros vac´ıa, para explicar su funcionamiento. 83 Secci´on G.2 G. Manual de usuario Figura G.4: Interfaz denoising par´ametros al inicio de ejecuci´on Punto 1. De nuevo se puede escoger el m´etodo que se desee. Adem´as, se debe elegir qu´e par´ametro se va a evaluar colocando un tick donde se desee. Una vez se elige el par´ametro, se asignan unos par´ametros por defecto en s´ı mismo y los dem´as. En el escogido se habilitan tanto el m´ın, m´ax y el paso. En los dem´as, s´olo se habilita el m´ın, donde se coloca el valor asignado al otro par´ametro. Se eval´ua el algoritmo para el valor de cada par´ametro seleccionado entre el m´ın y el m´ax, iterando seg´un el paso. La cantidad de veces que se ejecuta el algoritmo es por tanto pmax−pmin ppaso siempre redondeando hacia arriba para asegurar que se eval´uan los dos extremos. As´ı, conforme mayor sea la diferencia entre el min y el max y menor sea el paso, m´as tiempo tardar´a en ejecutarse la evaluaci´on de par´ametros. Tambi´en se debe considerar que en algunos algoritmos hay l´ımites de convergencia en algunos de ellos. Si los l´ımites son demasiado amplios podr´ıan originarse problemas de no convergencia en alg´un m´etodo. Punto 2. Tras la evaluaci´on, aparecer´a el ratio Se˜nal-Ruido si se escoge λoα y las iteraciones hasta convergencia para τyθ. La gr´afica representar´a el valor evaluado en el rango del par´ametro. Se muestra un ejemplo solucionado. Se trata de la evaluaci´on de αen el algoritmo de Huber ante un ruido gaussiano. 84 G. Manual de usuario Secci´on G.3 Figura G.5: Interfaz denoising par´ametros al final de ejecuci´on G.3. Ventana Zooming Una vez pulsada la opci´on de Zooming en la ventana principal, aparece la siguiente interfaz de usuario. Figura G.6: Interfaz zooming al inicio de la ejecuci´on Punto 1. La organizaci´on de la interfaz ha cambiado un poco. En este caso, aparece un hueco grande donde aparecer´a la imagen ampliada. Al ser negro en el inicio no se puede ver la zona donde aparecer´a la imagen peque˜na, que es la 85 Secci´on G.3 G. Manual de usuario esquina superior izquierda de la zona negra. Zooming consiste en una ampliaci´on de la imagen de partida. Si esta imagen es de 100x100 y se amplia con un factor de 5, la imagen resultante es de 500x500. En una interfaz compacta es imposible de modelar, pues la imagen resultante se podr´ıa solapar con la gr´afica o los par´ametros si hubiera un factor de escalado muy grande. La decisi´on que se tom´o para representar de una forma gr´afica la ampliaci´on es fijar el tama˜no de la grande en la interfaz, variando en tama˜no de la peque˜na. El algoritmo trabaja con los tama˜nos reales pero a la hora de mostrar por pantalla se ajusta a lo dicho anteriormente. As´ı se puede apreciar c´omo var´ıa la relaci´on de tama˜nos. Se debe introducir un factor de ampliaci´on lo suficientemente grande respecto a la imagen original para que supere el tama˜no de 400x400 que tiene el hueco para la imagen grande. Si es m´as peque˜na el algoritmo funciona igual y la imagen queda representada de la misma manera, pero para ajustarse al tama˜no se realiza una interpolaci´on lineal que provocar´a que el resultado final sea una mezcla del primal dual (lo que se realiza en el algoritmo) y una interpolaci´on lineal (lo que se trata de evitar al realizar el primal dual). Punto 2. El factor de escalado introducido por defecto es 2. Una vez que se edita, el tama˜no asignado a la imagen peque˜na var´ıa, aunque si no se ha ejecutado el programa no se aprecia al ser la zona de color negro. Se expone a continuaci´on la interfaz resultado con dos diferentes escalados: 2 y 10. Figura G.7: Interfaz zooming al final de la ejecuci´on, s = 2 86 G. Manual de usuario Secci´on G.3 Figura G.8: Interfaz zooming al final de la ejecuci´on, s = 10 Con los ejemplos se entiende algo mejor c´omo se debe interpretar la interfaz. Para distintos escalados, se observa la diferencia existente entre la imagen original y la imagen ampliada. Adem´as, se pueden hacer pruebas con el par´ametro λ, observando que aunque en la imagen de escalado 2 parece ´optimo, en la imagen de escala 10 el resultado est´a m´as difuminado. Tambi´en a igual valor de τyθ, se puede observar en la gr´afica y tiempo de c´omputo la diferencia en iteraciones hasta convergencia y en tiempo de c´omputo existente en funci´on del escalado. Se muestra a continuaci´on la interfaz de la evaluaci´on de par´ametros antes y despu´es de la ejecuci´on: Figura G.9: Interfaz zooming par´ametros al inicio de la ejecuci´on Punto 1. De nuevo se puede escoger qu´e par´ametros evaluar. En este caso el resultado no depende de αde forma que el par´ametro ha desaparecido. Los diferentes valores del resto son inicializados con par´ametros que convergen y, al igual 87 Secci´on G.3 G. Manual de usuario que en la interfaz denoising, pueden variar como desee el usuario. Es cr´ıtico tener cuidado con la introducci´on del par´ametro λ, que mediante una dependencia con τyspuede provocar la no convergencia del m´etodo y podr´ıa dar como resultado una imagen completamente blanca. Es habitual en evaluaci´on de par´ametros colocar l´ımites de evaluaci´on muy separados con un paso muy peque˜no con el fin de obtener una gr´afica precisa. Estos l´ımites deben introducirse con cuidado o la gr´afica de SnR podr´ıa carecer de sentido. Si se observa en dicha gr´afica algo il´ogico tan como una ca´ıda o una subida muy brusca en el valor de SnR, probablemente sea porque el m´etodo no ha convergido y se puede saber a partir de qu´e valor de λsucede este fen´omeno. Se muestra a continuaci´on el resultado del m´etodo convergiendo y no convergiendo: Figura G.10: Interfaz zooming par´ametros al final de la ejecuci´on que converge Figura G.11: Interfaz zooming par´ametros al final de la ejecuci´on que no converge Se observa para este segundo caso que el m´etodo a partir de λ= 3200 no converge. Seg´un la ecuaci´on: 88 G. Manual de usuario Secci´on G.4 λ < µs4/(s2−2) (G.1) Sabiendo que µ=1 τ, que τ= 0.01 y que s= 4, s´olo se asegura convergencia para valores de λmenores a 1828,57. Este valor es aproximadamente el punto en el que en el primer experimento la SnR empieza a decrecer. Se debe calcular este valor antes de lanzar ninguna simulaci´on, porque sino podr´ıa aparecer un resultado como el de la figura G.11. G.4. Ventana Deconvolution Una vez pulsada la opci´on de Deconvolution en la ventana principal, aparece la siguiente interfaz de usuario. Figura G.12: Interfaz deconvolution al inicio de la ejecuci´on La interfaz es calcada a la de denoising. La diferencia reside en que donde antes se generaba ruido, ahora se genera movimiento. El valor de la longitud en p´ıxeles de movimiento est´a acotado, ya que la imagen tiene una dimensi´on y en caso de exigir un movimiento demasiado grande, el algoritmo toma p´ıxeles exteriores a la imagen, volvi´endola completamente negra. Se sugiere no introducir un valor superior a 300. El ´angulo no tiene limitaciones. 89 Secci´on G.4 G. Manual de usuario Una vez ejecutado el programa se obtiene lo siguiente: Figura G.13: Interfaz deconvolution al final de la ejecuci´on La pesta˜na de los par´ametros sin ejecuci´on toma la siguiente forma: Figura G.14: Interfaz deconvolution par´ametros al inicio de la ejecuci´on Esta interfaz no introduce novedad, de nuevo no se deben introducir valores elevados de longitud. Si se ejecuta el algoritmo aparece un resultado de la siguiente manera: 90 G. Manual de usuario Secci´on G.4 Figura G.15: Interfaz deconvolution par´ametros al final de la ejecuci´on 91