Paralelización de un algoritmo de ray tracing para arrays de mililentes

Instituto Universitario de Sistemas Inteligentes y Aplicaciones Numéricas en Ingeniería TRABAJO FIN DE MÁSTER PARALELIZACIÓN DE UN ALGORITMO DE RAY TRACING PARA ARRAYS DE MILILENTES Titulación: Máster Oficial en Sistemas Inteligentes y Aplicaciones Numéricas en Ingeniería Autor: D. Antonio Cristo Suárez Rodríguez Tutores: Dr. Eduardo M. Rodríguez Barrera Dr. Domingo J. Benítez Díaz Fecha: Julio de 2015 "There are three kinds of lies: lies, damned lies, and benchmarks." Variación de una frase falsamente atribuida a Benjamin Disraeli por Mark Twain. Índice general I Memoria 1 1. Introducción 3 1.1. Motivación .................................. 3 1.2. Descripción del problema .......................... 4 1.3. Objetivos .................................. 6 1.4. Organización de la memoria ........................ 6 2. Antecedentes 9 2.1. Comunicaciones por luz visible ....................... 9 2.2. Visores tridimensionales .......................... 11 2.2.1. Barrera de paralaje ......................... 12 2.2.2. Lentes lenticulares ......................... 12 2.3. Tecnologías de paralelización ........................ 12 2.3.1. Memoria compartida ........................ 12 2.3.2. GPU ................................. 14 3. Modelo matemático 19 3.1. Justificación del uso de mililentes ..................... 19 3.2. Óptica geométrica .............................. 19 3.2.1. Reflexión .............................. 21 3.2.2. Refracción .............................. 21 3.2.3. Ángulo crítico ............................ 22 3.2.4. Ecuaciones de Fresnel ........................ 23 3.3. Fuentes ópticas ............................... 23 3.4. Lentes .................................... 25 3.5. Preprocesado ................................ 26 3.5.1. Preprocesado en azimut ...................... 26 3.5.2. Preprocesado en elevación ..................... 28 4. Implementación del sistema 31 4.1. Topología del problema ........................... 31 4.2. Estructuras de datos ............................ 31 4.3. Explotación de la simetría ......................... 33 4.4. Descripción del algoritmo ......................... 34 4.4.1. Definición del escenario ....................... 34 4.4.2. Cálculo de límites .......................... 36 4.4.3. Simulación .............................. 36 4.5. Código secuencial .............................. 39 4.6. Código memoria compartida ........................ 40 v vi ÍNDICE GENERAL 4.7. Código GPU ................................. 41 5. Resultados 43 5.1. Equipo empleado .............................. 43 5.2. Estudio de la precisión ........................... 43 5.3. Resultados de la paralelización ....................... 47 5.3.1. OpenMP ............................... 47 5.3.2. CUDA R .............................. 49 5.4. Resultados del simulador .......................... 56 6. Conclusiones 61 II Bibliografía 63 Bibliografía 65 Índice de figuras 1.1. Escenario del problema ........................... 4 1.2. Flujo del sistema .............................. 6 2.1. Esquemas DCO-OFDM y ACO-OFDM .................. 10 2.2. Comparación entre barrera de paralaje y lentes lenticulares ....... 11 2.3. Esquema de funcionamiento del modelo fork/join ............ 13 2.4. Comparación del espacio dedicado a cómputo entre una CPU y una GPU 15 2.5. Arquitectura del entorno CUDA R .................... 16 2.6. Arquitectura Nvidia Kepler ........................ 16 2.7. Transferencias de memoria entre el host y el device ........... 17 2.8. Jerarquía de bloques e hilos ........................ 18 3.1. RED ONE R ............................... 20 3.2. Ley de Snell ................................. 20 3.3. Diagrama de radiación lambertiano para m= 20 ............. 24 3.4. Lente planoconvexa ............................. 25 3.5. Imagen captada con un array de microlentes encima del sensor . . . . . 25 3.6. MLA150-5C ................................. 26 3.7. Situación de partida ............................ 27 3.8. Preprocesado en azimut .......................... 27 3.9. Preprocesado en elevación ......................... 29 3.10. Cálculo de los ángulos de elevación .................... 29 3.11. Función a resolver de forma numérica ................... 30 4.1. Diagrama del algorimo ........................... 32 4.2. Clúster de 9 mililentes y 4 ledes por mililente ............... 34 4.3. Malla de ángulos sólidos a simular sobre el plano XY .......... 36 4.4. Impactos sobre el plano imagen ...................... 37 4.5. Formato del vector de rayos ........................ 40 5.1. Precisión en función del ángulo mitad (índice lambertiano) ....... 45 5.2. Precisión en función de la separación de la fuente ............ 46 5.3. Precisión en función del nivel de muestreo angular ............ 47 5.4. Aceleraciones para memoria compartida ................. 48 5.5. Eficiencias para memoria compartida ................... 49 5.6. Aceleraciones para memoria compartida en doble precisión ....... 50 5.7. Eficiencias para memoria compartida en doble precisión ......... 50 5.8. Aceleraciones para GPU de 1 a 32 hilos por bloque ........... 52 5.9. Aceleraciones para GPU de 64 a 1024 hilos por bloque ......... 52 vii viii ÍNDICE DE FIGURAS 5.10. Aceleraciones para GPU contando sólo la simulación ........... 53 5.11. Aceleración media según el número de hilos por bloque ......... 54 5.12. Aceleraciones para GPU con memoria no paginada ........... 54 5.13. Ejecución solapada de kernel y transferencia de memoria ........ 55 5.14. Aceleraciones con ejecución solapada para distinto número de flujos . . 56 5.15. Renderizado de un array de 128 ×128 .................. 58 5.16. Clúster de 49 mililentes y 64 ledes ..................... 59 5.17. Renderizado de un array de 128 ×128 modificado ............ 59 Índice de tablas 4.1. Estructura de datos para un rayo ..................... 32 4.2. Estructura de datos para una lente .................... 32 4.3. Estructura de datos para una interfaz ................... 33 4.4. Resultados del profiler ........................... 40 5.1. Características del equipo empleado .................... 44 5.2. Error cuadrático medio según el tipo de coma flotante empleado . . . . 51 5.3. Resultados de nvprof ............................ 52 5.4. Resultados de nvprof para pinned memory ................ 54 5.5. Resultados de nvprof para ejecución solapada de 2 flujos ........ 55 ix xvi Lista de acrónimos Acrónimo Descripción PHY Physical Layer PN Positive-Negative PWM Pulse-Width Modulation RAM Random-Access Memory SDMA Spatial-Division Multiple Access SIMT Single Instruction, Multiple Threads SLM Spatial Light Modulator SM Streaming Multiprocessors TFM Trabajo Fin de Máster ULL Universidad de La Laguna VLC Visible Light Communications WDM Wavelength-Division Multiplexing WSN Wireless Sensor Networks Parte I Memoria 1 Capítulo 1 Introducción En el presente capítulo se expone la motivación del trabajo y se introduce una visión general del problema así como los objetivos considerados. Por último, se expone la distribución de los contenidos de esta memoria. 1.1. Motivación Para entender cuál es la motivación de este trabajo, es necesario citar al grupo CÁmara FAse DIStancia (CAFADIS) de la Universidad de La Laguna (ULL). En colaboración con el Instituto Astrofísico de Canarias (IAC), se ha desarrollado un prototipo que detecta el frente de onda óptico y estima la distancia [1]. Dicho prototipo consiste en un conjunto de lentes y microlentes que puede adosarse a cualquier cámara convencional, desde cámaras réflex digitales hasta cámaras de cine. Entre sus aplicaciones, en el año 2010 desarrollaron un algoritmo de 3DTV en tiempo real que alcanza hasta 200 puntos de vista en 24 planos focales distintos. En astronomía, su aplicación directa pasa por corregir las aberraciones producidas por las distintas capas de la atmósfera, produciendo imágenes de mayor calidad. Además, han implementado modelos tanto en software (en GPU, Graphics Processing Unit) como en hardware (en FPGA, Field Programmable Gate Array). Desde el punto de vista de las comunicaciones, fruto de la colaboración entre el IDeTIC (Instituto para el Desarrollo Tecnológico y la Innovación en Comunicaciones) y el CAFADIS, en [2] se propone el uso de cámaras plenópticas en redes de sensores ópticas no guiadas. De forma inherente, emplear una cámara convencional como receptor introduce el concepto de diversidad espacial en el receptor, sin embargo, en este trabajo se plantea la utilización de plenópticas en escenarios con alta densidad de sensores donde se pueden dar casos de problemas near-far. Actualmente, dicha línea de investigación sigue activa a la espera de obtener resultados. A raíz de esta última colaboración, surgió la posibilidad de emplear un array de mililentes en el lado transmisor con el fin de introducir técnicas de Spatial-Division Multiple Access (SDMA) en redes de sensores ópticas inalámbricas. Esta configuración 3 41.2. Descripción del problema es precisamente la que se emplea en imagen integral donde, al situar un array de mililentes sobre una pantalla convencional, se consigue reproducir una imagen tridimensional sin emplear gafas u otros dispositivos. A este tipo de sistemas que no hacen uso de elementos externos se los conoce como sistemas autoesteroscópicos, en los que se puede recrear efectos de profundidad y/o paralaje. No obstante, en nuestro caso, y centrándonos en las redes de sensores ópticas inalámbricas, no existe la necesidad de formar imágenes puesto que los canales son modulados en intensidad con detección directa (IM/DD, Intensity Modulation with Direct Detection). Por lo tanto, se trata de un subproblema de imagen integral adaptado a comunicaciones. Es aquí donde comienza este trabajo, intentando modelar dicho problema para adaptarlo a nuestras necesidades. 1.2. Descripción del problema El título de este Trabajo Fin de Máster (TFM) hace alusión a la paralelización de un algoritmo, en concreto, al empleo de trazado de rayos para arrays de mililentes. Sin embargo, una de las primeras tareas a realizar para modelar el problema fue describir el escenario y los actores principales del mismo. Para ello, se parte del escenario mostrado en la Figura 1.1 donde se encuentran tres elementos: un display, un array de mililentes y un plano imagen. Display Plano imagenArray de mililentes z1z2 (ξ,η) (x,y) (u,v) Figura 1.1: Escenario del problema Por un lado, el display es el elemento emisor de información que se ha decidido modelar como un array de ledes; por otro lado, las mililentes se ha decidido que 1. Introducción 5 aumenten un orden de magnitud respecto a las microlentes empleadas en una cámara plenóptica (de cientos de micras a unidades de milímetros, esto se discutirá en la Sección 3.1); y por último, el plano imagen que aunque no se corresponde con un elemento físico nos permite obtener resultados. No se ha encontrado ningún trabajo en la literatura que aborde el estudio de este escenario con fines de comunicaciones. Sí se ha encontrado trabajos relacionados con la iluminación que sintetizan patrones sobre un array de ledes como hizo Moreno en [3] o que homogeneízan el patrón de radiación de un led por medio de arrays de mililentes tal como recoge Schreiber [4]. Por consiguiente, antes del planteamiento de este TFM, se pensó en obtener de forma analítica las ecuaciones que rigen el escenario presentado en la Figura 1.1 con el fin de obtener un conjunto de funciones base con las que realizar un proceso de síntesis. Pronto se comprobó que dicha tarea entrañaba una dificultad excesiva y se replanteó el problema para obtener mediante simulación el funcional descrito en la Ecuación 1.1. S(R(θ, φ)) W/sr (1.1) donde Srepresenta la intensidad radiante (W/sr) de una mililente en función de la intensidad radiante de un led, R, que, a su vez, depende de las coordenadas esféricas θyφ. Si planteamos la ecuación que describe el escenario completo de la Figura 1.1, se obtiene la Ecuación 1.2. T=X iX j Si,j X uX v αu,vRu,v (θi,j, φi,j)!W/sr (1.2) en este caso, se ha generalizado la Ecuación 1.1 donde la intensidad radiante de una mililente viene conformada por la intensdiad radiante de muchas fuentes moduladas en intensidad por el parámetro α. Es por ello que para obtener dicho funcional se pensó en un simulador que recogiese un número significativo de casos, para posteriormente realizar un ajuste de funciones para un escenario dado. Planteando grosso modo el sistema a optimizar, este consiste en: un primer bloque de simulación, una obtención del funcional para ese escenario y una extracción de parámetros. El criterio de parada para la síntesis es un determinado patrón de radiación a una distancia dada, de ahí el plano imagen. El flujo del sistema se presenta en la Figura 1.2. Carece de sentido modelar el sistema desde el punto de vista de la óptica física puesto que no pueden sucederse fenómenos de difracción ni de interferencia debido a la diferencia relativa de dimensiones entre la longitud de onda (nanómetros) y el sistema óptico (milímetros). Además, el uso de fuentes no coherentes como los ledes hace que sólo se pueda trabajar en intensidad. La opción lógica pasa por modelar de forma geométrica. 61.3. Objetivos Simulación Parámetros Figura 1.2: Flujo del sistema 1.3. Objetivos Este TFM tiene como objetivo principal la implementación utilizando las arquitecturas paralelas de memoria compartida y de tipo GPU de una aplicación de comunicaciones ópticas basada en técnicas de ray tracing. A su vez, este objetivo se desglosa en los siguientes: Estudiar la aplicación de las diferentes arquitecturas de paralelización sobre el problema de ray tracing. Análisis y diseño de la aplicación de ray tracing en las arquitecturas paralelas de memoria compartida y de tipo GPU. Implementación del código en las arquitecturas paralelas de memoria compartida y de tipo GPU. Evaluación de las prestaciones de la aplicación de ray tracing en las arquitecturas paralelas elegidas. 1.4. Organización de la memoria El presente trabajo se ha dividido de la siguiente manera: 1. Introducción 7 Parte IMemoria Se presenta la memoria del trabajo, compuesta de seis capítulos detallados a continuación: •Capítulo 1Introducción Se describen aspectos generales del problema a abordar, los objetivos del presente TFM, así como la estructura del documento. •Capítulo 2Antecedentes Se realiza una descripción del estado del arte de las distintas disciplinas involucradas. •Capítulo 3Modelo matemático Se define el modelo matemático seguido así como el algoritmo resultante. •Capítulo 4Implementación Se explica el proceso seguido para implementar el código en las distintas arquitecturas. •Capítulo 5Resultados Se presenta y analiza los resultados obtenidos en distintas fases del trabajo. •Capítulo 6Conclusiones Se detalla las conclusiones obtenidas a raíz de los resultados conseguidos. Parte II Bibliografía Se indica la bibliografía consultada para la realización del TFM. Capítulo 2 Antecedentes En este capítulo se esboza el estado del arte de aquellas disciplinas que entroncan con la realización de este trabajo. En concreto, se abordan temas de comunicaciones por luz visible, visores tridimensionales y tecnologías de paralelización. 2.1. Comunicaciones por luz visible Dentro de las comunicaciones ópticas no guiadas, el uso del canal VLC (Visible Light Communications) como enlace de bajada no sólo ofrece un amplio ancho de banda sin regular y una alternativa a las bandas 2,4 GHz ya muy saturadas, sino que además usa una infraestructura ya existente en todas las viviendas u oficinas. A consecuencia de esto se ha producido la eclosión de distintas iniciativas normativas como el estándar IEEE (Institute of Electrical and Electronics Engineers) 802.15.7 [5]. Su PHY (Physical Layer) III denominada CSK (Color-Shift Keying) es similar a una FSK (Frequency-Shift Keying) aunque en longitud de onda, no confundir con la multiplexación por longitud de onda (WDM, Wavelength-Division Multiplexing) donde se empelan canales independientes. En concreto, la CSK se basa en la codificación de la información en las amplitudes relativas entre las portadoras ópticas [6]. Aunque en el estándar se presentan una serie de reglas para la generación de constelaciones, estas pueden cambiarse empleando algoritmos de optimización como en [7] [8] o para añadir capacidades multiusuario directamente en la capa física [9]. Aparte de las técnicas de codificación propuestas por el estándar, durante los últimos años multitud de grupos de investigación han estudiado el uso de otras modulaciones para implementar sistemas VLC [10]. Por ejemplo, se han explorado el uso de técnicas OFDM (Orthogonal Frequency-Division Multiplexing) con el fin de lograr sistemas de comunicaciones con una eficiencia espectral alta y elevada robustez frente a multipropagación [11]. Para poder emplear estas señales bipolares en el dominio óptico se recurre comúnmente a dos técnicas bien diferenciadas: DCO-OFDM (DC-Biased Optical OFDM) y ACO-OFDM (Asymmetrically Clipped Optical OFDM) [12]. La primera de ellas se basa en introducir un nivel de continua para polarizar al led en la región lineal, mientras que la segunda opción elimina las componentes pares de la IFFT (Inverse FFT), reduciendo a la mitad la eficiencia espectral de la modulación 9 16 2.3. Tecnologías de paralelización Figura 2.5: Arquitectura del entorno CUDA R  Figura 2.6: Arquitectura Nvidia Kepler dos extremos es necesario haber reservado memoria tanto en origen como en destino mediante las funciones: 1. Host. 2. Antecedentes 17 malloc: asigna el número especificado de bytes. calloc: como malloc, además de inicializar a cero. cudaMallocHost: reserva la memoria (pinned memory) directamente sobre la RAM (Random-Access Memory) de forma que no requiere del host para realizar la transferencia por DMA (ver Figura 2.7). Figura 2.7: Transferencias de memoria entre el host y el device 2. Device. cudaMalloc: reserva memoria en la GPU. Una vez se ha inicializado la memoria correctamente en la GPU, se invoca al kernel que no es más que una función en C que se ejecuta tantas veces como hilos se hayan reservado en la tarjeta gráfica, en lugar de una única vez como en el procesador. Para indicar que la función definida en el código es un kernel, se emplea el atributo __global__, mientras que para llamarlo la notación es la siguiente: kernel<<<numBlocks, threadsPerBlock>>> La jerarquía existente en la arquitectura de la GPU, tiene su análogo a nivel de CUDA R con los bloques y los hilos como se muestra en la Figura 2.8. Cada vez que se comienza un kernel, es necesario especificar el número de bloques y el número de hilos por bloque. En las arquitecturas actuales ambos números pueden tomar la forma de ternas de hasta 3 elementos. Un kernel no puede, a su vez, llamar a funciones declaradas en el espacio de host. Para ello existe otro atributo, __device__. En el caso de que una función deba ser llamada por ambos extremos se debe emplear dos atributos, __device__ y__host__. Los SMX ejecutan los hilos de forma SIMT (Single Instruction, Multiple Threads), donde todos los cores de un mismo grupo (warp, agrupación de 32 hilos) ejecutan la misma operación a la vez. En el caso de que existan hilos con saltos condicionales o que terminen antes que otros, algunos núcleos son desactivados, redundando de forma negativa en el rendimiento [22]. 18 2.3. Tecnologías de paralelización Grid Block (1, 1) Thread (0, 0) Thread (1, 0) Thread (2, 0) Thread (3, 0) Thread (0, 1) Thread (1, 1) Thread (2, 1 ) Thread (3, 1) Thread (0, 2) Thread (1, 2) Thread (2, 2 ) Thread (3, 2) Block (2, 1) Block (1, 1) Block (0, 1) Block (2, 0) Block (1, 0) Block (0, 0) Figura 2.8: Jerarquía de bloques e hilos Capítulo 3 Modelo matemático Este capítulo describe el algoritmo seguido para implementar el simulador objeto del TFM. El modelo matemático parte de un preprocesado que permite aumentar la densidad de rayos útiles para el cálculo del funcional presentado en la Ecuación 1.1. A partir de ahí, el modelado del resto de componentes del escenario sigue los cauces habituales. 3.1. Justificación del uso de mililentes En [2] la cámara empleada fue una RED ONE R [23] (ver Figura 3.1) con una resolución de 5120 (h) ×2700 (v) píxeles, donde el tamaño del sensor MysteriumTM empleado es de 27,7 mm ×14,6 mm. Esto unido a que el pitch del array de microlentes era de tan solo 130 µm, permitía que bajo cada microlente hubiese aproximadamente 24 píxeles por dimensión. Si utilizáramos el mismo array con el fin de transmitir, nos encontraríamos con el problema de que los ledes son un orden de magnitud mayores a los fotodiodos. Por ejemplo, si tomamos como referencia una de las mejores pantallas OLED (Organic Light-Emitting Diode) existentes en el mercado, como la que monta el Samsung Galaxy S6, nos encontramos que, con una densidad de 577 PPI, el tamaño de píxel es de 44 µm. Empleando el mismo array únicamente contaríamos con 3 píxeles por dimensión lo cual es insuficiente a todas luces. Por este motivo, se decidió emplear mililentes con un pitch de 1,5 mm para poder alojar un número mayor de elementos. 3.2. Óptica geométrica Como se avanzó en la Sección 1.2, el modelado físico del problema se corresponde con un problema de trazado de rayos. Esto significa que la propagación de la luz obedece a los fenómenos de reflexión y refracción. En la Figura 3.2 se muestra la interfaz entre dos materiales con índices de refracción η1yη2. La dirección del rayo incidente viene dada por el vector unitario, i; de la misma forma, las direcciones de los rayos reflejado y refractado son ryt, respectivamente. 19 20 3.2. Óptica geométrica Figura 3.1: RED ONE R  Figura 3.2: Ley de Snell Igualmente, se define el vector normal, n, perpendicular a la interfaz de los dos materiales. Por lo tanto, todos los vectores cumplen la condición de la Ecuación 3.1. ||i|| =||r|| =||t|| =||n|| = 1 (3.1) Cada vector, a su vez, puede descomponerse en sus componentes normal y tangencial. Se emplea la notación v⊥para denotar la componente normal del vector vy lo propio con vkpara la componente tangencial. La componente normal de un vector se obtiene a partir de la proyección de este sobre el vector normal, tal y como se muestra en la Ecuación 3.2. v⊥= (v·n)n(3.2) Por lo que la componente tangencial puede obtenerse restando al vector original la componente normal, como se muestra en la Ecuación 3.3. vk=v−v⊥(3.3) 3. Modelo matemático 21 Por definición, el producto escalar entre ambas componentes es nulo, lo que quiere decir que ambas componentes son ortogonales entre sí; por consiguiente, se cumple la condición de la Ecuación 3.4. ||v||2=||vk||2+||v⊥||2(3.4) Los ángulos de incidencia, reflexión y refracción se denominan θi,θryθt; y se definen siempre como el menor ángulo positivo entre la dirección del rayo y el vector normal. Mediante trigonometría se obtiene las relaciones de la Ecuación 3.5. cos θv=||v⊥|| =±v·n sin θv=||vk|| (3.5) 3.2.1. Reflexión De los dos fenómenos, reflexión y refracción, la reflexión es el más sencillo y modela el choque mecánico de un rayo sobre una superficie reflectora. La ley de la reflexión nos dice que el ángulo de incidencia es igual al ángulo de reflexión: θr=θi. Si se desarrolla dicha igualdad como en la Ecuación 3.6. ||r⊥|| = cos θr= cos θi=||i⊥|| ||rk|| = sin θr= sin θi=||ik|| (3.6) Por simple inspección de la Figura 3.2 se puede deducir la Ecuación 3.7. r⊥=−i⊥ rk=ik(3.7) Finalmente, la dirección del rayo reflejado dada por el vector rno será más que la Ecuación 3.8. r=ik−r⊥ = [i−(i·n)n]−(i·n)n =i−2(i·n)n(3.8) 3.2.2. Refracción La refracción se basa, casi en su totalidad, en la ley de Snell presentada en la Ecuación 3.9. Esta nos dice que el producto de los índices de refracción y los senos de los ángulos de ambos medios debe ser igual. 22 3.2. Óptica geométrica η1sin θi=η2sin θt⇒sin θt=η1 η2 sin θi(3.9) De la Ecuación 3.9 se deduce inmediatamente que cuando el sin θi>η2 η1, el sin θt debería ser mayor a la unidad; lo cual es imposible. En ese caso, se produce lo que se conoce como reflexión total interna. Al igual que con la reflexión, se descompone el rayo refractado en sus partes tangencial y normal. Aplicando las Ecuaciones 3.5 y3.9, obtenemos la Ecuación 3.10. ||tk|| =η1 η2 ||ik|| (3.10) Las componentes tangenciales de ambos rayos son paralelas por lo que si al rayo incidente le sustraemos la parte normal (Ecuaciones 3.3), nos queda la Ecuación 3.11. tk=η1 η2 [i+ cos θin](3.11) Teniendo en cuenta que siempre se está trabajando con vectores unitarios y aplicando el teorema de Pitágoras, la componente normal se obtiene en la Ecuación 3.12. t⊥=−q1− ||tk||2n(3.12) Uniendo las Ecuaciones 3.3 y3.12, además de aprovechar las relaciones entre los ángulos y los rayos; se puede omitir el cálculo de funciones trigonométricas de forma explícita, según puede verse en la Ecuación 3.13. t=ri+rc −p1−r2(1 −c2)n(3.13) donde r=η1 η2yc=−n·i. 3.2.3. Ángulo crítico En la Subsección anterior se nombró la reflexión total interna como un fenómeno que provoca que no exista rayo refractado, por lo tanto, no hay transmisión de potencia al segundo medio. La condición para que esto no ocurra está implícita en el cálculo del rayo refractado. Si el radicando de la Ecuación 3.13 es negativo, significa que el rayo no se propaga. Se obvia el cálculo del ángulo crítico para no generar rayos que no se transmitan porque en la Sección 3.5 se combina con otra condición propia del problema. 3. Modelo matemático 23 3.2.4. Ecuaciones de Fresnel De toda la luz que llega a una interfaz entre dos materiales no absorbentes: una parte se refleja de nuevo hacia el medio del que proviene, mientras que el resto se transmite hacia el segundo medio. De forma matemática esto se traduce en la Ecuación 3.14. T+R= 1 (3.14) donde TyRsignifican transmitancia y reflectividad, respectivamente. La cantidad de luz reflejada o transmitida depende de los índices de refracción y del ángulo de incidencia. Las ecuaciones de Fresnel describen las amplitudes de las ondas reflejada y refractada en función de la amplitud de la onda incidente. En general, la reflectividad para luz polarizada se presenta en la Ecuación 3.15. R⊥(θi) = η1cos θi−η2cos θt η1cos θi+η2cos θt2 Rk(θi) = η2cos θi−η1cos θt η2cos θi+η1cos θt2 (3.15) No obstante, al trabajar con luz no polarizada simplemente se promedia la reflectividad para ambas polarizaciones. Con todo, las expresiones finales para la reflectividad y la transmitancia pueden observarse en la Ecuación 3.16. R(θi) = (R⊥(θi)+Rk(θi) 2si no hay reflexión total interna 1si hay reflexión total interna T(θi) = 1 −R(θi)(3.16) 3.3. Fuentes ópticas Las principales fuentes empleadas en los sistemas de comunicaciones ópticas son el diodo láser y el led. Su estructura básica, en ambos casos, es la heterounión o lo que es lo mismo, la unión de dos semiconductores con energías de gap distintas. La región de emisión es una unión PN (Positive-Negative) de semiconductores III–V de gap directo que al ser polarizada en directa provoca que los portadores mayoritarios se difundan y recombinen, emitiendo energía en forma de luz (proceso radiativo) o disipándose en forma de calor (proceso no radiativo). La principal diferencia entre los diodos láser y los ledes es la coherencia de la luz emitida. La radiación producida por un diodo láser se genera en una cavidad resonante, lo que le confiere una alta coherencia tanto espacial como temporal. Esto se traduce en una gran monocromaticidad y una alta directividad, siendo especialmente indicados para los enlaces punto a punto. Por contra, en un led no hay tal cavidad y la radiación 24 3.3. Fuentes ópticas resultante tiene una anchura espectral considerable y además no coherente. La mayoría de las fuentes ópticas del mercado presentan diagramas de radiación lambertianos del estilo de la curva de la Figura 3.3 y responden a la Ecuación 3.17. 0 1 R(θ, m) Figura 3.3: Diagrama de radiación lambertiano para m= 20 R(θ, m) = m+ 1 2πPTcosm(θ) W/sr (3.17) En este trabajo únicamente se ha considerado fuentes tipo led monocromáticas. Debido a que se pretende obtener el funcional descrito en la Ecuación 1.1, se necesita generar rayos en todo el ángulo sólido radiado por el emisor. Es por ello que se ha decidido mallar el espacio de forma regular y determinista para obtener información espacial relevante al problema. La dirección de salida de un rayo puede ser determinada por dos ángulos definidos en coordenadas esféricas: θyφ. En un principio, se pensó barrer todo el ángulo sólido, sin embargo, este hecho incurre en una gran cantidad de rayos que sufren de reflexión total interna en la interfaz de salida de la lente; por lo tanto, en la Sección 3.5 se introduce el método optimizado de generación de rayos. En cualquier caso, un led queda definido por lo siguientes parámetros: Potencia total de emisión: PT. Número de niveles en elevación: 2T. Número de niveles en azimut: 2P. Índice del patrón lambertiano: m. 3. Modelo matemático 25 3.4. Lentes Una lente es cualquier objeto que sea capaz de desviar los rayos de luz mediante refracción. Además, se puede usar para enfocar la luz mientras que un prisma, a pesar de que también refracta la luz, no la enfoca. Los arrays de lentes vistos, tanto en la literatura como en el mercado, para aplicaciones relacionadas con imagen integral son mayoritariamente planoconvexos, donde una superficie de la lente es plana y la otra sobresale hacia afuera como la presentada en la Figura 3.4. La diferencia entre los distintos arrays viene dada por el perfil utilizado: cilíndrico, esférico, hexagonal, etc. Figura 3.4: Lente planoconvexa A la hora de modelarse el array de mililentes, se ha hecho uso de un modelo real como el MLA150-5C (perfil esférico) de Thorlabs [24] que es muy similar al empleado en el grupo CAFADIS. En particular, este modelo cuenta con una máscara que imposibilita el paso de la luz si no es a través de la mililente, aumentando así el contraste (ver Figura 3.5). Este hecho será clave en la Sección 3.5. Los parámetros necesarios en la definición de una lente se muestran en la Figura 3.6. Figura 3.5: Imagen captada con un array de microlentes encima del sensor Material (índice de refracción): n. Perfil convexo: esférico. 32 4.2. Estructuras de datos Generar rayo Snell y Fresnel Snell y Fresnel Impactar rayo Rayos generados Rayos refractados Rayos impactados Propagación Figura 4.1: Diagrama del algorimo Ray Campo Significado Tipo de dato Posición Coordenadas cartesianas de la posición del rayo. float/double Vector director Vector unitario que define la dirección de propagación del rayo. float/double Distancia Distancia recorrida desde la fuente. float/double Potencia Potencia por unidad de ángulo sólido. float/double Etiqueta Identifica la procedencia de cada rayo mediante su fuente de emisión. unsigned int Tabla 4.1: Estructura de datos para un rayo Mediante esta abstracción en el código, la información de los ledes se traslada a los rayos durante la inicialización de los mismos. A partir de ahí, la única alusión a los ledes estará en la etiqueta que forma parte de los campos del registro Ray. En el caso de las mililentes, a pesar de que se utiliza una estructura similar, registros, se empleó una división jerárquica con vistas a flexibilizar la ampliación del simulador a otros tipos de lentes, distintas a las plano-convexas o modificando el perfil. El registro Lens (ver Tabla 4.2) consta de dos campos que son otros dos registros del tipo Interface:Left yRight. A su vez, la estructura Interface posee los campos de la Tabla 4.3. Lens Campo Significado Tipo de dato Left Primera interfaz de la lente (medio-lente). struct Interface Right Segunda interfaz de la lente (lente-medio). struct Interface Tabla 4.2: Estructura de datos para una lente En este TFM, al emplear un modelo concreto de mililente, únicamente se han implementado los tipos plana y esférica. En estos casos, el parámetro indica el radio de la lente que puede ser infinito (interfaz plana), positivo (interfaz convexa) o negativo 4. Implementación del sistema 33 Interface Campo Significado Tipo de dato Posición Coordenadas cartesianas del centro óptico. float/double Tipo Especifica el perfil de la lente. unsigned int Parámetro En función del tipo de interfaz tiene un significado diferente. float/double Índice de refracción Material de fabricación. float/double Tabla 4.3: Estructura de datos para una interfaz (interfaz cóncava). Si se añadiesen perfiles ad-hoc, el parámetro podría tomar el significado que se desease. Para finalizar, un escenario de simulación se define a través de una serie de parámetros, los cuales son: Distancia a las mililentes y al plano de impacto. Índices de refracción de los medios. Número de ledes y mililentes. Separación, diámetro y espesor de las mililentes. Tamaño de la protuberación (cúpula) y radio de curvatura de las mililentes. Índice del emisor lambertiano y potencia total de un emisor. Número de niveles en azimut y elevación por pareja de led-mililente. 4.3. Explotación de la simetría Antes de abordar la descripción del algoritmo, se hace necesario nombrar la estrategia seguida para simular. A pesar de que no se ha nombrado en ningún momento cuántos ledes forman un display o cuántas mililentes hay por pulgada, el número de elementos se dispara a poco que se introduzcan valores comerciales. Asimismo, si se supone que los ledes y las mililentes están alineados, es decir, una fila de ledes es paralela a otra de mililentes; se puede comprobar que el problema es básicamente repetitivo. Por este motivo, se replanteó la idea inicial de simular toda la pantalla y, en su lugar, tomar sólo una parte representativa de la misma de forma que, posteriormente, por superposición pueda recomponerse el patrón de radiación completo. Para ello se planteó agrupar las mililentes en clústeres según su conectividad (vecindad), empleando 8-conectividad. De esta forma los clústeres pueden contener 1, 9, ..., (2·n+1)2mililentes, siendo n∈Nun parámetro de simulación. El tamaño óptimo del clúster para que la simulación sea precisa se podría determinar mediante un criterio geométrico, en función de la potencia emitida en los ángulos sólidos subtendidos bajo las mililentes, no obstante, en el presente TFM, se hace un pequeño estudio numérico de la precisión en la Sección 5.2. De igual forma, los ledes pueden agruparse sobre la mililente central del clúster, cubriendo todos los casos posibles. En este caso, el número de ledes sigue potencias pares de dos según 22m, con m∈N. En la Figura 4.2 se presenta la idea de agrupar las mililentes y los ledes para agilizar el cálculo. 34 4.4. Descripción del algoritmo −1,5 0 1,5 −1,5 0 1,5 Geometr´ıa Distancia (mm) Distancia (mm) Figura 4.2: Clúster de 9 mililentes y 4 ledes por mililente Siguiendo este razonamiento, se podría seguir explotando la simetría del problema y calcular un único tramo 0,π 4sobre la dimensión φy aplicar simetría de revolución para cubrir todos los ángulos. No obstante, se descartó llevar hasta el límite la aplicación de simetrías puesto que para síntesis no aporta ninguna ventaja. En cambio, manejar un clúster puede resultar conveniente llegado el caso. 4.4. Descripción del algoritmo En el Algoritmo 1se introduce el pseudo-código del simulador completo. De forma general, el código puede dividirse en las siguientes subsecciones: definición del escenario, cálculo de límites y simulación. 4.4.1. Definición del escenario En primer lugar, se leen los parámetros de la simulación a través de la función ReadConfiguration. A partir de los parámetros de entrada, se inicializan las estructuras Lens eInterface mediante las funciones InitGeometry eInitScenario. InitGeometry calcula los centros de las mililentes y los ledes en función del clúster y la densidad de fuentes deseada, mientras que InitScenario solicita la memoria para las estructuras que alojan la información de las mililentes y las inicializa. 4. Implementación del sistema 35 Algorithm 1 Simulador 1: procedure lenslet 2: ReadConfiguration (file in,params) 3: leds ←22m 4: lenses ←2·n+ 1 5: InitGeometry (ledsCenters,lensesCenters) 6: InitScenario (array) 7: RaysPropagation (rays) 8: WriteResults (file out,rays) 9: end procedure 10: procedure RaysPropagation(rays) 11: GetPhiLimits (φm´ın,φm´ax) 12: P←2p 13: T←2t 14: for i←1, leds do 15: for j←1, lenses do 16: dφ←φm´ax(i,j)−φm´ın(i,j) P−1 17: for k←1, P do 18: φ←φm´ın(i, j) + k·dφ 19: GetSecondOrderSolutions (sm´ın,sm´ax) 20: GetThetaFromDistance (θm´ın,θm´ax) 21: dθ←θm´ax−θm´ın T−1 22: for l←1, T do 23: θ←θm´ın +l·dθ 24: GenerateRay (rays(i, j, k, l),m,φ,θ,PT) 25: CalculateOutput (rays(i, j, k, l),array(j)) 26: CalculateImpact (rays(i, j, k, l),distance) 27: end for 28: end for 29: end for 30: end for 31: end procedure 32: procedure CalculateOutput(ray,lens) 33: CalculatePlane (ray,left) 34: SnellFresnel (ray,left) 35: CalculateSphere (ray,right) 36: SnellFresnel (ray,right) 37: end procedure 36 4.4. Descripción del algoritmo 4.4.2. Cálculo de límites Tras realizar las acciones anteriores, donde únicamente se han definido las condiciones de contorno del escenario de simulación, se procede al preprocesado que se introdujo en la Sección 3.5. Tanto el cálculo de límites como la simulación están dentro de una rutina mayor denominada RaysPropagation. Primero, para cada par led-mililente se calculan los ángulos azimutales límites mediante la rutina GetPhiLimits. Según se barre los ángulos φi, se calculan los ángulos de elevación máximo y mínimo en dos pasos con las rutinas GetSecondOrderSolutions y GetThetaFromDistance. Con esto se determinan los ángulos sólidos útiles desde el punto de vista de la simulación. En la Figura 4.3 se presenta un ejemplo para un clúster de 9 mililentes y un único led. −0.6 −0.4 0.4 0.5 0.6 0.7 ´ Angul o s´oli do entrante 1 -0.2 0 0.2 0.5 0.6 0.7 0.8 ´ Angul o s´oli do entrante 2 0.4 0.6 0.4 0.5 0.6 0.7 ´ Angul o s´oli do entrante 3 -0.8 -0.6 -0.2 0 0.2 ´ Angul o s´oli do entrante 4 -0.2 0 0.2 -0.2 0 0.2 ´ Angul o s´oli do entrante 5 0.6 0.8 -0.2 0 0.2 ´ Angulo s´olido e ntrante 6 -0.6 -0.4 -0.7 -0.6 -0.5 -0.4 ´ Angul o s´oli do entrante 7 -0.2 0 0.2 -0.8 -0.7 -0.6 -0.5 ´ Angul o s´oli do entrante 8 0.4 0.6 -0.7 -0.6 -0.5 -0.4 ´ Angul o s´oli do entrante 9 Figura 4.3: Malla de ángulos sólidos a simular sobre el plano XY 4.4.3. Simulación Una vez la dirección de cada rayo queda prefijada con los pasos anteriores, se pasa a generar cada uno siguiendo el patrón de radiación lambertiano (modificable a través del parámetro m). A pesar de que la potencia de cada led también es un parámetro de la simulación, el preprocesado introducido fuerza a que no se genere todo el diagrama de radiación completo y con ello, si se realizase la integral sobre los rayos generados, no se obtendría la potencia introducida. La diferencia entre el valor introducido y el de la integral se corresponde con las pérdidas por absorción en la máscara de las mililentes, mientras que las pérdidas por reflexión total interna se calculan mediante las ecuaciones de Fresnel. Si se cambiase el tipo de mililentes a uno sin máscara, bastaría con desechar el preprocesado y generar todos los rayos. El proceso siguiente se realiza dentro de un 4. Implementación del sistema 37 mismo bucle, excepto los rayos directos de cada fuente que se calculan en un bucle externo. La simulación comprende los siguientes pasos: Generación de los rayos en las direcciones calculadas. Rutina GenerateRay. Cálculo del camino óptico hasta la salida de las mililentes a través de la rutina CalculateOutput. Esta, a su vez, se subdivide en: •Primero, impacto con la interfaz plana de las mililentes. Rutina CalculatePlane. •Seguidamente, refracción del medio original a las mililentes. Rutinas SnellFresnel oSnellSchlick. •Posteriormente, cálculo del impacto con la cúpula de las mililentes mediante la rutina CalculateSphere. •Por último, segunda refracción con las mismas rutinas SnellFresnel o SnellSchlick. Cálculo del impacto con el plano imagen. Rutina CalculateImpact. Los resultados obtenidos se almacenan en un fichero externo para ser procesados en software externo al TFM con la rutina WriteResults. Un ejemplo de los rayos almacenados puede verse en la Figura 4.4 donde, para un clúster de 9 mililentes y 4 ledes, pueden verse los impactos sobre un plano situado a 3 metros del origen. −10 −5 0 5 10 −8 −6 −4 −2 0 2 4 6 8 Impactos Distancia (m) Distancia (m) Figura 4.4: Impactos sobre el plano imagen 38 4.4. Descripción del algoritmo Rutinas de cálculo de impactos En el escenario planteado para simular pueden ocurrir dos tipos de impactos: con un plano o con una esfera. El primero se da al impactar con la lente plano-convexa por su interfaz plana, así como con el plano imagen (rutina CalculatePlane). La segunda al salir de la mililente por la cúpula que no está cubierta por la máscara (rutina CalculateSphere). En ambos casos, se calcula el punto de impacto del rayo en coordenadas cartesianas y se actualiza tanto el punto de aplicación del rayo como la distancia recorrida. Estas distancias se calculan de la siguiente forma según las Ecuaciones 4.1 y4.2. dplane =zplane −zray vz (4.1) dsphere =η±pη2+ξ(4.2) En la Ecuación 4.1 se puede prescindir del valor absoluto puesto que la propagación de los rayos siempre se realiza sobre el eje Zpositivo. No obstante, la Ecuación 4.2 es exactamente la misma que la Ecuación 3.18 puesto que el problema del corte de una recta con un círculo es equivalente al corte de una recta con una esfera, añadiendo una dimensión extra. En este caso, no interesa obtener las dos soluciones que nos devuelve el radicando sino la que sea mayor que cero. El significado físico de esto es que la distancia positiva es el corte con la cúpula, mientras que la distancia negativa, aunque es también un corte con la esfera, no existe en la realidad. Esto es así porque el rayo está dentro de la esfera. En caso contrario, ambas soluciones serían positivas. Rutinas de cálculo de refracción Con vistas a reutilizar cálculos se decidió unir el cálculo del rayo refractado con las ecuaciones de Fresnel, de este modo se evita calcular las funciones trigonométricas de forma reiterada (rutina SnellFresnel). Como la refracción depende de la normal en el punto de impacto, se debe diferenciar cuándo se impacta contra un plano o contra una esfera. El cálculo de normales se realiza como en las Ecuaciones 4.3 y4.4. nplane =ez(4.3) nsphere =xray −lensC ||xray −lensC|| (4.4) A partir de aquí, se sigue con las ecuaciones presentadas en el Capítulo 3para el rayo refractado (Ecuación 3.13) y las pérdidas de Fresnel (Ecuación 3.16). No obstante, en la literatura se encontró un método optimizado de cálculo de las ecuaciones de Fresnel: la aproximación de Schlick [25]. En las Ecuaciones 4.5 y4.6 se resume sus expresiones. 4. Implementación del sistema 39 RSchlick (θi) =      R0+ (1 −R0) (1 −cos θi)5η1≤η2 R0+ (1 −R0) (1 −cos θt)5η1> η2si no hay reflexión total interna 1η1> η2si hay reflexión total interna (4.5) con R0: R0=η1−η2 η1+η22 (4.6) En [26] se asegura que el cálculo de las pérdidas es hasta un 30 % más rápido en comparación a la ecuación de Fresnel para luz no polarizada si se evita el uso de funciones de exponenciación del tipo pow, en caso contrario, se vuelve el doble de lenta. Esto ocurre cuando se compara el cálculo aislado de ecuaciones de Fresnel, sin embargo, en nuestro caso se han integrado en una función que engloba tanto el cálculo de la refracción como el de las pérdidas, por lo que esta aceleración ya no resulta tan evidente. 4.5. Código secuencial La situación de partida del TFM fue un prototipo secuencial desarrollado en MATLAB R . A partir de este prototipo, se reescribió en C el código secuencial. De partida, sabiendo que habría que portarlo en última instancia a CUDA R , se intentó, en la medida de lo posible, eliminar las posibles bifurcaciones en el código (estructuras if-then-else oswitch) para minimizar la warp divergence. Esta sucede cuando no todos los hilos de un warp ejecutan la misma instrucción, normalmente cuando hay ramas condicionales. En este caso, cada rama del código se ejecuta de forma secuencial y los hilos que no cumplen la condición se desactivan. Además, se parametrizó el código de forma que pudiese cambiarse la precisión entre simple y doble mediante directivas #define. De esta forma, la estructura Ray ocupa en memoria 36 bytes cuando se trata de float y 72 bytes en el caso de double. Además, se decidió cómo irían los rayos alojados en memoria para facilitar los accesos a la misma, esto es, que se hiciesen de forma secuencial y, en la medida de lo posible, alineados. En un primer momento, se crearon 3 arrays de rayos unidimensionales para los distintos resultados necesarios, aunque tras algunos cambios en el orden de los bucles buscando el mejor rendimiento posible se optó por crear un único array que englobase a todos. A esto hay que añadir que, dentro del preprocesado, por cada ángulo en azimut se generan varios ángulos en elevación por lo que esa jerarquía era recomendable mantenerla en memoria. Con todo, en la Figura 4.5 se muestra el formato del array de rayos. Así, de mayor a menor orden jerárquico: primero, se subdivide en fuentes (ledes); seguidamente, en lentes; posteriormente, en azimut; después, en elevación y; por último, 40 4.6. Código memoria compartida Led0Led1LedN Lente0Lente1LenteM ... ... ... ... ... G0G1GQ R0R1RQ I0I1IQ G: rayo generado R: rayo refractado I: rayo impactado Figura 4.5: Formato del vector de rayos el rayo en las tres posiciones en las que se desea obtener. De esta forma, según se recorre el array mediante un índice global se va poblando el mismo. 4.6. Código memoria compartida Una vez comprobado que los resultados obtenidos por la versión secuencial en C se correspondían con aquellos obtenidos en MATLAB R , se hizo uso de un profiler incluido en el IDE (Integrated Development Environment) Xcode para localizar los hot-spots del código. El resultado se presenta en la Tabla 4.4. Tiempo ( %) Tiempo Self Nombre 100 % 56477,1 ms 0ms Main thread 99,99 % 56477,0 ms 0ms main 99,99 % 56458,9 ms 6402 ms RaysPropagation 29,8 % 16861 ms 16861 ms SnellFresnel 28,9 % 16343,8 ms 514,1 ms CalculateOutput 14,9 % 8434,5 ms 7732,5 ms GenerateRay 9,5 % 5369,5 ms 5369,5 ms CalculatePlane 7,5 % 4269,3 ms 4269,3 ms CalculateSphere 2,5 % 1415 ms 1404,3 ms GetThetaFromDistance Tabla 4.4: Resultados del profiler En concreto, se observa una simulación consistente en 121 mililentes y 1024 leds con un nivel de muestreo de 32 niveles por coordenada angular (1024 puntos). De entre 4. Implementación del sistema 41 todos los procesos, RaysPropagation destaca sobre los demás en duración. Recordemos que dicha rutina aloja tanto el preprocesado como la simulación propiamente dicha. Sin embargo, puede comprobarse que acto y seguido las rutinas con mayor peso son SnellFresnel yCalculateOutput, mientras que las rutinas GetThetaFromDistance oGetPhiLimits apenas repercuten sobre el resultado final. Por este motivo, se decidió paralelizar el bucle de simulación presentado en la Subsección 4.4.3 que comprende las rutinas que más tiempo demandan dentro del código secuencial. Dicho bucle está compuesto, a su vez, por 4 bucles: ledes, mililentes, ángulo azimutal y ángulo de elevación. Los rayos, una vez creados, son independientes por lo que a priori la paralelización podría hacerse sobre cualquier bucle. No obstante, los dos bucles más internos, los angulares, dependen de datos calculados en los dos más externos. Por consiguiente, la decisión de paralelizar queda reducida a hacerlo por ledes o por mililentes. Debido al artefacto del clúster, el número de ledes, normalmente, es sensiblemente mayor al número de mililentes por lo que se decidió paralelizar el bucle más externo y distribuir las iteraciones entre el número de hilos. El número de hilos se decidió pasar por parámetro al código. Ya que OpenMP considera como variables compartidas todas aquellas variables existentes antes de la paralelización, cada hilo deberá poseer la información relativa a su rayo de forma privada. Esta información son los ángulos θyφ. Además, estos dependen de los diferenciales de ángulo, dθy dφ, y de las variables correspondientes al preprocesado en elevación: sm´ın,sm´ax,θm´ın yθm´ax. Aunque pueda parecer un número excesivo de variables de tipo private, hay que entender que cada hilo se ocupa de un único rayo al mismo tiempo y que estos pueden pertenecer a pares led-mililente diferentes, con lo que hay que preservar los valores geométricos. Los índices de los bucles también se han declarado como variables privadas para poder acceder correctamente al array de rayos. Con todo, el código se ha paralelizado utilizando la siguiente directiva sobre el bucle externo: #pragma omp parallel for private(i, j, k, l, diffP, phi, solMin, solMax, thetaMin, thetaMax, theta, diffT) 4.7. Código GPU La versión para CUDA R del código secuencial se basa en añadir a las rutinas secuenciales que se llamen desde dentro del bucle paralelizado en OpenMP el atributo __device__ para que se puedan ejecutar en la GPU. La única excepción es la función dot_product a la que se le ha añadido los atributos __host__ y__device__. De la misma forma, se ha declarado dos kernels: uno que engloba a todos los rayos del bucle, cudaRaysPropagation, y otro para los rayos directos, cudaDirectRays. En cuanto a la estructura de bloques e hilos a emplear, se decidió fijar el número de hilos por bloque mediante un parámetro del código y, en función de su valor, calcular el número de bloques necesario. Tanto el grid de bloques como los bloques de hilos son unidimensionales, acordes al vector de rayos que se emplea. El número de modificaciones 48 5.3. Resultados de la paralelización razón por la que se obtienen aceleraciones lineales para 2 y 4 hilos, mientras que para 8 y 16 la tasa disminuye. Destacar que para todos los hilos existe una serie de casos donde la aceleración es ligeramente inferior a la unidad (resaltados en rojo), esto es debido a que, a pesar de que el número de rayos aumenta, en el escenario, únicamente hay un led. Si recordamos de la Sección 4.6, la opción elegida fue paralelizar por fuentes por lo que no existe ninguna ventaja frente al caso secuencial. 103104105106107108 0.8 1 1.2 1.4 1.6 1.8 2 2.2 Dos hilos N´umero de rayos Aceleraci´on 103104105106107108 0.5 1 1.5 2 2.5 3 3.5 4Cuatro hilos N´umero de rayos Aceleraci´on 103104105106107108 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5Ocho hilos N´umero de rayos Aceleraci´on 103104105106107108 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5Diecis´eis hilos N´umero de rayos Aceleraci´on Figura 5.4: Aceleraciones para memoria compartida En la Figura 5.5 se presenta la eficiencia para los distintos casos, se comprueba a simple vista que la eficiencia no varía apenas a partir de 4 hilos puesto que no hay más procesadores y el tiempo es muy similar. Hay que recordar que se están empleando núcleos lógicos (software) que emulan el comportamiento de uno real y por lo tanto, son resultados subóptimos. De los resultados obtenidos, se observa que la paralelización es bastante eficiente puesto que la eficiencia para los casos de 2 y 4 es cercana a la unidad. Si se hubiese simulado en un procesador con mayor número de núcleos, se podría haber comprobado que el comportamiento debe ser el mismo. De la misma forma, se ha obviado el estudio de otras planificaciones, estáticas o dinámicas, en el código por una sencilla razón: todos los rayos ejecutan las mismas operaciones. Esto es así porque se ha empleado un muestreo determinista para cada par led-mililente y además porque un rayo a pesar de no refractarse, a efectos de cálculo, sí se tiene en cuenta por lo que se sigue operando sobre él. 5. Resultados 49 103104105106107108 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 Dos hilos N´umero de rayos Eficiencia 103104105106107108 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1Cuatro hilos N´umero de rayos Eficiencia 103104105106107108 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Ocho hilos N´umero de rayos Eficiencia 103104105106107108 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Diecis´eis hilos N´umero de rayos Eficiencia Figura 5.5: Eficiencias para memoria compartida Doble precisión Aprovechando que el código está parametrizado para trabajar tanto en coma flotante de simple y doble precisión, se repitieron las simulaciones utilizando en este caso aritmética de doble precisión. Los resultados pueden observarse en las Figuras 5.6 y5.7. En estos casos, salvo para 2 hilos, las aceleraciones no son lineales. Además, se observa el mismo comportamiento errático en 8 y 16 hilos cuando se supera el número de cores físicos disponibles en el procesador. Llegados a este punto, cabe preguntarse si el aumento de la precisión es asumible teniendo en cuenta la pérdida de rendimiento. Para ello se compara en las Tabla 5.2 el error cuadrático medio entre el prototipo MATLAB R y el código multihilo para 9 mililentes y 4 ledes. Los guarismos de error que se manejan para simple precisión están en el orden de la millonésima, mientras que para doble precisión estos disminuyen seis órdenes de magnitud aproximadamente. Teniendo en cuenta que la versión en simple precisión es más que suficiente para el objetivo propuesto y que la versión en doble precisión no mantiene el comportamiento lineal con el número de hilos, desde el punto de vista del autor, no compensa aumentar la precisión en ningún caso. 5.3.2. CUDA R  En la versión CUDA R , el número de bloques se calcula a partir del número de rayos a simular y el número de hilos por bloque, así que se ha barrido todos los escenarios 50 5.3. Resultados de la paralelización 103104105106107 1 1.5 2 2.5 Dos hilos N´umero de rayos Aceleraci´on 103104105106107 0.5 1 1.5 2 2.5 3 3.5 Cuatro hilos N´umero de rayos Aceleraci´on 103104105106107 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5Ocho hilos N´umero de rayos Aceleraci´on 103104105106107 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5Diecis´eis hilos N´umero de rayos Aceleraci´on Figura 5.6: Aceleraciones para memoria compartida en doble precisión 103104105106107 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 Dos hilos N´umero de rayos Eficiencia 103104105106107 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Cuatro hilos N´umero de rayos Eficiencia 103104105106107 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Ocho hilos N´umero de rayos Eficiencia 103104105106107 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Diecis´eis hilos N´umero de rayos Eficiencia Figura 5.7: Eficiencias para memoria compartida en doble precisión 5. Resultados 51 Error cuadrático medio en simple precisión Rayos generados Distancia 0 Potencia 2,9370e-07 Posición [5,4210e-20 5,4210e-20 0] Vector director [2,7908e-07 2,7784e-07 3,0087e-07] Rayos refractados Distancia 2,8972e-09 Potencia 3,2028e-07 Posición [2,1771e-09 2,1780e-09 2,8223e-09] Vector director [2,8213e-07 2,8156e-07 3,6580e-07] Rayos impactados Distancia 7,2630e-06 Potencia 3,2028e-07 Posición [5,0440e-06 5,0761e-06 1,5399e-16] Vector director [2,8213e-07 2,8156e-07 3,6580e-07] Error cuadrático medio en doble precisión Rayos generados Distancia 0 Potencia 5,1938e-15 Posición [5,4210e-20 5,4210e-20 0] Vector director [2,9964e-14 2,9972e-14 1,6284e-14] Rayos refractados Distancia 5,9922e-17 Potencia 3,6580e-14 Posición [7,9599e-17 7,9637e-17 2,6966e-17] Vector director [1,7801e-14 1,7802e-14 4,4949e-14] Rayos impactados Distancia 1,4574e-12 Potencia 3,6580e-14 Posición [1,0794e-12 1,0801e-12 1,5399e-16] Vector director [1,7801e-14 1,7802e-14 4,4949e-14] Tabla 5.2: Error cuadrático medio según el tipo de coma flotante empleado para un número de hilos tal que 2hcon h= 0,...,10. Para representar mejor los datos se ha decidido dividir los resultados en las Figuras 5.8 y5.9. En la Figura 5.8 se muestra los resultados desde 1 hasta 32 hilos por bloque. En todos los casos se consigue aceleraciones por encima de la unidad a partir de 4 millones de rayos. Además, el rendimiento aumenta según aumenta el número de hilos, lo que de forma implícita implica una reducción del número de bloques. Sin embargo, en la Figura 5.9 puede comprobarse como esto deja de ser cierto a partir de 64 hilos. En efecto, para 64 hilos el rendimiento disminuye, mientras que para el resto de casos vuelve a aumentar pero sin mejorar el caso de 32 hilos. A la luz de los resultados, el mejor caso se da cuando se eligen 32 hilos por bloque para todos los escenarios considerados. No obstante, las aceleraciones obtenidas son pobres en comparación con el paradigma de memoria compartida. Para buscar una explicación a este hecho, se debe tener en cuenta el paso de memoria entre la GPU y la CPU para recolectar los resultados de la simulación. Los resultados de la simulación son 3 arrays de rayos (generados, refractados e impactados) que en muchos casos tienen varios millones de elementos hasta alcanzar prácticamente la capacidad total de memoria de la GPU. Para comprobarlo, se ha hecho uso del profiler que proporciona 52 5.3. Resultados de la paralelización 0 2 4 6 8 10 12 14 16 18 x 106 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5De 1 a 32 hilos N´umero de rayos Aceleraci´on 1 2 4 8 16 32 Figura 5.8: Aceleraciones para GPU de 1 a 32 hilos por bloque 0 2 4 6 8 10 12 14 16 18 x 106 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5De 64 a 1024 hilos N´umero de rayos Aceleraci´on 64 128 256 512 1024 Figura 5.9: Aceleraciones para GPU de 64 a 1024 hilos por bloque NVIDIA, nvprof, para 1 mililente y 1024 ledes. Los resultados de los kernels ymemcpy, así como las llamadas principales a la API, se muestran en la Tabla 5.3. Llamadas a los kernels ymemcpy Tiempo ( %) Tiempo Llamadas Media Mínimo Máximo Nombre 59,41 % 380,25 ms 1 380,25 ms 380,25 ms 380,25 ms CUDA memcpy DtoH 40,59 % 259,79 ms 1 259,79 ms 259,79 ms 259,79 ms cudaRaysPropagation 0,01 % 46,752 µs1 46,752 µs46,752 µs46,752 µscudaDirectRays 0,00 % 4,7360 µs5 947 ns 576 ns 1,4400 µsCUDA memcpy HtoD Llamadas a la API Tiempo ( %) Tiempo Llamadas Media Mínimo Máximo Nombre 54,36 % 380,49 ms 6 63,414 ms 3,3810 µs380,45 ms cudaMemcpy 37,11 % 259,80 ms 1 259,80 ms 259,80 ms 259,80 ms cudaEventSynchronize 5,04 % 35,262 ms 6 5,8770 ms 4,3220 µs34,455 ms cudaMalloc 3,30 % 23,097 ms 1 23,097 ms 23,097 µs23,097 ms cudaDeviceReset Tabla 5.3: Resultados de nvprof 5. Resultados 53 De hecho, si se mide la fracción de código paralelizada, el bucle de rayos, los resultados son bastante mejores. Hay que tener en cuenta que cuando serializa la llamada a un kernel alojado en la GPU, el host no queda bloqueado, es decir, puede seguir ejecutando código en paralelo al procesador gráfico, de forma que únicamente quedará bloqueado en caso de hacer una petición de lectura de memoria de tipo memcpy entre device-host. Este hecho genera un problema a la hora de realizar mediciones de tiempo, ya que con instrucciones de procesador comunes como clock no es posible medir los tiempos de cómputo en la GPU. Para ello, en [27] se describe cómo hacerlo correctamente empleando cudaEventSynchronize. La Figuras 5.10 recoge los resultados exclusivamente de la simulación. 0 0.5 1 1.5 2 x 107 0 2 4 6 8 10 12 14 16 18 20 De 1 a 32 hilos N´umero de rayos Aceleraci´on 1 2 4 8 16 32 0 2 4 6 8 10 12 14 16 18 x 106 0 2 4 6 8 10 12 De 64 a 1024 hilos N´umero de rayos Aceleraci´on 64 128 256 512 1024 Figura 5.10: Aceleraciones para GPU contando sólo la simulación Cuando se tiene en cuenta la simulación los resultados son prácticamente planos para todos los escenarios, sin importar la configuración geométrica escogida. De la misma forma, las aceleraciones son mayores a la unidad para todos los casos independientemente del número de rayos. Sin embargo, la tendencia observada en las aceleraciones globales se mantiene; el caso con 32 hilos es el que mejor se comporta, consiguiendo aceleraciones alrededor de 18 para ciertos casos. Para mostrar mejor este hecho, en la Figura 5.11 se comprueba la aceleración media para los distintos números de hilos donde se muestra de una forma más clara las zonas con aumento o disminución del rendimiento. Cambiando el código paralelo en CUDA R se podría hacer un análisis forzando el número de bloques pero no se creyó necesario. De los resultados obtenidos, puede concluirse que el verdadero cuello de botella de la aplicación reside en el paso de memoria entre la GPU y la CPU. En [28] se explica el uso de memoria no paginada (pinned) para agilizar la transferencia por DMA. Grosso modo se sustituye la función malloc, propia de C, por cudaMallocHost. Los resultados lejos de mejorar empeoran como puede apreciarse en la Figura 5.12. Esto es debido a que la transferencia de memoria se ha reducido a costa de aumentar el tiempo de reserva de memoria. Esto es fácilmente comprobable gracias a nvprof con la misma configuración: 1024 ledes y una sola mililente. En la Tabla 5.4 se muestra tanto los resultados de las principales funciones como los de la API. La principal diferencia con la Tabla 5.3 es que el tiempo consumido por cudaMemcpy se ha visto reducido en un orden de magnitud, de cientos de milisegundos a decenas; pero, en cambio, las funciones cudaHostAlloc ycudaFreeHost consumen 54 5.3. Resultados de la paralelización 1 2 4 8 16 32 64 128 256 512 1024 0 2 4 6 8 10 12 14 16 18 20 N´umero de hilos Aceleraci´on Figura 5.11: Aceleración media según el número de hilos por bloque 0 0.5 1 1.5 2 x 107 0 0.5 1 1.5 2 2.5 3 3.5 4De 1 a 32 hilos N´umero de rayos Aceleraci´on 1 2 4 8 16 32 0 2 4 6 8 10 12 14 16 18 x 106 0 0.5 1 1.5 2 2.5 3 3.5 4De 64 a 1024 hilos N´umero de rayos Aceleraci´on 64 128 256 512 1024 Figura 5.12: Aceleraciones para GPU con memoria no paginada Llamadas a los kernels ymemcpy Tiempo ( %) Tiempo Llamadas Media Mínimo Máximo Nombre 64,92 % 262,37 ms 1 262,37 ms 262,37 ms 262,37 ms cudaRaysPropagation 35,07 % 141,75 ms 1 141,75 ms 141,75 ms 141,75 ms CUDA memcpy DtoH 0,01 % 46,879 µs1 46,879 µs46,879 µs46,879 µscudaDirectRays 0,00 % 4,6720 µs5 934 ns 576 ns 1,4080 µsCUDA memcpy HtoD Llamadas a la API Tiempo ( %) Tiempo Llamadas Media Mínimo Máximo Nombre 32,20 % 304,40 ms 1 304,40 ms 304,40 ms 304,40 ms cudaHostAlloc 27,75 % 262,37 ms 1 262,37 ms 262,37 ms 262,37 ms cudaEventSynchronize 18,78 % 177,58 ms 1 177,58 ms 177,58 ms 177,58 ms cudaFreeHost 15,00 % 141,83 ms 6 23,639 ms 3,0980 µs141,80 ms cudaMalloc 11,62 % 27,200 ms 6 4,5334 ms 3,9090 µs27,165 ms cudaMemcpy 9,68 % 22,647 ms 1 22,647 ms 22,647 ms 22,647 ms cudaDeviceReset Tabla 5.4: Resultados de nvprof para pinned memory incluso más (en torno a 100 milisegundos). El tiempo de simulación de los kernels, cudaRaysPropagation ycudaDirectRays, no se ve afectado como cabría esperar puesto que sólo se ha intentado acelerar la transferencia de datos. 5. Resultados 55 La última optimización que se aplicó para intentar ganar algo de tiempo en la transferencia se encuentra en [29]. La idea básica es solapar la transferencia de memoria con la ejecución del kernel a modo de pipeline. Por un lado, la GPU debe ser compatible: en nuestro caso lo es y cuenta con un copy engine. Por el otro, la ejecución de la transferencia de memoria así como el kernel debe hacerse en un stream distinto al de por defecto (null stream). Por último, la memoria debe ser reservada en formato pinned, lo cual ya está hecho del paso anterior. En la Figura 5.13 se muestra gráficamente lo que se intenta conseguir. Figura 5.13: Ejecución solapada de kernel y transferencia de memoria En GPU con capacidad de cómputo superior o igual a 3.5 ambas versiones son indistinguibles gracias al Hyper-Q, que permite ejecutar varios kernels de forma concurrente; sin embargo, la GPU sobre la que se realizaron las simulaciones posee una capacidad de cómputo 3.0 por lo que se ha implementado la versión 2 que es la más eficiente. Destacar que la tranferencia H2D, en nuestro caso, es prácticamente nula. Los resultados (ver Figura 5.14), de nuevo, no son destacables puesto que, para 32 hilos, el mejor caso, estamos dividiendo el tiempo de transferencia y de ejecución pero las funciones cudaHostAlloc ycudaFreeHost no han cambiado como se comprueba, a su vez, en la Tabla 5.5 para el caso anteriormente considerado. Llamadas a los kernels ymemcpy Tiempo ( %) Tiempo Llamadas Media Mínimo Máximo Nombre 73,48 % 392,88 ms 2 196,44 ms 127,02 ms 265,86 ms cudaRaysPropagation 26,51 % 141,75 ms 3 47,251 ms 9,7920 µs70,896 ms CUDA memcpy DtoH 0,01 % 46,847 µs1 46,847 µs46,847 µs46,847 µscudaDirectRays 0,00 % 4,7680 µs5 953 ns 544 ns 1,4720 µsCUDA memcpy HtoD Llamadas a la API Tiempo ( %) Tiempo Llamadas Media Mínimo Máximo Nombre 46,22 % 463,58 ms 1 463,58 ms 463,58 ms 463,58 ms cudaEventSynchronize 30,23 % 303,22 ms 1 303,22 ms 303,22 ms 303,22 ms cudaHostAlloc 17,66 % 177,08 ms 1 177,08 ms 177,08 ms 177,58 ms cudaFreeHost 3,47 % 34,769 ms 6 5,7949 ms 2,9660 µs33,946 ms cudaMalloc 2,28 % 22,860 ms 1 22,860 ms 22,860 ms 22,860 ms cudaDeviceReset Tabla 5.5: Resultados de nvprof para ejecución solapada de 2 flujos 56 5.4. Resultados del simulador 0 2 4 6 8 10 12 14 16 18 x 106 0 0.5 1 1.5 2 2.5 3 3.5 432 hilos N´umero de rayos Aceleraci´on 1 2 4 8 Figura 5.14: Aceleraciones con ejecución solapada para distinto número de flujos El resto de funciones llamadas está dentro del orden del tiempo consumido en la Tabla 5.4, pero, al haber partido la ejecución de la simulación en dos, resulta que el kernel principal tarda algo más, yendo en detrimento del desempeño del sistema. Definitivamente, el cuello de botella de traerse los resultados a CPU penaliza sobremanera el rendimiento global de la implementación en GPU, por lo tanto, si se desease realizar el sintetizador con esta implementación se deberá considerar realizarlo en kernel para así aumentar el número de operaciones por rayo e intentar que la aplicación se convierta en intensiva en cómputo. A tenor de los resultados obtenidos con ambas implementaciones, memoria compartida y GPU, parece razonable pensar que la opción correcta para continuar el trabajo iniciado en este TFM es la de OpenMP. Doble precisión Se ha omitido los resultados relativos a las simulaciones sobre la arquitectura de GPU en doble precisión puesto que la discusión sobre la precisión ya se hizo para el paradigma de multihilo y además, al igual que entonces, el rendimiento es inferior. 5.4. Resultados del simulador En todo momento, se ha hablado de la simulación de un clúster de mililentes pero no de una pantalla completa, esto es, recordemos, debido al principio de superposición subyacente al escenario planteado. Aunque con esto es más que suficiente para obtener el funcional descrito en la Ecuación 1.1, en esta etapa temprana del trabajo se creyó adecuado realizar una recreación de una pantalla completa sobre el plano imagen. Es por ello que se realizó un código auxiliar secuencial denominado render que, como su nombre indica, se encarga de renderizar en escala de grises la luz proyectada sobre el plano imagen. El código recibe, por un lado, los resultados del código paralelizado en la memoria 5. Resultados 57 más el tamaño de la pantalla en píxeles así como el área de integración en el plano imagen. Por otro lado, necesita una matriz de conmutación que indique la intensidad de cada píxel como un parámetro comprendido entre 0 y 1. De esta forma, se permite un cierto grado de flexibilidad, los mismos que se tendrían en el caso de un sintetizador. No se ha hecho alusión a este código en el resto de la memoria puesto que se trata de algo totalmente accesorio que únicamente sirve para completar los resultados obtenidos, parecido a los scripts que se han desarrollado para tratar los datos. Para ilustrar esto, se ha elegido un escenario con los siguientes datos: Distancia a las mililentes: 1 mm. Distancia al plano de impacto: 3 m. Índice de refracción del medio: 1 (aire). Índice de refracción de las mililentes: 1,46 (cuarzo). Número de ledes: 16 por mililente. Número de mililentes del clúster: 49 Separación de las mililentes: 1,5 mm. Diámetro de las mililentes: 1,46 mm. Espesor de las mililentes: 1,24 mm. Tamaño de la protuberancia (cúpula): 114,73 µm. Radio de curvatura de las mililentes: 2,38 mm. Índice del emisor lambertiano: 1. Potencia total de un emisor: 1 W. Número de niveles en azimut y elevación: 25. Dimensión del array: 128 ×128. Área de integración: 1 cm2. Matriz de conmutación: todos encendidos. El resultado de simular la pantalla completa sobre el plano de impacto se puede ver en la Figura 5.15. Se ha decidido generar una imagen en escala de grises para dar una idea de donde recae la energía. El código, render, toma a partir de los bordes físicos del array un ángulo arbitrario, 60◦, e integra donde está la mayor parte de la energía. En muchos casos, la ausencia de rayos en una determinada zona implica que esa zona es oscura y no se debe en ningún caso a un problema de muestreo. Hay que destacar que debido al tamaño del array, 128 ×128, y del número de mililentes por array, 64, el número total de mililentes que hacen falta son 16, formando Bibliografía [1] J. Rodriguez-Ramos, J. Marichal-Hernandez, J. Luke, J. Trujillo-Sevilla, M. Puga, M. Lopez, J. Fernandez-Valdivia, C. Dominguez-Conde, J. C. Sanluis, F. Rosa, V. Guadalupe, H. Quintero, C. Militello, L. Rodriguez-Ramos, R. Lopez, I. Montilla, and B. Femenia, “New developments at CAFADIS plenoptic camera,” in Information Optics (WIO), 2011 10th Euro-American Workshop on, pp. 1–3, June 2011. [2] V. Guerra, C. Suarez-Rodriguez, S. Rodriguez, R. Perez-Jimenez, and J. Rodriguez-Ramos, “Plenoptics for optical wireless sensor networks,” in Information Optics (WIO), 2013 12th Workshop on, pp. 1–3, July 2013. [3] Creating a desired lighting pattern with an LED array, vol. 7058, 2008. [4] Homogeneous LED-illumination using microlens arrays, vol. 5942, 2005. [5] “IEEE Standard for Local and Metropolitan Area Networks–Part 15.7: Short-Range Wireless Optical Communication Using Visible Light,” IEEE Std 802.15.7-2011, pp. 1–309, Sept 2011. [6] S. Rajagopal, R. Roberts, and S.-K. Lim, “IEEE 802.15.7 visible light communication: modulation schemes and dimming support,” Communications Magazine, IEEE, vol. 50, pp. 72–82, March 2012. [7] E. Monteiro and S. Hranilovic, “Constellation design for color-shift keying using interior point methods,” in Globecom Workshops (GC Wkshps), 2012 IEEE, pp. 1224–1228, Dec 2012. [8] R. Drost and B. Sadler, “Constellation design for color-shift keying using billiards algorithms,” in GLOBECOM Workshops (GC Wkshps), 2010 IEEE, pp. 980–984, Dec 2010. [9] J. Luna-Rivera, R. Perez-Jimenez, J. Rabadan-Borjes, J. Rufo-Torres, V. Guerra, and C. Suarez-Rodriguez, “Multiuser CSK scheme for indoor visible light communications,” Opt. Express, vol. 22, pp. 24256–24267, Oct 2014. [10] O. Gonzalez, R. Perez-Jimenez, S. Rodriguez, J. Rabadan, and A. Ayala, “OFDM over indoor wireless optical channel,” Optoelectronics, IEE Proceedings -, vol. 152, pp. 199–204, Aug 2005. [11] J. Armstrong, “OFDM for Optical Communications,” Lightwave Technology, Journal of, vol. 27, pp. 189–204, Feb 2009. 65 66 BIBLIOGRAFÍA [12] S. Dissanayake and J. Armstrong, “Comparison of ACO-OFDM, DCO-OFDM and ADO-OFDM in IM/DD Systems,” Lightwave Technology, Journal of, vol. 31, pp. 1063–1072, April 2013. [13] Y. Jun, “Modulation and demodulation apparatuses and methods for wired/wireless communication system,” July 9 2009. US Patent App. 11/989,620. [14] N. Fernando, Y. Hong, and E. Viterbo, “Flip-OFDM for optical wireless communications,” in Information Theory Workshop (ITW), 2011 IEEE, pp. 5–9, Oct 2011. [15] V. Guerra, C. Suarez-Rodriguez, O. El-Asmar, J. Rabadan, and R. Perez-Jimenez, “Pulse width modulated optical OFDM,” in IEEE ICC 2015 - First Workshop on Visible Light Communications and Networking (VLCN) (ICC’15 - Workshops 24), (London, United Kingdom), June 2015. [16] S.-M. Kim and S.-M. Kim, “Performance improvement of visible light communications using optical beamforming,” in Ubiquitous and Future Networks (ICUFN), 2013 Fifth International Conference on, pp. 362–365, July 2013. [17] F. Yaras, H. Kang, and L. Onural, “State of the Art in Holographic Displays: A Survey,” Display Technology, Journal of, vol. 6, pp. 443–454, Oct 2010. [18] “IHS Technology.” https://technology.ihs.com/389045/. Último acceso en junio de 2015. [19] “DigInfo TV.” http://www.diginfo.tv/v/10-0155-r-en.php. Último acceso en junio de 2015. [20] “OpenMP.” http://openmp.org/wp/. Último acceso en junio de 2015. [21] “GeForce GTX 680 Whitepaper.” http://international.download.nvidia. com/webassets/en_US/pdf/GeForce-GTX-680-Whitepaper-FINAL.pdf. Último acceso en junio de 2015. [22] “CUDA R C Programming Guide.” http://docs.nvidia.com/cuda/cuda-cprogramming-guide/index.html#simt-architecture. Último acceso en junio de 2015. [23] “RED ONE R .” http://www.red.com/products/red-one. Último acceso en marzo de 2015. [24] “Microlens Arrays.” http://www.thorlabs.de/newgrouppage9.cfm? objectgroup_id=2861. Último acceso en marzo de 2015. [25] C. Schlick, “An Inexpensive BRDF Model for Physically-based Rendering,” Computer Graphics Forum, vol. 13, no. 3, pp. 233–246, 1994. [26] B. De Greve, “Reflections and Refractions in Ray Tracing ,” 2007. [27] “CUDA R C Best Practices Guide.” http://docs.nvidia.com/cuda/pdf/CUDA_ C_Best_Practices_Guide.pdf. Último acceso en junio de 2015. BIBLIOGRAFÍA 67 [28] “How to Optimize Data Transfers in CUDA C/C++.” http://devblogs. nvidia.com/parallelforall/how-optimize-data-transfers-cuda-cc/. Último acceso en junio de 2015. [29] “How to Overlap Data Transfers in CUDA C/C++.” http://devblogs.nvidia. com/parallelforall/how-overlap-data-transfers-cuda-cc/. Último acceso en junio de 2015.