Repositorio Institucional de Documentos

Proyecto Final de Carrera Ingenier´ıa en Inform´atica Multi-view image rendering for holographic stereogram printing David ´ Angel Valle Badenas Enero 2013 Departamento de Inform´atica e Ingenier´ıa de Sistemas Escuela de Ingenier´ıa y Arquitectura Universidad de Zaragoza Director: Atanas Gotchev Ponente: Francisco Jos´e Ser´on ii RESUMEN Desde que el ser humano ha sido capaz de reproducir objetos del mundo real en material sint´etico, el objetivo ha estado determinado por concebir sistemas de experiencias virtuales m´as precisos, desde las primeras fotograf´ıas estereogr´aficas hasta las contempor´aneas producciones cinematogr´aficas en 3D. Los productos actualmente disponibles en el mercado son consecuencia de investigaciones iniciadas d´ecadas atr´as. Sin embargo, nuevas exploraciones en tecnolog´ıas m´as complejas han dado lugar a nuevas inmersiones de realidad virtual. Los estereogramas hologr´aficos, los cuales proveen una perspectiva m´as amplia y una mayor cantidad de puntos de vista, se suponen ser el siguiente paso en la evoluci´on de los sistemas visuales 3D. Esta tecnolog´ıa no ha sido introducida en el mercado todav´ıa, sin embargo, el progreso realizado en el campo del hardware que permite una mayor cantidad de almacenamiento y proceso de datos y, mejores dise˜nos y algoritmos de software que consiguen resultados de forma m´as eficiente, har´an posible la disponibilidad de pantallas hologr´aficas al consumidor en un futuro cercano. Actualmente, los objetos 3D pueden ser manipulados f´acilmente y la posibilidad de obtener su informaci´on digital lleva al desarrollo de hologramas generados por ordenador basados en sistemas multi-vista. La captura y post-procesado de las vistas necesarias para crear un holograma generado por ordenador conlleva un coste de almacenamiento y computaci´on alto. Esta tesis afronta el problema de optimizaci´on del n´umero de vistas o c´amaras requeridas para reconstruir una escena compuesta por un n´umero de vistas totalmente densa. El m´etodo desarrollado parte de una escena 3D y un set completo de vistas. Solo algunas de esas vistas son utilizadas para reproducir la escena completa. A trav´es del an´alisis de las im´agenes en el dominio espacial y de frecuencia de Fourier, se han utilizado t´ecnicas de procesado de se˜nal para interpolar y regenerar la escena original. Se han considerado un n´umero diferente de vistas disponibles inicialmente, las cuales son procesadas para reconstruir la escena a posteriori, de manera que se pueda establecer un ratio entre el n´umero m´ınimo de c´amaras necesarias y el n´umero total de c´amaras que la escena contiene. Los resultados de este estudio sugieren, comparando la semejanza entre las escenas originales y sus versiones optimizadas, que la reconstrucci´on de todas las vistas puede ser obtenida a partir de muchas menos c´amaras y sin una p´erdida de calidad considerable. El siguiente texto muestra el an´alisis del problema y las t´ecnicas utilizadas para resolverlo. La memoria concluye con una valoraci´on de las t´ecnicas y herramientas desarrolladas, as´ı como de los resultados obtenidos. iii iv AGRADECIMIENTOS Este proyecto ha sido realizado en la Tampereen Teknillinen Yliopisto (Tampere University of Technology) en el departamento de Signal Processing. El proyecto supone un Proyecto Final de Carrera por la Universidad de Zaragoza, dentro del programa Erasmus. Quiero dar las gracias a mi supervisor de proyecto Dr. Atanas Gotchev por toda su dedicaci´on y paciencia, especialmente en aquellos momentos en los que mis conocimientos previos no eran suficientes para entender algunos temas. Tambi´en quiero agradecer al profesor Dr. Robert Bregovic por su tiempo dedicado a resolver mis dudas. Quiero agradecer a mi supervisor en Espa˜na, Dr. Francisco J. Ser´on, por su apoyo y tiempo, especialmente cuando s´olo me era posible estar en Espa˜na por un tiempo limitado. Finalmente, quiero dar las gracias a mi familia, amigos y personas especiales que me han apoyado a lo largo de todo el proyecto. v vi The scientist is not a person who gives the right answers, he is one who asks the right questions. Claude L´evi-Strauss vii viii Tabla de contenidos 1 Introducci´on 1 1.1 Contextoymotivaci´on ............................. 1 1.2 Alcance y objetivo del proyecto . . . . . . . . . . . . . . . . . . . . . . . . 6 1.3 Estructura.................................... 6 2 Conceptos te´oricos 7 2.1 M´etodo de holograma generado por ordenador . . . . . . . . . . . . . . . . 7 2.2 Representaci´ondeEPI............................. 8 2.3 Requisitos de la propuesta . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 3 An´alisis 13 3.1 Relaciones entre las caracter´ısticas de la escena y la representaci´on en el dominiodeFourier ............................... 13 3.2 An´alisis de la transformaci´on de dominio de las im´agenes EPI . . . . . . . 14 3.3 Filtrado de EPI y optimizaci´on del n´umero de vistas . . . . . . . . . . . . 17 4 Dise˜no del sistema multi-vista 20 4.1 Construyendo el sistema multi-vista para escenas 3D . . . . . . . . . . . . 21 4.2 Creaci´on y pre-filtrado de EPIs . . . . . . . . . . . . . . . . . . . . . . . . 23 4.3 Optimizaci´on del n´umero de vistas . . . . . . . . . . . . . . . . . . . . . . 25 5 Resultados 29 5.1 AlgoritmoFastDCT .............................. 29 ix simulando la escena 3D. Como las hendiduras son grabadas horizontalmente, siguiendo el esquema HPO, el vidente ve siempre la misma perspectiva vertical independientemente de la posici´on vertical desde donde se observe. De hecho, la imagen se mover´a como si estuviera en el plano de proyecci´on vertical. De acuerdo a la grabaci´on hologr´afica presentada en [2], se utilizan dos haces de luz coherente. El haz de l´aser es dividido en dos, originando un haz de referencia y uno del objeto. Este haz del objeto es ampliado cuando pasa a trav´es de la c´amara y las lentes de proyecci´on. Ilumina el holograma proyectando en la pantalla de proyecci´on una imagen grabada en film de cine en la c´amara. Los fotogramas de proyecci´on representan un objeto 2D localizado en un plano de proyecci´on sin l´ımites, y al mismo tiempo la placa con el mecanismo de hendiduras est´an encarados uno al otro, como se puede ver en la Fig. 1.4 con m´as precisi´on. La placa hologr´afica esta expuesta a la imagen proyectada y el haz de referencia. Sin embargo, s´olo una hendidura del holograma est´a expuesta, con un pieza opaca perforada con agujeros cubriendo la placa del holograma. Figura 1.4: Disposibi´on de la mesa hologr´afica vista desde arriba (Fuente [2]). Cada hendidura tiene que ser iluminada con el haz de referencia en la misma direcci´on, independientemente de la posici´on lateral de la hendidura. Esto se consigue juntando el haz con una lente de colimaci´on o Fresnel, el mismo principio usado en los faros. Fig. 1.5 da una idea m´as detallada de c´omo un colimador funciona. (a) (b) (c) Figura 1.5: a) Lente Fresnel usada en faros. b) Lente plano-convexa equivalente. c) Esquema de un colimador que alinea los haces produciendo una salida paralela. Otra configuraci´on t´ıpica usada en grabaci´on hologr´afica consiste en lanzar dos haces 4 de l´aser a un objeto f´ısico y grabar la interferencias de longitud de onda que el objeto desprende junto con el haz de referencia. En la Fig. 1.6 se detalla el esquema. En vez de lentes de colimaci´on, se usan lentes y espejos para redirigir los haces hacia la placa hologr´afica. La grabaci´on es realizada en la oscuridad para evitar cualquier tipo de interferencia con la luz natural. (a) (b) Figura 1.6: a) Esquema de grabaci´on hologr´afico con un objeto expuesto a un haz de l´aseres. Bob Mellish. b) Fotograf´ıa del mismo escenario en un laboratorio. College of Optics & Photonics, University of Central Florida. Cuando una placa hologr´afica es iluminada por un haz de referencia equivalente, se puede ver una reconstrucci´on del objeto. Si la posici´on de la placa o el observador cambia, la apariencia del objeto cambia tambi´en como si fuera un objeto real. Un ejemplo de reconstrucci´on hologr´afica se muestra en la Fig. 1.7. Figura 1.7: El perro es reconstru´ıdo aplicando un haz l´aser de referencia al holograma. Holograma presentado por el Dr.Kenji Yamamoto en la 3D Media Training School, Tampere, 2012. 5 1.2 Alcance y objetivo del proyecto El objetivo de este proyecto es entender el proceso de creaci´on de hologramas para desarrollar nuevas t´ecnicas. Se desarolla el software necesario con el prop´osito de proveer la informaci´on necesaria para los siguientes etapas del proceso. Los hitos a alcanzar en este proyecto incluyen: •Dise˜no y desarrollo de un script que provee un set de im´agenes renderizadas correspondientes a un sistema horizontal-parallax only desde una escena 3D. •Creaci´on de im´agenes epipolares (EPI) a partir de las im´agenes renderizadas, permitiendo representar todas las vistas en una imagen. •An´alisis de las im´agenes EPI en el dominio de la frecuencia. Relaci´on entre caracter´ısticas de la escena y su representaci´on en el dominio de Fourier. •Dise˜no y desarrollo de un filtro para procesas las im´agenes EPI, optimizando el n´umero m´ınimo de c´amaras. 1.3 Estructura Este documento est´a organizado con la siguiente estructura: El cap´ıtulo 2 describe el estado del arte relacionado con el tema. El cap´ıtulo 3 presenta el an´alisis te´orico que es combinado con los experimentos. El cap´ıtulo 4 explica todas las herramientas desarrolladas para llevar a cabo los experimentos y justificar su prop´osito. El cap´ıtulo 5 muestra los resultados obtenidos y su representaci´on gr´afica. El cap´ıtulo 6 completa la investigaci´on proveyendo las conclusiones y el trabajo futuro que puede proseguir el presente estudio. El cap´ıtulo 7 contiene las referencias bibliogr´aficas usadas para el desarrollo de este estudio. Los cap´ıtulos siguientes contienen los ap´endices incluyendo el resto del desarrollo completado durante este trabajo. La tesis desarrollada originalmente en ingl´es conformar´a el primer anexo. 6 2. Conceptos te´oricos Este cap´ıtulo revisa el estado del arte de trabajos anteriores as´ı como todo el conocimiento necesario para entender la investigaci´on realizada. Tambi´en revisa la propuesta de impresi´on de estereogramas hologr´aficos. 2.1 M´etodo de holograma generado por ordenador Como se ha explicado en el cap´ıtulo anterior, un holograma contiene informaci´on del campo de luz. Un holograma estereogr´afico es una versi´on discretizada compuesta de varios elementos de holograma llamados hogel (hologram element). Cada hogel contiene informaci´on 3D desde varias perspectivas, la cual ha sido previamente grabada dependiendo del ´angulo de refracci´on (i.e: el punto de vista del observador). Un holograma completamente computado ofrece un infinito n´umero de perspectivas. Por otro lado, un estereograma hologr´afico ofrece s´olo un n´umero finito de perspectivas, usando una reconstrucci´on por longitud de onda. Los estereogramas hologr´aficos forman una aproximaci´on de la longitud de onda decrementando el tiempo de computaci´on. Esta aproximaci´on es hecha con un n´umero discreto de perspectivas. La propuesta desarrollada en este estudio persigue la optimizaci´on del n´umero de vistas necesarias [6], obteniendo una percepci´on hologr´afica similar. Un m´etodo de codificaci´on relacionado con la difracci´on [8] hace referencia al estereograma hologr´afico como la suma de las amplitudes moduladas de la se˜nal chirp. Una se˜nal chirp incrementa su frecuencia con el tiempo. Su amplitud modulada en el plano del holograma produce un set de emisores direccionales en el plano emisor. Cada chirp enfoca la luz para crear un punto emisor, mientras la luz de las vistas dependiente del ´angulo es codificada en su modulaci´on de amplitud. De acuerdo a Quinn [8], la computaci´on paralela vectorial de un holograma est´a dividida en tres pasos: 1. Pre-computaci´on de los vectores chirp en texturas chirp. 2. Rendering multi-vista de la escena directamente en los vectores de modulaci´on alma- 7 cenados en la textura de modulaci´on usando una c´amara de doble cono (capturando vista delantera y trasera). 3. Montaje del holograma juntando los vectores de modulaci´on y chirp mediante la obtenci´on de la textura y ejecutando un producto escalar. Este trabajo se enfocar´a en una porci´on que es desarrollada durante el segundo paso, con algunas variaciones, e.g: horizontal parallax-only en vez de usar c´amaras con doble captura. Las diferencias exactas se discutir´an en Secci´on 2.3 al final de este cap´ıtulo, donde se presentan los requisitos del m´etodo desarrollado. Figura 2.1: A la izquierda, el holograma y una imagen aumentada de un hogel. A la derecha, un p´ıxel almacenado en el sensor de la c´amara, representando el mismo p´ıxel de la escena, es grabado en el hogel. Los sensores de la c´amara (en color rojo) est´an situados en realidad en la c´amara. Un m´etodo m´as general de entender c´omo un sistema de rendering multi-vista funciona esta basado en los hogels introducidos previamente. Fig. 2.1 mostrar´a inteligiblemente la siguiente explicaci´on. Cada hogel contiene la onda de luz correspondiente a un p´ıxel en la escena, desde cada punto de vista. Esta informaci´on puede ser reconstru´ıda aplicando la luz apropiada al hogel, que emitir´a la informaci´on grabada dependiendo del punto de vista. Varias c´amaras graban desde diferentes posiciones a lo de la pista de c´amaras. Cada sensor est´a compuesto por muchos p´ıxeles. Un solo hogel contiene la informaci´on asociada al mismo p´ıxel en cada vista. Esto significa que contiene informaci´on de tantas direcciones diferentes como n´umero de vistas usadas en la escena. 2.2 Representaci´on de EPI La correlaci´on entre im´agenes desde diferentes posiciones en la misma escena est´atica connota la coherencia de perspectiva. La principal distinci´on entre la coherencia de perspectiva y coherencia temporal es que esta ´ultima es que la primera es m´as restrictiva 8 que la segunda, en relaci´on a los cambios en la imagen a lo largo del tiempo. Las caracter´ısticas de la escena, i.e: cambios en el sombreado y geometr´ıa, pueden ser analizados en el dominio epipolar. Una Epipolar Plane Image (EPI) es una representaci´on 2D del campo de luz. Constituye la transformaci´on de un set de im´agenes multi-vista en un corte multi-perspectiva del campo de luz. Para obtener esta representaci´on, primero las c´amaras tienen que ser dispuestas como se muestra en Fig. 2.2. Halle [5] describe estas dos diferentes disposiciones como c´amaras regular shearing (RS): planar regular shearing (PSR) y linear regular shearing (LSR). Figura 2.2: A la izquierda, la disposici´on LSR provee horizontal parallax-only mientras que a la derecha, PSR, captura informaci´on horizontal y vertical. La disposici´on LSR, la cual consiste de un array de c´amaras representa una configuraci´on HPO, proveyendo toda la informaci´on necesaria con una sola fila de c´amaras, pero removiendo el parallax vertical. El vidente recibir´a la misma imagen independientemente de posici´on vertical que tome. Por otra parte, la disposici´on PSR provee full parallax, pero requiere muchas m´as c´amaras. Comprender c´omo se forman las EPIs ayuda a entender el sistema de impresi´on hologr´afico, dado que las EPIs conforman un elemento esencial en el proceso. Imagina que todas las im´agenes renderizadas del sistema multi-vista son apiladas como una baraja de cartas. La imagen correspondiente a la vista m´as a la derecha se coloca al frente, y la m´as situada a la izquierda, atr´as. A A continuaci´on, el volumen espacio-perspectiva creado se filetea horizontalmente. El corte multi-perspectiva obtenido, i.e: la EPI, contiene la misma fila equivalente de cada vista y est´a tumbada, mirando hacia arriba. En otras palabras, la ´ultima fila de la EPI pertenece a la vista m´as a la derecha y la primera fila pertenece a la de m´as a la izquierda. Fig. 2.3 muestra gr´aficamente la t´ecnica usada para obtener la EPI. Aunque Fig. 2.3 muestra una sola EPI, todas ellas deber´ıan ser calculadas para construir el volumen espacio-perspectiva completo. Existen tantas EPIs como la resoluci´on vertical de las vistas. La conveniencia de usar EPIs es fuertemente notable en relaci´on con las posibilidades de usar algunas caracter´ısticas en futuros desarrollos. Por ejemplo, su estructura linear 9 Figura 2.3: A la izquierda, 5 vistas. En el medio, ambas perspectivas 3D muestran c´omo las vistas son apiladas. El marco rojo designa la fila elegida para esta EPI, que reposa horizontalmente. A la derecha, la EPI resultante contiene la vista m´as a la izquierda arriba, y la vista m´as a la derecha abajo. N´otese c´omo el ojo derecho del camale´on desaparece de la EPI seg´un la c´amara se mueve a la derecha. permite aplicar efectivamente algoritmos de interpolaci´on. Esto es totalmente importante porque uno de los objetivos de este proyecto (explicado en m´as detalle en el cap´ıtulo siguiente) responde a la optimizaci´on del n´umero de vistas. Adicionalmente, las EPIs re´unen caracter´ısticas extra, como detecci´on de caras, sombreado y oclusi´on de objetos. Este ´ultimo puede ser observado en Fig. 2.4 donde se ha seleccionado una fila diferente. Figura 2.4: A la izquierda, el camale´on 3D y la fila elegida marcada en rojo. A la derecha, la EPI resultante donde las oclusiones producidas por ambas piernas est´an encerradas en un c´ırculo rojo. 2.3 Requisitos de la propuesta Despu´es de presentar los conceptos te´oricos necesarios, este cap´ıtulo describir´a los requisitos adoptados por la propuesta realizada y las razones de por qu´e ciertos par´ametros son elegidos y por qu´e otros son descartados. Antes de todo, como el sistema de impresi´on hologr´afico est´a basado en im´agenes multi-vista, se tiene que elegir una configuraci´on determinada con referencia a las c´amaras: 10 especifica su posici´on y rotaci´on con relaci´on a la escena. Hay disponibles cuatro configuraciones: c´amaras apuntando al frente paralelamente al plano de proyecci´on y tres tipos de c´amaras recentradas; dos de ellas est´an autom´aticamente descartadas por las distorsiones que introducen. Fig. 2.5 muestra de una manera m´as inteligible la disposici´on de las c´amaras en estas cuatro configuraciones. El tama˜no del plano de proyecci´on y su distancia a la linea de c´amaras determinar´a el ´angulo de visi´on y distancia focal especificada en la lente de la c´amara. Adem´as, el plano de proyecci´on conforma la secci´on cruzada con la pir´amide proyectada desde la c´amara, haciendo capaz a la c´amara de capturar apropiadamente el ´area correspondiente. Figura 2.5: A la izquierda, c´amaras paralelas. En el medio, c´amaras convergentes enfocan al centro de la escena a lo largo de la pista y giran alrededor del centro de la escena. Estas dos configuraciones conllevan a distorsiones. A la derecha, c´amaras recentradas con los sensores desplazados lateralmente de acuerdo con la lente para mantener la imagen centrada. Por lo tanto, la segunda y la tercera configuraci´on deben ser descartadas. Violan las correspondencias en t´erminos de geometr´ıa entre la posici´on de las c´amaras y la posici´on del observador. El lugar en el ´area donde las im´agenes son tomadas debe corresponder con el punto de vista del estereograma. La pista a lo largo de la cual se mueven las c´amaras, debe ser recta, manteniendo la misma geometr´ıa que la placa hologr´afica plana. La c´amara toma una vista desde cada posici´on correspondiente a la posici´on de las hendiduras. Adem´as, cada punto debe permanecer en el mismo lugar en la escena independientemente desde d´onde es capturado el punto de vista. De esta manera, la profundidad se mantiene en cada vista. Escenas complejas con puntos muy cercanos y distantes deber´ıan conservar la misma profundidad para cada punto de vista. Estas dos condiciones hacen que las dos configuraciones mencionadas arriban sean incompatibles. Adem´as, la EPI obtenida, por ejemplo, con la segunda configuraci´on, carece de propiedades que la hagan ´util, como l´ıneas rectas. Con respecto a las c´amaras paralelas o recentradas, las primeras tienen una desventaja. La mayor´ıa de los puntos de vista contienen informaci´on que no puede ser usada. Por lo tanto, es informaci´on in´util que es procesada pero nunca grabada en el holograma final. Las c´amaras paralelas conllevan una configuraci´on m´as sencilla, pero su desventaja implica demasiada informaci´on fuera del plano de proyecci´on. La cuarta configuraci´on mostrada en Fig. 2.5 es la elegida. Una manera de conseguir esas im´agenes renderizadas, como es explicado m´as en profundidad en el siguiente cap´ıtulo, consiste en usar la primera configuraci´on y c´amaras con un campo de visi´on que abarque el plano de proyecci´on incluso 11 (a) (b) (c) Figura 2.6: a) Las c´amaras paralelas preservan las correspondencias geom´etricas, pero la mayor parte, en negro, suponen informaci´on in´util. b) La EPI presenta geometr´ıa deformada causada por una configuraci´on err´onea. c) Las c´amaras recentradas proveen informaci´on ´util independientemente del punto de vista y preservan las correspondencias geom´etricas. desde el punto de vista m´as extremo (lateralmente hablando). Fig. 2.6 presenta las EPIs finales que ser´ıan producidas por las diversas configuraciones de c´amara. Los par´ametros elegidos para el siguiente desarrollo son los siguientes. La resoluci´on de la c´amara es de 600 p´ıxeles y 3 capas dado el espacio de color RGB. El n´umero de vistas elegidas es 512, lo que simplificar´a post-procesados de los renders al trabajar con potencias de 2. Renderizar 512 vistas de una escena compleja puede resultar en un proceso de computaci´on largo y pesado. Optimizar el n´umero de c´amaras a trav´es de interpolaci´on ser´a una de las tareas a ser pulidas. 12 3. An´alisis Como ha sido expuesto en cap´ıtulos anteriores, optimizar el n´umero de c´amaras implica el an´alisis de interpolaci´on y re-muestreado. Con este prop´osito, el an´alisis de EPIs tiene que hacerse en el dominio espacial y de frecuencia para entender como manipular las im´agenes. 3.1 Relaciones entre las caracter´ısticas de la escena y la representaci´on en el dominio de Fourier Antes de afrontar el estudio del filtrado e interpolaci´on de im´agenes epipolares, es apropiado tener un primer contacto con las propiedades de la escena y su representaci´on en el dominio de Fourier, i.e: la profundidad. La profundidad juega un papel muy importante en estereogramas hologr´aficos. Como diferentes profundidades afectan la representaci´on de objetos en el dominio epipolar y de frecuencia se explicar´a en esta secci´on. Inicialmente, se considera un punto z(v,t) en una escena en particular (Fig. 3.1). El punto es visto por dos c´amaras, c1 yc2 con una distancia focal f, las cuales se mueven a lo largo de la pista de c´amara t, grabando la posici´on horizontal donde caer´ıan los dos puntos (el mismo punto visto desde dos posiciones diferentes). Estas dos posiciones donde el punto z(v,t) es visto desde las c´amaras son vpara c1 yv’ para c2. La funci´on de disparidad v−v0=ft/z describe la geometr´ıa de la escena. Fig. 3.2 muestra el modelo de una se˜nal con profundidad constante. Las l´ıneas rojas representan la profundidad m´ınima y m´axima. Tambi´en se observa el soporte espectral del campo de luz. La l´ınea azul, con un ´angulo de 45ºrespecto de los ejes representa la profundidad constante. Las l´ıneas rojas determinan los l´ımites de de profundidad del soporte espectral. La inclinaci´on considerada en este dominio puede ser traducida como la profundidad de los objetos en la escena. De acuerdo a Chai[6], ”Cualquier escena con una profundidad entre zmin yzmax tendr´a su soporte espectral continuo limitado en el dominio de frecuencia”. El uso de im´agenes epipolares conecta las caracter´ısticas de la escena, como la profundidad, con el dominio de Fourier, el cual puede ser adecuadamente manipulado. Finalmente, un par de ejemplo 13 4. Dise˜no del sistema multi-vista Este cap´ıtulo examina las t´ecnicas usadas para obtener los elementos necesarios para construir un sistema de holograma generado por ordenador. Las siguientes secciones explicar´an las herramientas desarrolladas para obtener apropiadamente las vistas de una escena virtual 3D, el post-proceso ejecutado a ellas y las t´ecnicas de re-muestreo que consiguen la deseada optimizaci´on en cuanto al n´umero de vistas. Figura 4.1: 1) Rendering y recortado de las im´agenes para obtener la parte ´util de ellas. 2) Creaci´on de EPIs de densidad completa y pre-filtradas para obtener EPIs anti-aliased y re-muestradas. 3) Interpolaci´on vertical de las EPIs para recuperar su resoluci´on original. 4) Reconstrucci´on de las vistas a partir de las EPIs interpoladas. A trav´es de las fases de desarrollo y testado, se han creado m´as de 32.000 ficheros y se han manipulado casi 3GB de informaci´on, lo que requiere una organizaci´on de directorios clara y estructurada: /. functions.blend/py/m (1) scene_name/epis_density (2) scene_name/epis_density_interpolation (3) scene_name/parallel_r_density (4) scene_name/parallel_rc_density (5) scene_name/parallel_density_interpolation (6) 1) El directorio representa todas las funciones scripts de Blender, Python y Matlab. 2) El directorio contiene las EPIs de densidad completa y pre-filtradas con 5 factores de interpolaci´on. 3) Cada directorio contiene las EPIs interpoladas con 5 factores de interpolaci´on. 4) El directorio contiene los renders originales. 5) El directorio contiene las vistas recortadas a la resoluci´on adecuada. 6) El directorio contiene las vistas reconstruidas para cada factor de interpolaci´on. En el ´arbol de directorios presentado, densidad es representado por un n´umero (512 en este caso) e interpolaci´on por uno de los 5 kernels de interpolaci´on usado. 20 4.1 Construyendo el sistema multi-vista para escenas 3D El primer paso antes de analizar las vistas de la escena implica obtener dichas vistas. El software utilizado es Blender como escena virtual 3D y configuraci´on de c´amaras. El script que coloca las c´amaras en su posici´on correcta est´a escrito en Python, e incluye una interfaz gr´afica (Fig. 4.2. (a) (b) Figura 4.2: a) Interfaz para el modo paralelo. b) Interfaz para el modo convergente. Inicialmente, el script desarrollado ofrece dos posibles configuraciones, incluyendo c´amaras paralelas y convergentes (Fig. 4.3. En cap´ıtulos anteriores ya se ha comentado la raz´on por la que las c´amaras convergentes no proveer´an las propiedades geom´etricas apropiadas. Sin embargo, se ha desarrollado este modo con el objetivo de obtener pruebas gr´aficas de la distorsi´on producida y estudiar sus consecuencias en el dominio epipolar. (a) (b) Figura 4.3: a) Posicionamiento de las c´amaras en el modo paralelo. b) Posicionamiento de las c´amaras en el modo convergente. Los atributos que definen la escena est´an relacionados con los siguientes par´ametros: •Mode: Cambia entre el modo paralelo y convergente. •Width Field: Longitud total de la pista de c´amara (Unidades Blender * 1000). 21 •Field Density: N´umero total de c´amaras. •Distance to scene (Disponible s´olo en el modo convergente): Distancia al centro de la escena donde las c´amaras enfocar´an. De acuerdo al modo paralelo, las c´amaras se mueven a lo largo de la pista de c´amara, capturando vistas paralelas a la escena. La c´amara en el centro de la pista tiene el campo de visi´on necesario para abarcar la escena entera. Sin embargo, cuando la c´amara se mueve de un lado a otro, la vista capturada contiene informaci´on no necesaria. Fig. 4.4 muestra c´omo informaci´on vac´ıa es grabada en vez de p´ıxeles pertenecientes a la escena. Figura 4.4: La c´amara n´umero 2, en el extremo de la pista, captura informaci´on no necesaria (en rojo). Una manera simple de resolver este problema es ensanchar el campo de visi´on de la c´amara, adquiriendo la porci´on entera del plano de proyecci´on desde cualquier punto de vista. Esta soluci´on requiere m´as informaci´on que renderizar, aunque las regiones adicionales renderizadas suelen contener pocos o ning´un objeto, lo que no implica apenas extra tiempo de rendering. Una vez todas las vistas han sido renderizadas, es el momento de post-procesarlas y convertirlas en vistas con la resoluci´on adecuada. Un script escrito en Matlab (desde ahora en adelante, todos los scripts usados para procesar im´agenes estar´an basados en Matlab) recorta la imagen manteniendo solo la regi´on correspondiente con el plano de proyecci´on. Las im´agenes procesadas tendr´an la misma resoluci´on despu´es de este recortado. Fig. 4.5 muestra c´omo el siguiente script (crop renders.m) transformar´a los renders iniciales en im´agenes adecuadas para procedimientos posteriores. function crop_renders(object,density,h_res) % Object: Name of the scene % Density: Number of views -> Vertical resolution of EPIs % h_res: Horizontal resolution of views % It crops renders to normal resolution ... for i=1:density img = imread(strcat(name_in,int2str(density-i+1)),’png’); cropped = img(:,(i-1)*2+1:(i-1)*2+h_res,:); imwrite(cropped,strcat(name_out,int2str(density-i+1),’.png’)); end 22 Figura 4.5: La c´amara se caracteriza por un campo de visi´on m´as ancho. El ´area innecesaria de la imagen es descartada en la imagen de la derecha. La escena trabajada en este trabajo contiene 512 vistas, y los scripts de rendering producen im´agenes de 1622x600 p´ıxeles. Esas im´agenes est´an recortadas a una resoluci´on final de 600x600 p´ıxeles. La siguiente Fig. (4.6) contiene algunos ejemplos producidos. Figura 4.6: De izquierda a derecha, muestras de las vistas n´umero 1, 64, 128, 192, 256, 320, 384, 448 y 512. 4.2 Creaci´on y pre-filtrado de EPIs Todas las vistas obtenidas en la secci´on anterior componen la representaci´on espacial. Para conseguir una representaci´on del campo de luz de la escena, las vistas tienen que ser transformadas en im´agenes epipolares. Los siguientes pasos incluyen la creaci´on de las ya mencionadas im´agenes epipolares y el pre-filtrado de esas im´agenes para remover el aliasing. El primer script se encarga de crear EPIs de densidad completa, las cuales tienen el mismo n´umero de columnas que n´umero de vistas renderizadas de la escena (512 en este caso). Un total de 600 diferentes EPIs ser´an creadas, dado que las vistas tienen 600 p´ıxeles de resoluci´on vertical. Esto es completado por el siguiente script, el cual toma la misma fila de cada vista para crear una ´unica EPI. Por ejemplo, para crear la primera EPI, asigna la primera fila de la primera vista a la primera fila de la EPI, la primera fila de la segunda vista a la segunda fila de la EPI, etc. function epi(object,density,v_res) % Object: Name of the scene % Density: Number of views -> Vertical resolution of EPIs % v_res: Vertical resolution of views -> Number of EPIs % It creates epipolar images with full density 23 Como resultado de este script en esta particular escena, las primeras y ´ultimas 30-40 filas son completamente negras, dado que ning´un objeto aparece al comienzo o final de las vistas. Fig. 4.7 muestra algunas muestras obtenidas por esta funci´on. Figura 4.7: Desde la izquierda a la derecha, muestras de las EPIs n´umero 1, 100, 200, 300, 400, 500 y 600. Como ha sido explicado, el filtro anti-aliasing es requerido antes de la decimaci´on. La funci´on presentada m´as abajo (una versi´on simplificada) aplica un filtro anti-aliasing. Funciona con 5 factores de interpolaci´on diferentes [2, 4, 8, 16, 32] todos potencia de 2, lo que facilita el procesado de transformadas de Fourier. El re-muestreado produce una salida mostrada en Fig. 4.8. function prefiltering(object,density,v_res) % Object: Name of the scene % Density: Number of views -> Vertical resolution of EPIs % v_res: Vertical resolution of views -> Number of EPIs % It filters EPIs with anti-aliasing filter % and down-samples vertically in 5 factors [2, 4, 8, 16, 32] %Preparation of extended signal epi = im2double(...); epi_ext(1+factor:end-factor,:,:) = epi; epi_ext(1:factor,:,:) = epi(1,:,:); epi_ext(end-factor+1:end,:,:) = epi(end,:,:); %Fourier transform, filtering and normalization epi_f = fft(epi_ext); normal = density/((density/factor) + 2); fraction = density/(2*factor); epi_fdec = [epi_f(1:fraction+2,:,:); epi_f(end-fraction+1:end,:,:)] / normal; %Inverse Fourier transform and cropping extended samples. epi_factor = ifft(epi_fdec); epi_final = epi_factor(2:end-1,:,:); 24 Figura 4.8: La EPI original y las 5 EPIs anti-aliased con sus 5 factores de interpolaci´on diferentes. 4.3 Optimizaci´on del n´umero de vistas En este punto del dise˜no, los siguientes desarrollos trabajar´an en las vistas convertidas en EPIs filtradas, para convertirlas de nuevo a su resoluci´on original. El re-muestreado que optimiza el n´umero de vistas para construir el sistema multi-vista. Los scripts responsables de esta fase son presentados a continuaci´on: function run_epi_interp(object) % Object: Name of the scene % It runs epi_interp() with 5 factors [2, 4, 8, 16, 32] and % 5 interpolation kernels [triangle, cubic, lanczos2, lanczos3, fastdct] function epi_interp(object,density,factor,v_res,interpolation) % Object: Name of the scene % Density: Number of views -> Vertical resolution of EPIs % Factor: Interpolation magnitude % v_res: Vertical resolution of views -> Number of EPIs % interpolation: Interpolation kernel used % It up-samples vertically EPIs La funci´on epi interp utiliza la funci´on de Matlab imresize() con los kernels de interpolaci´on bilinear, bicubic, lanczos2 and lanczos3, o la siguiente funci´on fdctEPI, la cual implementa el algoritmo de escalado fast DCT [10]. Este algoritmo ha sido implementado para escalar se˜nales 1D, y m´as tarde adaptado a dos dimensiones. Todos los scripts codificados que mejoran el c´odigo en cuanto a tiempo de computaci´on son presentados en los ap´endices. El ejemplo siguiente (versi´on simplificada) tiene como objetivo demostrar y presentar los fundamentos b´asicos de la interpolaci´on fast DCT. function img_out = fdctEPI(img_in,factor) % img_in: input image 25 % factor: interpolation magnitud % It up-samples vertically img_in using fastDCT interpolation No = rows_in; n = repmat(0:rows_in-1,rows_out,1); k = repmat((0:rows_out-1)’,1,rows_in); DCT = sincd(2*No-2, 2*rows_in, (k+1/2)/factor -n -1/2) .* ... cos( pi/(2*rows_in) * ( (k+1/2)/factor -n -1/2) ) + ... sincd((2*No-2), 2*rows_in, (k+1/2)/factor +n +1/2) .* ... cos( pi/(2*rows_in) * ( (k+1/2)/factor +n +1/2) ); % Multiplication by DCT-Matrix to all image component layers for i=1:components img_out(:,:,i) = DCT*img_in(:,:,i); end % Truncating negative and greater than 1 values img_out(img_out<0) = 0; img_out(img_out>1) = 1; function result = sincd(M,N,x) % sincd(M,N,x) = sin(piMx/N)/(Nsin(pix/N)) % Digital sinc function used in fastDCT Una de las secciones cr´ıticas en el c´odigo de fdctEPI() es la poblaci´on de nyk, dado que representan las matrices que substituyen los bucles iniciales, que operan por filas y columnas. Ejecutar esta operaci´on con dos bucles anidados hace la computaci´on ineficiente. En una versi´on preliminar, uno de los bucles fue substituido por una operaci´on vectorial, pero popular dos matrices y efectuar operaciones lineales con ellas es m´as coherente. Matriz kynresponden al siguiente patr´on: k=      0 0 · · · 0 1 1 · · · 1 . . .. . ..... . . n n · · · n       n=      0 1 · · · n 0 1 · · · n . . .. . ..... . . 0 1 · · · n       Tres m´etodos de poblaci´on de matrices fueron considerados. % Kronecker tensor product ns = kron([1:size_in]-1,ones(1,size_out)’); ks = kron(ones(1,size_in),[1:size_out]’-1); 26 % Replicating matrices ns = repmat([0:size_in-1],size_out,1); ks = repmat([0:size_out-1]’,1,size_in); % Manual with ones temp= [0:size_in-1]; ns= temp(ones(1,size_out),:); temp= [0:size_out-1]’; ks= temp(:,ones(1,size_in)); Los resultados para unas resoluciones de 300, 600 and 900 p´ıxeles y factores 2-10 se muestran en la Fig. 4.9. Figura 4.9: El m´etodo kron es notablemente peor que los otros dos. No hay grandes diferencias entre el manual y ones, pero al mismo tiempo, manual requiere menos memoria y n´umero de instrucciones. La reconstrucci´on de vistas desde las EPIs interpoladas usando los kernels explicados arriba es directa. Una vista es formaba apilando filas que provienen de la misma fila de cada EPI. Las siguientes funciones son las responsables de la fase final. function run_reconstruct_views(object) % object: Name of the scene % It runs reconstruct_views() with 5 factors [2, 4, 8, 16, 32] and % 5 interpolation kernels [triangle, cubic, lanczos2, lanczos3, fastdct] function reconstruct_views(object,density,factor,v_res,interpolation) % Object: Name of the scene % Density: Number of views -> Vertical resolution of EPIs 27 % Factor: Factor of interpolation applied vertically % v_res: Vertical resolution of views -> Number of EPIs % interpolation: Method/kernel used for re-sampling % It reconstructs all the views from the re-sampled EPIs La Fig. 4.10 presenta el proceso completo visto desde una perspectiva gr´afica. Figura 4.10: La imagen muestra todos los pasos seguidos en este trabajo, y c´omo las im´agenes son transformadas y manipuladas en cada fase. 28 5. Resultados En este cap´ıtulo se presentan e interpretan todo los resultados gr´aficos y num´ericos. Los tests principales ejecutados en este trabajo incluyen el algoritmo fast DCT [10] y otros kernels de interpolaci´on, as´ı como la comparaci´on final entre las vistas originales e interpoladas. 5.1 Algoritmo Fast DCT El algoritmo Fast DCT ha sido usado para re-muestrear verticalmente Epipolar Planar Images. Previo a los tests de im´agenes, un set de se˜nales 1D fueron testadas para un mejor entendimiento de los efectos producidos en ciertos casos. Fig. (5.1) muestra alguno de los efectos que el escalado produce en dichas se˜nales. Consecuencia de este efecto, los bordes de las im´agenes epipolares que contienen b´asicamente lineas y franjas, presentan desenfoque y artefactos. Otro ejemplo, testado para confrontar los resultados de [10] presenta iterativamente el escalado de una imagen de texto. Los resultados arrojan la impresi´on de que independientemente del n´umero de iteraciones ejecutadas, el texto permanece inteligible, a diferencia de lo que ocurre con otros m´etodos de interpolaci´on. Fig. 5.2 muestra claramente la diferencia entre dichos m´etodos. As it was mentioned above, the border artefacts present will have to be tackled down in posterior steps. The artefacts appearing beyond the borders will take its utility away. Figure 5.3 depicts why the displayed results are not useful for epipolar image interpolation. Como ha sido mencionado previamente, los artefactos presentes en los bordes tendr´an que ser abordados en pasos posteriores. Fig. 5.3 explica por qu´e los resultados mostrados no son ´utiles para interpolaci´on epipolar. Otro aspecto considerado en este texto es el tiempo de computaci´on empleado por cada uno de los m´etodos. Los resultados presentados en Fig. 5.4 muestran que la implementaci´on de Fast DCT realizada para Matlab es m´as lenta que las funciones built-in en Matlab, dado que Fast DCT comprende dos funciones externas. El hardware usado en la ejecuci´on de estos tests no influencia los resultados, dado que no se est´an comparando mediciones con informaci´on en otro hardware. Adem´as, cada test fue ejecutado 10 veces 29 6.2 Trabajo futuro Este proyecto estudia un intento inicial de establecer una metodolog´ıa para un rendering multi-vista. Usa c´amaras virtuales debido a la libertad de manipulaci´on de sus par´ametros. Tambi´en toma ventaja del entorno 3D para renderizar objetos y escenas que pueden ser adaptados de mejor manera a los requisitos establecidos, a diferencia de los objetos en el mundo real. El siguiente paso a seguir implicar´ıa el uso de c´amaras f´ısicas para adquirir las vistas, nuevos m´etodos de interpolaci´on incluyendo t´ecnicas mixtas. Varias propuestas avanzadas con respecto la configuraci´on de c´amaras, siendo estas re-centradas m´as complejas que lleven a una configuraci´on full-parallax. Tambi´en, par´ametros intr´ınsecos del las lentes en las c´amaras reales jugar´an un papel importante en los temas relacionados con la profundidad. Llevar este estudio a posible aplicaciones es el siguiente paso a dar. Por ejemplo, juegos gr´aficos por ordenador con interacci´on en tiempo real requieren m´as y m´as simplificados modelos de pol´ıgonos, dado que toman la ventaja de mapping de texturas avanzados. Los modelos, m´as simples, pueden ser manipulados mejor en el dominio epipolar, y la calidad del producto final ser´a m´as dependiente a nivel de textura. El trabajo de este proyecto simboliza un n´umero independiente de fases para obtener un resultado final. Sin embargo, un estudio futuro podr´ıa desarrollar todos estos pasos en una ´unica unidad de producci´on. Las especificaciones de dicho software de producci´on variar´ıa considerablemente, dependiente de las diferentes plataformas, herramientas y software de desarrollo usado durante el proceso. Si una herramienta de producci´on tiene que ser desarrollada, algunos aspectos tendr´ıan que ser tenidos en cuenta. Por ejemplo, todos los algoritmos deber´ıan de ser escritos en lenguajes compilados, donde una verdadera optimizaci´on podr´ıa tener lugar. Mex o C podr´ıan ser un buen punto de partida como alternativa a Matlab. Algunas partes del c´odigo que requieren una inmensa transferencia de entrada/salida podr´ıa ser paralelizados. El aspecto m´as importante en la optimizaci´on de c´odigo reside en la Ley de Amdahl: la optimizaci´on de un bloque de c´odigo optimizar´a la aplicaci´on general en la misma proporci´on que ese c´odigo representa el proceso general. Algunas veces un gran esfuerzo es puesto en optimizar algo que, como ha ocurrido en el estudio de poblaci´on de matrices, no refleja un beneficio digno en comparaci´on con el tiempo invertido. 6.3 Valoraci´on personal Esta tesis ha significado un proyecto completamente especial en el per´ıodo que ha tenido lugar, por varias razones. Primero, afront´e este proyecto en un ambiente internacional, como parte de mi estancia Erasmus en Finlandia en la Tampere university of Technology. En segundo lugar, ser investigador en el departamento de Signal Processing me ha 36 mostrado una cara de la universidad en la que no hab´ıa estado involucrado antes. Tambi´en me ha dado la posibilidad de acudir a la 3D Media Training School en el verano de 2012 en Tampere. Finalmente, aunque mis conocimientos previos son enteramente basados en Computer Science, el reto de adaptar mi conocimiento y comprender algoritmos y propiedades relacionados con el procesado de se˜nal (la mayor´ıa de ellos no estudiados en mi programa con anterioridad) me ha impulsado a afrontar los problemas desde una perspectiva diferente. Adem´as, hacer uso de mi experiencia en software y computaci´on ha sido de ayuda a la hora de desarrollar ciertos algoritmos. Mi visi´on de c´omo las im´agenes son procesadas y c´omo esa tecnolog´ıa es usada ha cambiado despu´es de completar esta tesis. 37 38 Bibliography [1] E. N. Leith. White light holograms. Scientific American, (1976). [2] Michal W. Halle. The Generalized Holographic Stereogram. Massachusetts Institute of Technology, (1991). [3] Michal W. Halle, Adam B. Kropp. Fast computer graphics rendering for full parallax spatial displays. Brigham and Women’s Hospital, Massachusetts Institute of Technology, (1997). [4] Ravikanth Pappu et al.A Generalized Pipeline for Preview and Rendering of Synthetic Holograms. Massachusetts Institute of Technology, Interval Research Corporation, (1997). [5] Michal W. Halle. Multiple Viewpoint Rendering. Brigham and Women’s Hospital, (1998). [6] Jin-Xiang Chai et al.Plenoptic Sampling. Microsoft Research China, (2000). [7] Wendy Plesniak et al.Reconfigurable Image Projection Holograms. Brigham and Women’s Hospital, MIT Media Laboratory, ThingMagic, Inc., (2006). [8] Smithwick, Quinn Y. J. et al.Real-time Shader Rendering of Holographic Stereograms. Society of Photo-Optical Instrumentation Engineers, (2009). [9] Melania Paturzo et al.Holographic Display of synthetic 3D dynamic scene. 3D Research Center and Springer, (2010). [10] Leonid Yaroslavsky, Leonid Bilevich. Fast DCT-based Algorithm for Signal and Accurate Scaling . Tel Aviv University, (2012). 39 40 7. Ap´endices Las secciones siguiente contienen todo el material que no pod´ıa ser incluido en el texto principal debido a su extensi´on o relevancia. 7.1 Original MSc Thesis at Tampere University of Technology El primer anexo incluye el proyecto realizado originalmente en ingl´es, en la Tampere University of Technology, Tampere, Finlandia. La numeraci´on de dichas p´aginas no se ha cambiado para mantener una coherencia con la tabla de contenidos e ´ındice de figuras con las que el documento comienza. 41 Table of Contents 1 Introduction 46 1.1 Context and motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 1.2 Scope and objective of the project . . . . . . . . . . . . . . . . . . . . . . . 50 1.3 Structure .................................... 51 2 Background 52 2.1 Computer generated hologram method . . . . . . . . . . . . . . . . . . . . 52 2.2 EPIrepresentation ............................... 53 2.3 Requirements of the approach . . . . . . . . . . . . . . . . . . . . . . . . . 55 3 Analysis 58 3.1 Relations between scene characteristics and Fourier domain representation 58 3.2 Transform-domain analysis of EPI images . . . . . . . . . . . . . . . . . . 60 3.3 Filtering of EPI and optimizing the number of views . . . . . . . . . . . . 61 4 Design of the multi-view system 64 4.1 Building a multi-view system from 3D scenes . . . . . . . . . . . . . . . . . 65 4.2 Creation and pre-filtering of EPIs . . . . . . . . . . . . . . . . . . . . . . . 67 4.3 Optimizing number of views . . . . . . . . . . . . . . . . . . . . . . . . . . 68 5 Results 73 5.1 FastDCTalgorithm .............................. 73 5.2 Experimentalresults .............................. 75 i 6 Conclusions 78 6.1 Conclusions about the results . . . . . . . . . . . . . . . . . . . . . . . . . 78 6.2 FutureWork................................... 80 6.3 PersonalAssessment .............................. 80 ii List of Figures 1.1 HolograminStarWars............................. 47 1.2 Horizontal Parallax Only vs. Full-parallax stereograms . . . . . . . . . . . 48 1.3 Cross-section of a camera with projection plane . . . . . . . . . . . . . . . 48 1.4 Holographic table layout . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 1.5 Collimatinglens................................. 49 1.6 Typical holographic recording scheme . . . . . . . . . . . . . . . . . . . . . 50 1.7 Example of holography reconstruction . . . . . . . . . . . . . . . . . . . . . 50 2.1 Representation of the hogel (hologram element) . . . . . . . . . . . . . . . 53 2.2 Camera arrangement for full or horizontal only parallax . . . . . . . . . . . 54 2.3 EPI construction process . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 2.4 Example of occlusion in EPI . . . . . . . . . . . . . . . . . . . . . . . . . . 55 2.5 Four possible configurations for camera arrangement . . . . . . . . . . . . . 56 2.6 EPIs obtained with different camera configurations . . . . . . . . . . . . . 56 3.1 Depth function represented in spatial and epipolar domain . . . . . . . . . 58 3.2 Spectral support of light field bounded in frequency domain. . . . . . . . . 59 3.3 EPI representation and 2D Fourier transform . . . . . . . . . . . . . . . . . 59 3.4 Anti-aliasing steps for a one-dimensional signal. . . . . . . . . . . . . . . . 61 3.5 Comparison between anti-aliasing filtered and non filtered EPIs. . . . . . . 62 4.1 Scheme of developing process . . . . . . . . . . . . . . . . . . . . . . . . . 64 iii 4.2 Interfaces for parallel and convergent scripts in Blender . . . . . . . . . . . 65 4.3 3D example of parallel and convergent cameras in Blender . . . . . . . . . 65 4.4 Cropping problem of views with parallel cameras . . . . . . . . . . . . . . 66 4.5 Solution for cropping renders. . . . . . . . . . . . . . . . . . . . . . . . . . 67 4.6 Samplesofviews................................. 67 4.7 SamplesofEPIs. ................................ 68 4.8 Samples of EPIs filtered with anti-aliasing. . . . . . . . . . . . . . . . . . . 69 4.9 Populating matrices with kron, repmat and ones. . . . . . . . . . . . . . . 71 4.10 Graphical perspective of the whole process. . . . . . . . . . . . . . . . . . . 72 5.1 Sample signals scaled with Fast DCT algorithm . . . . . . . . . . . . . . . 73 5.2 Comparison iterative scaling between bilinear, bicubic and fast DCT. . . . 74 5.3 Border effects in the DCT method . . . . . . . . . . . . . . . . . . . . . . . 74 5.4 Time comparison between fast DCT implementation and rest of interpolationkernels.................................... 75 5.5 PSNR of reconstructed views for every interpolation method. . . . . . . . . 76 5.6 Average of PSNR for every factor against aliased EPIs-based views. . . . . 77 6.1 Average of PSNR for every factor and range of acceptable PSNR values. . . 79 1 2. Background This chapter reviews the state of art in the area as well as all the necessary knowledge to understand the topic. It overviews the holographic stereogram printing approach. 2.1 Computer generated hologram method As explained in the previous chapter, a hologram contains light field information. A holographic stereogram is a discretised version composed of several hologram elements called hogels. Every hogel contains 3D information from various perspectives, which was previously recorded depending on the the angle of refraction (i.e: viewer’s point of view). A fully computed hologram offers an infinite number of perspectives. On the other hand, a holographic stereogram offers only a finite number of perspectives instead, using wavefront reconstruction. Holographic stereograms form an approximation of the wavefront, for decreased computation time. This approximation is done with a discrete number of perspectives. The approach developed in this thesis will pursuit the optimization of number the views needed [6], obtaining a similar holographic perception. One method of diffraction-specific encoding [8] refers to the holographic stereogram as a sum of overlapping amplitude modulated chirp gratings. A chirp signal increases its frequency over time. This amplitude modulated chirped grating on the hologram plane produces a set of view-directional emitters on an emitter plane. Each chirp focuses light to create a point emitter, while the angle-dependent brightnesses of the views are encoded in the amplitude modulation. According to Quinn [8], the parallel vector computation of a hologram is divided into three steps: 1. Pre-computing the chirp vectors into a chirp texture. 2. Multi-view rendering of the scene directly into the modulation vectors stored in the modulation texture using a double frustum camera. 3. Assembling the hologram by gathering chirp and modulation vectors via texture fetches and then performing a dot product. 52 This thesis will focus on a portion that is developed during the second, with some variations, e.g., horizontal parallax-only cameras instead of using double frustum. The exact differences will be discussed in Section 2.3 at the end of this chapter, where the requirements of the developed approach are presented. Figure 2.1: On the left, the hologram and an magnified picture of a single hogel. On the right, the pixel stored in the sensor, representing the same pixel of the scene, is recorded in the hogel. The camera sensors (in red color), are placed in the actual camera. A more general approach to understand how the multi-view rendering system works is based on the hogels introduced previously. Figure 2.1 will picture intelligibly the following explanation. Every hogel contains the light wavefront corresponding to one pixel in the scene, coming from every point of view. This information can be reconstructed applying the proper light to the hogel, then emitting the data recorded depending on the point of view. Several cameras record from different positions along the camera track line. Each sensor is composed of many pixels. A single hogel contains the information associate to the same pixel of every view. It means that it holds data from as many different directions as number of views used in the scene. 2.2 EPI representation The correlation between images from different positions of the same static scene connotes the perspective coherence. The main distinction between perspective and temporal coherence is that the former is more restricted than the latter, related to changes in images along time. Scene characteristics, i.e: shading and geometric changes, can be analysed in the epipolar domain. An Epipolar Plane Image (EPI) is a 2D representation of the light field. It constitutes the transformation of a set of multi-view images into a multi-perspective light field cut. In order to obtain this representation, first cameras should be arranged as shown in Figure 53 2.2 exhibits. Halle [5] describes these two different arrangements as regular shearing (RS) cameras: planar regular shearing (PSR) and linear regular shearing (LSR). Figure 2.2: On the left, the LSR disposition provides only horizontal parallax while the right dispotition, PSR, captures horizontal and vertical information. The LSR disposition, which consists of a single array of cameras represents a HPO setting, providing all the necessary information with only one row of cameras, but removing vertical parallax. Viewer will receive the same image as he sees from a different vertical position. On the other hand, the PSR disposition supplies full parallax, but, it requires many more cameras. Apprehending how EPIs are formed helps to understand the holographic printing system, since EPIs conform an essential element within the process. Imagine that all images rendered from the multi-view system are stacked as deck of cards. The image corresponding to the rightmost view is placed on the front, and the leftmost view on the back. Then, the created spatio-perspective volume is sliced horizontally. The multiperspective cut obtained, i.e: the EPI, holds the same equivalent row of every view and it is lying down, pointing at the top. In other words, the bottom row of the EPI belongs to the rightmost view and the top row belongs to the leftmost view. Figure 2.3 shows graphically the technique used to obtain the EPI. Although Figure 2.3 displays a single EPI, all of them should be calculated in order to build the entire spatio-perspective volume. There are as many EPIs as height resolution of the camera plane. The convenience of using EPIS is strongly notable with regard to possibilities of using some of the features in further developments. For instance, their linear structure allows to effectively apply interpolation algorithms. This is utterly important as one of the goals of this project (explained in more detail in the next chapter) responds to optimizing the number of views. Additionally, EPIs gather extra features, such as side detection, shading and occlusion of the objects. The last one can be observed in the Figure 2.4 where a different row has been chosen. 54 Figure 2.3: On the left, five views representing from the leftmost to the rightmost views, with some intermediate ones in between. In the middle, both 3D perspectives show the views stacked. The red frame designates the row chosen for this EPI, which reposes horizontally. On the right, the resultant EPI contains the leftmost views on the top and the rightmost views on the bottom. Notice how the right eye of the chameleon disappear from the EPI as the camera goes to the right. Figure 2.4: On the left, our 3D chameleon and the row chosen marked in red. On the right, the resultant EPI where occlusions produced by both pair of legs are enclosed in a red circle. 2.3 Requirements of the approach After presenting the necessary background, this chapter will describe the requirements adopted for the approach undertaken and the reasons of why certain parameters are chosen and why others are dismissed. First of all, as the holographic printing system is multi-view image-based, a determined setting with reference to cameras has to be chosen: it specifies their placement and position in relation to the scene. Four configurations are available: cameras pointing ahead parallel to the projection plane and three sorts of re-centering cameras; two of them are automatically discarded because of distortions introduced. Figure 2.5 displays in a more understandable manner the arrangement of the cameras in those four configurations. The size of the projection plane and its distance to the camera track, or slit, will determine the angle of view and focal length specified in the lens of the camera. Furthermore, the projection plane accommodates the cross-section with the pyramid projected from the 55 camera, thus making the camera able to image properly the corresponding area. Figure 2.5: On the left, parallel cameras. In the middle, convergent cameras heading the center of the scene along the track and spinning around the center. Theses two configurations lead to distortions. On the right, recentering cameras with their rear films shifted according to the lens to keep the image centered. Then, the second and third configurations must be discarded. They violate the correspondences in terms of geometry between the position of the cameras and the position of the viewer. The place in the area where images are taken must correspond with the stereogram point of view. The track where the camera moves along must to be straight, remaining with the same geometry as the holographic plane plate. The camera takes view from every position corresponding to the position of the slits. Besides, every point must remain in the same place in the scene regardless where the viewpoint is capturing from. This way, depth is maintained in every view. Complex scenes with too distant and too close points should remain at the same depth from each viewpoint. With a tilted film respecting the projection plane, its contained image would be scaled vertically, known as keystone effect. These two conditions make the two configurations mentioned above incompatible. Furthermore, the EPI obtained, for instance, with the second configuration, lacks of properties that make it useful such as linear tracks. (a) (b) (c) Figure 2.6: a) Parallel cameras preserve geometry correspondences, but most of the black region means useless data. b) EPI presents deformed geometry caused by wrong camera configuration. c) Recentered cameras provides useful data regardless the point of view and preserve geometry correspondences. Regarding the parallel or recentering cameras, the former have a disadvantage. Most of the viewpoints contain image data that could never be seen. Therefore, it is useless 56 information that is processed but never recorded in the final hologram. Parallel cameras brings an easier setting, but its disadvantage is that it implies too much data out of the projection plane. The fourth configuration showed in Figure 2.5 is considered. A way to achieve these rendered images, as it is explained in the chapter below, consists of using the first configuration and cameras with a wider field of view that accommodate the projection plane even from the most extreme points of views. Figure 2.6 presents the final EPIs that will be produced by the diverse camera configurations. The parameters chosen for later development are as follows. The resolution of the camera is 600 square pixels and 3 layers for RGB color space. The number of views chosen is 512, what will simplify some post-processing of the renders working with some power of 2. Rendering 512 views from a complex scene might result in a heavy and long computation process. Optimizing the number of cameras through interpolation will be one of the tasks to be polished. 57 3. Analysis As it was stated in previous chapters, optimizing the number of cameras involves the analysis of interpolation and re-sampling. With this purpose, the analysis of EPIs has to done in both spacial and frequency domain to understand how to manipulate the images. 3.1 Relations between scene characteristics and Fourier domain representation Before approaching the study of filtering and interpolation of epipolar images, it may be appropriate a first contact with a property of scenes and its Fourier domain representation, i.e: depth. Depth plays an important role in holographic stereograms. How different depths affect the representation of object in epipolar and frequency domain will be explained in this section. Initially, one point z(v,t) in a particular scene is considered (see Figure 3.1). It is seen by two cameras, c1 and c2 with focal length f, which go alone the camera track t, recording the horizontal position where the two points would fall onto. The two positions where the point z(v,t) fall onto the cameras are vfor c1 and v’ for c2. The disparity function v−v0=ft/z describes the geometry of the scene. Figure 3.1: On the left, the point placed in the scene captured by both cameras. On the right, line formed by stacking pixel captured along the camera track. Figure 3.2 depicts the model of a signal with constant depth. The red lines represent 58 the minimum and maximum depth. The same figure also shows the spectral support of light field. The blue line, with angle of 45ºrespect to the axis represents the constant depth. The red lines set the spectral support boundaries of depth. The tilt pictured in this domain can be translated into the depth of the signals in the scene. Figure 3.2: On the left, a model of constant depth. On the right, the frequency domain representation with boundaries for minimum and maximum depth. According to Chai[6], ”Any scene with a depth between zmin and zmax will have its continuous spectral support bounded in the frequency domain”. The use of epipolar images connects scene characteristics such as depth with the Fourier domain, which can be adequately manipulated. Finally, a couple of examples of EPIs and their 2D fast Fourier Transform (FFT2D) is shown in Figure 3.3. Figure 3.3: The image above displays an example of EPI and its Fourier transform. The image below, with an object behind the first shown above, is further in the scene, and its Fourier transform shows it in its tilt. 59 3.2 Transform-domain analysis of EPI images Epipolar plane images reconstructed with all the available views represent the perfect horizontal parallax-only scenario. No information is missing. However, one of the goals of this thesis is to reconstruct the original views with an smaller number of cameras. A straightforward way to generate this scenario is decimation, which directly discards explicit rows. For instance, a decimation of factor 2 would dismiss even rows obtaining a final image with half vertical resolution, a decimation of factor 4 would dismiss three rows for every row that remains and so on. Although this method to manipulate images is clear and it might sound correct initially, it introduces aliasing. Therefore, an anti-aliasing filter was required. The filter is a low-pass type and it discriminates the high-frequency information from the image. It will produce some visible changes in the plane color regions close to abrupt color changes. Aside from this collateral effect, it softens the sharpness on the borders. The loss of information in the pre-filtering is worth respect to the benefits provided by the posterior interpolations, when applied on anti-aliased images, as the final results will prove. The filter is applied column-wise given an image (img). With the intention of describing the essence of this technique, a one-dimensional signal (sgn) is considered, representing one column in the image. Figure 3.4 shows this process in more details. Firstly, the image needs to be enlarged vertically by double of the factor applied, and the first and last value of the signal extended to the beginning and end of the signal correspondingly. The reason of this arrangement resides in the circular property of Fourier transform. If this is not taken into account, the first pixels would interpolate with the last pixels, producing undesirable effects such as mirror at the extremes. Once the signal sgn is prepared, the Discrete Fourier Transform is applied. sgn fourier =DFT(sgn) The signal is now handled in the frequency domain. In this example an interpolation factor of 2 is considered, hence, half of the signal has to be filtered. The filter removes the middle half samples, corresponding to high-frequency. Besides this filtered, the signal has to be normalized, because it should not contain the same energy. The next formulas represent the theory behind this method for this example, not the actual development. normalization =N.ofcameras N.ofcameras 2+ 2 sgn decimated =sgn fourier(1 : N.ofcameras 4+ 2,3 4N.ofcameras + 1 : end) normalization 60 (a) (b) (c) (d) Figure 3.4: a) Original signal of 512 samples. b) Fourier transform of the signal. Highlighted area contains high-frequency region to be filtered. c) Filtered and normalized Fourier transform. d) Final signal, which softened and interpolated after the low-pass filter is applied. After this step, the signal is filtered in the frequency domain. The final phase include the Inverse Fourier transform, and discard the samples added at the extremes. sgn final extended =IDFTsng decimated sgn final =sgn final extended(2 : end −2) Figure 3.5 compares two epipolar images with and without pre-filtering. Nevertheless, simple decimation has been applied with the aim of comparing the results, obtained at the point of reconstructing the initially discriminated views. 3.3 Filtering of EPI and optimizing the number of views At this point of the process, the views are converted into EPIs in five available factors of interpolation. The actual optimization step takes place now, where the pre-filtered 61 completely black, by reason of no objects appearing in the top and bottom rows of the views. Figure 4.7 shows some output samples of this function. Figure 4.7: From left to right, samples of EPIs number 1, 100, 200, 300, 400, 500 and 600. As it was explained, the anti-aliasing filter is required before decimation. In the function presented below (a simplified version of the actual one) applies an anti-aliasing filter. It works with 5 different interpolation factors [2, 4, 8, 16, 32] power of 2, which facilitates the processing for Fourier transforms. The down-sampling produces an output shown in the figure 4.8. function prefiltering(object,density,v_res) % Object: Name of the scene % Density: Number of views -> Vertical resolution of EPIs % v_res: Vertical resolution of views -> Number of EPIs % It filters EPIs with anti-aliasing filter % and down-samples vertically in 5 factors [2, 4, 8, 16, 32] %Preparation of extended signal epi = im2double(...); epi_ext(1+factor:end-factor,:,:) = epi; epi_ext(1:factor,:,:) = epi(1,:,:); epi_ext(end-factor+1:end,:,:) = epi(end,:,:); %Fourier transform, filtering and normalization epi_f = fft(epi_ext); normal = density/((density/factor) + 2); fraction = density/(2*factor); epi_fdec = [epi_f(1:fraction+2,:,:); epi_f(end-fraction+1:end,:,:)] / normal; %Inverse Fourier transform and cropping extended samples. epi_factor = ifft(epi_fdec); epi_final = epi_factor(2:end-1,:,:); 4.3 Optimizing number of views At this moment in the design, the following developments will work on the views converted into filtered EPIs to bring them back to their original resolution. The up-sampling that optimizes the number of views necessary to build the multi-view system. 68 Figure 4.8: The original EPI and the 5 anti-aliased EPIs with 5 different interpolation factors. The scripts responsible for this step are presented below: function run_epi_interp(object) % Object: Name of the scene % It runs epi_interp() with 5 factors [2, 4, 8, 16, 32] and % 5 interpolation kernels [triangle, cubic, lanczos2, lanczos3, fastdct] function epi_interp(object,density,factor,v_res,interpolation) % Object: Name of the scene % Density: Number of views -> Vertical resolution of EPIs % Factor: Interpolation magnitude % v_res: Vertical resolution of views -> Number of EPIs % interpolation: Interpolation kernel used % It up-samples vertically EPIs The function epi interp uses either imresize() Matlab function for bilinear, bicubic, lanczos2 and lanczos3 kernels interpolation, or the next function fdctEPI, which implements the fast DCT-based scaling algorithm [10]. This algorithm was implemented to scaling 1D signals, and later it was adapted to two dimensions. All the scripts coded that improved that code time-wise in every version are shown in the appendices. The next sample of code (simplified version) aims to demonstrate and present the basics of the fast DCT interpolation. function img_out = fdctEPI(img_in,factor) % img_in: input image % factor: interpolation magnitud % It up-samples vertically img_in using fastDCT interpolation No = rows_in; n = repmat(0:rows_in-1,rows_out,1); k = repmat((0:rows_out-1)’,1,rows_in); DCT = sincd(2*No-2, 2*rows_in, (k+1/2)/factor -n -1/2) .* ... 69 cos( pi/(2*rows_in) * ( (k+1/2)/factor -n -1/2) ) + ... sincd((2*No-2), 2*rows_in, (k+1/2)/factor +n +1/2) .* ... cos( pi/(2*rows_in) * ( (k+1/2)/factor +n +1/2) ); % Multiplication by DCT-Matrix to all image component layers for i=1:components img_out(:,:,i) = DCT*img_in(:,:,i); end % Truncating negative and greater than 1 values img_out(img_out<0) = 0; img_out(img_out>1) = 1; function result = sincd(M,N,x) % sincd(M,N,x) = sin(piMx/N)/(Nsin(pix/N)) % Digital sinc function used in fastDCT One of the critic code sections in fdctEPI() is the population of nand k, as they represent the matrices which substitute the initial loops, operated by rows and columns. Executing this operation with two nested loops makes its computation inefficient. In a preliminary version one of the loops was substituted by a vectorized operation, but, populating two matrices and carrying linear operation is more coherent. Matrix kand nrespond to the following pattern: k=      0 0 · · · 0 1 1 · · · 1 . . .. . ..... . . n n · · · n       n=      0 1 · · · n 0 1 · · · n . . .. . ..... . . 0 1 · · · n       Three methods of populating matrices might be considered. % Kronecker tensor product ns = kron([1:size_in]-1,ones(1,size_out)’); ks = kron(ones(1,size_in),[1:size_out]’-1); % Replicating matrices ns = repmat([0:size_in-1],size_out,1); ks = repmat([0:size_out-1]’,1,size_in); % Manual with ones temp= [0:size_in-1]; ns= temp(ones(1,size_out),:); temp= [0:size_out-1]’; ks= temp(:,ones(1,size_in)); 70 The results for resolution 300, 600 and 900 pixels and factors 2-10 are displayed in Figure 4.9. Figure 4.9: Kron method is notably worse than the other two. There is no difference between repmat and manual, at the same time, manual requires less memory and number of instructions. The reconstruction of the views from the interpolated EPIs using the kernels explained above is straightforward. A view is formed by stacking rows coming from the same row of every EPI. The following functions are responsible for the final phase. function run_reconstruct_views(object) % object: Name of the scene % It runs reconstruct_views() with 5 factors [2, 4, 8, 16, 32] and % 5 interpolation kernels [triangle, cubic, lanczos2, lanczos3, fastdct] function reconstruct_views(object,density,factor,v_res,interpolation) % Object: Name of the scene % Density: Number of views -> Vertical resolution of EPIs % Factor: Factor of interpolation applied vertically % v_res: Vertical resolution of views -> Number of EPIs % interpolation: Method/kernel used for re-sampling % It reconstructs all the views from the re-sampled EPIs The whole process seen from a graphical perspective is pictured in Figure 4.10. 71 Figure 4.10: The figure shows all the steps taken in this work, and how images are transformed and manipulated in every step. 72 5. Results In this chapter all graphic and numeric results will be presented and interpreted. The main tests executed in this work involved the fast DCT algorithm [10] and other interpolation kernels, as well as the final comparison tests between original and interpolated views. 5.1 Fast DCT algorithm Fast DCT scaling algorithm has been used to up-sample vertically epipolar planar images. Before the image tests, a set of peculiar one-dimensional signals was tested towards a better understanding of its effects in the particular cases. The next figure (5.1) shows the effects of scaling in some peculiar signals and simple images based on those signals. Figure 5.1: Scaling produced over the sharp signal causes sinusoidal effect. 73 Consequence of this effect, the borders of epipolar images, which contain basically lines and stripes, present blurring. Another example, tested to confront the outcome of the paper [10] against this work, presents iteratively scaling of a text image. The results throw the impression that regardless how many iterations are executed, the text remains readable unlike with other interpolation kernels. Figure 5.2 show the difference between interpolation kernels clearly. (a) (b) (c) (d) Figure 5.2: a) Original text image. 75 iterative zoom-in&zoom-out scaling, factor sqrt(2) b) by bilinear algorithm, c) by bicubic algorithm, d) by Fast DCT algorithm. As it was mentioned above, the border artefacts present will have to be tackled down in posterior steps. The artefacts appearing beyond the borders will take its utility away. Figure 5.3 depicts why the displayed results are not useful for epipolar image interpolation. (a) (b) Figure 5.3: Text image scaled 1000 times with factor sqrt(2). a) Region of output text image displayed alike [10]. b) Complete text image, contains the artifacts on the borders. Another aspect considered in this text is the computing time spent by all these methods. The results presented in Figure 5.4 prove that the implementation of Fast DCT coded in Matlab is slower than built-in functions Matlab, since fast DCT is comprised of two external functions. The hardware used in the execution of these tests does not influence the results, as we are not comparing time-wise data with any other results obtained in other hardware. Furthermore, every test was executed 10 times in the same conditions. Figure 5.4 show the average of all the executions. 74 Figure 5.4: The implementation of the fast DCT algorithm is considerably slower than built-in algorithms. 5.2 Experimental results The most important result to be discussed involves original and interpolated EPIs. This comparison was accomplished through PSNR (Peak Signal-to-Noise Ratio) measurement, as explained in the analysis chapter. The first experimental result with PSNR checks every reconstructed view for every interpolation method. The following script was written to obtain the results: function test_psnr(object,density,factor) % Object: Name of the scene % Density: Number of views % Factor: Interpolation magnitud % It calculates the PSNR for every factor and view The visible peaks that appear in Figure 5.5 represent the true views, non-interpolated remaining in the filtered EPI. A view closer to a true view will get, therefore, a better PSNR result. According to this, the extreme views might show the worst values, since they do not have neighbour views to interpolate with. Figure 5.5 belongs to the test with factor 8. The tests with the rest of interpolation factors are available in the appendices. 75 Figure 5.5: Peaks represent original and therefore non-interpolated views. Views on the extreme tend to generate worse result because of their isolation. Also, another script was written to compare all the interpolation factors in the same figure. The complete script is available in the appendix: function PSNR = test_psnr_average(object,density) % Object: Name of the scene % Density: Number of views % Factor: Interpolation magnitud % It calculates the PSNR for every factor and its average between all the % views The same graph is used to contrast the PSNR values obtained in the tests with prefiltered EPIs confronted to the same test with only decimated EPIs. In the case of views based on aliased epipolar images, the maximum factor applied was 12, since it was expected that greater factors would provide PSNR values utterly unacceptable. Pre-filtered EPIs-based views are presented as points, and non-pre-filtered as lines. The first aspect to notice is that the difference between interpolation kernels is not sufficiente in such magnitude, hence, the interpolation factor will be the decisive parameter. Unlike it could have seemed in the analysis of Fourier domain of EPIs, dismissing certain information related to high-frequency, it is proved now that an anti-alising filter is sufficient useful. In the next chapter, this graph 5.6 will be discussed further. The minimum number of cameras needed to reconstruct an acceptable multi-view system will be determined accordingly to this last figure. 76 Figure 5.6: Views based on pre-filtered EPIs obtain much better results for the same interpolation factor than views based on simple decimated EPIs. 77 84 7.2 Metodolog´ıa de trabajo La metodolog´ıa utilizada en este proyecto se caracteriza por circunstancias particulares. Previa a la presentaci´on del proyecto como una MSc Thesis, se solicit´o el dise˜no de una versi´on preliminar de la aplicaci´on Blender, la cual renderiza la escena desde diferentes puntos de vista, como una prueba de conocimientos de programaci´on para continuar con el resto del desarrollo. M´as adelante, despu´es del an´alisis de trabajos anteriores, cada fase fue analizada, dise˜nada, implementada y finalmente testada. Si los tests muestran que los requisitos se han cumplido, la fase de desarrollo se considera finalizada. La escritura de la memoria comenz´o simult´aneamente con las fases de dise˜no y testado. 7.3 Diagrama de Gantt Figura 7.1: Diagrama de Gantt. El diagrama muestra un periodo en verano donde el proyecto estuvo paralizado. La raz´on es un trabajo de verano en Tampere, Finlandia, que me permiti´o extender mi estancia en este pa´ıs desde Septiembre en adelante, cuando las ayudas econ´omicas de la beca Erasmus hab´ıan finalizado. 7.4 Herramientas utilizadas A lo largo del an´alisis, dise˜no y escritura de la tesis se utilizaron varias herramientas. Blender es una plataforma de c´odigo abierto para gr´aficos por ordenador 3D, utilizada para configurar el escenario con la escena virtual y c´amaras, as´ı como herramienta de rendering. Los scripts escritos en Blender est´an basados en Python. La herramienta principal usada a lo largo del proyecto ha sido Matlab, para el procesado de im´agenes y 85 el desarrollo y testado de algoritmos. En la escritura de la presente memoria, el kernel LaTeX y el editor TeX Maker han sido los elegidos, junto con GIMP para manipular las im´agenes presentadas en este trabajo. 7.5 C´odigo - Sistema de rendering multi-vista Fichero holographic camera 1 2 0.py # Author: David Valle # Departament of Signal Processing (SGN) # David Valle (230113) [email protected] # bl_info = { ’name’: "Holographic Camera", ’author’: "David Valle <[email protected]>", ’version’: (1, 0, 0), ’blender’: (2, 6, 0), ’api’: 41098, ’location’: "Select a Camera > Properties Panel > Camera Panel > Holographic Camera", ’description’: "Render images for a future holographic model", ’warning’: "", ’wiki_url’: "", ’tracker_url’: "", ’category’: "Object"} import bpy import mathutils from math import * from bpy.props import * # GUI (Panel) class OBJECT_PT_holographic_camera(bpy.types.Panel): bl_label = "Holographic Camera" bl_space_type = "PROPERTIES" bl_region_type = "WINDOW" bl_context = "data" # show this add-on only in the Camera-Data-Panel @classmethod 86 def poll(self, context): return context.active_object.type == ’CAMERA’ bpy.types.Object.holographic_width_field = bpy.props.FloatProperty( attr="holographic_width_field", name=’holographic_width_field’, description=’Camera Width Field’, min=0.0, soft_min=0.0, max=1000000, soft_max=1000000, default=100) bpy.types.Object.holographic_field_density = bpy.props.IntProperty( attr="holographic_field_density", name=’holographic_field_density’, description=’Camera Field Density’, min=0, soft_min=0, max=1000000, soft_max=1000000, default=5) bpy.types.Object.holographic_distance_to_scene = bpy.props.FloatProperty( attr="holographic_distance_to_scene", name=’holographic_distance_to_scene’, description=’Distance to Scene’, min=0.000, soft_min=0.000, max=1000000, soft_max=1000000, default=300) bpy.types.Object.show_mesh_distance = bpy.props.BoolProperty( attr="show_mesh_distance", name="show_mesh_distance", default=False) bpy.types.Object.holocamera_type = bpy.props.EnumProperty( attr="holocamera_type", items=( ("PARALLEL", "Parallel", "Parallel camera distribution"), ("CONVERGE", "Converge", "Camera distribution with convergency")), name="holocamera_type", description="", default="PARALLEL") # Draw def draw(self, context): layout = self.layout holocamera = context.scene.camera row = layout.row() row.prop(holocamera, "holocamera_type", text="Holographic Camera Type", expand=True) row = layout.row() row.prop(holocamera, "holographic_width_field", text="Width Field") 87 row = layout.row() row.prop(holocamera, "holographic_field_density", text="Field Density") if holocamera.holocamera_type == "CONVERGE": row = layout.row() row.prop(holocamera, "holographic_distance_to_scene", text="Distance to Scene") if holocamera.holocamera_type == "PARALLEL": pass row = layout.row() row.separator() row = layout.row() row.operator(’holocamera.render_holographic_cameras’) row = layout.row() row.operator(’holocamera.clean_holographic_cameras’) # Operator ’Render Holographic Cameras’ class OBJECT_OT_render_holographic_cameras(bpy.types.Operator): bl_label = ’Render Holocameras’ bl_idname = ’holocamera.render_holographic_cameras’ bl_description = ’Render the Holographic Cameras’ bl_options = {’REGISTER’, ’UNDO’} #Call the operator ’Render Holographic Cameras’ def execute(self, context): cam = bpy.data.objects["Camera"] density = cam.holographic_field_density width = (cam.holographic_width_field/1000)/(density-1) if cam.holocamera_type == "CONVERGE" : #Distance to focus, 1/100 blender units distance= cam.holographic_distance_to_scene #Save original location of camera to restore it at the end orig_loc = cam.location.copy() #Save original rotation of camera to restore it at the end orig_rot = cam.rotation_euler.copy() local_pos = cam.location.copy() 88 local_rot = cam.rotation_euler.copy() print(’local_rot(x,y,z)= ’+str(local_rot)) #Hint: Way to bring local coord. to world coord. (world = local*matrix) matrix = bpy.data.objects[’Camera’].matrix_world.copy() sceneKey = bpy.data.scenes.keys()[0] scene = bpy.data.scenes[sceneKey] scene.camera = cam for object in scene.objects: if object.name[:7] == "Camera.": scene.objects.unlink(object) i = floor(density / 2) if 0 == 0 : #for i in range(density): for j in range(density): #Currently location of camera print() print(’(’+str(i)+’,’+str(j)+’)’) local_pos.x = (j - (density-1)/2)*width local_pos.y = ((density-1)/2 - i)*width cam.location = matrix*local_pos.copy() print(’local_pos(x,y,z)= ’+str(local_pos)) print(’world_pos(x,y,z)= ’+str(cam.location)) #Rotate cameras in converge method if cam.holocamera_type == "CONVERGE" : #Currently rotation of camera cam.rotation_euler.x -= atan( ((density - 1) / 2 - i)*width / distance) cam.rotation_euler.y -= atan( ((density - 1) / 2 - j)*width / distance) cam.rotation_euler.z = 0.0 #Set output file path #scene.render.filepath = ’//chameleon/converge_600/’+str(i)+’_’+str(j) scene.render.filepath = ’//chameleon/converge_600/’+str(j) else: 89 #Parallel #scene.render.filepath = ’//chameleon/parallel_300/’+str(i)+’_’+str(j) scene.render.filepath = ’//chameleon/parallel_r_512/’+str(j+1) #Create multiple cameras #bpy.ops.object.camera_add(location=(cam.location.x, cam.location.y, cam.location.z),rotation=(cam.rotation_euler.x, cam.rotation_euler.y, cam.rotation_euler.z)) #Render camera bpy.ops.render.render( write_still=True ) #Restoring initial position and rotation of camera cam.location = orig_loc.copy() cam.rotation_euler = orig_rot.copy() cam = bpy.data.objects["Camera"] bpy.context.scene.objects.active = cam return {’FINISHED’} # Operator ’Render Holographic Cameras’ class OBJECT_OT_render_holographic_cameras(bpy.types.Operator): bl_label = ’Clean Holocameras’ bl_idname = ’holocamera.clean_holographic_cameras’ bl_description = ’Clean the Holographic Cameras’ bl_options = {’REGISTER’, ’UNDO’} #Call the operator ’Render Holographic Cameras’ def execute(self, context): sceneKey = bpy.data.scenes.keys()[0] scene = bpy.data.scenes[sceneKey] for object in scene.objects: if object.name[:7] == "Camera.": scene.objects.unlink(object) cam = bpy.data.objects["Camera"] bpy.context.scene.objects.active= cam 90 return {’FINISHED’} # Register def register(): bpy.utils.register_module(__name__) def unregister(): bpy.utils.unregister_module(__name__) if __name__ == "__main__": register() 91 7.6 C´odigo - Recortado de renders Fichero crop renders.m % Tampere University of Technology (TUT) % Departament of Signal Processing (SGN) % David Valle (230113) [email protected] function crop_renders(object,density,h_res) % Object: Name of the scene % Density: Number of views -> Vertical resolution of EPIs % h_res: Horizontal resolution of views % It crops renders to normal resolution name_in = strcat(object,’/parallel_r_’,int2str(density),’/’); name_out = strcat(object,’/parallel_rc_’,int2str(density),’/’); for i=1:density i img = imread(strcat(name_in,int2str(density-i+1)),’png’); cropped = img(:,(i-1)*2+1:(i-1)*2+h_res,:); imwrite(cropped,strcat(name_out,int2str(density-i+1),’.png’)); end end 92 7.7 C´odigo - Creaci´on de EPIs con densidad completa Fichero epi.m % Tampere University of Technology (TUT) % Departament of Signal Processing (SGN) % David Valle (230113) [email protected] % Creation of EPI’s (full density) function epi(object,density,v_res) % Object: Name of the scene % Density: Number of views -> Vertical resolution of EPIs % v_res: Vertical resolution of views -> Number of EPIs % It creates epipolar images with full density % Concatenation of input/output path name_in = strcat(object,’/parallel_rc_’,int2str(density),’/’); name_out = strcat(object,’/epis_’,int2str(density),’/’); % Construction of EPI, one per row, same amount as vertical % resolution of views for i=1:v_res i for j=1:density img = imread(strcat(name_in,int2str(j)),’png’); img_epi(j,:,:) = img(i,:,:); end imwrite(img_epi,strcat(name_out,int2str(i),’.png’)); end end 93 7.11 C´odigo - Algoritmos Fast DCT Fichero fast dct1D.m con sus 3 versiones mejoradas. % Tampere University of Technology (TUT) % Departament of Signal Processing (SGN) % David Valle (230113) [email protected] % Implementation of a fast DCT-based interpolation % a = SUM(n=0,N-1)an( A + B ) % A = sincd(2No-2, 2N, (k+1/2)/d -n -1/2 ) * cos(pi/2N((k+1/2)/d -n -1/2)) % B = sincd(2No-2, 2N, (k+1/2)/d +n +1/2 ) * cos(pi/2N((k+1/2)/d +n +1/2)) % No = min(N,|dN|) % sincd(M,N,x) = sin(piMx/N)/(Nsin(pix/N)) % Explanation of variables % a = Output Signal % an = Input Signal % N = signal length % k = current output signal index % d = scaling factor function output = fast_dct1D(sample,factor) if( round( log2(size(sample,2)) ) ˜= log2(size(sample,2)) ) error(’Lenght of input signal should be power of 2’); end if( round( log2(factor) ) ˜= log2(factor) ) error(’Factor should produce an output signal power of 2’); end %1D Sample size_in= size(sample,2); size_out= size_in*factor; output= zeros(1,size_out); if factor >= 1 No= min(size_in,ceil(factor*size_in)); else No= min(size_in,floor(factor*size_in)); end %Version 1 100 for k= 0:size_out-1 for n= 0:size_in-1 part1= sincd(2*No-2, 2*size_in, (k+1/2)/factor -n -1/2); part3= sincd(2*No-2, 2*size_in, (k+1/2)/factor +n +1/2); part2= cos( (pi/(2*size_in) ) * ( (k+1/2)/factor -n -1/2) ); part4= cos( (pi/(2*size_in) ) * ( (k+1/2)/factor +n +1/2) ); output(k+1)= output(k+1) + sample(n+1)*(part1*part2+part3*part4); end end %Version 2 DCT= zeros(size_out,size_in); %Calc of DCT matrix for k= 0:size_out-1 DCT(k+1,1:size_in)= ... sincd(2*No-2, 2*size_in, (k+1/2)/factor -(0:size_in-1) -1/2) .* ... cos( pi/(2*size_in) * ( (k+1/2)/factor -(0:size_in-1) -1/2) ) + ... sincd(2*No-2, 2*size_in, (k+1/2)/factor +(0:size_in-1) +1/2) .* ... cos( pi/(2*size_in) * ( (k+1/2)/factor +(0:size_in-1) +1/2) ); end output= (DCT*sample’)’; %Final version %Constructing DCT matrix n = repmat(0:size_in-1,size_out,1); k = repmat((0:size_out-1)’,1,size_in); DCT=sincd(2*No-2, 2*size_in, (k+1/2)/factor -n -1/2) .* ... cos( pi/(2*size_in) * ( (k+1/2)/factor -n -1/2) ) + ... sincd((2*No-2), 2*size_in, (k+1/2)/factor +n +1/2) .* ... cos( pi/(2*size_in) * ( (k+1/2)/factor +n +1/2) ); output= (DCT*sample’)’; end Fichero fast dct2D.m con sus 3 versiones mejoradas. % Tampere University of Technology (TUT) % Departament of Signal Processing (SGN) % David Valle (230113) [email protected] % Implementation of a fast DCT-based interpolation in 2 dimensions %function output = fast_dct2D(image,factor_x,factor_y) 101 function e = fast_dct2D t = cputime; %Important to convert to double, imread returns uint8. img= im2double(imread(’chrome.png’)); img_in= im2double(imread(’chrome.png’)); rows= size(img,1); cols= size(img,2); factor= 1; cols_out= cols*factor; img_rows= zeros(rows,cols_out); rows_out= rows*factor; img_out= zeros(rows_out,cols_out); % Version 1 for n=1:1 %Zoom-in for i=1:rows img_rows(i,:)= fast_dct1D(img_in(i,:),factor); end for i=1:cols_out img_out(:,i)= fast_dct1D(img_rows(:,i).’,factor).’; end %Zoom-out for i=1:cols_out img_rows(:,i)= fast_dct1D(img_out(:,i).’,1/factor).’; end for i=1:rows img_in(i,:)= fast_dct1D(img_rows(i,:),1/factor); end end % Version 2 %Pre-allocating matrices img_out= zeros(size_out,size_out); DCT= zeros(size_out,size_in); for k= 0:size_out-1 DCT(k+1,1:size_in)= ... sincd(2*No-2, 2*size_in, (k+1/2)/factor -(0:size_in-1) -1/2 ) .* ... cos( pi/(2*size_in) * ( (k+1/2)/factor -(0:size_in-1) -1/2 ) ) + ... sincd(2*No-2, 2*size_in, (k+1/2)/factor +(0:size_in-1) +1/2 ) .* ... cos( pi/(2*size_in) * ( (k+1/2)/factor +(0:size_in-1) +1/2 ) ); 102 end % Final version n = repmat(0:size_in-1,size_out,1); k = repmat((0:size_out-1)’,1,size_in); DCT=sincd(2*No-2, 2*size_in, (k+1/2)/factor -n -1/2) .* ... cos( pi/(2*size_in) * ( (k+1/2)/factor -n -1/2) ) + ... sincd((2*No-2), 2*size_in, (k+1/2)/factor +n +1/2) .* ... cos( pi/(2*size_in) * ( (k+1/2)/factor +n +1/2) ); img_out= DCT*img_in*DCT’; % Truncating negative and greater than 1 values img_out(img_out<0)= 0; img_out(img_out>1)= 1; end 103 7.12 C´odigo - Tests Fast DCT Fichero run fdct1D.m % Tampere University of Technology (TUT) % Departament of Signal Processing (SGN) % David Valle (230113) [email protected] % Testing Fast DCT-based interpolation %Ramp signal signal= 1:128; factor= 3; figure output= fastdct1D(signal,factor); subplot(1,2,1), plot(1:size(signal,2),signal,’b’,1:size(output,2),output,’r’); title(sprintf(’Fast-DCT interpolation factor %i’,factor)); output= fdct1D(signal,1/factor); subplot(1,2,2), plot(1:size(signal,2),signal,’b’,1:size(output,2),output,’r’); title(sprintf(’Fast-DCT interpolation factor 1/%i’,factor)); legend(’input’,’output’); clear all %Ramp with step signal signal= 1:128; for i= 64:96 signal(i)= i-20; end factor= 2; figure output= fastdct1D(signal,factor); subplot(1,2,1), plot(1:size(signal,2),signal,’b’,1:size(output,2),output,’r’); title(sprintf(’Fast-DCT interpolation factor %i’,factor)); output= fdct1D(signal,1/factor); subplot(1,2,2), plot(1:size(signal,2),signal,’b’,1:size(output,2),output,’r’); title(sprintf(’Fast-DCT interpolation factor 1/%i’,factor)); legend(’input’,’output’); clear all %Sinusoidal signal signal= sind(1:512); factor= 2; figure output= fastdct1D(signal,factor); 104 subplot(1,2,1), plot(1:size(signal,2),signal,’b’,1:size(output,2),output,’r’); title(sprintf(’Fast-DCT interpolation factor %i’,factor)); output= fdct1D(signal,1/factor); subplot(1,2,2), plot(1:size(signal,2),signal,’b’,1:size(output,2),output,’r’); title(sprintf(’Fast-DCT interpolation factor 1/%i’,factor)); legend(’input’,’output’); clear all %Pulse signal for i=1:64 signal(i)= 20; end for i=65:192 signal(i)= 200; end for i=193:256 signal(i)= 20; end factor= 2; figure output= fastdct1D(signal,factor); subplot(1,2,1), plot(1:size(signal,2),signal,’b’,1:size(output,2),output,’r’); title(sprintf(’Fast-DCT interpolation factor %i’,factor)); output= fdct1D(signal,1/factor); subplot(1,2,2), plot(1:size(signal,2),signal,’b’,1:size(output,2),output,’r’); title(sprintf(’Fast-DCT interpolation factor 1/%i’,factor)); legend(’input’,’output’); clear all Fichero run fdct2D.m % Tampere University of Technology (TUT) % Departament of Signal Processing (SGN) % David Valle (230113) [email protected] % Testing Fast DCT-based interpolation clear all; tic; img_in=im2double(imread(’text.png’)); out2=img_in; for i=1:1000 out=fdct2D(out2,2); out2=fdct2D(out,2); end 105 toc; figure; imshow(img_in); figure; imshow(out); imwrite(out,’text75_out.png’); figure; imshow(out2); imwrite(out2,’text75_out2.png’); Fichero population test.m % Tampere University of Technology (TUT) % Departament of Signal Processing (SGN) % David Valle (230113) [email protected] % Testing fastest way to populate matrix with given parameters % Parameters are according renders we should take % Around 600x600pixels and factor [2:10] clear all times= []; tests= [2:10; repmat([300; 600; 900],1,9)]; for j=1:size(tests,2) for i=2:size(tests,1) factor= tests(1,j); size_in= tests(i,j); size_out=factor*size_in; clear ns ks tic; ns = kron([1:size_in]-1,ones(1,size_out)’); ks = kron(ones(1,size_in),[1:size_out]’-1); times(j,i-1) = toc; clear ns ks tic; ns = repmat([0:size_in-1],size_out,1); ks = repmat([0:size_out-1]’,1,size_in); times(j,i+2) = toc; clear ns ks tic; temp= [0:size_in-1]; ns= temp(ones(1,size_out),:); temp= [0:size_out-1]’; ks= temp(:,ones(1,size_in)); 106 times(j,i+5) = toc; end end figure; bar(times,’grouped’); title(’Comparison population of matrix’,’FontWeight’,’bold’); times_leg = legend(’kron300’,’kron600’,’kron900’, ... ’repmat300’,’repmat600’,’repmat900’, ... ’ones300’,’ones600’,’ones900’); set(times_leg,’Location’,’NorthWest’) xlabel(’Cases: 2 to 10 factor, 300 - 600 - 900pixels’); ylabel(’Seconds’); 107 7.13 C´odigo - Tests Fichero run epi interp.m % Tampere University of Technology (TUT) % Departament of Signal Processing (SGN) % David Valle (230113) [email protected] function run_epi_interp(object) % Object: Name of the scene % It runs epi_interp() with 5 factors [2, 4, 8, 16, 32] and % 5 interpolation kernels [triangle, cubic, lanczos2, lanczos3, fastdct] disp(’Processing epi_interp factor 2 density 256’); epi_interp(object,512,2,600,’triangle’); epi_interp(object,512,2,600,’cubic’); epi_interp(object,512,2,600,’lanczos2’); epi_interp(object,512,2,600,’lanczos3’); epi_interp(object,512,2,600,’fastdct’); disp(’Processing epi_interp factor 4 density 128’); epi_interp(object,512,4,600,’triangle’); epi_interp(object,512,4,600,’cubic’); epi_interp(object,512,4,600,’lanczos2’); epi_interp(object,512,4,600,’lanczos3’); epi_interp(object,512,4,600,’fastdct’); disp(’Processing epi_interp factor 8 density 64’); epi_interp(object,512,8,600,’triangle’); epi_interp(object,512,8,600,’cubic’); epi_interp(object,512,8,600,’lanczos2’); epi_interp(object,512,8,600,’lanczos3’); epi_interp(object,512,8,600,’fastdct’); disp(’Processing epi_interp factor 16 density 32’); epi_interp(object,512,16,600,’triangle’); epi_interp(object,512,16,600,’cubic’); epi_interp(object,512,16,600,’lanczos2’); epi_interp(object,512,16,600,’lanczos3’); epi_interp(object,512,16,600,’fastdct’); disp(’Processing epi_interp factor 32 density 16’); epi_interp(object,512,32,600,’triangle’); 108 epi_interp(object,512,32,600,’cubic’); epi_interp(object,512,32,600,’lanczos2’); epi_interp(object,512,32,600,’lanczos3’); epi_interp(object,512,32,600,’fastdct’); figure; timing= reshape(timing,5,6)’; plot(timing); title(’Interpolation time comparison’,’FontWeight’,’bold’); timing_leg = legend(’triangle’,’cubic’,’lanczos2’, ’lanczos3’,’fastdct’); set(timing_leg,’Location’,’NorthWest’) xlabel(’Interpolation factors: 32, 16, 8, 6, 4, 2’); ylabel(’Seconds’); end Fichero run reconstruct views.m % Tampere University of Technology (TUT) % Departament of Signal Processing (SGN) % David Valle (230113) [email protected] % Script that runs reconstruct_view with all possible configurations function run_reconstruct_views(object) % object: Name of the scene % It runs reconstruct_views() with 5 factors [2, 4, 8, 16, 32] and % 5 interpolation kernels [triangle, cubic, lanczos2, lanczos3, fastdct] disp(’Processing reconstruct_views factor 2 density 256’); reconstruct_views(object,512,2,600,’triangle’); reconstruct_views(object,512,2,600,’cubic’); reconstruct_views(object,512,2,600,’lanczos2’); reconstruct_views(object,512,2,600,’lanczos3’); reconstruct_views(object,512,2,600,’fastdct’); disp(’Processing reconstruct_views factor 4 density 128’); reconstruct_views(object,512,4,600,’triangle’); reconstruct_views(object,512,4,600,’cubic’); reconstruct_views(object,512,4,600,’lanczos2’); reconstruct_views(object,512,4,600,’lanczos3’); reconstruct_views(object,512,4,600,’fastdct’); disp(’Processing reconstruct_views factor 8 density 64’); reconstruct_views(object,512,8,600,’triangle’); 109 Figura 7.6: PSNR de las vistas reconstruidas con cada m´etodo de interpolaci´on, factor 32. 116