Contributions to non-iterative super-resolution algorithms enabling efficient real-time FPGA implementation for resolution enhancement of video sequences

ContributionstoNon-IterativeSuper-ResolutionAlgorithmsEnablingEfficient Real-TimeFPGAImplementationforResolutionEnhancementofVideoSequences TESIS DOCTORAL IUMA TomaszSzydzik LasPalmasdeGranCanaria,Noviembre2015 InstitutoUniversitariodeMicroelectrónicaAplicada Tomasz Szydzik TESISDOCTORAL TESISDOCTORAL ContributionstoNon-Iterative Super-ResolutionAlgorithmsEnabling EficientReal-TimeFPGAImplementationfor ResolutionEnhancementofVideoSequences 2015 D. PEDRO PÉREZ CARBALLO SECRETARIO DEL INSTITUTO UNIVERSITARIO DE MICROELECTRÓNICA APLICADA DE LA UNIVERSIDAD DE LAS PALMAS DE GRAN CANARIA, CERTIFICA, Que el Consejo de Doctores del Departamento en su sesión de fecha 13 de noviembre de 2015 tomó el acuerdo de dar el consentimiento para su tramitación, a la tesis doctoral titulada “Contributions to Non-Iterative Super-Resolution Algorithms Enabling Efficient Real-Time FPGA Implementation for Resolution Enhancement of Video Sequences” presentada por el doctorando Tomasz Szydzik y dirigida por los Doctores D. Gustavo I. Marrero Callicó, y D. Antonio Nuñez Ordoñez. Y para que así conste, y a efectos de lo previsto en el Artº 6 del Reglamento para la elaboración, defensa, tribunal y evaluación de tesis doctorales de la Universidad de Las Palmas de Gran Canaria, firmo la presente en Las Palmas de Gran Canaria, a13 de noviembre de 2015. Instituto: Instituto Universitario de Microelectrónica Aplicada Programa de doctorado: Ingeniería de Telecomunicación Avanzada Título de la Tesis “CONTRIBUCIONES AL PROCESO DE SÚPER-RESOLUCIÓN MEDIANTE TÉCNICAS DE FILTROS SELECTIVOS, TOPOLOGÍA DE MACROBLOQUES ADAPTABLE Y SISTEMAS MULTI-CÁMARA” Tesis Doctoral presentada por D. Eduardo G. Quevedo Gutiérrez Dirigida por el Dr. D. Gustavo I. Marrero Callicó Codirigida por el Dr. D. Félix B. Tobajas Guerrero El Director, El Codirector El Doctorando, Las Palmas de Gran Canaria, a 14 de abril de 2015 Instituto: INSTITUTO UNIVERSITARIO DE MICROELECTRÓNICA APLICADA Programa de doctorado: INGENIERÍA DE TELECOMUNICACIÓN AVANZADA Título de la Tesis CONTRIBUTIONS TO NON-ITERATIVE SUPER-RESOLUTION ALGORITHMS ENABLING EFFICIENT REAL-TIME FPGA IMPLEMENTATION FOR RESOLUTION ENHANCEMENT OF VIDEO SEQUENCES Tesis Doctoral presentada por D. TOMASZ SZYDZIK Dirigida por Dr. D. GUSTAVO I. MARRERO CALLICÓ Codirigida por Dr. D. ANTONIO NUÑEZ ORDOÑEZ El Director, (firma) El Codirector (firma) El Doctorando, (firma) Las Palmas de Gran Canaria, a 13 de noviembre 2015 D. PEDRO PÉREZ CARBALLO SECRETARIO DEL INSTITUTO UNIVERSITARIO DE MICROELECTRÓNICA APLICADA DE LA UNIVERSIDAD DE LAS PALMAS DE GRAN CANARIA, CERTIFICA, Que el Consejo de Doctores del Departamento en su sesión de fecha 13 de noviembre de 2015 ha acordado que la tesis doctoral titulada “Contributions to Non-Iterative Super-Resolution Algorithms Enabling Efficient Real-Time FPGA Implementation for Resolution Enhancement of Video Sequences” presentada por el doctorando Tomasz Szydzik y dirigida por los Doctores D. Gustavo I. Marrero Callicó, y D. Antonio Nuñez Ordoñez reúne todo los requisitos para optar a la acreditación de DOCTORADO EUROPEO. Y para que así conste, firmo la presente en Las Palmas de Gran Canaria, a 13 de noviembre de 2015. Instituto: Instituto Universitario de Microelectrónica Aplicada Programa de doctorado: Ingeniería de Telecomunicación Avanzada Título de la Tesis “CONTRIBUCIONES AL PROCESO DE SÚPER-RESOLUCIÓN MEDIANTE TÉCNICAS DE FILTROS SELECTIVOS, TOPOLOGÍA DE MACROBLOQUES ADAPTABLE Y SISTEMAS MULTI-CÁMARA” Tesis Doctoral presentada por D. Eduardo G. Quevedo Gutiérrez Dirigida por el Dr. D. Gustavo I. Marrero Callicó Codirigida por el Dr. D. Félix B. Tobajas Guerrero El Director, El Codirector El Doctorando, Las Palmas de Gran Canaria, a 14 de abril de 2015 CONTRIBUTIONS TO NON-ITERATIVE SUPER-RESOLUTION ALGORITHMS ENABLING EFFICIENT REAL-TIME FPGA IMPLEMENTATION FOR RESOLUTION ENHANCEMENT OF VIDEO SEQUENCES Tomasz Szydzik Institute for Applied Microelectronics University of Las Palmas of Gran Canaria This dissertation is submitted for the degree of Doctor of Philosophy Las Palmas de Gran Canaria, November 2015 calidad de salida sub-óptima. La presente tesis doctoral aborda los principales desafíos encontrados actualmente a la hora de desarrollar implementaciones hardware de las técnicas de SR para secuencias de vídeo, dónde cabe destacar las siguientes: la capacidad de ejecución en tiempo real, una alta eficiencia de la implementación y la preservación de la misma calidad de salida de la imagen súper-resuelta que la obtenida con el software. El objetivo final de esta tesis ha sido proporcionar una implementación hardware caracterizada por las propiedades antes mencionadas. La implementación del algoritmo base en software no consideró en ningún momento desafíos cruciales de la implementación hardware, caracterizándose por una alta utilización de recursos de memoria que impiden su aplicación directa en tecnologías basadas en FPGAs utilizando solamente la memoria disponible en el dispositivo. Con el objetivo de superar estas limitaciones, se ha propuesto un flujo de ejecución modificado que funciona con elementos de grano más fino, aplicando la súper-resolución en el contexto de una ejecución centrada sólo en macro-bloques. Los resultados de las modificaciones propuestas conducen a una reducción significativa de la ocupación de memoria a expensas de un aumento del tráfico de memoria. El valor mínimo y máximo calculado del factor de reducción en la ocupación de la memoria asociada con el cambio del flujo de ejecución a nivel de macrobloques está entre 3,5 y 16, dependiendo de los valores de los parámetros del algoritmo. Al mismo tiempo, el valor mínimo y máximo calculado del factor de aumento en el tráfico de memoria asociado con dicho cambio de flujo está entre 1,1 y 16,9. Para cubrir las necesidades de una implementación hardware planificada se ha desarrollado una metodología de diseño de alto nivel. La metodología establecida define una jerarquía de niveles de abstracción, en la que cada uno de los modelos de la jerarquía se obtiene basándose en el modelo de nivel superior de dicha jerarquía. Esto hace que los modelos estén fuertemente ligados entre sí, lo que facilita la propagación de las modificaciones entre modelos de la jerarquía. El uso de una representación intermedia codificada usando SystemC genérico aumenta la portabilidad del diseño al permitir la síntesis de alto nivel dirigida a una serie de lenguajes HDL y de tecnologías a partir de la misma descripción de alto nivel. El nivel de detalles proporcionado por la descripción usando esta metodología facilita la reutilización de otros flujos establecido en las implementaciones de algoritmos similares. La arquitectura final alcanzó las prestaciones deseadas de 24 fps con una frecuencia de operación de 109 MHz utilizando el dispositivo FPGA xc5vj70t-l (Xilinx: tecnología Virtex5). La comparación realizada con el estado del arte ha resultado satisfactoria, dado que la ocupación de los recursos lógicos del sistema propuesto es hasta 5 veces menor que iv la reportada para el estado del arte usando el mismo tipo de tecnología FPGA. Los resultados de síntesis de la implementación hardware: (i) han demostrado la capacidad de alcanzar prestaciones en tiempo real usando tecnología FPGA, preservando al mismo tiempo la calidad de las imágenes de salida súper-resueltas al mismo nivel que el ofrecido por la referencia software, y (ii) han demostrado la fidelidad de los cambios realizados a nivel algorítmico y la validez de la metodología de implementación establecida. Los resultados de síntesis obtenidos suponen una contribución al estado del arte gracias a la implementación del flujo de procesamiento a nivel de MB (flujo típico de los algoritmos de compresión de imagen/vídeo) lo que supone una mayor eficiencia de la implementación. v Acknowledgements I would not have been able to complete this journey without the aid and support of countless people over the past four years. Foremost, I would like to express my gratitude to my supervisors, Prof. Gustavo I. Marrero Callicó and Prof. Antonio Nuñez Ordoñez, who have been greatly supportive and have guided me during this research and while writing this dissertation, offering constructive comments and warm encouragement. Over the years, I have received funding from several entities, which have supported me while I completed my research. I would like to thank Barco Electronic Systems S.A., Movidius Ltd., the European Network of Excellence on High Performance and Embedded Architecture and Compilation (HiPEAC) and the Institute for Applied Microelectronics (IUMA) for their financial support. In particular I would like to thank David Moloney and Luigi Albani for their generous support and for sharing their immense knowledge. I’d like to thank also my fellow labmates, for all the stimulating discussions, the fun we have had, the coffees and their patience. Last, but not least, I am deeply grateful to my parents for generously offering me support and the education that has made it possible for me to get here. Thanks also to my sister, Alina, and the rest of my family and for believing in me and instilling confidence in me. vii Contents Abstract i Resumen iii Acknowledgements vii Contents ix List of Figures xv List of Tables xxi List of Acronyms xxiii 1 Introduction 1 1.1 Challenges................................... 1 1.2 Researchmotivation.............................. 6 1.3 Objectives................................... 9 1.4 Organization of this thesis . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2 Image enhancement using super-resolution 15 2.1 Introduction.................................. 15 2.2 Super-resolution basics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.2.1 Applications of super-resolution . . . . . . . . . . . . . . . . . . . 16 2.2.2 Describing the super-resolution process . . . . . . . . . . . . . . . 18 2.2.3 Super-resolution in spatial domain . . . . . . . . . . . . . . . . . . 24 2.2.4 Super-resolution of video sequences . . . . . . . . . . . . . . . . . 30 2.3 Imageregistration............................... 34 2.3.1 Classification of image registration techniques . . . . . . . . . . . 36 Contents 2.3.2 Fundamental components and steps of image registration . . . . . . 37 2.3.3 Featuredomain............................ 39 2.3.4 Transformations ........................... 42 2.3.5 Searchspace ............................. 46 2.3.6 Searchstrategies ........................... 48 2.4 Super-resolution hardware implementations . . . . . . . . . . . . . . . . . 54 2.4.1 Multi-frame SRIR implementations . . . . . . . . . . . . . . . . . 55 2.4.2 Single-image SRIR implementations . . . . . . . . . . . . . . . . . 58 2.5 Conclusions.................................. 61 3 Reference super-resolution algorithm 63 3.1 Introduction.................................. 63 3.2 The non-uniform grid projection algorithm . . . . . . . . . . . . . . . . . . 64 3.2.1 Image registration . . . . . . . . . . . . . . . . . . . . . . . . . . 64 3.2.2 Fusion-based reconstruction . . . . . . . . . . . . . . . . . . . . . 65 3.2.3 Non-uniform interpolation . . . . . . . . . . . . . . . . . . . . . . 67 3.3 Mathematical model of the non-iterative classical super-resolution . . . . . 67 3.3.1 Overview of the non-iterative restoration-interpolation classical multiimage super-resolution . . . . . . . . . . . . . . . . . . . . . . . . 69 3.3.2 Image registration for the non-uniform grid projection algorithm . . 72 3.3.3 Mathematical model of the image registration stage . . . . . . . . . 75 3.3.4 Mathematical model of the restoration-interpolation stage . . . . . 80 3.4 Quantitative evaluation of the reference software implementation . . . . . . 82 3.4.1 Image quality assessment . . . . . . . . . . . . . . . . . . . . . . . 83 3.4.2 Output quality as a function of SR parameters . . . . . . . . . . . . 88 3.5 Conclusions.................................. 104 4 Proposed super-resolution algorithm 107 4.1 Introduction.................................. 107 4.2 Execution flows of the NUGPA . . . . . . . . . . . . . . . . . . . . . . . . 109 4.2.1 Reference coarse grain execution flow . . . . . . . . . . . . . . . . 109 4.2.2 Proposed finer grain execution flow . . . . . . . . . . . . . . . . . 111 4.3 Evaluation of memory occupancy . . . . . . . . . . . . . . . . . . . . . . 112 4.3.1 Identification of memory storage requirements . . . . . . . . . . . 112 4.3.2 Quantitative analysis of memory occupancy . . . . . . . . . . . . . 118 x Contents 4.4 Evaluation of memory accesses count carried out by the macroblockand frame-level NUGPA execution flows . . . . . . . . . . . . . . . . . . . . . 125 4.4.1 Modeling memory accesses . . . . . . . . . . . . . . . . . . . . . 125 4.4.2 Quantitative evaluation of memory traffic . . . . . . . . . . . . . . 132 4.5 Evaluation of the trade-offs of the proposed execution flow . . . . . . . . . 141 4.6 Optimization of memory traffic . . . . . . . . . . . . . . . . . . . . . . . . 143 4.6.1 Optimization technique . . . . . . . . . . . . . . . . . . . . . . . . 143 4.6.2 Optimization results . . . . . . . . . . . . . . . . . . . . . . . . . 144 4.7 Conclusions.................................. 145 5 Methodology 149 5.1 Introduction.................................. 149 5.2 Modeling of hardware-targeted systems . . . . . . . . . . . . . . . . . . . 150 5.2.1 High-level synthesis . . . . . . . . . . . . . . . . . . . . . . . . . 150 5.2.2 Electronic System-Level design . . . . . . . . . . . . . . . . . . . 153 5.2.3 The SystemC language . . . . . . . . . . . . . . . . . . . . . . . . 159 5.3 Design and implementation methodology . . . . . . . . . . . . . . . . . . 162 5.3.1 Established ESL design flow . . . . . . . . . . . . . . . . . . . . . 165 5.3.2 Structural decomposition and encapsulation . . . . . . . . . . . . . 167 5.3.3 Untimed/loosely timed TLM model . . . . . . . . . . . . . . . . . 173 5.3.4 ESL refinement methodology . . . . . . . . . . . . . . . . . . . . 174 5.4 Verification methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . 183 5.4.1 Algorithmic model validation . . . . . . . . . . . . . . . . . . . . 184 5.4.2 ESLverification ........................... 185 5.4.3 RTLverification ........................... 188 5.5 Conclusions.................................. 191 6 Hardware implementation 195 6.1 Introduction.................................. 195 6.2 Hardware implementation of the NUGP algorithm using FPGAs . . . . . . 196 6.2.1 Super-resolution system overview . . . . . . . . . . . . . . . . . . 197 6.2.2 Super-resolution core: architecture evolution . . . . . . . . . . . . 198 6.2.3 Super-resolution core: implemented/timed model architecture andorganization ........................... 201 6.3 Implementation challenges . . . . . . . . . . . . . . . . . . . . . . . . . . 204 6.3.1 Provision of synthesis-time customization . . . . . . . . . . . . . . 204 xi Contents 6.3.2 Memory related challenges . . . . . . . . . . . . . . . . . . . . . . 209 6.3.3 Data dependency in holes filling ................... 216 6.3.4 Variable size of search area . . . . . . . . . . . . . . . . . . . . . . 224 6.3.5 Grids construction and management . . . . . . . . . . . . . . . . . 228 6.3.6 Frame window size determination . . . . . . . . . . . . . . . . . . 232 6.3.7 Division operation . . . . . . . . . . . . . . . . . . . . . . . . . . 234 6.4 Results..................................... 240 6.4.1 Core-level implementation evaluation . . . . . . . . . . . . . . . . 241 6.4.2 Implementation results for each module . . . . . . . . . . . . . . . 243 6.5 Conclusions.................................. 245 7 Conclusions and future work 247 7.1 Introduction.................................. 247 7.2 Conclusions.................................. 248 7.3 Futurework.................................. 251 Appendix A Publications 255 Appendix B Resumen en castellano 259 B.1 Introducción.................................. 259 B.2 Imágenes de súper-resolución y secuencias de vídeo . . . . . . . . . . . . . 261 B.2.1 Aplicaciones de la súper resolución . . . . . . . . . . . . . . . . . 261 B.2.2 Clasificación de las técnicas de súper-resolución . . . . . . . . . . 262 B.2.3 Contextualización del algoritmo NUGP . . . . . . . . . . . . . . . 262 B.3 El algoritmo de proyección sobre la cuadrícula no-uniforme . . . . . . . . 266 B.3.1 Flujo de ejecución del algoritmo NUGP . . . . . . . . . . . . . . . 266 B.3.2 Evaluación cuantitativa de la calidad de la imagen súper-resuelta del algoritmoNUGP........................... 269 B.3.3 Flujo de ejecución con granularidad fina . . . . . . . . . . . . . . . 271 B.3.4 Evaluación del flujo de ejecución a nivel de macro bloque . . . . . 272 B.4 Metodología de diseño y verificación . . . . . . . . . . . . . . . . . . . . . 277 B.4.1 Visión general de la metodología de diseño y verificación . . . . . . 277 B.5 Implementación en FPGA . . . . . . . . . . . . . . . . . . . . . . . . . . 282 B.5.1 Visión general del sistema de súper-resolución . . . . . . . . . . . 282 B.5.2 Desafíos en la implementación . . . . . . . . . . . . . . . . . . . 284 xii Contents B.5.3 Arquitectura del núcleo de súper-resolución: arquitectura implementada y organización . . . . . . . . . . . . . . . . . . . . . . . 286 B.5.4 Resultados .............................. 289 B.6 Conclusiones ................................. 292 B.7 Líneas futuras de investigación . . . . . . . . . . . . . . . . . . . . . . . . 294 References 295 xiii List of Tables 2.1 Frequency vs. spatial domain SR . . . . . . . . . . . . . . . . . . . . . . . 20 2.2 Hierarchy and properties of 2D coordinate transformations. . . . . . . . . 44 2.3 Summary of hardware implementations executing in multi-frame context presentedinSection2.4. ........................... 60 3.1 Algorithm parameters and their values used in the study on image quality. . 90 3.2 Test sequences used in the experiment. . . . . . . . . . . . . . . . . . . . . 95 4.1 Equations for memory occupancy estimation of the identified memory types. 114 4.2 Equations for memory requirements of the algorithm steps defined for the analyzed software implementations. . . . . . . . . . . . . . . . . . . . . . 117 4.3 Computed minimal and maximal memory occupancy of the super-resolution kernel for frameand MB-level flows and QCIF input. . . . . . . . . . . . 118 4.4 Memory occupancy of the frame-level super-resolution kernel for investigated algorithm parameters values. . . . . . . . . . . . . . . . . . . . . . . 119 4.5 Computed memory storage requirements of defined entities of the analyzed software implementations . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 4.6 Equations modeling the count of accesses to memories carried out by the algorithmsteps................................. 129 4.7 Equations modeling the count of accesses of the identified algorithm steps required for processing of one frame. . . . . . . . . . . . . . . . . . . . . . 131 4.8 Noticed minimal, maximal and average values of probability figures. . . . . 132 4.9 Computed average values of probability figures for groups sharing algorithmparametervalue. ............................ 133 4.10 Computed minimal and maximal memory access counts of frameand MBlevel super-resolution kernels. . . . . . . . . . . . . . . . . . . . . . . . . 135 xxi List of Tables 4.11 Computed memory access counts of defined entities of the analyzed softwareimplementations. ............................ 139 5.1 Tools and their role in the Agility Complier centered deployment. . . . . . 166 5.2 Tools and their role in the C-to-Silicon centered deployment. . . . . . . . . 166 6.1 Expected savings gained by allowing variable search area size vs padding with zeros for a QCIF reference frame. . . . . . . . . . . . . . . . . . . . . 229 6.2 An example of the SFW cardinality computations, WLR flag and u(i)values. 234 6.3 Resulting hardware performance for two FPGA devices in comparison with [BB08]. .................................... 241 6.4 PSNR observed at the output (average over 90 initial frames; CIF; luma component)................................... 242 6.5 Objective quality boost over interpolation cost in LUTs. . . . . . . . . . . . 242 6.6 Observed resource occupancy for some SRK parameters values combinations.243 B.1 Secuencias de test utilizadas en los experimentos. . . . . . . . . . . . . . . 270 B.2 PSNR observado de la secuencia de salida. . . . . . . . . . . . . . . . . . . 270 B.3 Resultados de las prestaciones hardware para dos dispositivos FPGA en comparacióncon[BB08]............................ 290 B.4 Figura de mérito boostCostSR frente al coste de interpolación en términos deLUTs. ................................... 291 xxii List of Acronyms ANSI American National Standards Institute ASIC Application Specific Integrated Circuit AT Approximately Timed BCA Bus Cycle Accurate BFM Bus Functional Model BM Block Matching BRAM Block Random Access Memory BSynth Behavioral Synthesizable CA Cycle Accurate CCF Cross-Correlation Function CFT Continuous Fourier Transform CGI Computer-Generated Imagery CIF Common Intermediate Format CMOS Complementary Metal Oxide Semiconductor COTS Commercial Off-The-Shelf CPU Central Processing Unit CT Computed Tomography CUDA Computer Unified Device Architecture xxiii List of Acronyms DFT Discrete Fourier Transform DoF Degree of Freedom DSP Digital Signal Processor DUT Design Under Test EDIF Electronic Design Interchange Format ESA European Space Agency FL Frame-Level FLM Functional Level Model FPGA Field Programmable Gate Array FR FRame GPU Graphics Processing Unit HD High Definition HDL Hardware Description Language HDTV High Definition TeleVision HEVC High Efficiency Video Coding HL High Level HLA High Level (of) Abstraction HLS High Level Synthesis HR High Resolution I/O Input/Output IBP Iterative Back Projection IP Intellectual Property IR Image Registration xxiv List of Acronyms IUMA Institute for Applied Microelectronics LR Low Resolution LT Loosely Timed LUT Look-Up Table MAD Mean Absolute Difference MAD Multiply and ADd MAE Mean Absolute Error MAP Maximum A Posteriori MB Macroblock MBL Macroblock-Level ME Motion Estimation ML Maximum Likelihood MPSoC MultiProcessor System on Chip MRI Magnetic Resonance Imaging MSE Mean Square Error MSSIM Multi-scale Structural SIMilarity N/A Not Available NTSC National Television System Committee NUGP Non-Uniform Grid Projection NUGPA Non-Uniform Grid Projection Algorithm OSCI Open SystemC Initiative PCA Pin Cycle Accurate PDC Pixel Difference Classification xxv List of Acronyms PLD Programmable Logic Device POCS Projection Onto Convex Sets PSNR Peak Signal-to-Noise Ratio QCIF Quarter Common Intermediate Format QVGA Quarter Video Graphics Array RAM Random Access Memory RB Reorder Buffer RF Reference Frame RGB Red, Green, Blue ROI Region Of Interest ROM Read Only Memory ROV Remotely Operated Vehicle RS Reference Structure RT Register Transfer (accurate) RTL Register Transfer Level SA Search Area SAM System Architectural Model SAR Search Area Radius SFW Sliding Frame Window SLM System Level Model SoC System-On-Chip SPM System Performance Model SR Super-Resolution xxvi List of Acronyms SRIR Super-Resolution Image Reconstruction SSDA Sequential Similarity Detection Algorithm SSIM Structural SIMilarity STORM Stochastic Optical Reconstruction Microscopy SVC Scalable Video Coding SXGA Super Extended Graphics Array TF Timed Functional TLM Transaction Level Modeling TMAC Total Memory Access Count TMR Total Memory storage Requirements TVM Transaction Verification Model ULPGC University of Las Palmas de Gran Canaria UT UnTimed VGA Video Graphics Array VHDL VHSIC Hardware Description Language VHSIC Very High Speed Integrated Circuit VLIW Very Long Instruction Word WVGA Wide Video Graphics Array WXGA Wide Extended Graphics Array XGA Extended Graphics Array XST Xilinx Synthesis Technology xxvii Chapter 1 Introduction 1.1 Challenges We live in a reality dominated by visual content, more specifically, by High Resolution (HR) visual content. High resolution images are of crucial importance in a number of areas centered on two main applications, namely the improvement of pictorial information for human interpretation, and robust automatic machine vision. In order to fit with the new level of expected visual quality, the Low Resolution (LR) contents quality must be augmented in a process called resolution enhancement. Image resolution describes the details contained in an image: the higher the resolution of an image, the higher the amount of details it contains. The resolution of a digital image can refer to: pixel resolution, spatial resolution, spectral resolution, temporal resolution, and/or radiometric resolution. In this context, this PhD Thesis focuses on contributing to the improvement of spatial resolution. (A) Considering a digital image as composed of picture elements (pels), spatial resolution of digital images is expressed (or measured) as the total number of these elements and their density understood as the number of the elements per area unit. In the context of digital images, a singular discrete pel is usually referenced as a pixel. In practice, the spatial resolution of most of the imaging systems is limited by the properties of the used sensor rather than the optical diffraction. The most intuitive solution to increase spatial resolution is to reduce the pixel size (i.e. increase the number of pixels per area unit) by sensor manufacturing techniques. However, as the pixel size decreases so does the number of photons captured and the amount of light being integrated by each pixel per time unit. In consequence the sensor is more susceptible to shot noise, that can degrade the image quality severely. The acceptable value of the shot noise 1 CHAPTER1.–INTRODUCTION First SR algorithm [Ger74] First multiple image SR in frequency domain [TH84] First hallucination (neural network) [Mjo85] Non-uniform grid projection algorithm [MC03] Shift and add [FM06] First non-parameteric [PC09] 1974 1984 1985 2003 2006 2009 First multiple image SR in spatial domain [PKS87] 1987 Fig. 1.2 Excerpt of the time line of proposals of SR algorithms. parallelism all of the data would be required to be readily available from memories, the execution units have to be replicated/pipelined, etc. On one hand the processing has to be parallelized, on the other hand achievable parallelism is limited by the memory required for its exploitation. All of the above has led to a point where it is not guaranteed that the proposed contributions can be easily propagated onto hardware. Hence, there is a growing need of provision of know-how or procedures and recommendations on how to consider the implications of possible hardware implementations at the time of algorithm definition in order to help to close the ever-growing gap between software and hardware implementations. As pointed out, software based implementations are not capable of providing HR realtime performance when mapped on available technology. The currently available hardware implementations with real-time performance provide an output quality lower than their software implementation counterparts. Even the state-of-the-art implementations using FieldProgrammable Gate Array (FPGA) technology have their performance limited by memory storage requirements and/or memory access latency. Thus, only a limited number of iterations or reference frames are implemented, leading to sub-optimal output quality. Based on our experience with compression algorithms we find much similarity in the nature of processing carried out by SRIR and compression algorithms and the problems faced by their implementations. Both types of processing deploy complex algorithms that work with large amounts of data and require hardware implementations to process them in real-time. Moreover, many of the compression algorithms pose hard real-time constraints on the execution time. Hardware implementations of compression algorithms can provide solutions for how to implement algorithms that process images so that their hardware implementations can achieve real-time performance. 8 1.3 OBJECTIVES This doctoral thesis provides an example of the required modifications that alleviate the aforementioned problems and therefore could result in the future hardware implementations of super-resolution based resolution enhancement, becoming efficient enough to allow their use in mainstream consumer electronics. In this line, we propose an SRIR implementation based on a non-iterative version of the Non-Uniform Grid Projection Algorithm (hereafter the NUGP algorithm or NUGPA). Direct implementation of the reference version of this algorithm is considered unviable due to high memory requirements. In this work we evaluate and implement modifications to the NUGPA execution flow so that its processing can be applied in a way that only requires a singular macroblock execution context to be readily available from memories. This approach has achieved so much success in the context of compression algorithms implementation that we hope to replicate it for the needs of SRIR. This change has some significant implications on the overall algorithm execution flow that have to be considered at the time of implementation. Thus, the algorithm is modified in such ways that not only provide real-time performance and observed quality of the super-resolved output that matches the level of the software state-of-the-art, but also facilitate efficient hardware implementation. Only by providing all of these characteristics the high fidelity resolution enhancement can find its way out of research facilities to everyone’s home. The lessons learned during the process of NUGPA implementation are expected to contribute to the understanding of the hardware implementation challenges and limitations knowledge by the software/algorithmic designers. In the long term, this could contribute to the process of providing more hardware-aware implementations, simplifying the process of hardware implementation and most likely leading to increased efficiency and performance. 1.3 Objectives The ultimate goal of this thesis is to provide solutions at the algorithmic and architectural level to the fusion-based SR algorithm developed by the Institute for Applied Microelectronics (IUMA) of the University of Las Palmas de Gran Canaria (ULPGC) that would allow reaching real-time execution when implemented in FPGA devices. In order to achieve this goal, this doctoral thesis is focused on reaching the following objectives: 1. Contextualize this work by exploring the state-of-the-art of super-resolution methods and their hardware implementations. In this process, the state-of-the-art of SR algorithms and a complete classification of the available methods will be described. The study is carried out in order to: (i) highlight strengths and weaknesses of the NUGPA approach in comparison with other approaches, and (ii) justify the decision 9 CHAPTER1.–INTRODUCTION of using NUGPA as the base algorithm for hardware implementation capable of real-time execution. Following, the known approaches to hardware implementations of SR in FPGAs will be presented. At this point it is important to pay attention: (i) to the parallelization strategies exploited, (ii) techniques applied for improving the final performance of the system, (iii) avoidance of system bottlenecks, and (iv) required trade-offs. Based on the obtained knowledge, an adequate development strategy will be developed. 2. Provide algorithmic contributions to facilitate hardware implementation of the NUGPA super-resolution kernel (SRK) using only on-chip memory. The available version of the NUPGA does not contemplate hardware implementation in programmable logic devices. Similarly, the available baseline software was not implemented with the view of limiting the use of resources. Thus, in order to increase the chances of a successful and efficient implementation using FPGA devices the algorithm data flow should be revised and modified if needed. 3. Create a design flow and implementation methodology that increases the chances of obtaining an efficient implementation. The methodology should consider the fact that the SR algorithm and its software implementation are constantly under development by facilitating rapid and robust propagation of the modifications made at the algorithmic level to the hardware implementation. The created design and implementation methodology should be general enough to be able to be applied to other similar algorithms. At the same time, the description should be specific enough to remove any doubts on the steps required to provide high quality results. All these requirements lead to a combination of Electronic System Level (ESL) and Register Transfer Level (RTL) synthesis flow. 4. Demonstrate the viability of reaching real-time performance by providing a baseline implementation of the NUGPA SR kernel in a FPGA device. Apart from the performance goal, the implemented design must be able to fulfill additional expectations of which the most relevant are listed below: (a) Efficiency: in this specific case understood as minimization of resources (logical and memory) occupancy while providing a solution that meets the targeted performance in terms of execution speed. The hardware implementation is expected to be successful in preserving the observed super-resolved image quality at the level offered by software implementations. 10 1.4 ORGANIZATION OF THIS THESIS (b) Portability: this characteristic is related to the fact that we expect the implementation to implement solutions and mechanisms that facilitate migration to a different technology (different FPGA family/vendor) without having to redesign it completely from scratch. Having in mind that FPGA implementations are frequently a required prototyping step for custom hardware, this migration should also help in the road towards ASIC and SoC implementation targets. 1.4 Organization of this thesis The present document is structured into seven chapters that somewhat follow the order of the objectives given above. The first three chapters (chapters 1–3) are dedicated to presenting the problematic, the state of the art and the motivations for developing this Doctoral thesis. These chapters provide information required to fully understand the issues being raised in this thesis by readers who might be not familiar with the topic of super-resolution, its current state-of-the-art, and the NUGP algorithm in particular. The introductory nature of these chapters serves to set the scene for presenting the main contributions of this work. Each of the chapters of this work is considered to be self-confined, that is, understanding of the main contributions presented in a chapter does not require previous reading of any of the chapters preceding it. Thus, readers interested in particular objectives of this thesis (modeling, methodology, or hardware implementation) may start reading directly from the chapter of their choosing. Chapter 2: IMAGE ENHANCEMENT USING SUPER-RESOLUTION This chapter introduces the basic concepts of the super-resolution process and the state-of-the-art of the SR methods. The presented data allow contextualizing the noniterative non-uniform grid projection algorithm within the state-of-the-art, show its weaknesses and strengths, and justify its choice as a base for hardware implementation. The final objective of this chapter is to describe the available approaches to hardware implementations of SR in FPGAs, exposing their weaknesses, highlighting the scarce number of such implementations and the need of higher efficiency in terms of resources occupancy in order to overcome the observed limitations. Chapter 3: THE NON-UNIFORM GRID PROJECTION ALGORITHM (NUGPA) Provision of efficient hardware implementation of any algorithm requires in-detail knowledge on the algorithm control and data flow, as well as their dependency on the 11 CHAPTER1.–INTRODUCTION algorithm parameters configuration. This chapter focuses on presentation and evaluation of the base SR algorithm and its baseline software implementation. The description emphasizes the importance of using a single-pass (non-iterative) approach in order to facilitate hardware implementation and identifies memory occupancy as the likely factor precluding hardware implementation of the reference implementation. This chapter finalizes presenting results of a study on the quality of the superresolved image produced by the software implementation. The presented quantitative data proves that, when appropriately configured, the base implementation is capable of providing output whose quality of results is not only better than interpolation but also on-par or better than the one provided by state-of-the-art hardware SR implementations. Chapter 4: PROPOSED SUPER-RESOLUTION ALGORITHM The NUGPA algorithm was expected to pose memory requirements that were prohibitively high for a hardware implementation meeting the proposed objectives. In order to tackle this problem the algorithm dependency on Frame-Level (FL) buffers had to be eliminated. The focus of this chapter is the set of modifications of the NUGPA execution flow proposed in order to facilitate its hardware implementation. The chapter presents models and quantitatively evaluates memory occupancy and traffic of both the reference and the proposed execution flows. The obtained data show that the proposed modifications are expected to lead to reduction in memory occupancy at the cost of increased memory traffic. The chapter concludes presenting an example of system bottleneck identification data and then algorithmic-level transformation in order to optimize the system. Chapter 5: IMPLEMENTATION METHODOLOGY Most of the contributions to the state-of-the-art of SR methods focus on introducing algorithmic-level contributions that are aimed at improving the observable output quality. The effectiveness of these proposals is usually proved mathematically or by means of software implementations. The latter, in most cases, are not capable of applying the processing in real-time. Meeting real-time constraints usually requires the algorithm to be mapped onto hardware. With the ever-growing complexity of the algorithms the organization of the process of mapping becomes a challenge itself. In our implementation we have opted for using a tiered hierarchy of abstraction models. Each of the defined tiers of abstraction is accompanied by a verification set-up. The use of multiple models of abstraction that are progressively more accurate allows re12 1.4 ORGANIZATION OF THIS THESIS ducing the gap separating the created system descriptions and greatly facilitates the propagation of changes from one tier to another. This chapter focuses on the presentation of the methodology established to carry out the process of mapping the NUGP algorithm onto FPGA technology. The used models and refinement steps carried out at each stage of mapping are organized and presented in the order that follows the design and implementation flow stages. The presented description provides details on how to set-up the verification environment in a way that leverages mixed-language simulation allowing the same environment to be used for both ESL and RTL simulations. Chapter 6: HARDWARE IMPLEMENTATION Definition of an implementation methodology is the first step required for successful implementation. A well-defined methodology is the stepping stone of any hardware implementation and significantly increases the probability of successful mapping. Even if provided with a well-defined implementation methodology, meeting the implementation goals is not guaranteed as the designer has yet to tackle architecturelevel implementation challenges. Solutions for these challenges are constrained by the targeted system characteristics (performance, output quality, efficiency, etc.) whose implications have to be taken into consideration at the moment of balancing the implementation trade-offs. Delivery of a hardware implementation that reaches the targeted performance while maximizing the efficiency of resources utilization is not a trivial task. This chapter presents the process and the results of implementation of the NUGP algorithm operating at MB-level in FPGA technology. The description of the implementation process provides information on the system architecture model evolution, final organization, and the used abstraction levels. The tackled implementation challenges are presented in great detail along with the implemented solutions. The observed synthesis results prove that the MB-level processing contributes towards real-time implementations of the NUGPA by significantly reducing the implementations memory occupancy while preserving software-level output quality. The comparison with the state-of-the-art yields satisfactory results in terms of implementation efficiency considered as the observed gain in the quality of the super-resolved output (vs interpolation) per resource utilized. Chapter 7: CONCLUSIONS AND FUTURE WORK This chapter presents a recapitulation of the contributions provided in this doctoral dissertation and highlights their relevance to the super-resolution field. This document 13 CHAPTER1.–INTRODUCTION concludes with a presentation of further research lines, which might complement and enhance some of the aspects developed in this thesis. 14 Chapter 2 Image enhancement using super-resolution 2.1 Introduction Super-resolution comprises algorithms that attempt to reconstruct the high resolution image corrupted due to the limitations of the imaging system. Typically these limitations arise in the system’s sensor resolution and/or the optics causing aliasing to occur as the sampling process does not meet the requirements defined by the Nyquist sampling theorem. Aliasing in digital images is often considered as a nuisance and (both optical and digital) filters are designed to avoid aliasing in digital cameras. However, aliasing also contains extra high-frequency information with additional details about the scene. Super-resolution algorithms extract the information present in the aliasing to reconstruct a higher-resolution image [Mil10]. Due to its ill-posed nature and vast applications, super-resolution has been a very active area of research since its introduction in 1974. Over the last 40 years several methods have been proposed under the umbrella of super-resolution. Most super-resolution algorithms consist of two main stages: (i) meta data registration, where the images are precisely aligned and analyzed, and (ii) image reconstruction, where the input data and meta data are used to estimate a higher resolution image. The most recent and complete taxonomy of super-resolution methods, presented in [NM14], classifies the introduced SR techniques into four main families: (i) frequency based approaches, (ii) interpolation-based approaches, (iii) regularization-based approaches, (iv) example-based approaches. The first three categories reconstruct (or super-resolve) the image from a set of lower resolution input images (multiple-image context), while the last example-based approaches are capable of achiev15 CHAPTER2.–IMAGE ENHANCEMENT USING SUPER-RESOLUTION ing the same objective by exploiting the information provided by an image and a database with a prori attained knowledge. Recently, super-resolution has gain a lot of attention from industry which resulted in multiple software and hardware super-resolution-based solutions destined for the medical, military and even consumer electronics market. 2.2 Super-resolution basics The general definition of the super-resolution process given in the introduction will now be refined. In particular: (i) the main characteristics of the super-resolution process will be described, and (ii) a classification of the super-resolution methods will be introduced. The Image Registration (IR) stage that forms a critical part of the multi-image approaches is presented in details in Section 2.3. 2.2.1 Applications of super-resolution In this work we define the objective of super-resolution in reference to the limitations of the imaging system that the SR technique is designed to overcome. The limitations tackled by the SR methods can be classified into the following three categories: optical,geometric, and temporal. 2.2.1.1 Optical limitations The diffraction-limited resolution theory was formulated by Ernst Abbe in 1873 and later refined by Lord Rayleigh in 1896. These researchers contemplated the separation necessary between two airy patterns in order to distinguish them as separate entities. They have observed and specified an upper limit —the cut-off spatial-frequency— below which the two point sources are considered to be resolved and can readily be distinguished. Beyond the cut-off spatial-frequency structural details fail to be correctly transferred into the optical image and cannot be resolved. In order to transcend the limitations of optical imaging systems a sequence of images which meet the Rayleigh criterion can be used. By manipulating the imaging conditions between the subsequent captures (i.e. using spatially-variant excitation light, gathering light over a larger set of angles around the specimen, swapping spatial frequency bands beyond the cut-off spatial frequency for one inside it (activators), changing the excitation light, etc.) it is possible to generate and transfer complementary subsets of information to different observations. Given that these observations are stochastically unique, then the extracted information can be fused, using an SR method, in order to form the final 16 2.2 SUPER-RESOLUTION BASICS higher-resolution image. This method allows two particles that appear to be merged, i.e. form a single point in a particular image, to be discriminated into two independent particles. Even in the case of flawless fabrication, glass-based systems resolution is still hampered by an ultimate limit in optical resolution that is imposed by the diffraction of visible light wavefronts. The use of complementary information from the context is considered to determine the exact position of the two adjacent particles even below the Rayleigh limit, effectively breaking the diffraction barrier limiting the optical imaging systems. Hence, SR methods from this category have found great success in glass-based microscopy. 2.2.1.2 Geometric limitations In practice, the spatial resolution of most of the imaging systems is limited by the properties of the used sensor rather than the optical diffraction. The most intuitive solution to increase spatial resolution is to reduce the pixel size (i.e. increase the number of pixels per unit area) by sensor manufacturing techniques. However, as the pixel size decreases so does the number of photons captured and the amount of light being integrated by each pixel per time unit. In consequence the sensor is more susceptible to shot noise that can degrade the image quality severely. The acceptable value of the shot noise imposes a technology dependent limitation on the minimal pixel size. In case of 0.35 µmCMOS process the optimal pixel size is estimated at about 40 µm2. This size has been already reached by the current image sensor technology. On the other hand, increasing the chip size to accommodate more pixels, and thus augment the spatial resolution, has its own limitations, being the most important the corresponding increase in sensor area and capacitance. Increased capacitance results in slower transfer rates limiting the temporal resolution and system’s applications. Increasing area size results in higher fabrication costs. Along with the arising requirement of higher precision optics, these properties significantly increase the production costs. The abovementioned ‘technology’ based solutions are considered to be not cost-efficient. Therefore, a new approach toward increasing spatial resolution is required to overcome these limitations of the sensors and manufacturing technology. With the rise of easily accessible computation power, spatial resolution enhancement using soft-processing became a promising alternative. In the context of limited spatial resolution, super-resolution comprises the methods that allow increasing spatial resolution beyond the native resolution of the sensor. 17 CHAPTER2.–IMAGE ENHANCEMENT USING SUPER-RESOLUTION the best match(es) are identified using the method’s match criteria (e.g. similarity, spatial compatibility, smoothness, etc.). One or more of these matches are used to compute the reconstruction weights that minimize the reconstruction error for the input LR patch. The weights are used in combinations with the HR data associated with the matched entries to hallucinate the HR patch. The estimated HR representation is projected onto the HR grid and optionally back-projected to feed the error estimation loop and adapt (if necessary) the processing of the subsequent patches. The initial HR representation of the entire frame is further regularized and refined using the reconstruction constraint. The reconstruction constraint’s penalty terms are also derived from the training set. 2.2.2.2.3 Combined or multi-patch approaches. The advantages of the described approaches can be combined. The combined framework is based on an observation (justified statistically) that patches in a single natural image tend to redundantly recur many times inside the image, both within the same scale, as well as across different scales. Given an input image, the combined framework creates a pyramid of its scaled representations and looks for recurrence of patches within and across the scales. Patches recurring within the same image scale (with sub-pixel misalignments) can be regarded as if extracted from different observations of the same scene and used to impose constraints to the HR information. This forms the basis for applying the multi-patch counterpart of the classical SR approach as illustrated in Fig. 2.4. Recurrence of patches across different image scales, as illustrated in Fig. 2.5, implicitly provides examples of low-resolution/high-resolution pairs of patches, thus, giving rise to learning-based super-resolution from a single image. The great benefits of the combined framework are the capability to obtain super-resolved observations from as little as a single low resolution image, the eliminations of the external a priori-created database and the off-line training phase. 2.2.3 Super-resolution in spatial domain Spatial domain approaches offer powerful methods to deal with the ill-posed problem of image restoration. The most recent and complete taxonomy of SR methods has been proposed in [NM14] and is presented in Fig. 2.6. This figure uses the execution context as the main classification criteria. The singleand multiple-images approaches are further sub-divided into families based on the used SR technique. 24 2.2 SUPER-RESOLUTION BASICS patches are used, e.g., 5×5, such patch repetitions occur abundantly within and across image scales, even when we do not visually perceive any obvious repetitive structure in the image. This is due to the fact that very small patches often contain only an edge, a corner, etc. such patches are found abundantly in multiple image scales of almost any natural image. Moreover, due to the perspective projection of cameras, images tend to contain scene-specific information in diminishing sizes (diminishing toward the horizon), thus recurring in multiple scales of the same image. We statistically tested this observation on the Berkeley Segmentation Database1(Fig. 2). More specifically, we tested the hypothesis that small 5×5patches in a single natural grayscale image, when removing their DC (their average grayscale), tend to recur many times within and across scales of the same image. The test was performed as follows: Each image Iin the Berkeley database was first converted to a grayscale image. We then generated from Ia cascade of images of decreasing resolutions {Is}, scaled (down) by scale factors of 1.25sfor s=0,−1, .., −6 (I0=I). The size of the smallest resolution image was 1.25−6=0.26 of the size of the source image I(in each dimension). Each 5×5patch in the source image Iwas compared against the 5×5patches in all the images {Is}(without their DC), measuring how many similar2patches it has in each image scale. This intra-image patch statistics was computed separately for each image. The resulting independent statistics were then averaged across all the images in the database (300 images), and are shown in Fig. 2a. Note that, on the average, more than 90% of the patches in an image have 9or more other similar patches in the same image at the original image scale (‘within scale’). Moreover, more than 80% of the input patches have 9or more similar patches in 0.41=1.25−4of the input scale, and 70% of them have 9or more similar patches in 0.26=1.25−6of the input scale. Recurrence of patches forms the basis for our singleimage super-resolution approach. Since the impact of super-resolution is expressed mostly in highly detailed image regions (edges, corners, texture, etc.), we wish to eliminate the effect of uniform patches on the above statistics. Therefore, we repeated the same experiment using only 25% of the source patches with the highest intensity variance. This excludes the uniform and low-frequency patches, 1www.eecs.berkeley.edu/Research/Projects/CS/vision/grouping/segbench 2Distances between patches were measured using gaussian-weighted SSD. Note that textured patches tend to have much larger SSD errors than smooth (low-variance) patches when compared to other very similarlooking patches (especially in the presence of inevitable sub-pixel misalignments). Thus, for each patch we compute a patch-specific ‘good distance’, by measuring its (gaussian-weighted) SSD with a slightlymisaligned copy of itself (by 0.5 pixel). This forms our distance threshold for each patch: Patches with distance below this threshold are considered similar to the source patch. (a) Classical Multi-Image SR (b) Single-Image Multi-Patch SR Figure 3: (a) Low-res pixels in multiple low-res images impose multiple linear constraints on the high-res unknowns within the support of their blur kernels. (b) Recurring patches within a single low-res image can be regarded as if extracted from multiple different low-res images of the same high resolution scene, thus inducing multiple linear constraints on the high-res unknowns. maintaining mostly patches of edges, corners, and texture. The resulting graphs are displayed in Fig. 2b. Although there is a slight drop in patch recurrence, the basic observation still holds even for the high-frequency patches: Most of them recur several times within and across scales of the same image (more than 80% of the patches recur 9or more times in the original image scale; more than 70% recur 9 or more times at 0.41 of the input scale, and 60% of them recur 9or more times in 0.26 of the input scale.) In principle, the lowest image scale in which we can still find recurrence of a source patch, provides an indication of its maximal potential resolution increase using our approach (when the only available information is the image itself). This is pixel-dependent, and can be estimated at every pixel in the image. 3. Single Image SR – A Unified Framework Recurrence of patches within the same image scale forms the basis for applying the Classical SR constraints to information from a single image (Sec. 3.1). Recurrence of patches across different scales gives rise to ExampleBased SR from a single image, with no prior examples (Sec. 3.2). Moreover, these two different approaches to SR can be combined into a single unified computational framework (Sec. 3.3). 3.1. Employing in-scale patch redundancy In the classical Multi-Image Super-resolution (e.g., [12, 5, 8]), a set of low-resolution images {L1, ..., Ln}of the same scene (at subpixel misalignments) is given, and the goal is to recover their mutual high-resolution source image H. Each low resolution image Lj(j=1,...,n)is assumed to have been generated from Hby a blur and subsampling process: Lj=H∗Bj↓sj, where ↓denotes a subsampling operation, sjis the scale reduction factor (the subsampling rate) between Hand Lj, and Bj(q)is the corresponding blur kernel (the Point Spread Function – PSF), represented in the high-resolution coordinate system – see Fig. 3a. Thus, each low-resolution pixel p=(x, y)in each low-resolution image Ljinduces one linear constraint on 351 (a) Classical multiple image SR. patches are used, e.g., 5×5, such patch repetitions occur abundantly within and across image scales, even when we do not visually perceive any obvious repetitive structure in the image. This is due to the fact that very small patches often contain only an edge, a corner, etc. such patches are found abundantly in multiple image scales of almost any natural image. Moreover, due to the perspective projection of cameras, images tend to contain scene-specific information in diminishing sizes (diminishing toward the horizon), thus recurring in multiple scales of the same image. We statistically tested this observation on the Berkeley Segmentation Database1(Fig. 2). More specifically, we tested the hypothesis that small 5×5patches in a single natural grayscale image, when removing their DC (their average grayscale), tend to recur many times within and across scales of the same image. The test was performed as follows: Each image Iin the Berkeley database was first converted to a grayscale image. We then generated from Ia cascade of images of decreasing resolutions {Is}, scaled (down) by scale factors of 1.25sfor s=0,−1, .., −6 (I0=I). The size of the smallest resolution image was 1.25−6=0.26 of the size of the source image I(in each dimension). Each 5×5patch in the source image Iwas compared against the 5×5patches in all the images {Is}(without their DC), measuring how many similar2patches it has in each image scale. This intra-image patch statistics was computed separately for each image. The resulting independent statistics were then averaged across all the images in the database (300 images), and are shown in Fig. 2a. Note that, on the average, more than 90% of the patches in an image have 9or more other similar patches in the same image at the original image scale (‘within scale’). Moreover, more than 80% of the input patches have 9or more similar patches in 0.41=1.25−4of the input scale, and 70% of them have 9or more similar patches in 0.26=1.25−6of the input scale. Recurrence of patches forms the basis for our singleimage super-resolution approach. Since the impact of super-resolution is expressed mostly in highly detailed image regions (edges, corners, texture, etc.), we wish to eliminate the effect of uniform patches on the above statistics. Therefore, we repeated the same experiment using only 25% of the source patches with the highest intensity variance. This excludes the uniform and low-frequency patches, 1www.eecs.berkeley.edu/Research/Projects/CS/vision/grouping/segbench 2Distances between patches were measured using gaussian-weighted SSD. Note that textured patches tend to have much larger SSD errors than smooth (low-variance) patches when compared to other very similarlooking patches (especially in the presence of inevitable sub-pixel misalignments). Thus, for each patch we compute a patch-specific ‘good distance’, by measuring its (gaussian-weighted) SSD with a slightlymisaligned copy of itself (by 0.5 pixel). This forms our distance threshold for each patch: Patches with distance below this threshold are considered similar to the source patch. (a) Classical Multi-Image SR (b) Single-Image Multi-Patch SR Figure 3: (a) Low-res pixels in multiple low-res images impose multiple linear constraints on the high-res unknowns within the support of their blur kernels. (b) Recurring patches within a single low-res image can be regarded as if extracted from multiple different low-res images of the same high resolution scene, thus inducing multiple linear constraints on the high-res unknowns. maintaining mostly patches of edges, corners, and texture. The resulting graphs are displayed in Fig. 2b. Although there is a slight drop in patch recurrence, the basic observation still holds even for the high-frequency patches: Most of them recur several times within and across scales of the same image (more than 80% of the patches recur 9or more times in the original image scale; more than 70% recur 9 or more times at 0.41 of the input scale, and 60% of them recur 9or more times in 0.26 of the input scale.) In principle, the lowest image scale in which we can still find recurrence of a source patch, provides an indication of its maximal potential resolution increase using our approach (when the only available information is the image itself). This is pixel-dependent, and can be estimated at every pixel in the image. 3. Single Image SR – A Unified Framework Recurrence of patches within the same image scale forms the basis for applying the Classical SR constraints to information from a single image (Sec. 3.1). Recurrence of patches across different scales gives rise to ExampleBased SR from a single image, with no prior examples (Sec. 3.2). Moreover, these two different approaches to SR can be combined into a single unified computational framework (Sec. 3.3). 3.1. Employing in-scale patch redundancy In the classical Multi-Image Super-resolution (e.g., [12, 5, 8]), a set of low-resolution images {L1, ..., Ln}of the same scene (at subpixel misalignments) is given, and the goal is to recover their mutual high-resolution source image H. Each low resolution image Lj(j=1,...,n)is assumed to have been generated from Hby a blur and subsampling process: Lj=H∗Bj↓sj, where ↓denotes a subsampling operation, sjis the scale reduction factor (the subsampling rate) between Hand Lj, and Bj(q)is the corresponding blur kernel (the Point Spread Function – PSF), represented in the high-resolution coordinate system – see Fig. 3a. Thus, each low-resolution pixel p=(x, y)in each low-resolution image Ljinduces one linear constraint on 351 (b) Single-image multi-patch SR. Fig. 2.4 Contributions of patches (Bi) recurring in LR images (Li) to the process of HR image (H) reconstruction [GBI09]. Fig. 2.5 Example of patch recurrence across scales resulting in implicit LR-HR pairs identification. 2.2.3.1 (A) Multiple-images super-resolution in spatial domain The multiple-images-based approaches comprise five families of techniques: (A.1) Direct. The theoretical basis of these techniques is the non-uniform sampling theory which allows for the reconstruction of functions from samples taken at nonuniformly distributed locations.These methods follow the pure classical SR approach presented in Fig. 1.1. Given a set of LR observations, one of the LR images is chosen as the target image and the others are registered against it. Following, the target image is scaled up by a specific scaling factor and the other LR images are projected (that is scaled and warped or shifted) into the reference grid using the data obtained from the registration stage. Then, the HR image is generated by fusing all the images together. The missing data are created by means of non-uniform interpolation. Finally, an optional deblurring/refinement kernel may be applied to the result. The fusing process is usually implemented 25 CHAPTER2.–IMAGE ENHANCEMENT USING SUPER-RESOLUTION Spatial Domain (A) Multiple Images (A.1) Direct Shift and Add Non-parameteric (A.2) Stochastic Maximum Likelihood Maximum a Posteriori (A.3) Set Theoretic Projection onto Convex Sets Ellipsoid Method (A.4) Iterative Back Projection (A.5) (Iterative) Adaptive Filtering (B) Single Image(B.1) Learning-based Feature Pyramid Belief Networks Projection Neaural Networks Manifold Tensor Compresive Sensing (B.2) Reconstruction Primal Sketches Gradient Profiles Fields of Experts Frequency Domain Fourier Wavelet Super-resolution Fig. 2.6 Excerpt from the SR taxonomy proposed recently in [NM14]. as a weighted sum, thus, this approach is usually referenced to as shift and add [FREM03, FREM04b, FREM04a, FM06, AIAOM13]. More recent direct methods carry out the registration and fusion at patch-level. That is, the LR images are first divided into patches. Then every patch of the reference image is compared to a set of patches, including the corresponding patch and some patches in its vicinity, in the remaining LR images. Based on the similarity of these patches and their distances from the current patch, a weight is associated to every patch indicating the contribution of this patch in the process of producing the output HR patch. This approach helps to handle occlusion, local and complex motions (e.g. facial expressions change). The family of methods that deploy this type of processing is referenced to as non-parametric direct SR [PE09, PETM09, TMPE09, BBV10, CCL10, HS12, ZLRG12]. (A.2) Stochastic. Methods from this family consider HR image and motions among low resolution observations as stochastic variables. These variables are used to construct a Bayesian framework which is used to estimate the unknown values by minimization of a cost function. The statistical approaches vary in the ways: they treat the degradation process, the LR/HR observations priors (defined as our beliefs about the sensitivity and specificity of a process/variable) they use and the inference method used by the framework. The two most widely used statistical approaches to SR are the maximal likelihood (ML) [CKK+96, SS94, SS96, LSZ97, EHO01, HYCC07] and the maximum a posteriori (MAP) [GH97, HBA97, BS99, RC01, TM10, CNYA12, KI12, YZS12, ZZLH12]. The ML ap26 2.2 SUPER-RESOLUTION BASICS proach only considers the relationship among the low resolution observations and the original high resolution. In contrast, the MAP approach incorporates also the prior image model to reflect the expectation of the unknown high resolution image. In most cases, the use of MAP-based approaches is encouraged as the preferred one. In mathematics, a process of introducing additional information in order to solve an ill-posed problem is referred to as regularization. Thus, many of the stochastic approaches are also referred to as regularization-based SR. More on stochastic approaches to SR can be found in [Mil10]. (A.3) Set theoretic reconstruction. Set theoretic reconstruction assumes the SR process to be a feasible linear optimization problem with rational data that can be solved by generating a sequence of sets of possible solutions whose cardinality uniformly decreases. The subsequent sets are subsets of the previous set containing the desired images, that is the images that fulfill (or intersect with) a subsequent constraint convex set. Defining constraints in terms of convex sets is flexible and allows incorporating even non-linear and non-parametric constraints or priors. In the context of SR, the most important set theoretic reconstruction algorithms have been the ellipsoid method and the projections onto convex sets (POCS) [SO89, EF97, HKK97, PA98, BS00, CTH03]. In the POCS method, the subsequent sets are determined by finding the intersection points of the constraint and the current solution set. In the ellipsoid method the subsets are ellipsoid bound by the constraints set. More on set-theoretic approaches to SR can be found in [Mil10]. (A.4) Iterative back projection (IBP). IBP is similar in its approach to the back-projection used in computer tomography image reconstruction. This method creates an initial guess of the high resolution image which is refined by iteratively projecting the difference (the residual) between the observed low resolution images and the simulated (back-projected) low resolution images [IP92, CD00, ZP00, ZRAP01, JWB03, BEZN04, BEZN05, YB06, DHWG07, GDAG09]. (A.5) (Iterative) adaptive filtering. This SR approach is most appropriate when the motion model and the point spread function model commute. The filters usually deployed for SR are the Kalman and Wiener adaptive filters. These methods are in effect linear minimum mean square error estimators and do not allow inclusion of non-linear a priori constraints. Although limited in terms of performance 27 CHAPTER2.–IMAGE ENHANCEMENT USING SUPER-RESOLUTION understood as the output quality, these approaches offer very fast execution and have been mainly used in SR of video sequences [EF99b, EF99a, EF99c, CB07, CB06, CB08, KKC+10]. Each of the above-presented approaches has its advantages and drawbacks. The direct approach is basic and intuitive method of super-resolution with its implementations having relatively low computational complexity. On the other hand, the direct approach assumes that the blur and noise characteristics are identical across all low resolution images. These methods rely on very accurate registration between images and thus are very sensitive to any errors produced in the image registration stage which can be easily propagated to the SR kernel. Moreover, the step-by-step forward approach with separation of the motion estimation and image fusion steps does not guarantee optimality of the restoration. The IBP is also understood intuitively and easily. However, this approach preserves the ill-posed nature of the inverse problem. Thus, there are multiple possible solutions and the algorithm either converges to one or it may oscillate between them. Although, the latter can be dealt with by incorporating a priori knowledge about the solution, the inclusion of the a priori constraints is not easily achieved in the IBP method and usually requires additional regularization terms. The possibility to conveniently include the a priori information, along with their theoretical simplicity, is the biggest advantage of the set theoretic methods. Nevertheless, these methods have high computational cost, slow convergence and lead to a solution that is not guaranteed to be optimal. Moreover, the use of the POCS method has the disadvantage of producing non-unique solution that depends on the initial guess. Practical implementations of POCS also have problems in dealing with the reconstruction of discontinuities of edges and details. Compared with other approaches, stochastic SR offer better robustness, flexibility in modeling noise characteristics and easiness of including the a priori knowledge of the solution in the probability priors. In certain cases (e.g. the noise process is white Gaussian, and MAP estimation with convex energy functions), this method ensures the uniqueness of the solution. The uniqueness of the solution allows the use of efficient gradient descent methods and is necessary for optimal reconstruction. Recently, approaches looking to combine the properties of the stochastic approaches with other families have been undertaken. The ML/MAP-POCS joint SR estimates the HR image by minimizing the ML/MAP cost functional while enforcing the containment of the solution within the intersection of the convex constraint. The advantage of the combined approach is that it ensures a single optimal solution, which is not the typical case when POCS approach is used on its own, while allowing any kinds of constraints and priors, even the ones that may present as impossible for purely stochastic approaches. 28 2.2 SUPER-RESOLUTION BASICS 2.2.3.2 (B) Single-image super-resolution in spatial domain Single-image SR is based on interpolation and hallucination of the high-fidelity details (usually high frequency part of the image) of the HR image. Single-image SR approaches comprise two families of techniques: (B.1) Learning based. Learning-based approaches [BK00a, CZ01, BK00b, PRZ03, WT03, JG05, ZPJ08, KCR09, ZH12, ZHD12] use pairs of LR and HR training inputs to create a dictionary that relates the low-resolution data (or its features) with their corresponding high-resolution representation. During run-time, the given LR test input is used to traverse the dictionary in order to find the closest match (or a set of closest matches) and retrieve the associated HR-data to be used in the hallucination. The learning-based approaches are further divided based on the way the learning is carried out, the structure of the dictionary, the way the search is carried out and how the matches contribute to the hallucination process. The most recent work [NM14] defines seven sub-families of example-based SR approaches, namely: features pyramids, belief networks, projection, neural networks, manifold, tensors, and compressive sensing. Details on each of these approaches can be found in [NM14]. (B.2) Reconstruction based. The a priori knowledge is generalized in order to encapsulate techniques for reconstruction of primitives (and/or shapes), or statistics that can be used to carry out the reconstruction. The SR is carried out by applying the adequate reconstruction technique to each of the identified primitives (like edges, ridges, corners, etc.) or using the provided statistics to form a gradient based constraint to be applied to the reconstruction process [SZTS03, RB05, Fat07, SLJT08, SSXS08, XSW09, HFL10, XSW10, SSXS11]. Some of the single-image methods allow their execution in a multi-images context. These methods jointly exploit the information learned from a given high resolution training data set, as well as that provided by multiple low resolution observations. Also, the statistical approach can be incorporated in the example-based approach where the found patches are used as the prior image model and then merged into a MAP framework cost function in order to arrive at the closed form solution of the desired high-resolution image. 29 CHAPTER2.–IMAGE ENHANCEMENT USING SUPER-RESOLUTION 2.2.4 Super-resolution of video sequences The SR processing of video sequences is usually required to be dynamic, that is to produce an output sequence that (at least) maintains the number of frames and the frame rate of the input sequence. The scenario of dynamic SR imposes high constraints on the execution time. In practice, the choice of algorithms for application of SR on video sequences is limited to approaches that offer a favorable trade-off between performance, execution time and resources requirements. The existing approaches to SR of video sequences can be classified into the following four categories: (i) the sliding-window-based approach [SKR07, NHBS07, NSLZ07, PJ07], (ii) the sequential approach [EF99c, EF99b, FEM06, CB07], (iii) the simultaneous approach [BS99, ZM07, AMK03], and (iv) the learning-based approach [BBM03, DKA04, KHX+06]. The first three approaches require multi-images/frames context, while the learning-based SR is capable of executing in a single-image/frame context. The learning-based approach to SR of video sequences is straightforward: subsequently apply learning-based SR of still images on each of the frames of the sequence. In the case of execution in multipleimages context, this straightforward approach is not the preferred one. Its use is not efficient as it leads to reloading of data and does not (even try to) take advantage of the relationships between the subsequent frames of the sequence. As for now, software implementations offer performance that is not sufficient for real-time execution, even though the single-image context shows inherent parallelism that allows simultaneous multi-frame SR. The available implementations of learning-based SR of video sequences have not been implemented in hardware. Considering that in the case of video sequences availability of multiple LR images/frames is assured, the multi-image SR techniques constitute the preferred approach, especially in the context of hardware implementation. 2.2.4.1 Multi-frame approaches to SR of video sequences As aforementioned, the multi-image approaches require a number of frames to execute. This set of frames, called the frame window (or the working window), forms the execution context. The frame window comprises two types of images: the target image that is being super-resolved and a number of reference images (or reference frames (RF or RFnr)). Frame window usually comprises the immediate neighbors of the to-be-processed frame as illustrated in Fig. 2.7. A consequence of using a multi-image context of execution is the requirement of registering the frames belonging to the current working window. In the context of SR of video 30 2.2 SUPER-RESOLUTION BASICS Target image Reference image 1 Reference image 2 Reference image 3 Reference image 4 Set of #RF reference images Frame window (set of #RF + 1 images) Fig. 2.7 Example of frame window comprising 5 images (4 reference images or reference frames (RF); RF =4) and one target image). sequences there are two main approaches to image registration and (effectively) context updating being deployed, namely: the anchored and the progressive approach. In the anchored approach, one of the frames of the context is chosen as the target frame and the other frames misalignments are registered in relation to this frame. In case of the progressive registration, the current frame is registered in relational to its immediate temporal neighborhood made up by the frames that precede it. In practice, multi-image SR for video sequences is implemented using one of the three following approaches illustrated in Fig. 2.8: Sliding-window-based SR approach. Deploys anchored registration over a set of consecutive low resolution frames which are later combined producing one high resolution image frame. The execution context is moved across the input frames to produce successive high resolution frames sequentially. The major drawback of this approach is that the temporal correlations among the consecutively reconstructed high resolution images are largely unexploited. Also, memory occupancy associated with this approach is high as all required data from the context need to be readily available from the memory. Sequential SR approach. The sequential SR approach tries to exploit the temporally correlated information provided by the established high resolution images. Correct identification and exploitation of temporal correlation of HR images using their LR observations is a challenging problem. In an effort to do so, progressive registration and previous HR frame can be used in tandem during the process of estimating the current HR frame. This allows to reduce the cardinality of the execution context while still including (by means of propagation) information from multiple frames. Reduction of the execution context usually leads to significant reduction of memory occupancy and computational complexity. On the other hand, the propagation may lead to er31 CHAPTER2.–IMAGE ENHANCEMENT USING SUPER-RESOLUTION High resolution sequence Low resolution sequence (a) Sliding frame window approach. High resolution sequence Low resolution sequence (b) Sequential approach. High resolution sequence Low resolution sequence (c) Simultaneous approach. Fig. 2.8 Approaches to super-resolution of video sequences executing in multiframe context. Based on [TM11]. 32 2.2 SUPER-RESOLUTION BASICS rors escalation and drift if robust fusion and outliers elimination mechanisms are not deployed. Simultaneous SR video approach. This approach tackles the problem of simultaneous reconstruction of multiple high resolution images. The typical approach is to impose temporal smoothness constraint on the prior image model and produce multiple images from the context at once. Alternatively, multiple contexts could be processed in parallel i.e. multiple sliding windows could be selected and used to reconstruct multiple high resolution frames simultaneously. In both cases, simultaneous SR poses very high memory occupancy requirements, since high resolution and low resolution images from multiple contexts are required to be readily available from memory at the same time throughout the reconstruction process. It should be noticed that on-line video SR is significantly more demanding in terms of performance, execution speed and memory occupancy than SR of still images [Goh13]. The highest requirements are presented by the simultaneous SR where frames from multiple contexts need to be stored and/or hardware resources have to be duplicated (or pipelined) if the computations are to be parallelized. On the other end of the requirements spectrum is the sequential approach, whose requirements in terms of resources and computational power seem to be the lowest. Nevertheless, in order to be successful, this approach requires algorithms that allow inclusion of a priori knowledge to regularize the SR process. Thus, stochastic/regularization-based approaches are preferred over direct/IBP-based ones. The former are known to be mathematically involved and suffer from high execution times. In the context of real-time on-line dynamic SR of video sequences less involved (and thus faster) algorithms are preferred. The interpolation-based, namely, the direct and IBP methods have found most success in meeting real-time constraints. 2.2.4.2 Dynamic super-resolution using sliding frame window Dynamic SR using the sliding-frame-window approach is provided by moving the frame window across the video sequence and performing the SR process once for each frame of the sequence. After the processing of one frame is terminated, the least recently loaded frame is discarded, the remaining frames are shifted, and a temporarily subsequent (next) frame is loaded. For a sequence composed of nframes, processing of the whole sequence is defined as carrying out the super-resolution process for each frame t:t≤nwith the set κof frames contributing in the fusion process being limited to a size κ≤n(in practice 3≤κ≤17) and updated for each of the frames being processed. The set κrepresents the 33 CHAPTER2.–IMAGE ENHANCEMENT USING SUPER-RESOLUTION Image Registration Domain Area-Based Correlation-Like Methods Mutual Information Methods Fourier Methods Phase-Correlation Matching in DCT Domain Matching in Wavelet Domain Feature-Based Methods Using Spatial Relations Methods Using Invariant Descriptors Relaxation Methods Pyramids and Wavelets Acquisition Registration Basis Extrinsic (introduced foreign object) Intrinsic (image information) Non-Image Based (sensor coordinates) Modality Mono-Modal Multi-Modal Other Application Multi-View Multi-Temporal Multi-Modal Scene-to-Model Transformation Dimensionality 2D 3D Model Translantional Rigid Affine Projective Random Interaction (user role) Interactive Semiautomatic Automatic Mapping Global Local Hybrid Fig. 2.11 Classification of the most popular approaches to image registration. allowing efficient implementation. If uncorrected variations have not been eliminated by the feature space and similarity metric, then the search for the optimal match is also made more difficult, since there are more likely to be several local optima and a less monotonic space [Bro92]. This problem is solved to some extent by feature detection step and similarity metric. Feature-based class encapsulates the algorithms that work by extracting a sparse set of features in the images which are then matched against each other. To provide efficient registration these methods rely on the existence of salient points in the input. The problem with feature-based matching is that, typically, good, matchable features such as corner points are sparse while poor easily mismatched features such as edges, are denser. Even when reasonably unique features are available, establishing the correct correspondences can be problematic, especially for cases when occlusion occurs. Image registration of images with features occlusion is known to be the most challenging scenario for feature-based methods, and is likely to lead to matching errors [BB95]. 40 2.3 IMAGE REGISTRATION Area-based approaches, sometimes referenced as correlation-like or template matching, are less sensitive to these problems. These techniques do not rely on the presence of salient points (called control points), rather, they consider areas of the image as features used in matching—shifting the emphasis from feature detection to the feature matching step. The matching step is carried out by directly minimizing point-to-point dissimilarities for predefined sized regions/patches. The variable window sizes can be used near occlusion boundaries to handle multiple motions. This matching and minimization in most cases are carried out in the spatial domain using directly the illumination values (so-called direct methods). However, it is possible to carry out the process using representations in the frequency domain which offers advantages in noise sensitivity and computational complexity. The limitations of the area-based methods originate from their basic idea: the predefined size/shape window (rectangular windows are most often used) suits the registration of images which locally differ only by a translation. If images are deformed by more complex transformations, this type of window is not able to cover the same parts of the scene in the target and sensed images. Feature detection is not implemented by the area-based methods rending them more sensitive to acquisition conditions (and noise). Hence, a more robust and noise tolerant similarity metric are required to be used in these algorithms. Some of the aforementioned limitations are alleviated by carrying out the processing in the frequency domain. Frequency domain is less susceptible to differing conditions of illumination since illumination changes are usually slow varying and therefore concentrated at low-spatial frequencies. Similarly, the techniques using frequency domain are relatively scene independent and useful for images acquired from different sensors since it is insensitive to changes in spectral energy. The scene independence is further strengthened in the case of phase-correlation methods that use only the phase information making the correlation measure invariant to linear changes in brightness. By using the frequency domain, the Fourier methods achieve excellent robustness against correlated and frequency-dependent noise. On the other hand, if the images have significant white noise, noise which is spread across all frequencies, then the location of the peak will be inaccurate since the phase difference at each frequency is corrupted. In this case, use of the spatial cross-correlation is better. Also the Fourier methods are applicable only for images which have been at most rigidly misaligned. To summarize, the use of feature-based methods is recommended if the images contain enough distinctive and easily detectable objects. This is usually the case of applications in remote sensing and computer vision. Moreover, feature-based approaches have also the advantage of being more robust against scene movement, and are potentially faster, if imple41 CHAPTER2.–IMAGE ENHANCEMENT USING SUPER-RESOLUTION mented in the right way. It should be noted that registration methods using simultaneously both area-based and feature-based approaches have recently started to appear [HZ07]. 2.3.4 Transformations In the context of image registration a transformation can be defined as the mathematical model that maps pixel coordinates from one image to another. Most registration techniques involve searching over the space of transformations of a certain type to find the optimal transformation for a particular problem. Thus, before the registration and alignment of images happens, it is necessary to establish the basis of the possible search spaces and applicability of these transformations. These transformations can be expressed by a variety of models ranging from simple 2D translations to complex 3D movements of elastic body models. This section will focus only on the models and mappings which are of relevance for the task of 2D image registration. 2.3.4.1 Transformations primitives describing movement in 2 dimensions Let x2ddenote a geometric primitive of a point in 2-D space x2ddefined using a pair of values (x,y)as x2d= (x,y)∈R2. In computer graphics a point in 2-D space is usually represented in homogeneous coordinates as ˜x2d= (˜x,˜y,˜w)∈P2, where P2=R3−(0,0,0) is called the 2D projective space. The conversion between the Cartesian and homogeneous coordinates is described as ˜x2d= (˜x,˜y,˜w) = ˜w·(x,y,1) = ˜w·x,(2.1) where x= (x,y,1), created by extension of the coordinates number with third coordinate equal to ‘1’, is called the augmented vector. Homogeneous coordinates have the advantage of allowing all of the transformations to be expressed using multi-dimensional matrices with mapping parameters. Having established the representations of the 2D point, it is necessary to define and describe the basic set of primitives that carry out the transformations listed in Fig. 2.12. This description is based on the contents from [Sze10] covering only the transformations that are used most extensively. A graphical illustration of these primitives is presented in Fig. 2.12. Basic properties of the introduced below primitives are presented in Table 2.2. Translation. Translation is a pure shift in 2D space which preserves the orientation (size, 42 2.3 IMAGE REGISTRATION y x translation Euclidean similarity projective affine Fig. 2.12 Basic set of 2-D planar transformation [Sze10]. scale and shape) of the object. Having δ2drepresent the translation distance, translation can be written as x′ 2d=x2d+δ2d. Euclidean. Also known as 2D rigid body motion. It can be seen as a sequence of two transformations: rotation and translation. It preserves the inter-pel distances, shape and scales but allows changes of its orientation. Denoting rotation as R, rigid body motion is modeled as x′ 2d=R·x2d+δ2d. Similarity. Extends the Euclidean model by allowing the scale of the object to be changed. This process is modeled as x′ 2d=s·R·x2d+δ2d, where sis an arbitrary scale factor. One thing to note is that the similarity transform still preserves angles between lines. This translation is sometimes referenced as scaled rotation. Affine. Affine transformation is modeled as x′ 2d=A2×3·x2d, where x2dis such a representation of x2din the homogeneous coordinates space that x2d= (x,y,1)and A2×3is an arbitrary 2×3 matrix with mapping parameters. Affine transformation preserves the parallelism of lines. Projective. Also known as a perspective transform or homography. Allows even higher degree of freedom (DoF) than the affine transformation. Denoting ˜ H3×3an arbitrary 3×3 matrix with mapping parameters, the projective transformation is modeled as ˜x′ 2d=˜ H3×3·˜x2d. Straight lines remain straight after the transformation. 43 CHAPTER2.–IMAGE ENHANCEMENT USING SUPER-RESOLUTION TABLE 2.2 Hierarchy and properties of 2D coordinate transformations. An extract from [Sze10]. Transformation Matrix dim # DoF Preserves Translation 2×3 2 orientation Euclidean(rigid) 2×3 3 lengths Similarity 2×3 4 angles Affine 2×3 6 parallelism Projective 3×3 8 straight lines 2.3.4.2 Global and local transformations mapping Each of the presented transformations representing motion models can be constructed or mapped by considering different extent of the available support (control points (CP; featurebased methods) or pels (area-based)). The extent of support considered in the process of mapping parameters estimate determines whether each individual technique is classified as local or global. Global models use all available support for estimating one set of the mapping function parameters which accounts for all the variations between the images. In other words, a single equation is considered valid for the entire image, and it is used to describe the motion over the entire visual field. As these variances become more local, it will become progressively more difficult for a global point-mapping method to model all of the changes using one global transformation. In presence of significant and multiple local geometric variances or local 3D features observed from different viewpoints (resulting in different 3D-to-2D projections) global methods may fail to account for the misalignment between the images in a satisfactory way. In these cases, more than one transformation, which limit their support by using only a local neighbourhood, would be preferable. To accomplish that, local methods use not one, but a set of mapping parameters that vary across the different pieces of the support in order to account for different models of (local) variations. In other words, the mapping transformation is no longer a single mapping with one set of parameters independent of position. This allows local methods to be more powerful and handle many (combinations of) distortions that global methods cannot. Several authors have shown the superiority of the local or at least locally sensitive registration methods above the global ones in situations where local variations are present [YMB11, BAA05, CLT+08]. In all cases, there is a trade-off between the power of these methods and their corresponding computational and implementation cost. As illustrated in Fig. 2.13 and Fig. 2.14 the choice of mapping of the transformations directly influences the size and complexity of the search space. In practice, this manifests itself as execution time and memory requirements. Local methods are well-known to have the largest and most complex search 44 2.3 IMAGE REGISTRATION 1 DoF 2 DoF 4 DoF 8 DoF 16 DoF 1 DoF 2 DoF 4 DoF 8 DoF 16 DoF 1 DoF 2 DoF 4 DoF 8 DoF 16 DoF 4 combinations 8 combinations 32 combinations 1 10 100 1000 10000 100000 1000000 optical flow (1x1) patch (2x2) patch (4x4) patch (8x8) patch (16x16) global (64x64) Comparisions in matching [operations] Fig. 2.13 Number of matchings as a function of template size, degrees of freedom and candidate combinations. Optical flow uses only 1 DoF. 1 DoF 2 DoF 4 DoF 8 DoF 16 DoF 1 DoF 2 DoF 4 DoF 8 DoF 16 DoF 1 DoF 2 DoF 4 DoF 8 DoF 16 DoF 4 candidates 8 candidates 32 candidates 10000 100000 1000000 10000000 optical flow (1x1) patch (2x2) patch (4x4) patch (8x8) patch (16x16) global (64x64) Similarity evaluation [operations] Fig. 2.14 Number of similarity index evaluations as a function of template size, degrees of freedom and candidates number. Optical flow uses only 1 DoF. 45 CHAPTER2.–IMAGE ENHANCEMENT USING SUPER-RESOLUTION spaces. On the other hand, is many cases local registration allow pieces of the input to be registered without influencing other portions which have already been matched leaving more room for parallelization. Moreover, for many registration problems, both local and global distortions exist, and it is sometimes useful to take a hierarchical approach in finding the optimal transformation. In this case, the regions with local variations are registered using local methods with the rest of the support being described using one global transformation mapping. This may lead to significant reduction in computational complexity and implementation cost [BAA05, CLT+08]. 2.3.5 Search space Let search space denote the set of created by superposition of all allowable transformations with all their respective allowable parameters values in the defined feature space. Then, the aim of feature matching is to find the combination of transformation and its parameters from the search space which carries out the registration in an optimal way. To this end it is necessary to devise methods that (i) identify the transformation (and its parameters), and (ii) measure the quality of results on the preselected features in a quantitative way. The former and the latter methods are called the search strategy and the similarity criterion, respectively. The matching operation in spatial domain that traverses the set of possible transformations computing the corresponding similarity test results for a set of allowed templates is known as ‘template matching’ [TK08]. A search iteration is carried out by placing the template over the reference at the test location and computing the similarity metric. The similarity metric is constructed in a way that it returns higher values for templates placed over a reference that is more similar (has smaller differences between the corresponding intensities). The optimal transformation is the one for which the similarity measure reaches its global optimum (maximum in case of similarities). Finding the global optimum of similarity index is a multi-dimensional optimization problem, where the number of dimensions corresponds to the degrees of freedom of the class of the expected geometrical transformation. For example, in the case when the only allowable transformation is translation, then the search space is the set of all translations over the range of parameters (in this case displacement/horizontal and vertical shifts). Under the assumption of global mapping, the matching process would be carried out by iteratively translating the sensed image by a value from the defined range and quantifying how the template and the reference pels correlate. 46 2.3 IMAGE REGISTRATION The only way to assure that the displacement that results in the optimal transformation can be found, that is the global maximum of correlation, is for the matching to be carried out exhaustively over all possible displacements. 2.3.5.1 Computational complexity as a function of search space The computational intensity associated with template matching is determined by the size of the search space, type of transformations mapping and the computational intensity of the similarity measure estimation. The size of the search space is determined by the number of degrees of freedom of the (expected) transformation model and the range and grain of the allowable values of the parameters. The simplest solution of trying out all allowed combinations from the search space results in the template being transformed (translated, rotated, scaled, etc.) for each possible transformation of interest. This affects the number of combinations that need to be tested, skyrocketing with the increase of the number of degrees of freedom and/or search range. The type of mapping determines the number of sets of parameters that have to be computed to align the images being registered. This translates to the number of templates for which the matching process has to be carried out. In global mapping only one template is used. The number of templates used by local methods corresponds to the number of regions that the target image has been divided into. In case of area-based correlation-like image registration three types of templates can be defined based on the number of pels that they contain: a singular pel, a full frame and a block of pels. The choice of granularity of the template has a significant impact not only on the dimensionality of the search space but also on the quality of the results. Using a coarser-grain template reduces the number of templates (towards the limit of one) at the cost of potential results quality degradation should local distortions be present. The use of a finer-grain template has the disadvantage of increasing the number of templates used in matching introducing additional operations in the similarity measure computation. Piece-wise (local) search of transformations may result in a blocking effect, characterized by the appearance of artifacts at patch boundaries. Effective treatment of these variations would requires additional processing [Ric04]. Consequently, methods using fine-grain templates tend to have the largest and most complex search spaces. An exception to that rule is the case of using only one pel. This finest-grain template is an interesting case as it allows to reduce the complexity of the motion model to purely translational. Significant reduction in computational complexity, relatively low memory requirements, and intrinsic parallelism have made optical flow one of the preferred ways for hardware implementations in low cost systems. 47 CHAPTER2.–IMAGE ENHANCEMENT USING SUPER-RESOLUTION On the other hand, accurate estimation of the true transformations becomes more and more challenging for fine-grain templates. The change from coarse-to-fine template is carried out by limiting the extent of support used in matching. A well known trade-off of reduction in support is the loss of ability to accurately represent structural information. This makes fine-grain templates more susceptible to producing false motion estimation, especially when applied on homogeneous regions of the image [CLT+08]. To some extent this problem can be mitigated by using more robust (and more computationally intensive) similarity measures. 2.3.6 Search strategies In most cases, trying out all of the combinations from search space with high cardinality is too time-consuming to be practical. Due to its high computational intensity, the exhaustive search is rarely implemented for registration of images using models with more than 2 degrees of freedom. In practice, exhaustive search is used as the baseline to evaluate the quality of results of other algorithms with reduced computational intensity. The approaches to reduce the computational complexity of template matching aim at lowering the number of tests or making the test less expensive. The latter is implemented by using a mathematically less involved similarity metric and/or allowing early termination at various stages of computation should a threshold value be reached. The former can be enforced by limiting the set of allowable transformations (arbitrarily or based on an a priori knowledge) and/or using heuristic methods for matching. Reduction of transformations model limits the allowable transformations to those with lower number of degrees of freedom. Heuristic algorithms encapsulate a set of methods that limit the number of carried out tests by constraining the search range of parameters for which the matching is carried out. In case of the area-based matching, the space described by the set of combinations of parameter values of the search range is referenced as the search area. The search area in spatial domain is described by its (relative) span and the size of the template. An example of such a search area used in piecewise area-based correlation-like registration (block-matching) is presented in Fig. 2.15. This particular case assumes translational model over a range of (search radiusx,search radiusy) displacements and a square template composed of template sizex×template sizeypels. The pattern used to traverse the search range/space is known as the search strategy. Heuristic algorithms can be classified based on the deployed search strategy. It is difficult to present a complete classification of search strategies without sacrificing presentation clarity. Each search strategy has its advantages, disadvantages, sometimes limited domains 48 2.3 IMAGE REGISTRATION search radiusy search radiusy template heighty search area heighty search radiusxsearch radiusx template widthx search area widthx Fig. 2.15 Search area structure used in block matching of a template assuming translational model over a range of displacements. Block-Matching Single-Scale Scene-Adaptive Adaptive Search Area Dynamic Search Range Algorithms Dynamic Search Window Algorithms Block-Based Gradiend Descent Search Predicitive Adaptive Rood Pattern Search (New) Predictive Search Area Adaptive Template Size Variable Block Size Algorithms Scene-Fixed Full/Exhaustive Search Spiral Search 2D Logarithmic Search (New) Three-Step Search Four-Step Search One at a Time Search Binary Search Orthogonal Search Cross Search Conjugate Direction Search Parallel Hierarchical One-dimensional Search (New) Diamond Search Hexagon Search Simple and Efficient Search Multi-Scale Hierarchical Search Hierarchical Partial Distortion Search Hierarchical Block Matching Algorithm Pel Decimation Technique Adaptive Pel Decimation Technique Fig. 2.16 Classification of block-matching implementations based on the used search strategy. 49 CHAPTER2.–IMAGE ENHANCEMENT USING SUPER-RESOLUTION consecutive LR frames are assigned different HR samples values. XSR applied on T regions is able to undo smart down-sampling, as it knows the inverse of sample selection scheme. Flat regions are upscaled using interpolation and MT regions are super-resolved in the regular manner. The algorithm was mapped onto the Integrator CM922T-XA10 development board. The board includes an Excalibur XA10 device and off-chip memory. The former device comprises a programmable logic device, used to accelerate inverse discrete cosine transform (IDCT) and handle off-board communication, and an ARM9 embedded processor that manages the rest of the tasks. The described system was capable of super-resolving QCIF (176x144 pixels) resolution frames to CIF (352x288 pixels) resolution at the rate of 15 frames per second. The memory access latency and bus contention have been identified as the system bottlenecks. In [ABCC09] Angelopoulou et al. present a FPGA implementation of a super-resolution image reconstruction based on iterative back projection that executes in a multi-frame context. In this approach, additional details are reconstructed based on exploitation of sub-pixel shifts caused by warping. In order to facilitate parallelization and minimize the execution context memory occupancy, the processing is done at pixel level by means of weighted mean optical flow, referred to as weight based (picture elements) merging. Weights estimation is based on the inter-frames motion estimations for pel matching. The implementation reaches an operating frequency of 80 MHz, which the authors claim is sufficient to allow real-time execution outputting 25 VGA 2x super resolved frames. Nevertheless, the output quality is compromised due to the low number of implemented iterations (up to 10) limited by the available resources. The main weaknesses of this approach are its high memory requirements and the fact that, in order to output a super resolved image, multiple passes through the hardware are required. The design bottleneck was identified to be the triple buffering memory access scheme. In [BB08] Bowen et al. present a similar implementation of IBP algorithm targeting a development board hosting a FPGA device. The architecture onto which the algorithm was mapped is shown in Fig. 2.19. Motion estimation data, pixel values, and weights of aforementioned LR frames are loaded and used to fill-in the HR grids. Those two grids are merged with previous frame initial approximation to form initial approximation for the current frame. Missing pixels are reconstructed by means of modified nearest neighbor interpolation forming an initial approximation of the HR image. This approximation is later fed to the first iterative stage modules. Iterative stage modules carry out the refinement process using values from previous iterations, original pixel, and weight related parameter in similar fashion to the aforementioned weight based merging. Output of the last iterative 56 2.4 SUPER-RESOLUTION HARDWARE IMPLEMENTATIONS Fig. 2.19 System architecture deployed by [BB08]; source [BB08]. stage is considered the super-resolved pixel. The authors claim, that, with a loaded pipeline, their approach is capable of producing one super-resolved pel per cycle. The described architecture, implementing 10 iteration stages, was mapped onto a Xilinx XC2V6000 FPGA device reaching a frequency of 58 MHz. In this configuration the system was capable of super-resolving 61 CIF formatted LR images to 1280x720 pels per second. Nevertheless, a satisfying quality level requires at least 20 iteration stages, being out of reach for the target device. The number of implemented iteration stages is said to be limited only by the available on-chip memory. In [STdA13] Singla et al. implement a modified version of the NISR on a low cost NoC-based MPSoC platform comprising up to 4 Xilinx MicroBlaze soft-core processors. The processors are equipped with a 64 kB of local memory and access to DDR3 controller over an Arteris FlexNoc 2D-Mesh NoC. Each core executes its own copy of the algorithm that applies SR processing on a (statically assigned) chunk of the input frame loaded from off-device random access memory (RAM). Mapping onto Xilinx Spartan-6 LX45T FPGA results in utilizing 88% of the available slices (∼6003 slices or ∼38420 LUTs) and 97% of BRAMs (∼2025 Kb) while reaching operating frequency of almost 53 MHz. For this configuration a complete SR (including ME) of a QCIF frame to CIF (2x SR) takes around 60 seconds. 57 CHAPTER2.–IMAGE ENHANCEMENT USING SUPER-RESOLUTION 2.4.2 Single-image SRIR implementations In [Mal06] Mallat et al. propose a novel approach for carrying out SR in the frequency domain based on bandlets [MP07, PM08]. The algorithm does not rely on motion estimation but rather on a dictionary-based bandlet matching executing in a single-frame context. The meta data describing the patches of the LR image (frequency representation) is produced by means of geometric total variation estimation. Based on these metrics most suitable geometry reconstruction method, out of a provided a priori set, is found and applied. The processing differentiates the processing (and the reconstruction sets) between fine grain patches (areas rich in details) and coarse grain patches (plain/derived of details areas). The above-described super-resolution image reconstruction method (with noise removal) based on geometric spatio-temporal bandlets transformation is implemented by Let it wave LB101. The LB-101 is a technology converter for Broadcast/ProAV and Home Cinema applications capable of SD (Standard-Definition) to HD (High-Definition) upconversion as well as a cross-conversion from 1080i to 720p and 1080i to 1080p format. The LB-101 is available as a standalone integrated circuit with reference design, a soft macro IP implementable on Altera FPGAs and/or HardCopy Structured ASICs, or as a complete, easy to integrate module mezzanine board. LB-101M FPGA implementation requires 70000 logic elements (when mapped onto Altera Cyclone-II 70 device) and can be found in high-end electronics devices, i.e. in Analog Way’s HD Optimizer. In [NEC09] NEC has introduced its approach to single-image SRIR, namely the NEC µPD9245GJEC SoC. The device implements NEC Electronics proprietary single-frame super-resolution algorithm with blur reduction. Details on the algorithm have not been disclosed, however, it is most likely to be a dictionary based hallucination in frequency domain similar to the one implemented by the LB-101. NEC claims that its design is capable of upscaling: (i) quarter VGA (QVGA) resolution (320x240 pixels) to wide VGA (WVGA) resolution (800x480 pixels), (ii) NTSC format resolution (720x480 pixels) to wide extended graphics array (WXGA) resolution (1366x768 pixels) or super XGA (SXGA) (1280x1024 pixels), at the rate of 60 frames per second (fps). The µPD9245GJEC is also available as a soft macro destined for NEC’s cell-based ASIC (CB-90 and CB-12) and gate array (CMOS-12M) libraries. Okuhata et al. in [OIOS13] present an image up-converter that claims to deploy ‘superresolution’ to provide a more accurate depiction of edge and details than those of conventional interpolation algorithms. The algorithm checks each region in a given frame for presence of edges and then applies distinct interpolation functions to pixels which are part of an edge as well as pixels located in ‘smooth’ regions of the image. The fact that this algorithm 58 2.4 SUPER-RESOLUTION HARDWARE IMPLEMENTATIONS does not use a dictionary but a static set of coefficients suggests that the algorithm belongs to the reconstruction based single-image family or SR methods. Nevertheless, based on the given details this algorithm can also be classified as an adaptive form of linear interpolation which is capable of selecting a suitable set of convolution coefficients according to the edge orientation in the vicinity of the interpolated pel. When implemented using Verilog and mapped onto an Altera Arria II GX EP2 AGX125 EF35C4 device the implementation supports a maximum resolution of 1920x1080 pixels at a maximum frame rate of 60 fps at a 148.5 Mhz operating frequency. In [Goh14] Goshi proposes a novel SR method based on non-linear signal processing. The algorithm detects edges in the input frame using a high pass filter and uses them as the input to a non-linear function in order to create harmonic waves that have higher frequency than the LR input. Once saturated by a limiter, the created high frequency details are added to the upscalled version of the original input enhancing its quality. The quality enhancement claim is justified by presentation of 2DFFT results that show additional high frequencies in the output image. No complementary objective quality assessment values are presented. The author claims successful implementations in FPGAs but does not provide any implementation data justifying the claim. 59 CHAPTER2.–IMAGE ENHANCEMENT USING SUPER-RESOLUTION TABLE 2.3 Summary of hardware implementations executing in multi-frame context presented in Section 2.4. Contributor Year Method Platform Input SR fps LUTs BRAM Limitations Callico [CLL+06] 2006 NISR ISR Picasso COTS* (ARM+VLIW) N/A N/A N/A N/A N/A Memory occupancy and access latency. Implements ME. Callico [CN07] 2007 XSR (NISR) Excalibur XA10 (ARM+PLD) QCIF 2x 15 N/A N/A Memory access latency and bus contention. Angelopoulou [ABCC09] 2008 IBP Celoxica ADMXRC4SX (Virtex-4 FPGA) VGA 2x 25 5000+ 150+ Limited SR quality, memory access buffering scheme, multi-pass nature. Bowen [BB08] 2008 IBP Xilinx XC2V6000 (FPGA) CIF 4x 61 35707 ∼2400 Kb Limited SR quality due to limited resoureces (memory). Singla [STdA13] 2013 ESR (NISR) Xilinx S6LX45T (FPGA) CIF(/HD) 2x(/4x) 0.02 ∼38420 ∼2025 Kb Software running on soft-core processors. 60 2.5 CONCLUSIONS 2.5 Conclusions This chapter focused on the presentation of the basic concepts and the state-of-the-art of the super-resolution process. Over the last 40 years several approaches to the super-resolution problem have been developed. The most important approaches were described and organized in a complete taxonomy presented in this chapter. Super-resolution of video sequence introduces an additional level of complexity to the already complex super-resolution problem. The challenges and approaches to super-resolution of video sequences have been introduced. The non-uniform grid projection algorithm that is of relevance to this thesis has been contextualized as belonging to the interpolation-based family of the direct classical multi-frame super-resolution algorithms. The main advantage of this class of SR algorithms is their relatively low computational load, which is essential in making them suitable candidates for real-time implementations in hardware. Finally, the chapter concluded with the presentation of the state-of-the-art of the superresolution hardware implementations in FPGAs. The availability of such implementations is scarce. The main factors limiting the success of FPGA-targeted implementations have been identified as being resource-related. In particular: (i) Most implementations have been found to be limited by the available device memory. (ii) Iterative algorithms implementation have been found to be additionally limited by logical resources and tend to offer lower-then-the-reference output image quality due to the limited number of implemented iterations. In spite of the above drawbacks, FPGA devices are still the best platform for multi-frame SR systems prototyping in hardware. 61 Chapter 3 The non-uniform grid projection algorithm 3.1 Introduction In this work we tackle the challenge of providing real-time SR of video sequences by using the non-uniform grid projection algorithm that has been proposed in [MC03]. NUGPA is a fusion algorithm that looks for and exploits the non-redundant data encountered in a set of images due to the warping associated with movement and aliasing caused by the bandlimited registration sensors. Using the classification presented in Section 2.2.3, NUGPA falls into the interpolation-based fusion techniques of the direct multi-image family of spatial domain methods. The main advantage of this class of SR algorithms — the relatively low computational load, which is essential in making real-time applications possible — comes at the price of a limited degradation model. Moreover, the optimality of the reconstruction algorithm is not guaranteed, since the reconstruction step ignores the errors that occur in the interpolation stage. In this work we assume that the blur, noise characteristics, point spread function and decimation factor is common and space invariant in all LR images. Additionally, we will consider the relative motion to be purely translational. Over the years two versions of NUGPA have been developed by IUMA in cooperation with Philips research, namely, the iterative and non-iterative (or single-pass) versions. Based on the comparison of both versions presented in [MC03] we have opted for using the non-iterative version due to its deterministic execution time and better quality of the superresolved image. Additionally, implementations of iterative algorithms (e.g. iterated back projection) have been reported to not be able to implement sufficient number of iterations in order to provide satisfactory quality of the output image [ABCC08, BB08]. The single-pass 63 CHAPTER3.–REFERENCE SUPER-RESOLUTION ALGORITHM nature of the non-iterative NUGPA approach not only facilitates hardware implementations but also increases the chances of providing software-level quality of the output. The other main issue of the FPGA-targeted SR implementations, namely high memory requirements, will be tackled in the following chapter. 3.2 The non-uniform grid projection algorithm NUGPA carries out super-resolution image reconstruction in accord with the three-stages classical flow of fusion-based algorithms presented in Fig. 1.1. The stages defined by the classical approach, namely (i) pre-processing, (ii) multi-image fusion, and (iii) post-processing, correspond to the motion estimation, fusion-based reconstruction and non-uniform interpolation stages of the NUGPA, respectively. Due to this mapping, NUGPA is an example of the interpolation-restoration direct methods. Detailed description of execution flow is the focus of this section. 3.2.1 Image registration In the first stage the warping function and its metrics are estimated. The warping function usually is not known a priori, hence, regions which could contain supplementary information have to be found at run-time. What is known is that these regions are believed to differ only slightly from regions that they could enrich. In order to register images and find regions that are most probable of containing additional information the NUGPA super-resolution algorithms deploys a variant of block matching motion estimation. During image registration a set of neighboring pels (hereafter a macro-block (MB)) from the processed frame is matched against pels from other frames from the sliding frame window. For each MB only a limited set of pels (called search area, SA) confined within certain spatial vicinity (defined by the so called search area radius, SAR; expressed in number of pels) participate in the process of candidate set creation during motion estimation. The motion estimation process calculates the so called hyperdata, which in our case comprise the motion vector (MV) and the associated similarity criteria values. The motion vector is a vector that identifies the MB for which the similarity indices have the value closest to the optimum being sought for. An example of block-matching for a SFW comprising two reference frames, search area radius equal to MB width (MBwidth ), comprising (2×MBwidth +1)2candidates, is illustrated in Fig. 3.1 (for clarity only nine candidates are shown). 64 3.2 THE NON-UNIFORM GRID PROJECTION ALGORITHM SAD SAD Sliding Frame Window Search Area Motion Vector Patch being processed Best Match Candidate Fig. 3.1 Image registration using block matching motion estimation with sum of absolute differences (SAD) used as the similarity criterion. 3.2.2 Fusion-based reconstruction The hyperdata produced during the motion estimation stage are passed on to the SR kernel. The restoration processes starts with HR grid creation. The dimensions of this grid are determined by the used image registration precision (hereafter precisionir or precisionme). First, the so called up-holes transformation is carried out. During this transformation a HR grid gets filled with LR pels of the frame being super-resolved. Having g(x,y,t)representing pixels of the LR frame captured at time instance t, with spatial coordinates (x,y), the new HR spatial coordinates (x,y)are computed by multiplying the LR coordinates by the motion estimation precision, in accord with (3.1). (x,y) = (x∗precisionir,y∗precisionir)(3.1) 65 CHAPTER3.–REFERENCE SUPER-RESOLUTION ALGORITHM 3.3.2 Image registration for the non-uniform grid projection algorithm The choice of the appropriate image registration method is made per application. The selected image registration has to meet the registration requirements in terms of accuracy and execution time, especially if real-time execution is being targeted. Our implementation of the non-uniform grid projection algorithm carries out image registration by deploying block matching—a direct area-based image registration in spatial domain with piece-wise mapping. The used block-matching implementation limits the transformations model to purely translational ones. This is a common design choice for systems aimed at real-time processing. More complex transformations models (e.g. affine) prove to be simply beyond capabilities of today’s customer grade hardware. In fact, even real-time registration of transformations with 2 DoF can be a very demanding task, especially in cases of large search areas traversed using not efficient search strategies. 3.3.2.1 Provision of sub-pixel accuracy The non-uniform grid projection algorithm falls into the category of classical multi-frame super-resolution. As aforementioned, these algorithms are based on exploiting aliasing and/or variations created due to sub-pixel variations between the captures. Thus, in order to be able to provide any enhancement these algorithms require that the registration be carried out with higher-than-full-pixel accuracy. The techniques described up till now have been designed to carry out the registration without considering the sub-pixel domain. Even though, a sub-pixel estimate could be created based on the results of the full-pixel ones (e.g. weighted average of a set of ‘best’ results), the accuracy offered by this approach is in most cases insufficient. There are several techniques that allow to obtain a better accuracy [TH86, GSP86]. One common approach, and the one used by the authors of this work, is to modify a well known full-pixel technique and perform (some of) its steps at a finer, subpixel level. This approach operates on a pyramid of images where the search is cascaded through different accuracy levels. The representations of each level of the pyramid are created by re-sampling the previous coarser-level representation. Alternatively, interpolation can be avoided by using re-sampling based on Taylor series approximation. However, this approach is too complex to allow real-time execution. In our implementation, sub-pixel accuracy has been provided by preforming the search in three steps cascaded over three levels of accuracy. The estimations are cascaded from coarser-to-finer. Before being used, the results from the coarser-level are projected onto the finer-grain coordinates space. These coordinates after the projection are used as the refer72 3.3 MATHEMATICAL MODEL OF THE NON-ITERATIVE CLASSICAL SUPER-RESOLUTION 1 1 1 1 1 1 1 222 2 2 222 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 Fig. 3.5 Image registration using iterative refinement across three scales. ence point of finer-level search area. The search area encompasses coordinates that correspond to fractional coordinates at the coarser-level. The missing values of these coordinates, the so-called holes, are estimated using bilinear interpolation. Once the re-sampling is done the search continues. The coarsest-level corresponds to the full-pixel granularity search algorithm as the first stage of motion estimation, followed by a sub-pixel motion refinement. An example of a particular search carried out using the implemented block-matching algorithm is illustrated in Fig. 3.5. In order to estimate the hyperdata (results of the registration process) with quarter-pixel accuracy the searches are cascaded over three levels. Each finer-granularity level corresponds to re-sampling with a magnification factor of 2. 3.3.2.2 Similarity criterion The last component of the image registration puzzle is the selection of the similarity criterion. In this work, we carry out the similarity criterion selection under the assumption that the images being registered have been captured with insignificant variances in intensity. This is the typical case. This assumption allows to avoid normalization which is required if the above is assumption is not made. The most widely used similarity criterion for correlation-based image registration in spatial domain, among others, are the: Cross73 CHAPTER3.–REFERENCE SUPER-RESOLUTION ALGORITHM Correlation Function (CCF), Mean Square Error (MSE), Mean Absolute Error (MAE), Sum of Absolute Differences (SAD) and the Pixel Difference Classification (PDC). Cross-correlation is the basic criterion. It has proved itself to be useful for images which are misaligned by small rigid or even affine transformations. The problem with CCF and MSE is that, both of these criteria require square root computations, which implementation in hardware is costly. The criteria based on l2norm penalize heavily large projection errors. These errors are usually caused by variations that should be eliminated (e.g. outliers), but are also related to aliasing in higher frequencies. By suppressing the use of these data we may effectively limit the availability of the complementary information on which SRIR operates. The MAE, SAD and the PDC criteria are based on computing the absolute differences, making them computationally much simpler and easier to implement in hardware. Moreover, these methods deploy l1norm which can handle outliers much better than l2and is more robust. The possible disadvantage of these functions is the fact that they are not differentiable at the origin. This makes them not well suited for gradient-descent approaches. This is not the case of the modified three-step search. A crucial disadvantage of the PDC is that its use results in many ‘ties’ that require additional computations to be broken. Also, when compared with SAD, it requires additional comparisons even in the tie-less scenario. The fact that SAD neither requires the final division nor has to deal with fractional values (and their representation) makes it the preferred similarity criteria for hardware implementations. 3.3.2.3 Compilation determinable parameters Current version of the implemented block-matching defines a set of customizable parameters. These are referenced as the (template) Macro-Block width (MB width or MBwidth), Search Area Radius (SAR), the registration precision (precisionir) and the number of sensed images used in registration (simply, reference frames (RF)). Registration precision defines the finest-level accuracy, and indirectly the number of levels in the search hierarchy/pyramid. The reference software operated at quarter-pixel accuracy, corresponding to registration precision of 4 (precisionir =4). Setting the value of the macro-block width allows to choose the size used to tessellate the image (and matching). The size of MB significantly influences registration accuracy, as well as the computational and memory requirements. Smaller MB size relaxes memory requirements, at the cost of higher computational complexity and higher probability of ‘false’ motion vectors [BAA05]. This macro-block width is defined at compilation time and can be set to one of the supported values (4, 8 or 16). The chosen value is interpreted as the 74 3.3 MATHEMATICAL MODEL OF THE NON-ITERATIVE CLASSICAL SUPER-RESOLUTION one dimensional span of the template expressed in (full) pels, that is, as the number of pels in horizontal/vertical direction that make up the template. Current version supports only square-shaped templates (block). Run-time variable template size is not yet supported. NUGPA specifies two dimensions of search areas: one with fixed and with one parametrizable number of candidates. Both of the templates include all of the pels of the candidate template, but (can) differ in the number of included pels from its immediate vicinity. The fixed template includes only the pels which are the immediate neighbours of the candidate template. This effectively fixes the maximal number of candidates to 4×MBwidth +4. The customizable template allows inclusion of a larger range of the candidate template neighbors. The maximal distance at which a pel will be included in the search area is defined by the value of the search area radius parameter. The customizable template is used to define the search area used for matching at the full-pixel (coarsest) level. The candidates set is built by picking up every pel of a spiral-shaped path around the initial guess pixel (upper left corner of the block), turning SAR times, and starting with the top-left pixel in every turn. The cardinality of the set of candidates is modeled as card(Cp) = 1+2×precisionir ×(1+2+...+SAR) = 1+2×precisionir ×SAR×(SAR+1), (3.11) where Cprepresents the set of candidates of the parametrizable search area. The fixed template is used for performing the second and thirds step of search. Both of the defined templates are illustrated in Fig. 3.6 for the case of using a square-shaped template of MBwidth width and a search area radius of SAR pels. As aforementioned, the number of sensed frames is specified by setting the value of the RF parameter. Our algorithm uses the sliding-window-frame approach in which the registration process is carried out for each possible combination of the source and sensed images from the frame window. Thus, the time required for image registration grows proportionally with the number of used sensed images. The same frame window is used by the super-resolution process that follows image registration. Thus, this value has a significant impact not only on the execution time of image registration but also on the quality of the super-resolution process as a whole. 3.3.3 Mathematical model of the image registration stage As aforementioned, area-based image registration in spatial domain is carried out by means of template matching. In this work only the case of purely translational model in 2D space with piece-wise (block) mapping is considered. In this case, image registration is modeled 75 CHAPTER3.–REFERENCE SUPER-RESOLUTION ALGORITHM SAR SAR MBwidth search area height (step 1) SAR SAR MBwidth search area width (step 1) (a) Step 1 (customizable; full-pixel). 1 1 MBwidth search area height (step 2 & 3) 1 1 MBwidth search area width (step 2 & 3) (b) Step 2 & 3 (fixed; sub-pixel). Fig. 3.6 Search area templates used by different steps of the implemented matching. as a function of a set Tof Nblocks τi,i∈[1,N]∧∈N, a search space described as candidate set Cwith Melements, and a similarity function ε. The outcome of image registration, called the hyperdata, contains (i) a set of the minimal value of the dissimilarity criterion e(τ→C), (ii) a set of displacement vectors −→ ∆(τ→C)and (iii) auxiliary data. The hyperdata uniquely identify, for each of the defined template blocks, the candidate for which the minimal dissimilarity score has been observed. In classical SR the sets T(templates) and C (candidates) are exclusive and refer to regions of images of the same natural scene captured observed with different viewing conditions (temporal or spatial). Before the matching takes place, image registration process starts off by creating higher resolution representations of the input images. To this end, all the images undergo the upsampling transformation with the fiencapsulating an interpolation algorithm and scale s=rcard(X2) card(X2): g(x,y,t) = upsamples(g(x,y,t)).(3.12) Regions of these representations obtained through tessellation will be used as the templates τbeing registered. 3.3.3.1 Modeling exhaustive optical flow matching The modeling will start off with presentation of the model of exhaustive template matching of a single pel. This model will be then gradually extended in order to model the version 76 3.3 MATHEMATICAL MODEL OF THE NON-ITERATIVE CLASSICAL SUPER-RESOLUTION of block-matching used by the non-uniform grid projection algorithm implementation. The exhaustive template matching of a single pel describes the case in which the template is a singular pixel and the candidate set is the whole sensed frame (effectively N=M). Let us define εas a function that given intensities of a pair of pels returns a value that measures their similarity, being the pair a template-pel pτ:τ∈[1,N]located at (x,y), and a candidatepel pc:c∈[1,M]from the candidate set Clocated at coordinates (x+δxc,y+δyc). Then, this value, denoted as e(τ→c), for the assumed purely translational model, is modeled as e(τ→c)=ε(τ→c)(pτ(x,y),pc(x+δxc,y+δyc)).(3.13) The dissimilarity criterion is constructed in such a way that it assumes lower values for ‘better’ matches. Thus, the ‘best’ match for the template from the candidate set of pels is found by minimizing the dissimilarity function computed over all the pels of the candidate set (exhaustive search). Considering (3.13) in the context of all c∈[1,M], it can been seen that in the investigated case the candidate set corresponds to the whole search space and is exhaustively traversed by using all possible displacements between the template and the candidate pels. Let us define a set of all possible displacements ∆τ Cthat uniquely identifies each and every candidate from the candidate set Cfor a template τ. Then, the search for the best match of a single pel as the search for a displacement vector −→ ∆(τ→C)for which the dissimilarity function reaches its minimum can be modeled as: −→ ∆(τ→C),arg min ∆=(δxc,δyc)∈∆τ C ετ→c(pτ(x,y),pc(x+δxc,y+δyc)).(3.14) 3.3.3.2 Modeling exhaustive region-based matching Now, let us relax the assumption that the template contains only one pixel. This requires changes in the similarity criterion, so that it would allow all the pels contained within the limits of the pattern to contribute to the overall similarity score (sub-sampling is not considered here). Let us define a region W(x,y)as a local neighborhood centered at coordinates (x,y)and comprising pels whose coordinates are within the distance or radius of mpels [±m,±m]from the center of the set. Coordinates that point outside of the image boundaries form a corner case of the process. In this description it is assumed that these pels are detected handled by the function that computes the dissimilarity score. Now let us define ε(Wτ→Wc)as a function that computes the dissimilarity score given two regions Wτ(x,y)and Wc(xc,yc)centered, respectively, at coordinates (x,y)and (xc,yc):xc=x+δxc;yc=y+δyc. 77 CHAPTER3.–REFERENCE SUPER-RESOLUTION ALGORITHM Then the e(Wτ→Wc)is computed as e(Wτ→Wc)=ε(Wτ→Wc)(Wτ(x,y),Wc(xc,yc)).(3.15) Assuming that ε(Wτ→Wc)can be modeled as a combination, denoted as ϒ, of similarity measures computed element-wise over the pels of the two regions leads to e(Wτ→Wc)= m ∑ mx=−m m ∑ my=−m ϒ(ε(τ→c)(pτ(x+mx,y+my),pc(x+mx+δxc,y+my+δyc))). (3.16) Considering the template matching for sets of pels and assuming that all candidate regions Wcfrom the candidate set Cto be matched with the template region are uniquely identified by a set of possible displacements ∆Wτ C, the minimization problem (3.14) becomes −→ ∆(Wτ→C),arg min ∆=(δxc,δyc)∈∆Wτ C ε(Wτ→Wc)(Wτ(x,y),Wc(x+δxc,y+δyc)).(3.17) 3.3.3.3 Sum of absolute errors as a robust matching criterion The template matching implemented by NUGPA deploys the sum of absolute differences as the matching criterion. The value of SAD is computed between the template and the candidate patch located at coordinates identified by a vector representing relative displacement between the regions. Let us define a region as a series of n=mx×mypels luminance values represented as a lexicographically ordered vector pi= [pi1,pi2,...,pin−1,pin]T. Then, the dissimilarity e(Wτ→Wc)is determined by computing the sum of differences SADτ→cbetween the representations of the template pτiand candidate pci. That is, e(Wτ→Wc)=SADτ→c= n ∑ i=1|pci−pτi|.(3.18) Block matching is considered an effective scheme for textured regions, but exhibits performance degradation when carried out for regions with high homogeneity [BAA05]. In order to prevent quality degradation, homogeneous regions are identified and handled differently. For this purpose, additional SAD between pels belonging to a candidate MB ci and average pel value for this MB (pcicomputed as in (3.19)) is estimated, in accord with (3.20). This value, represented as SADintra ci, forms a part of the hyper data and is passed on to the SR kernel. pci=∑n j=1pcj n(3.19) 78 3.3 MATHEMATICAL MODEL OF THE NON-ITERATIVE CLASSICAL SUPER-RESOLUTION SADintra ci= n ∑ j=1|pci−pcj|(3.20) 3.3.3.4 Modeling of heuristic region-based matching As aforementioned, due to its computational cost the exhaustive (or full) search is not used in practice. Practical implementations of block-matching rely heavily on heuristic search strategies, which operate only on a subset of the possible candidate set Cdefined for a particular search space. In order to model heuristic search (or search strategies) let us define a function Φthat given the set of possible candidates C, the search strategy αand the reference coordinates xr,yrproduces a set of allowable candidatesCαon which the heuristic search operates. Cα=Φ(C,α,(xr,yr)) (3.21) 3.3.3.5 Modeling hierarchical sub-pixel matching In order to be able to achieve sub-pixel accuracy NUGPA carries out image registration in three steps, each executing at different accuracy-level. The results of registration are cascaded between the steps. The cascade is formed in the direction of the finer steps. The results of template matching at a coarser-level are used as the initial point to guide the matching in the subsequent, finer-accuracy space. Before the matching takes place the results have to be projected onto the finer-step’s grid and the vicinity of the initial guess has to be re-sampled in order to create values included in the subsequent-step’s search area. Re-sampling is carried out using the already defined upsampling operator with fibeing an bilinear interpolation kernel and the projection scale sK→Ffrom coarser Kto finer Fstep sK→F=scard(F) card(K).(3.22) The projection operation is carried out by multiplying the values of the coarser-accuracy results by the projection scale. Denoting the displacement returned from template matching carried out at the coarser as −→ ∆K (Wτ→C)= (δK x(Wτ→C),δK y(Wτ→C)), the projection scale sK→F, then corresponding displacement at a finer-level is modeled as −→ ∆F (Wτ→C)= (δF x(Wτ→C),δF y(Wτ→C)) = (sK→F·δK x(Wτ→C),sK→F·δK y(Wτ→C)) = sK→F·−→ ∆K (Wτ→C) (3.23) 79 CHAPTER3.–REFERENCE SUPER-RESOLUTION ALGORITHM Let’s label the accuracy levels from the coarsest (step 1) to the finest (step 3), respectively as full-pixel (FP), half-pixel (HP) and quarter-pixel (QP). Then, using the new terms and the defined operations the hierarchical matching of a single regionWτcarried out by the NUGPA is modeled as a series of three subsequent minimizations, that is −→ ∆FP (Wτ→CαFP ),arg min ∆=(δxc,δyc)∈∆Wτ CαFP SAD(Wτ→Wc)(Wτ(x,y),Wc(xc,yc)),(3.24) −→ ∆HP (Wτ→CαHP ),arg min ∆=(δxc,δyc)∈∆Wτ CαHP SAD(Wτ→Wc)(Wτ(x,y),Wc(˜xc,˜yc)),(3.25) −→ ∆QP (Wτ→CαQP ),arg min ∆=(δxc,δyc)∈∆Wτ CαQP SAD(Wτ→Wc)(Wτ(x,y),Wc(xc,yc)),(3.26) where, the candidate sets used in the sub-pixel matching are created as CαHP =Φ(C,αHP,(sFP→HP ·x+δHP x(Wτ→CαFP ),sFP→HP ·y+δHP y(Wτ→CαFP ))),(3.27) CαQP =Φ(C,αQP,(sHP→QP ·x+δQP x(Wτ→CαHP ),sHP→QP ·y+δQP y(Wτ→CαHP ))).(3.28) Registration of the whole image is carried out by repeating the above process for all regions of the source image. In case of performing the registration for multiple sensed images, the whole registration process is repeated for each combination of the source-sensed images. In the case of a set κcontaining kimages of which each is divided into wtemplate regions containing up to m×mpels, then, using the defined terms, the image registration process for image t:t∈κis executed M×(k−1)times for each template region of each frame i∈κ∧i=t. The result of the image registration process (the hyperdata) of the image tover a set of k−1 sensed images (all images except the image t) are w×(k−1) displacements in X2space and (corresponding) computed sums of absolute differences, forming a set {−→ ∆QP (Wτ→Ci),SAD(Wτ→Ci),SADintra ci}for τ: 1 ≤τ≤w∧τ∈N;i:i∈κ∧i=t and quarter-pixel accuracy of the X2space. 3.3.4 Mathematical model of the restoration-interpolation stage In case of the non-uniform grid projection algorithm the second stage of the SR process is the restoration-interpolation. As the name suggests, at the highest level of abstraction, the restoration-interpolation stage can be seen as composed of two stages: (i) the reconstruction 80 3.3 MATHEMATICAL MODEL OF THE NON-ITERATIVE CLASSICAL SUPER-RESOLUTION or fusion of the information extracted from the input images, and (ii) interpolation used to compute the data that are still missing. The reconstruction process carries out the fusion of the image to be super-resolved with the non-redundant data found in the HR projections of the remaining κimages. The non-redundant data are extracted from the locations identified by the hyperdata. In that sense, this part of the reconstruction process can be seen as an inverse of the image registration step. As it was the case with image restoration, the first step of the reconstruction process is the up-sampling of the input low resolution frames g(x,y,i):i∈κto the resolution of image registration. The resulting HR representations g(x,y,i):i∈κare created using the up-sampling operator (3.9) with the interpolation kernel fidefined to return a constant value of /0. This is modeled as g(x,y,i) = upsamplercard(X2) card(X2)(g(x,y,i)).(3.29) The super-resolution kernel k(x,y)processing is carried out by applying the fusion kernel FK on the frames from the set κ, carrying out interpolation process INT on the output of the fusion and downsampling the obtained representation. Let us defined the hyperdata obtained from the image registration of an image tfrom a set κcontaining kimages, over all of the remaining images of this set as Ωκ t. Then, k(x,y) = downsampleqprecisionir scalesr (INT (FK(g(x,y,κ),Ωκ t),int(x,y)),(3.30) where int(x,y)represents the interpolation kernel used in the interpolation process. Let us define the fusion operator ⊕that carries out the operation of projecting pels of one region onto another. Using this operator the fusion of an image over a set of images is modeled as (3.31) where ˆs(x,y,t)represents the outcome of the fusion process for frame t. ˆs(x,y,t) = k ∑ i=1,i=t g(x,y,t)⊕g(x,y,i)(3.31) Then the fusion process of the image captured at the time twith another image from κ can be modeled as an act of performing a series of fusions of regions of these images. Let ⊕operate in a way that allows it to carry out the projection of the region that is its righthand argument onto a grid defined for region which is its left hand argument. Assuming, cardinality of τto be wand that the fusions of each region result in an immediate update of 81 CHAPTER3.–REFERENCE SUPER-RESOLUTION ALGORITHM been exploited intensively in Computer-Generated Imagery (CGI) applications where the level of details of a texture or a model varies with the viewing conditions — depending on the distance from the observer and display density different sets of textures/models are being displayed [Dur96, DH97]. Multi-scale SSIM tackles the issues of incorporating the influence of these conditions on image quality by creating multiple scaled-down versions of the compared images and carrying the assessment not only for the original images but for all created versions. The scaled-down versions used in the assessment are created by iteratively applying a low-pass filter and downsampling the filtered image by a factor of 2 as illustrated in Fig. 3.7. This effectively removes the limitation of the evaluation being appropriate for only one specific viewing condition. Let us label the original image as scale 1, and the highest scale M as corresponding to the most downsampled image (obtained after M−1 iterations). At the i-th scale, the contrast comparison and the structure comparison are calculated and denoted as ci(x,y)and si(x,y), respectively. The luminance comparison is computed only for the highest scale (scale M) and is denoted as lM(x,y). The overall evaluation is obtained by combining the measurements at different scales using MSSIM(x,y) = [lM(x,y)]αM· M ∏ i=1 [ci(x,y)]βi·[si(x,y)]γi,(3.47) where, as in (3.45), αM,βi,γiare used to adjust the relative importance of different components. The values of these cross-scale parameters need to be evaluated empirically during the calibration process. The number of iterations required for the MSSIM to give relevant values is reasonably small. The authors of this method use five levels in their work [WSB03a]. For bigger images (in terms of number of pels) it may be justified to use more levels. 3.4.2 Output quality as a function of SR parameters The goal of the study presented in this section has been to provide a quantitative evaluation of the reference software implementation of NUGPA (hereafter SRiuma or SRiuma) and to determine the impact that the SR algorithm parameters, tabulated and briefly described in Table 3.1, have on the quality of the obtained super-resolved images. For the purpose of algorithm performance quantification a set of 120 combinations of the algorithm parameters values—MBwidth (4, 8 and 16), SAR (2, 4, 8 and 16), number of RFs (2, 4, 8, 12 and 16), and the scale (2 or 4)—has been defined to be used in simulations. 88 3.4 QUANTITATIVE EVALUATION OF THE REFERENCE SOFTWARE IMPLEMENTATION reference image test image F2F2F2 c(x,y) s(x,y) 1 1 c(x,y) s(x,y) M M c(x,y) s(x,y) 2 2 l(x,y) M similarity measure F2F2F2 Fig. 3.7 Multi-scale structural similarity measurement system. F: low-pass filter; 2↓downsampling by 2. Modified version of [WSB03a]. 3.4.2.1 Experimental set-up The reference image (and thus the output) format has been chosen to be the YUV 4:2:0p 8bit CIF (288x352 pels) sequence. This resulted in the LR input being a QCIF (144x176 pels) or a 72x88 pels frame format for the tested values of the super-resolution zoom (scale) of 2 and 4, respectively. Even though, due to high computational complexity and the fact of proven better correlation between the objective proximity measures (MSE, PSNR) with the results of subjective evaluation, it is common to use images of lower resolutions (32x32, 64x64, 128x128) [ESHEK12, chapter 2], in this study higher resolution inputs have been used in order to facilitate better quality of details inspection and more accurate MSSIM index assessment. Results of the study presented in [WBS05] suggest that the use of the color components does not significantly change the performance of the used quality measurement model. This was expected as it is a well-known fact that most of the energy is express by the luminance component [IP91]. These carried out experiments considered only the quality of luma component, which is the only one for which super-resolution is carried out. The simulations were carried out in accord with the flow presented in Fig. 3.8: 1. First the YUV 4:2:0p 8bit low resolution (QCIF or 72 ×88 pels) sequences were obtained from the CIF reference sequence. In order to do so, the CIF sequence was loaded, and decimated, forming a sequence that was then stored. 2. The LR sequences were used as input for tested algorithms (the choice of the benchmark is the focus of the next subsection), resulting in CIF YUV 4:2:0p 8bit output. 89 CHAPTER3.–REFERENCE SUPER-RESOLUTION ALGORITHM TABLE 3.1 Algorithm parameters and their values used in the study on image quality. Parameter Value Description MBwidth 4, 8, 16 Number of pels that a MB contains in each dimension. Determines the memory used for storage of a MB and the number of MBs in a frame. FRcols 176 Number of pel columns in a LR frame. FRrows 144 Number of pel rows in a LR frame. precisionir 4 The ratio between the 1D dimensions of the finest accuracy grid used in image registration and the fullpixel grid used for input representation. scale 2, 4 The ratio between the 1D dimensions of the outcome accuracy grid and the full pixel grid used for input representation; SAR 2, 4, 8, 16 The absolute distance (in full-pixel coordinates) used to determine if a pel is contained within the limits of the SA; With MBwidth determines the SAsize. intwindow 2 The absolute distance (in finest IR accuracy coordinates) used for the determination of the neighboring MB’s pels that can contribute to the interpolation process. 3. The output luma pels ( ˆ pelSR) and the original CIF sequence pels pelre f were postprocessed by the IQA that quantified the grade of similarity between these two produced representations by calculating, among others, the PSNR and SSIM indices. In calculations both images are consider a lexicographically ordered sets containing N pels (total number of pels in the image). (a) MSE has been computed using (3.35). (b) PSNR was then computed based on the MSE value using (3.36) which for an 8 bit pel representation becomes PSNR =10∗log10       2552 1 N N ∑ i=1 (ˆ pelSR(i)−pelre f (i))2       .(3.48) (c) Single-scale structural similarity index has been calculated as follows: i. First, the Cconstants were computed using (3.44) with K1=0.01, K1= 90 3.4 QUANTITATIVE EVALUATION OF THE REFERENCE SOFTWARE IMPLEMENTATION Decimation SRiuma Benchmark Image Quality Assessment Reference Input Output CIF YUV 4:2:0 Decimated YUV 4:2:0 Original Fig. 3.8 Data flow of the experiments used to assess the output quality of SRIR. 0.03, and dynamic range (L) of 255 (8-bit grayscale). ii. Then, the mean intensity µXwere computed. The luminance comparison was then a function of µXand µY. µX=1 N N ∑ i=1 xi(3.49) iii. Following, the mean intensity was removed from the signal. In discrete form, the resulting signal was the x−µX. Standard deviation (the square root of variance) was used as an estimate of the signal contrast. An unbiased estimate in discrete form was given by σX= 1 N−1 N ∑ i=1 (xi−µX)2!1 2 .(3.50) The contrast comparison was then a function of σXand σY. iv. Third, the signal was normalized (divided) by its own standard deviation, so that the two signals being compared have unit standard deviation. The structure comparison s(x,y)was conducted on these normalized signals (x− µX)σXand (x−µY)σY. v. Finally, the three components were combined as in (3.46). (d) Multi-scale similarity index has been computed by using i. First, the Cconstants were computed using (3.44) with K1=0.01, K1= 0.03, and dynamic range (L) of 255 (8-bit grayscale). ii. Low pass filter, implemented as convolution (Matlab filter2 function) of 91 CHAPTER3.–REFERENCE SUPER-RESOLUTION ALGORITHM a 11-by-11 window with Gaussian distribution (standard deviation =1.5) with the image being filtered, was applied and its results were downsampled by 2. iii. Contrast and structure similarities were computed as in the aforementioned single-scale SSIM flow. iv. The iteration counter was incremented. If the final iteration has been reached the luminance index lM(x,y)was computed. Otherwise, steps 3(d)ii–3(d)iii were repeated. v. The computed SSIM indices were combined using parameters values from [WSB03a] (M=5, β1=γ1=0.0448, β2=γ2=0.2856, β3=γ3=0.3001, β4=γ4=0.2363 and , α5=β5=γ5=0.1333) and (3.47) which became MSSIM(x,y) = [l5(x,y)]0.1333 · 5 ∏ i=1 [ci(x,y)]βi·[si(x,y)]γi.(3.51) 3.4.2.2 Determination of the image quality benchmark This initial aim of this study was to compare the results of our algorithm with two competitive solutions: one interpolation and one SRIR implementation. Before the assessment of the proposed SRIR implementation was to be carried out the image quality benchmark needed to be chosen. In order to do so the available competitive solutions presented in the state of the art Section 2.4 have been looked at. A natural competitor and the one that had been successfully implemented in hardware [ABCC08, BB08] would have been the IBP algorithms. Unfortunately, to the best knowledge of the authors, software implementations of these implementations were either not available [ABCC08, BB08] or the available software [LCA07] was found to be unable to process real-life video sequences without severely corrupting the output. In order to determine the interpolation algorithm to be used as the quality benchmark simulations for three test sequences using six interpolation algorithms have been carried out as presented in Fig. 3.9. The sequence used in the experiments are briefly described in Table 3.2. These sequence form part of the Xiph.org Video Test Media [derf’s collection] available online [Xip15]. The produced images allowed computing the PSNR and MSSIM indices. The computed PSNR and MSSIM indices are presented in Fig. 3.10 and Fig. 3.11, respectively. The results for the nearest neighbour, bilinear, bicubic, lanczos2 and lanczos3 interpolation [BB09, Key81, Get11, MN13] have been obtained by means of invocation of the 92 3.4 QUANTITATIVE EVALUATION OF THE REFERENCE SOFTWARE IMPLEMENTATION Decimation nearest Image Quality Assessment Reference Input Output CIF YUV 4:2:0 Decimated YUV 4:2:0 Original bilinear bicubic lanczos2 lanczos3 mean nearest Fig. 3.9 Data flow of the experiments used to choose the benchmark algorithm. imresize() function from the Image Processing Toolbox of the MatLab programming environment [Mat96]. The non-uniform mean nearest neighbours (neighbourhood of 5x5 pels) algorithm have been implemented in software using ANSI C. From the resulting MSSIM values it can be seen that the mean nearest neighbors implementation has consistently outperformed the other five implementations reporting the highest index value in all test cases. In terms of the reported PSNR, the use of the mean nearest neighbors interpolation leads to superior results in 5 out of 6 performed tests, being outperformed only in one case by the bilinear interpolation. Based on these results the decision to use the non-uniform mean nearest neighbors as the benchmark for image quality assessment of our algorithm has been made. 3.4.2.3 PSNR and MSSIM as a function of algorithm parameters Once the benchmark algorithm has been chosen, the flow presented in the previous section has been used to carry out simulations for all of the aforementioned 120 combinations. Test sequences used in experiments are the same that have been used in the benchmark determination tests. The experiments were conducted on a Sun Ultra 24 Workstation, hosting one Intel Core 2 Quad cpu [email protected] GHz, 3.25GB of RAM memory, and using Windows XP Professional operating system with Service Pack 3. The produced output has been postprocessed using functions from the Matlab MeTriX MuX Visual Quality Assessment Package library [Gau07]. Let us remind that the version of the SR software used in this study is not robust, thus it is not capable of eliminating outliers and it does not include mechanisms to 93 CHAPTER3.–REFERENCE SUPER-RESOLUTION ALGORITHM gfx Page 1 foreman mobile paris foreman mobile paris SCALE 2 SCALE 4 0 5 10 15 20 25 30 35 PSNR observed for tested interpolation methods nearest bilinear bicubic lanczos2 lanczos3 mean nearest PSNR [dB] Fig. 3.10 Average PSNR (300 initial frames) observed for the interpolated foreman, mobile and paris sequences. RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 SAR2 SAR4 SAR8 SAR12 SAR16 SAR2 SAR4 SAR8 SAR12 SAR16 SAR2 SAR4 SAR8 SAR12 SAR16 MB16 MB8 MB4 0,945 0,95 0,955 0,96 0,965 0,97 0,975 0,98 0,985 0,99 0,995 MSSIM INT vs SR Qcif->Cif (scale=2) foreman (INT x2) mobile (INT x2) paris (INT x2) foreman (SR x2) mobile (SR x2) paris (SR x2) MSSIM [ ] gfx Page 1 foreman mobile paris foreman mobile paris SCALE 2 SCALE 4 0 0.2 0.4 0.6 0.8 1 1.2 MSSIM observed for tested interpolation methods nearest bilinear bicubic lanczos2 lanczos3 mean nearest MSSIM [1] Fig. 3.11 Average output MSSIM (300 initial frames) observed the interpolated foreman, mobile and paris sequence. 94 3.4 QUANTITATIVE EVALUATION OF THE REFERENCE SOFTWARE IMPLEMENTATION TABLE 3.2 Test sequences used in the experiment. Sequence Visualization Description Foreman Presents a worker performing rapid random head movements, while talking to the camera. The sequence then changes to present the working site. This sequence presents global and rapid local movement, along with context change. Mobile An electrical train moving on the table, a spinning planetary system and a calendar hanging on the wall. The train is moving horizontally, whereas the calendar is pulled to present vertical movement. The ball rotates representing a mixture of both types of movement. Global movement is present as the camera is moving to the left. This sequence contains rich and fine grain textures. Paris Presents a couple sitting by a table and talking. While talking, the woman is juggling a ball, and the men is playing with a pen. The sequence presents rapid local movements. Global movement is absent. correctly handle context changes, nor capable of correctly handling non-translational movement (panning, rotations, etc.). Also, noise attenuation by averaging is significantly limited as no proximity-based fusion weights computation is carried out (weights are set equal to 1). The focus of this section is the presentation of the PSNR and MSSIM indices obtained for the foreman, mobile and paris sequences. These three sequences have been found to be representative for all of the carried out test. The PSNR values computed for the foreman, mobile and paris sequences for SR scales of 2 and 4 are presented in Fig. 3.12 and Fig. 3.14, respectively. The corresponding MSSIM values are shown in Fig. 3.13 and Fig. 3.15, respectively. From these figures one can see that there is a strong correlation between the trends observed for both scales. In other words, a change of the SRIR parameters has a similar effect on the way the observed outcome quality is changed for a given sequence for both scale values. Even though, there is a strong correlation between the observed outcome quality for results obtained for different scale values, that can be clearly seen when comparing both of the aforementioned figures side by side, the way the algorithm parameters impact the outcome quality is not as straightforward. The observed correlations have been 95 CHAPTER3.–REFERENCE SUPER-RESOLUTION ALGORITHM found to be heavily dependent on the used test sequence. In order to analyze the way that each of the parameters impact the output quality the change in the outcome quality indices observed for the combinations for which other-than-the-analyzed parameters are fixed has been investigated. In the case of the foreman sequence, the highest PSNR and MSSIM values have been observed for the macro-block comprising 8x8 pels, SAR of 16 pels, and 4 reference frames (2 in each temporal direction; (SFW comprising 5 frames)) for both tested scales values. For this sequence the change of the MB size parameter manifests itself as a (<0.5dB) translation towards lower PSNR values. The magnitude of the observed translation is about two times bigger in the case of the switch to MB size of 4 than for the case of switch to the value of 16. In the case of the foreman sequence, increase in search radius has always led to a higher PSNR, with the noticed gain being higher for combinations with higher number of available reference frames. The observed impact of the change in the number of used reference frames on the output quality has been the greatest from all of the tested parameters (for a fixed scale). For foreman, better results have been reported for smaller number of available reference frames (2 and 4) which outperform interpolation in most cases. For the case with 8 reference frames SRIR, managed to surpass interpolation in some cases (3) when provided with large search areas. From the carried out tests none of the combinations with 16 reference frames managed to match the quality of the benchmark interpolation. The gap between the computed values for different RFs numbers is reduced with the increase of the search area radius. The noticed values converged the fastest for the medium MB size (8x8) and the slowest for the largest tested MBs (16x16). The described type of correlation between the output’s quality observed for scale value of 2 holds true also for scale value of 4. The difference between the two diagrams are the overall lower PSNR and MSSIM values and the fact that, in case of PSNR, the use of SRIR resulted in significantly better (>2.0 dB) quality than the benchmark interpolation for all 60 tested parameters combinations. As far as MSSIM is considered, SRIR managed to outperform the interpolation only in some cases (≈20%) with medium/large MBs and relatively large search areas. In the case of the mobile sequence, the highest PSNR and MSSIM values have been observed for the macro-block comprising 16x16 pels, SAR of 2 pels, and 8 (scale = 2) or 16 reference frames (scale = 4). The outcome quality has been observed to be better for larger macroblocks with the exception of PSNR observed for the cases with 2 reference frames and scale of 4. In contrast with the trends observed for foreman, for the mobile test sequence use of more RFs resulted, in most cases, in higher quality index—showing that the algorithm was able of taking advantage of additional information contained in these frames. 96 3.4 QUANTITATIVE EVALUATION OF THE REFERENCE SOFTWARE IMPLEMENTATION RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 SAR2 SAR4 SAR8 SAR12 SAR16 SAR2 SAR4 SAR8 SAR12 SAR16 SAR2 SAR4 SAR8 SAR12 SAR16 MB16 MB8 MB4 22 23 24 25 26 27 28 29 30 PSNR INT vs SR Qcif->Cif (scale=2) foreman (INT x2) mobile (INT x2) paris (INT x2) foreman (SR x2) mobile (SR x2) paris (SR x2) PSNR [dB] SAR2 SAR4 SAR8 SAR12 SAR16 SAR2 SAR4 SAR8 SAR12 SAR16 SAR2 SAR4 SAR8 SAR12 SAR16 MB16 MB8 MB4 22,5 22,55 22,6 22,65 22,7 22,75 22,8 Paris Qcif->Cif (scale=2) INT RF2 RF4 RF8 RF16 PSNR[dB] Fig. 3.12 Average PSNR observed for the super-resolved and interpolated foreman, mobile and paris sequences (initial 300 frames; scale=2). RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 SAR2 SAR4 SAR8 SAR12 SAR16 SAR2 SAR4 SAR8 SAR12 SAR16 SAR2 SAR4 SAR8 SAR12 SAR16 MB16 MB8 MB4 0,945 0,95 0,955 0,96 0,965 0,97 0,975 0,98 0,985 0,99 0,995 MSSIM INT vs SR Qcif->Cif (scale=2) foreman (INT x2) mobile (INT x2) paris (INT x2) foreman (SR x2) mobile (SR x2) paris (SR x2) MSSIM [ ] SAR2 SAR4 SAR8 SAR12 SAR16 SAR2 SAR4 SAR8 SAR12 SAR16 SAR2 SAR4 SAR8 SAR12 SAR16 MB16 MB8 MB4 22,5 22,55 22,6 22,65 22,7 22,75 22,8 Paris Qcif->Cif (scale=2) INT RF2 RF4 RF8 RF16 PSNR[dB] Fig. 3.13 Average MSSIM observed for the super-resolved and interpolated foreman, mobile and paris sequences (initial 300 frames; scale=2). 97 CHAPTER3.–REFERENCE SUPER-RESOLUTION ALGORITHM large set of reference images. •Loss in MSSIM for combinations with gain in PSNR. This was expected, as the core version of the SR software does not have the property of being robust, thus many outliers take part in the fusion process, with the same impact as the non-outliers. This is almost guaranteed to lead to compromised image structure (and decreased MSSIM), but can result (on average) in only small MSE changes in cases where the average difference in luma values of outliers and the reference is low enough. Presented figures prove that, even though, for some test (algorithm values) combinations SRiuma has not been capable of outperforming the baseline interpolation on a regular basis, there is at least one combination that provides a better-than-interpolation quality of results for each of the tested sequences (with the exception of MSSIM for paris). Experiments have shown that it is quite difficult to indicate one configuration that would guarantee the best performance over a set of sequences. Depending on the test sequence optimal value of macro-block size, search area, and reference frames varies. Thereby, when designing a hardware implementation a large set of super resolution parameters values should be taken into account, and the resulting organization/architecture should be made as customizable as possible. 3.5 Conclusions The focus of this chapter was the non-uniform grid projection algorithm. The chapter started with a brief description of the algorithm flow, the used motion estimation and the fusionbased SR kernel. Next, the mathematical models of the algorithm and the aforementioned stages were developed. Following the theoretical introduction, the objective quality of the super-resolved image produced by the reference software implementation was quantitatively and qualitatively evaluated. For the needs of this evaluation, average PSNR and MSSIM values have been computed for three representative test sequences and 120 configurations of the software implementations. The obtained results have shown that the analyzed NUGPA software implementation: (i) can provide better than interpolation quality of results expressed as PSNR or MSSIM, (ii) the quality of the output super-resolved image is very sensitive to the changes in values of algorithm parameters, 104 3.5 CONCLUSIONS RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 SAR2 SAR4 SAR8 SAR12 SAR16 SAR2 SAR4 SAR8 SAR12 SAR16 SAR2 SAR4 SAR8 SAR12 SAR16 MB16 MB8 MB4 0 1 Correlation/agreement between PSNR and MSSIM foreman (scale 2) mobile (scale 2) paris (scale 2) Parameters values Agreement (>0) or difference (0) [ ] SAR2 SAR4 SAR8 SAR12 SAR16 SAR2 SAR4 SAR8 SAR12 SAR16 SAR2 SAR4 SAR8 SAR12 SAR16 MB16 MB8 MB4 22,5 22,55 22,6 22,65 22,7 22,75 22,8 Paris Qcif->Cif (scale=2) INT RF2 RF4 RF8 RF16 PSNR[dB] Fig. 3.20 Correlation between PSNR and MSSIM indices (scale=2). RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 RF4 RF16 SAR2 SAR4 SAR8 SAR12 SAR16 SAR2 SAR4 SAR8 SAR12 SAR16 SAR2 SAR4 SAR8 SAR12 SAR16 MB16 MB8 MB4 0 1 Correlation/agreement between PSNR and MSSIM foreman (scale 4) mobile (scale 4) paris (scale 4) Parameters values Agreement (>0) or difference (0) [ ] SAR2 SAR4 SAR8 SAR12 SAR16 SAR2 SAR4 SAR8 SAR12 SAR16 SAR2 SAR4 SAR8 SAR12 SAR16 MB16 MB8 MB4 22,5 22,55 22,6 22,65 22,7 22,75 22,8 Paris Qcif->Cif (scale=2) INT RF2 RF4 RF8 RF16 PSNR[dB] Fig. 3.21 Correlation between PSNR and MSSIM indices (scale=4). 105 CHAPTER3.–REFERENCE SUPER-RESOLUTION ALGORITHM (iii) obtaining satisfactory results requires the configuration to be fine-tuned for the particular input sequence, (iv) does not implement sufficient mechanisms to limit the impact of outliers on the observable output quality, and (v) presents less variance in observed quality of results for smaller macroblock sizes. 106 Chapter 4 Proposed super-resolution algorithm 4.1 Introduction As seen in chapter 2, the main two factors limiting the success of hardware implementations are the iterative nature of processing and high resource occupancy. The choice of the noniterative non-uniform projection algorithm described in chapter 3 solves the former issue, promising efficient hardware implementation. Unfortunately, the memory occupancy of the non-iterative algorithm presented in the previous chapter remains prohibitively high for an FPGA implementation that intends to avoid using any off-device memory. A quick evaluation of the reference software implementation and the algorithm data flow identified the frame buffers used by the fusion kernel as the main factor contributing to the total memory occupancy. The implemented execution flow required at least two of the frames being processed to be readily available from memory throughout the fusion process. Thus, in order to eliminate the frame-level buffers, the algorithm’s execution flow needs to be modified. In particular, the data path’s flow granularity has to be refined in order to reduce the algorithms memory requirements. Execution flow granularity is one of the main characteristics of any algorithm. The spectrum of possible granularity choices spans from one pel (finest granularity) up to a (set of) frame(s) (coarsest granularity). The choice of the execution flow has a significant impact on the implementation’s characteristics, among others, memory/logic occupancy, system’s throughput and scalability. The deciding factor that determines the choice of the granularity of the execution flow is the type of targeted implementation and its bottlenecks. Software (super-resolution) implementations are (usually) not limited by the available memory but rather by the memory hierarchy and access latency. Thus, these implementations tend to apply processing on the coarsest type of elements possible, usually a (multiple of a) whole 107 CHAPTER4.–PROPOSED SUPER-RESOLUTION ALGORITHM frame. By doing so, a higher amount of data can be loaded at once using coalesced accesses, leading to the minimization of the number of memory accesses, the delay caused by of memory access latency, and optimization of the overall observed implementation’s performance. On the other hand, processing of coarse grain structures like a whole frame requires significant amount of memory. In most cases, the required amount of memory is too big to fit into the internal device memory forcing the use of external memory. Moreover, when targeting a PLD/FPGA device the use of a coarser granularity leads to higher entropy of the system that increase the complexity of logic required to model the system at the hardware level. That may result in higher fan-out, increase the occupancy of resources and extend the critical route latency [Koc13], meaning that significant trade-offs may have to be made in order to achieve performance that allows the SR processing to be applied in real-time. Thus, hardware implementations do not process coarse grain structures, but rather carry out the SR transformation for finer granularity elements. The most successful implementations available up-to-the-date apply SR on the finest grain possible — a single pel. This approach, known as the optical flow, allows to tailor the processing by breaking it into stages (iteratively executed for subsets of the data) that fit within a defined memory budget, and then use streaming to meet the target performance (throughput) [BB08, ABCC09]. Nevertheless, hardware implementation of optical flow faces problems of its own. Due to the extremely fine granularity of the flow, its implementation requires a very deep pipeline. Implementation of deep pipelines has to deal with the overhead introduced by inter-stages glue logic and/or latches, high logic occupancy (less room for optimization — low reuse of resources), high power dissipation due to repeated transfers of pels belonging to the reference structures (RS; corresponds to RF or SA). All of the above have a significant impact on the overall system efficiency in terms of resource occupancy. In many cases, these requirements are so high that only a sub-optimal version of the algorithm can fit in the device, leading to a sub-optimal quality of the super-resolved image. In this work we propose a novel flow for the non-iterative NUGPA. The proposed flow internally applies full SR processing on a small set of pels (called macro-block). By effectively forming a cross-over point between the flows that use coarse or very fine granularity elements, the proposed flow is expected to allow to combine their strengths while alleviating their weaknesses. 108 4.2 EXECUTION FLOWS OF THE NUGPA 4.2 Execution flows of the NUGPA As aforementioned, SRIR is typically carried out in two steps: extraction of data that characterizes the input and reconstruction (where possible) of missing data. The former task encapsulates a series of computations needed to estimate values of metrics characterizing the input data. These metrics are used in the latter step (called the super-resolution kernel, SRK) to identify pels from reference structures that when merged with the data being upscaled can lead to additional details reconstruction. The SRK processing ends producing a higher resolution image (called the super-resolved image). In order register the images and extract the hyperdata the NUGPA uses block matching (BM) motion estimation. Detailed description and evaluation of motion estimation is presented in [BAA05] and [CLT+08]. This section focuses solely on the SRK, presenting and comparing in detail the reference and the proposed execution flows. 4.2.1 Reference coarse grain execution flow The base for the study has been a software implementation of a complete SR system, comprising ME and NUGPA SR kernel coded using ANSI C. The SRK execution flow implemented by this software is presented in Fig. 4.1. The SRK flow comprises two internal subflows: one for the luma component and another for the chroma components processing. Chroma picture elements carry a significantly lower energy concentration than luma pels, and thus have a significantly lower impact on the subjective quality of the super-resolved image [IP91]. Therefore, super-resolution is applied only on the luma pels. Spatial resolution of chroma components is augmented using bilinear interpolation [GW08], meaning that no additional details are being reconstructed for these components. The luma pels undergo full NUGPA processing presented in Fig. 4.1: 1. The frame to be processed is loaded from memory. As the luma and chroma components are to undergo different processing, thus, when presented with a YUV formatted LR input, the software starts the processing by separating one type of the picture elements from the other. 2. Value ‘0’ is used in step 7 to identify the pels whose value has to be created using interpolation. This allows to eliminate the need of an extra memory to store the map of pels to be interpolated. In order to distinguish the input LR pels with value equal to ‘0’ (called zeroes) from the pels marked for interpolation (called holes) the original 109 CHAPTER4.–PROPOSED SUPER-RESOLUTION ALGORITHM Fig. 4.1 The frame-level execution flow of the NUGPA super-resolution kernel. value of the zero is changed to be equal to ‘1’. This process will be referenced to as the zeroes2ones (or zeros2ones) transformation. 3. Two HR representations of the transformed LR frame, the holesGrid and the srGrid, are created. The HR representations dimensions are several times bigger than the LR ones. For quarter pixel ME the growth is of a factor of 4 for each dimension (resulting in total size growth of 16 times). The srGrid is constructed directly from the output of the zeroes2ones transformation. For holesGrid, pels belonging to the borders (outer rows and columns) of the zeroes2ones outcome are duplicated before the spatial augmentation takes place. Border duplication is carried out in order to (somewhat) solve the aperture problem that could arise in ME. 4. Steps 1–3 are repeated until all frames from the sliding frame window have their HR representations. 110 4.2 EXECUTION FLOWS OF THE NUGPA 5. Candidate pels specified by the ME metrics are extracted from all holesGrid of frames belonging to the SFW. Extracted data are fused with the srGrid. 6. Each element of the srGrid is divided by the sum of weights associated with it. Steps 5–6 make up the shiftAdd transformation. 7. The outcome of the shiftAdd transformation is usually composed of many pels with value equal to zero. These pels values get estimated by computing the mean value of pels belonging to a square-shaped neighborhood centered at the hole’s position. The process in which a hole is assigned a value is called holesFilling. 8. The srGrid after holesFilling is the super-resolved frame, with the scale factor equal to precision of used in image registration (equal to 4 for quarter pixel precision). When a different scale factor is specified, srGrid has to be adjusted to the expected dimensions. This task is carried out by the scale function. This function discards certain pels, producing the final super-resolved representation of the luma component. 9. Finally, the synchronization of the luma and chroma flows takes place. This step is needed in order to provide the super-resolved frame in YUV format. 10. The SFW update policy is applied. Typically, one of the reference frames is discarded (in the FIFO order) and a new frame is loaded in its place. 11. Steps 1–10 are repeated for each processed frame. The presented processing is carried out at frame-level. Transitions from one step to another take place only after step’s transformation has been applied on all pels fo the frame being super-resolved. The above-presented processing requires that HR representations of all frames from SFW be available from memory throughout the complete duration of the super-resolution process. 4.2.2 Proposed finer grain execution flow The used SRiuma reference software originally operated at frame-level, posing memory storage requirements precluding the planned FPGA implementation using only the on-device available Block RAM (hereafter BRAM) memory. In this work we propose a novel flow for a non-iterative fusion algorithms that mitigates the problem of high memory occupancy associated with the presented coarse-granularity flow. In contrast with the frame-level processing, which requires that HR representations of all frames from SFW be available from 111 CHAPTER4.–PROPOSED SUPER-RESOLUTION ALGORITHM memory throughout the complete duration of the super-resolution process, the proposed approach requires that only search areas, instead of whole frames, be available from memory. As a result, the requirements for memory storage are reduced, at the expense of increasing the number of memory accesses. Processing of a patch instead of only a singular pel allows more optimized logic synthesis which is expected to result in moderate (lower than for the optical flow) logic occupancy allowing full algorithm implementation on the targeted FPGA device. The proposed execution flow of the NUGPA is presented in Fig. 4.2. This flow is similar to the FL flow (compare with Fig. 4.1). The main differences between the flows are: 1. the size of the basic element being stored in local memories, 2. substitution of reference frames by search areas, 3. the fact that transition to the subsequent step of the algorithm is carried out after one macro-block (not one frame) has been processed, and 4. the necessity of managing inter-macroblock dependencies, not present in frame-level processing, that significantly complicate the holesFilling process. These changes result increase the complexity of memory management and call for different buffering schemes in order to mitigate the expected increase in memory traffic. 4.3 Evaluation of memory occupancy of the reference and the proposed NUGPA execution flows In order to quantify the memory storage requirements and their reduction associated with the change to the proposed macroblock-level (MBL) flow a theoretical study based on the available software implementation has been carried out. First, the memories used by the implementations had been identified. Then, their size and its dependency on the SR parameters have been evaluated. These steps have led to the creation of equations determining memory size for each type of identified memories. Based on these equations the memory storage requirements of each step of the investigated execution flows have been quantified and compared. 4.3.1 Identification of memory storage requirements In order to provide a quantitative evaluation of the memory storage requirements of the reference and the proposed NUGPA execution flows: (i) the memory storage requirements of the flows’ steps had to be identified and modeled as a function of algorithm’s parameters, 112 4.3 EVALUATION OF MEMORY OCCUPANCY Fig. 4.2 The macroblock-level execution flow of the NUGPA super-resolution kernel. (ii) the memories instantiated by each step had to be identified and their size modeled using the equations established in step (i). 4.3.1.1 Memory types and their dependency on algorithm parameters From an analysis of the NUGPA flow and the structures defined in the code of the software implementation we have inferred that the parameters that determine the size of memories are: (i) the number of frame rows FRrows and columns FRcols from which the size of the LR representations is derived, (ii) the macroblock width MBwidth which determines the size of the macroblock (we assume MBs are square shaped), (iii) the ME precision precisionme which determines the size of the HR representations, (iv) the SR scale scale which determines the size of the super resolved representations, (v) search area radius (SAR) which determines the size of the search area, and (vi) the interpolation window intwindow which determines the size of the memories used to solve data dependencies in holesFilling. Short 113 CHAPTER6.–HARDWARE IMPLEMENTATION of memory instances if the total size of the memory after consolidation is smaller than the minimal (vendor and technology dependent) size of a BRAM instance. Data coalescing plays an important role in provision of sufficient memory throughput where it was most required. An example of this is the implementation of the entity responsible of carrying out the interpolation process (holes filling), presented in Section 6.3.3.2. This process required up to 24 pel values to be loaded for each encountered hole. In order to meet the performance threshold pels were concatenated into bundles of four 8 bit values (32 bit data word) and stored in four memories (four-times replication). This organization allowed loading of all the necessary pels in time of four memory accesses, leaving more time for processing and facilitating meeting of the performance requirements. The introduction of an indirect communication scheme, using an inter-module shared memory as shown in Fig. 6.8(d), solved the problem of long synchronization periods encountered in direct communication. As a result, the time spent on data unpacking became independent of the number of accesses required for (coalesced) data word creation, as long as these words were readily available from the shared memory. Relaxed performance requirements let the initiator module to carry out data word creation at its own pace and, only after all data have been processed, synchronize with the target. A side effect of this is a seamless coalescing scheme in which (sub)modules dedicated to data coalescing could be introduced into the system, in order to carry out the communication and data coalescing in parallel with computations. In cases when the memory throughput offered by this solution is not sufficient, the effective memory throughput could be further increased using memory replication —so that multiple concurrent accesses to the shared memory can be carried out in parallel. 6.3.3 Data dependency in holes filling The amount of the additional information possible to be extracted and fused with the to-besuper-resolved LR pels is determined by the acquisition setup, deployed image registration (motion estimation) performance and the super-resolution kernel parameters values. Some of the crucial parameters are the search area size and number of frames in the sliding frames window (SFW). With the increase of these parameter values, the possibility of encountering additional information increases as more candidate-pels are being used in motion estimation. However, even in case of a favorable setup, the usable information tends to be scarce and its distribution is non-uniform. A common scenario is the one in which for most of the coordinates of the HR representation of the macroblock there’s no additional information. 216 6.3 IMPLEMENTATION CHALLENGES In other words, the post-fusion pel value associated with the HR coordinate is zero. For the sake of satisfactory output quality, these coordinates cannot be left with this value. Thus, for all coordinates with zero value (holes) new values need to be estimated. The estimation of these values, in the case of the considered (SRiuma) implementation of the NUGP algorithm is carried out by means of the non-uniform mean-nearest neighbor interpolation. The new value to substitute the one held by a ‘hole’ pel is estimated based on its neighbor pels’ values. 6.3.3.1 Problem statement and implications The interpolation process requires that a square-shaped neighborhood of the coordinate which value is being computed, be readily available for access. The boundaries of this neighborhood are determined by the value assigned to the intwindow parameter and frame borders (which effectively crop the neighborhood). The intwindow parameter specifies the maximal distance (computed as the absolute difference of HR coordinates) from the tobe-filled coordinate that a pel coordinate can have in order for its value to be allowed to contribute to the process of new value creation. A graphical interpretation of the interpolation window and the way it is used to specify the neighborhood spread is presented in Fig. 6.9. In order to determine the set of pels needed for a MB interpolation it is necessary to consider the interpolation neighborhood of all pels that constitute a MB. Thus, the set that contains all pels that can contribute to interpolation of all pels of a MB is defined as a superposition of pels contained in neighborhoods of all pels contained within the MB. As shown in Fig. 6.9, this set comprises all pels of the MB currently being processed and a set of pels belonging to adjacent MBs that are within the distance of intwindow from the MB borders. When super-resolution is carried out at frame-level the necessary neighborhood pels are always readily available from the memory as the HR grids represent a complete frame. The switch to the macroblock-level processing leads to elimination of the frame-level memories, introduces data partitioning into MBs and significant changes in the execution flow. All of the aforementioned changes effect in a situation in which only a subset of pels making up the frame is available at a given moment. Thus, for some coordinates not all of the pels required in the interpolation process are available. Namely, a subset of pels from the MB borders that fall within the so called interpolation window requires pels belonging to the adjacent borders of adjacent MBs. For most MBs, the information necessary to carry out the processing for these pels may not be available at the time of reception of this macroblock. As a consequence of not all frame pels being readily available from the beginning 217 CHAPTER6.–HARDWARE IMPLEMENTATION Frame height Frame width Interpolation window intwindow window Fig. 6.9 Interpolation window placement on a MB grid with MB boundaries shown. of the process, data dependencies arise and cause complications in the implementation of the interpolation process. In particular, the rightside columns and the downside rows of the macroblock being processed, that fall inside the interpolation window, cannot be processed as not all of the pels that make up their neighborhoods are readily available. When processing is carried out at MB-level, the value of the intwindow and MBsize parameters, apart from determining the overall amount of data for which the interpolation process has to be carried out, also determine the current workload, defined as the data for which the processing can be started at the moment of MB reception. When looked at from a different perspective, these parameters specify the pels whose processing does require pels from different macroblocks, and which effectively, present data dependencies and cannot be processed at reception time. For constant (non-variable) size MBs, a MB can have up to 8 immediate neighbors. The actual number depends on MB’s location within the frame. The data that presents MB dependencies can be further divided internally into 8 regions with different data dependencies. A simplified representation of this division is presented in Fig. 6.10(a), where the dependencies are labeled as: U, L, R and B, meaning that they dependent on upper, left, bottom or right neighbor’s pels, respectively. Regions marked with labels consisting of two letters present dependencies on pels of two neighboring MBs. The dependencies are mutual, and none of the coordinates that present data dependency can be processed until both pels values are readily available. Considering that MBs are processed in raster scan order, those 8 dependencies (no dependencies being the 9th type of dependency) can be generalized as resolvable or non-resolvable at the time of reception. The former type comprises data de218 6.3 IMPLEMENTATION CHALLENGES UR R BR Upper neighbor Left neigh bor BR BBR L UUR R BR B BL UL Upper left neighbor Upper right neighbor BL BL No aparent dependency (a) Data dependencies for available pels. Resolvable dependency Possibly unresolvable dependency L U R BR B BL UL UR (b) Division into defined regions of dependencies. Fig. 6.10 Data dependencies regions and types within the interpolation window at the time of a macroblock reception. pendencies that depend only on previously processed MBs received and can be processed if the data were stored and are available. The latter type comprises dependencies that cannot be resolved at the time of MB reception as they reference data that will be provided in the future and cannot be made available beforehand. One should note that the dependencies are not equally present for all MBs of a frame. Some of the dependencies do not arise for MBs that belong to the frame borders as some of their immediate neighbors do not exist. Assuming raster scan order of execution and that the left and upper neighbors pels can be made readily available, the dependencies of the U, L and UL type are resolvable at the time of MB reception. On the other hand, regions with the B and R dependencies are non-resolvable at that time. When new data are introduced into the system and processed, the non-resolvable, at the time of reception, dependencies of the previously processed MBs may become resolvable. The R dependencies become resolvable as soon as the next MB (right neighbor) is to be interpolated (with the exception of the MBs that do not have neighbor to the right). The B type dependencies hold until a complete row of MBs is processed (with the exception of the MBs that do not have neighbors below, for which the B type dependencies are not present). A generalized view on the division of the interpolation neighborhood into regions based on dependency resolution at the moment of MB reception is presented in Fig. 6.10(b). An additional side-effect of data dependencies is the necessity of supporting processing workload which cardinality varies along the frame coordinates. The workload increases along with the horizontal and vertical coordinates as the macroblocks from second row/col219 CHAPTER6.–HARDWARE IMPLEMENTATION umn process some data of the upper/left neighbors. Macroblocks in the last row/column, due to the lack of the bottom neighbors are responsible of processing more of their data. The least and the most workload has been observed for the first and the last macroblock of the frame, respectively. 6.3.3.2 Implemented solution During the development the following issues had to be tackled: (1) the value of the intwindow parameter had to be determined, (2) the actual workload for each MB had to be determined, (3) future processing of the temporarily unresolvable dependencies had to be allowed, and (4) sufficient memory throughput had to be assured. 6.3.3.2.1 Interpolation window size specification. When processing is carried out at MB-level, the value of the intwindow parameter, apart from determining the amount of the workload of the interpolation process, also determines the set of pels that present data dependencies and the size of interpolation buffers needed for data dependencies management. For our implementation, we have considered the use of three values of the interpolation window, namely the window radius of 1, 2 or 3 pels was contemplated. Based on the observed experimental results for these configurations the interpolation window value has been fixed to 2. Use of this value has an advantage of preventing the interpolation windows of the nonimmediate neighbors to overlap due to its value being smaller than the half of the spread of the minimal MB width/height. 6.3.3.2.2 Workload determination. As aforementioned, the amount of data with resolvable dependencies is directly related with MBs placement within the frame. Macroblocks belonging to the frame borders present different data dependencies that macroblocks from further within the frame, and effectively required different amount of workload to be processed. Going into details, MB placement determines the existence (or nonexistence) of its immediate neighbors, which specifies which data dependencies hold. In our implementation, the actual workload is being determined based on the distance (in HR coordinates) of the MB borders from the frame borders as this determines the existence of MB’s neighbors and the limits of the neighborhood (that cannot extend beyond frame borders). Considering this criterion, 9 possible MB positions within the frame can be distinguished. The distinguished 9 MB positions and the corresponding codifications of them as MB placement flags are presented in Fig. 6.11. Gray squares without a number, centered on the squares that contain a number, represent the immediate neighborhood of the latter. The 220 6.3 IMPLEMENTATION CHALLENGES Frame height Frame width 1 5 2 8 96 4 7 3 MB position bottom_border_flag right_border_flag left_border_flag top_border_flag 1001 1 201 1 0 3110 0 4 1 0 0 1 50010 6010 0 7 1 0 0 0 80001 90 0 0 0 Fig. 6.11 Placements of MBs within a frame that change the number of pels within interpolation window. MBs located in each of frame corners (marked as 1–4 in Fig. 6.11), belonging to a frame border but not located in the corners (5–8), and not belonging to frame borders (9) are assigned a different type. The distinguished MB positions are codified as a set of four bits, referred to as MB placement flags, namely: bottom_border_flag (or f lagsMB bottom), left_border_flag (f lagsMB le ft), right_border_flag (f lagsMB right) and top_border_flag (f lagsMB top ). Each of the flags was set if the current MB’s border also belonged, respectively, to the bottom, left, right and top border of the frame. The f lagsMB top and f lagsMB bottom flags, are set if MB’s upper or bottom border corresponded to frame upper or bottom border (the MB is located in the first, or last row), respectively. In a similar fashion, the f lagsMB le ft and f lagsMB right flags are set if the MB was located, respectively, in the first, or last column of the frame. Based on these flags the actual workload was determined and enforced by means of limiting the range of addresses accessed during the execution of interpolation. In our implementation the MB placement is expected to be determined and codified outside of the SRK module and then passed to the super-resolution kernel as a part of the input. The MB placement within the frame was set by the testbench and then propagated along the data path being passed from one module on to next one along with the (corresponding) MB pels. 221 CHAPTER6.–HARDWARE IMPLEMENTATION 6.3.3.2.3 Data buffering. Data dependencies resulted in the necessity of buffering the data that cannot be processed at the moment of MB reception until the data dependencies can be resolved. The data that have to be buffered are the ones belonging to the regions with unresolvable dependencies. As presented in Fig. 6.10(b), with the regard to the received MBs, these regions included the intwindow right-most columns and down-most rows. The right-most columns’ pels could be processed with the following macroblock and, thus, are stored in local registers. The down-most rows’ pels have to be stored until the pels of the macroblocks below them are available. This means that intwindow rows of a whole row (plus 1, due to dependency on bottom-right neighbor) of macroblocks have to be buffered. This set proved to be too numerous to be stored in registers and its storage was implemented using standalone module encapsulating RAM instance referenced in this work as the intbu f f er. For the used interpolation method, the pels that make up the execution context are the ones presented in Fig. 6.10(a). In our implementation, at MB-level the actual execution of holes filling is carried out over a set that comprises all of pels belonging to the execution context that are readily available. This set corresponds to a subset of the set presented in Fig. 6.10(a) which comprises the MB pels plus the pels of the left and upper neighbor that are stored in the buffers. In our implementation this set is stored in the so called pel map. A high-level view of the pel map organization is presented in Fig. 6.12. Once available, the MB pels are stored in a dedicated pel map area (defined address range). The regions with unresolvable dependencies (labeled as B/BL/BR and R/UR/BR) are stored in registers and RAMs until the dependencies are resolved. At that time, these pels correspond to the upper and left neighbor of the currently processed MB. Thus, when loaded from the buffer they are stored in the pel map in regions mapped to the corresponding regions of these neighbors. The current content of the pel map and the buffers loads/stores are managed based on the MB placement. Actual workload is specified based on the range (or number) of addresses for which holes filling is to be executed. The used address range is limited based on the received MB placement flags. A thing to notice is that the actual pel map implementation uses 1D linear addressing in order to meet the requirements of RAM instantiation. 6.3.3.2.4 Memory throughput assurance. For implemented interpolation window parameter value of 2 (intwindow =2), up to 24 values can contribute to the process of one hole value estimation. For high number of holes the sheer number of accesses becomes a likely system bottleneck. In order to tackle high performance requirements that the holes filling process poses and to meet the target performance this process has been decomposed into two stages, one that 222 6.3 IMPLEMENTATION CHALLENGES Current MB BR R BR UR BL B BR MB height Upper neighbor pels Left neighbor Upper left neighbor pels MB width int intwindow window HR HR Store in registers Load from memory Load from registers Store in memory Fig. 6.12 Pel map organization with labels of regions loaded/stored in the buffers. carries out memory management and the other that implements the interpolation. Each of these stages was encapsulated as a separate SystemC module, the HolesFillingPrep and the MeanNearest module, respectively. In accord with the scheme presented in Section 6.3.2.2 data passing between these modules was carried out using a shared memory with doubled memory address space. This memory represented the pel map structure and acted as the output/input of the respective modules. The HolesFillingPrep module was responsible for loading the appropriate MB context pels, management of the interpolation buffers and data dependencies. There were many reasons justifying the carried out decomposition into submodules. (i) Decomposition into two submodules effectively (almost) doubles the cycle budget for the holes filling transformation, allowing to use more cycles for each of the tasks. (ii) Smaller modules result in better performance of the tools in terms of optimization of the critical path latency (at the cost of overall area growth) (iii) Decoupling of data preparation from computations lowers the cost (in terms of area) of the computation units replication as replication of data preparation logic is limited. 223 CHAPTER6.–HARDWARE IMPLEMENTATION Holes filling Normalize MB MBHR Holes filling prep Intbuffer Pel map Interpolation engine Pel map MB∗ HR Pel map Interpolation engine Pel map MB∗ HR MB collector / 8 / 8 /8 / 32 / 32 / 32 / 32 / 32 / 32 / 32 / 32 / 8 / 8 / 8 / 8 Memory/RAM Module Data / Fig. 6.13 Holes filling implementation after decomposition into submodules, performance and memory throughput optimization. (iv) Finally, the remaining idle cycles leave some room for customization, in terms of broader ranges of algorithm parameters and/or data, without the need of module’s redesign. Even the decomposition the holes filling transformation did not assure that the worst case scenario can be processed within the (now doubled) cycle budget. Additional steps, facilitated by the decomposition were introduced. The MeanNearest module was replicated, allowing for the workload to be split between two computational units. In order to increase the memory throughput, additional cycles were used on implementation of data packing. Before being stored in pel map, the received context pels (8 bit values) are packed into one 32 bit word. The pel map memory used for data passing was also replicated, resulting in total of 4 pel maps being instantiated (2 for each MeanNearest module). Memory duplication requires additional resources but effectively doubles the overall memory throughput allowing for the MeanNearest module to issue up to two accesses simultaneously loading quadrupled elements, effectively accessing up to 8 pels in one cycle. All of the above modifications and tweaks reduced the number of required memory accesses while also reducing the amount of processing that needs to be carried out by a singular MeanNearest module instance. The data flow of the holes filling implementation after the described modifications is presented in Fig. 6.13. 6.3.4 Variable size of search area In frame-level processing flow, all pels belonging to a search area (SA) are readily available from the frame buffer of the super-resolution kernel. After switching to MB-level processing, the frame buffer was removed. This resulted in the necessity of search area pels 224 6.3 IMPLEMENTATION CHALLENGES extraction from a different source and their propagation to the super-resolution kernel. 6.3.4.1 Problem statement In our implementation, the search area pels are extracted outside of the super-resolution kernel and fed to it as input. Only the SA pels which coordinates are valid frame coordinates are used in the SR processing. The SA pels which coordinates point outside of frame borders are not available and cannot be loaded. The pels that are not loaded can be, either, (i) assigned a value equal to zero (hereinafter zero padded) and sent to the kernel, or (ii) not sent at all. The former scheme is easier to be managed as the number of pels that are sent to the kernel, search area grid (storage structure), and the MB upper left corner coordinates within the search area grid are constant throughout the execution, for a given SRK parameters configuration. On the other hand, zero padding can be considered a waste of processing power and cycles. This waste could be avoided by sending a set comprising only the available pels, which cardinality is not constant. A side effect of not sending the padding pels is the decoupling of the geometrical relations of pels being send from the actual search area. This leads to the requirement of the actual search area to be reconstructed from the pels bundle being received prior to its use. In both cases, the way in which the frame borders clip the search area and effectively limit the regions of valid data has to be somehow signaled to the core. 6.3.4.2 Implemented solution In order to avoid the necessity of padding and only transfer the relevant (non-padded) pels, our implementation allows receiving SAs that contain variable number of pels. In order to do so, two issues had to be solved: (i) the structure of the received pels bundle and the limits of valid data had to be signaled to the kernel, and (ii) the placement of the MV reference point, the upper left corner of the MB, within the received search area to be determined in order to allow correct MV processing and data extraction further in the data path. Both of the above are determined by the MB and search area location within the frame and the search area radius. The number of pels in search area depends the SAR and MBwidth SRK parameters and the SA placement within the frame. This relation is similar in nature to the problem of determination of the number of neighbors that a MB has which solution has been presented in Section 6.3.3.2.2. The SA placement determination has been solved in the same manner, that is by providing a description of the SA placement codified using four values which are 225 CHAPTER6.–HARDWARE IMPLEMENTATION reduction in memory occupancy. Moreover, implementation of the vendor-specific RAM IPs would have to be enforced by manual modifications of the resulting RTL description using vendor specific libraries. This would significantly limit the portability of the design. Thus, a different solution was sought for. 6.3.6 Frame window size determination The concept of sliding frame window was used in order to provide support for dynamic super-resolution. The exact number of frames that are included in the SFW varies throughout the execution. This behavior requires some mechanism to be handled properly. 6.3.6.1 Problem statement As stated in Section 2.2.4.1, the SFW can be seen as composed of two sets of frames: the frame being processed and a neighborhood comprising its preceding and succeeding frames. For the first frame of the sequence there are no previous reference frames available, thus, the SFW only contains the frames succeeding the one being processed. Having processed each frame the SFW has to be updated to include the processed frame in the preceding frames set. Analogously, when reaching the end of the sequence the number of SFW frames decreases, as there are no further frames (past the sequence limits) to be added to succeeding frames subset of the SFW. 6.3.6.2 Implemented solution In our implementation, a variable number of frames in the SFW is managed by means of internal counters and external flags signaling state transitions. At any given moment, after being initialized, the system can be in one of three states, corresponding to the scenarios in which: (i) the preceding frames number has to be incremented, (ii) the SFW is fully occupied, and (iii) the succeeding frames number has to be decremented. The distinction between the preceding and succeeding frames set is not reflected in the implementation, as processing of these sets is the same, and the system needs to know only the number of reference structures to access. For the first frame there are no previous reference frames and more than the maximal allowable succeeding frames available. Thus, initialization is carried out by setting the number of reference frames to the number of maximal value of the following reference frames (set at synthesis time). This value is incremented with each following SFW update, until the maximal allowed number is reached. This corresponds to the SFW being fully 232 6.3 IMPLEMENTATION CHALLENGES occupied. Once the last frame of the sequence is contained within the limits of the SFW, in order not to reference frames beyond the limits of the sequence, the following update has to result in decrementing the number of reference frames. For the same reason, the number of reference frames is decremented with each following SFW update, reaching the maximal number of preceding reference frames at the time of processing of the last frame of the sequence. The frames number update conditions are summarized in (6.9), where, nis the number of frames in the sequence, ni RF represents the count of the reference frames in the SFW for ith frame of the sequence (from 1 to n), maxPFS and maxSFS, refer to the maximal cardinality of the preceding and succeeding frames set. ni RF =                      maxSFS i f i =1, ni−1 RF +1i f 1<i≤maxSFS +1 ni−1 RF i f maxPFS +1<i≤n−maxSFS ni−1 RF −1i f n −maxSFS <i≤n maxPFS i f i =n. (6.9) In the presented implementation, the presence of the first frame of a sequence and the last frame within the SFW limits is signaled by setting or clearing a corresponding flag. These flags are referenced as the new sequence (NS) and window limit reached (WLR) flag, respectively. These signals generation logic does not form a part of the SR Kernel, and their values are considered one of the input parameters of the core. In our case these signals are set by the testbench module. After being received from the input, these flags travel along the data path with location flags. These flags and an internally stored frame counter are used in the process of determination of the SFW cardinality carried out by the ShadPrep module. The frame count, WLR and NS flags are checked for the first MB of each frame. If the WLR signal is set, the end of the sequence is within the forward frame window limits and the number of frames comprising SFW is to be decreased. Otherwise, the frame count is checked. If the locally stored frame number is lower than the number of frames to be included in the preceding frames set, the number of frames is increased. The count of frames is incremented with first MB of each frame and zeroed when the new NS flag is set. The number of reference frames currently in the SFW is propagated to ReapSteps with first MB of each frame. This value is checked in order to determine the ShadSteps with which the ReapSteps should synchronize and from which it will request data. The implemented process of SFW cardinality determination based on flags values is the one from (6.10) and (6.11), where irepresents the locally 233 CHAPTER6.–HARDWARE IMPLEMENTATION TABLE 6.2 An example of the SFW cardinality computations, WLR flag and u(i)values for a sequence comprising 300 frames (n=300) and the maximal number of frames allowed in the proceeding and succeeding frames sets equal to 2 (maxSFS =2and maxPFS =2). Frame Flag Update ith ni RF WLR NS u(i) 1 2 0 1 N/A 2 3 0 0 +1 3 4 0 0 +1 4···297 4 0 0 0 298 4 0 0 0 299 3 1 0 -1 300 2 1 0 -1 stored frame number and u(i)represents the current iteration update. An example of the SFW cardinality computations, NS/WLR flags and u(i)values for a sequence comprising 300 frames ( n=300) and the maximal allowed number of the preceding and succeeding frames set equal to 2 (maxSFS =2 and maxPFS =2) are presented in Table 6.2. ni RF =   maxSFS i f i =1, ni−1 RF +u(i)i f 1<i≤n (6.10) u(i) =                −1i f NSi=0and WLRi=1, 0i f ni−1 RF <maxSFS +maxPFS and WLRi=0, 1otherwise. (6.11) 6.3.7 Division operation Division operation is not a construct synthesizable by the used HL-synthesis tools. Thus, this operation had to be either emulated by other operations (i.e. bit shifts (+ addition) or multiplication) or implemented using a custom hardware module. 6.3.7.1 Problem statement The typical software implementations use an algorithm based on counting the number of subtractions carried out until the reminder becomes negative (e.g. [HN10]). Algorithm 1 234 6.3 IMPLEMENTATION CHALLENGES Algorithm 1 Division by successive substractions algorithm. 1: procedure DIVISION(dividend,divisor)◃Division of dividend by divisor 2: reminder ←dividend 3: quotient ←0 4: if quotient ≥divisor then 5: while reminder >=0do ◃While the reminder is positive 6: reminder ←reminder −divisor 7: quotient ←quotient +1 8: end while 9: else 10: ... ◃ Handle division by zero exception. 11: end if 12: return quotient ◃Return the result 13: end procedure presents pseudocode of such an approach. This algorithm is known in literature as the generalized division. The main drawback of this algorithm is its variable execution time. Additionally, the execution time is dependent on the run-time values being used in the process. When a dividend >> divisor the number of iterations (∼cycles) required for execution is increased significantly. These properties prohibit efficient hardware implementation. In order to facilitate meeting the targeted cycle constraints, division was decided to be implemented as a separate module. This solution was expected to allow data load and store operations to be executed in parallel with the division, and facilitate performance customization based on division module replication or pipelining. In order to alleviate the aforementioned issues we considered two options: (i) use of a vendor specific divider IP or (ii) provision of a custom implementation. In most cases the former is the better choice, as it offers the best performance and its introduction requires the least effort. Nevertheless, in our implementation we have decided to tackle the problem of division operation implementation by provision of a custom divider module. There were two main reasons for that choice, namely customization possibilities and portability provision. The range of operations supported by the vendor-specific IP were considered to reach beyond our needs, leading to waste of resources allocated for execution of operations never used in our implementation. Moreover, use of vendor-specific solution leads to limited portability of the design and increases the effort required for migration to other FPGA devices. Apart from facilitating the possible migration, a custom implementation allows to optimize for the particular use scenarios encountered in the system. 235 CHAPTER6.–HARDWARE IMPLEMENTATION 6.3.7.2 Implemented solution During the development the following issues had to be tackled: (1) division algorithm suitable for hardware implementation in FPGAs had to be chosen, (2) the chosen algorithm had to be customized and prepared for implementation in SystemC, and (3) implementation had to assure meeting the targeted cycle/latency budget. 6.3.7.2.1 Suitable algorithm determination. In the process of implementation several division algorithms were taken into consideration, among others, the general,non-restoring, converge,pre-inverted divisor and reciprocal division algorithms. All of these algorithms are presented and evaluated, also in terms of possible hardware implementation, in [DBS06]. Generalized division algorithm was not implemented due to variable number of cycles that it required, that became prohibitively high for dividends much greater than divisors. Convergence algorithms were not used due to the large precision of the intermediate representations, required by these algorithms to converge, that required a number of bits surpassing the one supported by the hardwired multiplication blocks. The pre-inverted divisor methods were discarded due to the wide range of possible divisor values that would require large look up tables. 6.3.7.2.2 Base 2 non-restoring division of integers. The implementation presented in this work uses the base 2 non-restoring division algorithm. The used algorithm is a variation of the one presented in [DBS06]. Having Xrepresent the dividend, Dthe divisor, Qthe quotient and Rthe remainder such that X=Q·D+R,|R|<|D|,X,Q,D,R∈N(6.12) the base 2 non-restoring division algorithm can be executed as follows: (1) Scale the divisor, so that it is greater than the dividend. (2) Check the result of X·Dand determine first reminder (r0) and quotient digit value (q0). (3) Compute the subsequent preminders and quotient digits value, where pspecifies the required precision. (4) Correct the final reminder and quotient value if needed. Invert the scaling of the quotient. 236 6.3 IMPLEMENTATION CHALLENGES Scaling of the divisor allows simplifying hardware implementation of the division. Correct scaling should result in the divisor having additional ‘leading’ digit over the dividend. The only side effect of (over)shifting too much to the right could be the generation of ‘redundant’ sign extension digit. Let mdefine the scaling implemented as shift by mpositions to the left. The result of X·2m·Dcan be determined by looking at sign bits of the arguments. If X·2m·D< 0, then the first reminder (r0) is incorrect, thus, the corresponding digit of the quotient (bit in case of base 2 system) (q0) is cleared. Otherwise, the reminder is correct and the corresponding digit is set as in (6.13). r0=   X+2m·D i f X ·D<0, X−2m·D i f X ·D≥0q0=   0i f r0·D<0, 1i f r0·D>0. (6.13) The following reminder is computed based on the correctness of the current reminder. Updating the remainder can result in: (i) correct remainder ri, when ri·D>0 — corresponding to quotient digit qi=1, or (ii) incorrect remainder ri, when ri·D<0 — corresponding to quotient digit qi=0. In case of the reminder ribeing incorrect, then: (i) the correct subsequent remainder is the ri=2·ri−1, (ii) the base for next quotient digit determination is ri+1=2·(2·ri−1)−D, and (iii) the value of base for next quotient digit determination could be obtained as a result of delayed correction by adding Das in ri+1=2·(2·ri−1−D)+D. Defining the current reminder as ri=2·ri−1−D, the subsequent p<mreminders (ri+1) are derived from the preceding one (ri−1) and the current quotient digit qi, as follows: qi=   0i f ri<0, 1i f ri≥0ri+1=   2·(2·ri−1−D)+D=2·ri+D i f qi=0, 2·(2·ri−1−D)−D=2·ri−D i f qi=1. (6.14) The remainder obtained in the last step (i=p) could also be incorrect. In the case of rp·D<0, remainder restoration, carried out by addition of D·2m, is required. Finally, the scaling applied to the divisor has to be inverted in order to obtain the final quotient value. The quotient scaling is carried out as in (6.15). Q=q0·2m+q1·2m−1+q2·2m−2+···+qp−1·2m−p+1+qp·2m−p.(6.15) The above-described algorithm forms the base and reference used for SystemC code generation and validation. The control flow of the resulting implementation of the base 2 237 CHAPTER6.–HARDWARE IMPLEMENTATION Algorithm 2 Base 2 non-restoring division algorithm with delayed correction. 1: procedure DIVISON(dividend,divisor,m,p) 2: ◃Division of dividend by divisor with precision p 3: divisor2m←divisor << m◃Prescale the divisor by the factor m>p 4: {quotient0,quotient1,...,quotientp−1,quotientp}←0 5: if dividend ·divisor ≥0then 6: reminder0←dividend −divisor 7: else 8: reminder0←dividend +divisor 9: end if 10: i←0 11: while i≤pdo ◃While the target precision pis not reached 12: if reminderi·divisor ≥0then 13: quotienti←1 14: reminderi+1←2·reminderi−divisor2m 15: else ◃In case of reminderi·divisor <0 16: quotienti←0 17: reminderi+1←2·reminderi+divisor2m 18: end if 19: i←i+1 20: end while 21: quotient = {quotient0,quotient1,...,quotientp−1,quotientp} 22: if reminderp·divisor <0then 23: reminderp←reminderp+divisor2m 24: end if 25: return (quotient,reminderp>> m)◃Return the result 26: end procedure non-restoring division is presented in Algorithm 2. 6.3.7.2.3 Division implementation. The execution time of the non-restoring division depends on the required precision pas this parameter determines the number of main loop iterations. In our implementation required precision is determined by the number of bits of the greater of the dividend and divisor. This relation between the execution time tdiv of the division and the bit width of the operands can be approximated as in (6.16). tdiv ∼max(dividendbitwidth,divisorbitwidth)≈log2(max(divident,divisor)) (6.16) In our implementation, division is required during the fusion and interpolation stages. For the former, the required number of bits satisfactory for divisor representation depends on 238 6.3 IMPLEMENTATION CHALLENGES the maximal number of frames in the reference frames in the SFW (nmax RF ) and maximal (allowed) value of the fusion weight (weightmax f usion). In case of the latter stage, the maximal divisor bit width is derived from the cardinality of the pel neighborhood specified by the interpolation window and the maximal value of the contribution weight (weightmax int ). For both stages, the dividend adds additional dependency on the maximal pel value (pelmax val ) with respect to the divisor dependencies. Thus, in both cases the dividend bit width dominates the overall requirements. These relations are summarized in (6.17) and (6.18), and, (6.19) and (6.20), respectively, for the fusion and interpolation stages. divisorbitwidth f usion ∼log2⌈nmax RF ⌉×⌈weightmax f usion⌉(6.17) dividendbitwidth f usion ∼log2⌈nmax RF ⌉×⌈weightmax f usion⌉×⌈pelmax val ⌉(6.18) divisorbitwidth interpolation ∼log2(int2 window +1)2×⌈weightmax int ⌉(6.19) dividendbitwidth interpolation ∼log2(int2 window +1)2×⌈weightmax int ⌉×⌈pelmax val ⌉(6.20) Substituting the parameters with their respective values we obtain the same requirement of 13 bits. Thus, implementation arguments bitwidth can be limited to 13 bits. dividendbitwidth f usion ≈log2(16×1×256)<log2(25+8)<13 (6.21) dividendbitwidth interpolation ≈log2(25×1×256)<log2(25+8)<13 (6.22) The use of 13 bits to represent the value of dividend and 12 to represent the divisor resulted in the necessity of 12 rounds and a latency of 13 cycles for internal computations [DBS06]. Assuming that capturing the input and latching the output require 1 cycle each, the total number of cycles can be estimated as Nrounds +3. For the mentioned configuration this gives 18 cycles meaning that up to 160 divisions could be carried out during the targeted 2867 cycles limit. The NormalizeMB had to be capable of performing up to 240 (holes) divisions and 256 load/store operations in the time reserved for one 4x4 MB processing. Thus, the limit in the number of cycles for a division was estimated at around 9 cycles. This was more than the maximal number of cycles required for division to be carried out. In order to meet the targeted performance, the division module could have been replicated with the workload being distributed to multiple instances. Nevertheless, this solution would further complicate the execution flow and require additional logic for available division input/output ports selection. Thus, in the final implementation the divider module was not replicated. Instead, the division was decomposed in to two stages that were later pipelined. This allowed the module to start and terminate a new division every 9 cycle (with 239 CHAPTER6.–HARDWARE IMPLEMENTATION full pipeline). Finally, a pipelined version of the divider module was implemented. In this version the division loop was split into two stages. The internal loops of the two stages are basically equivalent, with the exception of the starting values and ending condition. The second stage includes additional code for final reminder and quotient computation resulting in additional cycle of latency. This latency is matched by an additional 1 cycle stall in the first stage of the division resulting in Nrounds +4 being required for the division result to be readily available at the output. Nevertheless, by using a pipelined architecture the effective throughput (with loaded pipeline) is almost doubled, as the input/output data can be set every 9 cycles. This performance has been found sufficient for our implementation needs. Once the division implementation met the performance target, the existing execution flow of the modules using division had to be modified in order to allow efficient outsourcing of computations and concurrent execution. These modules shared a bidirectional relationship with the newly created division modules. As aforementioned, this type of relationship requires tight synchronization scheme, in order to maximize periods of non-blocking execution. The processing carried out by the divider module has been pipelined allowing one division to be issued every 9 cycles. To maximize performance the pipeline should be kept full and two division operations should be issued before the first outcome is received. If the processing carried out in the internal loop of the outsourcing module is organized in a way that allows synchronizing with the divider every 9 cycles (the number of cycles required for a pipeline stage execution) then the maximal pipeline utilization can be achieved reducing the division operation processing time by a factor of (nearly) 2, while avoiding the consequences of replication of the division module. 6.4 Results The results of the hardware implementation of the MB-level NUGPA core in and advanced FPGA technology are presented in Table 6.3. Implementation in the xc5vf70t-1 [Xil09] device (labeled as Tech1in TABLE 6.3) for MB size of 4x4 pels, scale equal to 2, and SFW containing 3 frames (RF=2) occupied 10291 LUTs, 16 DSP blocks, and 66 Block RAMs. This corresponds to device resource occupancy of, respectively, 22%, 12% and 44.6%. Total memory usage was 2.376 MB (44%). No off-chip memory was required. The final clock frequency is of 109 MHz, resulting in the capability to apply 2 x super-resolutions on 25 QCIF frames per second. 240 6.4 RESULTS TABLE 6.3 Resulting hardware performance for two FPGA devices in comparison with [BB08] for the particular case of RF=2, scale=2, SA=2, and MB=4. N/A stands for not available. SRiuma Reference [BB08] Technology Iteration stages [Xil07] Tech1[Xil09] Tech2[Xil07] 10 20 Occupancy [LUT] 10291 13031 35707 68317 BRAM 66 132 134 234 DSP blocks 16 2 N/A N/A Frequency [MHz] 109 68 58 58 Frame rate [fps] 25 16 61 61 6.4.1 Core-level implementation evaluation Apart from presenting the observed implementation results, Table 6.3 contains data that allow a direct comparison with the state-of-the-art implementation presented in [BB08]. In terms of raw performance the MBL level implementation, offers lower frame rate and supports lower range of spatial resolutions. Nevertheless, the competitive implementation only features 10 iterative stages, that is, roughly a half of the number of 20 stages assigned by the authors as the threshold for producing satisfying outcome quality. Thus, in order to offer a better comparison of implementation efficiency, an approximation of the resources required to implement a solution featuring 20 iteration stages is also presented. The logic and memory requirements for this configuration were estimated by addition of the cost of implementation of additional 10 iteration stages to the cost of implementation of the system already sporting 10 iteration stages. On the other hand, the implementation presented in this work was also mapped and routed for the device [Xil07] targeted in [BB08]. The results for this implementation are labeled as Tech2in Table 6.3, and can be used for direct comparison. The implementation presented in this work was optimized in order to offer base software quality while minimizing resources utilization. In order to measure the efficiency of the FPGA implementations, the LUTs per SR boost cost (boostCostSR) figure of merit is proposed. This figure represents the number of look-up tables the implementation required in order to offer an increase of 1% in measured average (objective) quality of the super-resolved image over the image produced using bilinear interpolation. It is desired to keep this figure as low as possible. The boostCostSR is estimated as in (6.23), where the LUTnr stands for the number of LUTs required for implementation and the PSNRSRi and the PSNRintivalues represent, respectively, the PSNR value of the super-resolved, and 241 CHAPTER7.–CONCLUSIONS AND FUTURE WORK 7.2 Conclusions 1. The review of the SR algorithms taxonomy presented in this work has identified the non-iterative direct and iterative back projection as the algorithms of the multi-frame SR family best suited for hardware implementation. A closer look at the state-of-theart of hardware implementations in FPGA has highlighted: (i) the trend of using execution flows that apply SR on fine-grain elements in order to increase exploitable parallelism and reduce memory utilization, (ii) the fact that the iterative approaches have to settle for sub-optimal output quality due to being limited by the amount of available resources, and (iii) the predominant impact of on-device memory occupancy and off-device memory access latency on the system’s capabilities. This review has demonstrated that there are still many challenges within the field in terms of hardware implementations efficiency and quality, as well as facilitation of hardware implementation by considering the context of limited resources typical of FPGA targets. 2. This doctoral thesis has tackled key contemporary challenges faced at the time of developing hardware implementations of the SR techniques for video sequences, namely: the delivery of real-time execution capability, provision of high implementation efficiency, and preservation of software-level quality of the super-resolved image observed at the output. The ultimate result of this thesis has been to provide a hardware implementation that is characterized by the above-mentioned properties. 3. We have created a complete taxonomy of the different approaches and algorithms developed in the literature to carry out super-resolution. 4. We have introduced mathematical models for the image registration stage and for the reconstruction stage of the non-iterative restoration-interpolation process. In an attempt of complementing the existing algorithmic descriptions with a consistent closed form analytical model that captures the layered computation scheme of the SR process. The different transformations operating over different image spaces are clarified and underlined. 5. The starting point of this research work was a software implementation of the noniterative version of the non-uniform grid projection algorithm coded in ANSI C. 248 7.2 CONCLUSIONS The review of the state-of-the-art has contextualized the NUGP algorithm as belonging to the interpolation-based family of the direct classical multi-frame superresolution algorithms, which is considered suitable for efficient hardware implementations. The presented quantitative results of the study on the objective quality of the super-resolved images produced by the reference software implementation prove the capability of the algorithm to provide quality better than other algorithms used by state-of-the-art hardware implementations. Considering the single-pass nature of the non-iterative version of the NUGP algorithm, we have concluded that this version of the algorithm meets the prerequisites for delivering a hardware implementation that meets real-time performance while preserving high quality of the super-resolved images. 6. The base software implementation did not consider hardware implementation challenges and was characterized by high use of memory resources precluding its direct implementation in FPGA technology using only on-device memory. With the goal of overcoming these weaknesses, we have proposed a modified execution flow that applies SR on a finer-grain element When implemented, this flow has removed the algorithm’s dependency on frame-level buffers requiring only single MB execution context to be readily available from on-device memory. The theoretical models, developed based on the software implementation, of the proposed and reference flows have been used to identify these flows’ bottlenecks, and quantitatively evaluate the gains and expenses associated with the change to the proposed finer-grain flow. 7. The results of the proposed modifications lead to significant memory occupancy reduction at the expense of increased memory traffic. The computed minimal and maximal value of the expected factor of reduction in memory occupancy associated with the switch to the MB-level flow was within the range of 3.5 to 16, depending on the SRK parameter values. The computed minimal and maximal value of the expected factor of increase in memory traffic associated with the change was within the range of 1.1 to 16.9. The exact value of these factors depends on the parameters value combination used to configure the SR kernel. For 87 out of the 96 investigated configurations (more than 90%) the factor of memory occupancy reduction has been greater than the factor of increase of memory traffic, proving the effectiveness of the proposed approach. 8. For the needs of the planned hardware implementation we have developed a highlevel implementation methodology. The established methodology defines a hierarchy 249 CHAPTER7.–CONCLUSIONS AND FUTURE WORK of levels of abstractions, in which each of the models is obtained from the model being one level higher in the hierarchy. This renders the models tightly-coupled, facilitating inter-model propagation of modifications. By using SystemC to code the intermediate models we have been able to reuse most of the ANSI C code. This facilitates rapid propagation of modifications made at the functional level through the intermediate models down to the HDL description. The use of intermediate representation increases design portability by allowing high-level synthesis to target a range of HDL languages (VHDL, Verilog, etc.) from the same higher-level description. Furthermore, by using C-based languages throughout the abstraction hierarchy we leave open the possibility of rapid migration of the models to other C-based languages, most notably the OpenCL standard. The level of details of the provided methodology description facilitates the reuse of the established flow in implementations of other similar algorithms. 9. We have provided a hardware implementation that is characterized by real-time performance and efficiency, while preserving portability and software-level quality of the output super-resolved images. The final architecture has met the targeted performance of 24 fps with operating frequency of 109 MHz using the xc5vf70t-1 device (Xilinx Virtex5 technology). The carried out comparison with the state-of-the-art has yielded satisfactory results as the observed logical resources occupancy for the proposed system was up to 5 times lower than the one reported for the state-of-the-art mapped using the same target FPGA technology [BB08]. The hardware implementation synthesis results: (i) have demonstrated the capability of reaching real-time performance in FPGA technology, while preserving the quality of the super-resolved images at the level offered by software implementation, and (ii) have proved the correctness of the presented algorithmic level changes and the established implementation methodology. 10. The portability of the design has been preserved by using several tiers of abstraction code in generic SystemC which allow the high-level description to be synthesized for a range of different languages (VHDL, Verilog, etc.) and targets (Altera, Actel, Xilinx, etc.) without (or with) vendor specific optimization enabled. The portability of the design is further reinforced by not using (explicitly) any vendor-specific accelerators (MACs, adders, dividers, etc.). To recap, this thesis: (i) has successfully carried out knowledge transfer from the domain of compression algorithms to the domain of SR algorithms. (ii) has quantitatively evaluated 250 7.3 FUTURE WORK the impact of algorithm modifications, (iii) has provided a high-level methodology of implementation that facilitates rapid propagation of modifications, and (iv) has delivered an efficient hardware implementation in FPGAs that reaches real-time capabilities while preserving the quality of the output super-resolved images. The observed synthesis results justify the claim that the proposed implementation contributes to the advancement of the state-of-the-art by implementing the MB-level processing flow (typical flow of image/video compression algorithms) and by delivering higher implementation efficiency. 7.3 Future work This work accomplishments leave room for new research lines focused on pushing further the state-of-the-art of SR implementations in hardware: • The proposed implementation could be explored for further improvement, in particular: (i) Implementation of the algorithmic level improvements proposed by Quevedo in [Gut15]. These enhancements would improve the robustness of the fusion kernel by detecting and eliminating outliers, leading to a further increase in the observable output quality. (ii) Extension of the range of supported SRK parameters values. Additionally, explore the possibilities of supporting variable macroblock size used in recent compression algorithms. (iii) Mapping onto heterogeneous platforms and encapsulation as an IP. This task would require an exploration of changes that could lead to a better hardware/- software partitioning, motion estimation methods, buffering scheme, elimination of adapters, and appropriate interconnections interface selection (e.g. AHB, AXI, etc.). • The work of Singla et al. [STdA13] presents an approach that uses multiple kernels of NUGPA applying SR in parallel on exclusive chunks of the input frame in order to speed-up the processing time. The presented approach is limited by memory accesses latency. As our approach uses only on-device memory its extension to such a configuration could allow alleviating these issues and is expected to result in significant execution speed-up. To provide performance scalability by means of replication, the architecture proposed in this work would require some modifications. In this line, it 251 CHAPTER7.–CONCLUSIONS AND FUTURE WORK would be interesting to explore: (i) further possibilities of reduction of memory occupancy, (ii) the possibility of elimination of the frame-level adapters from the SR kernel (in order to prevent their replication), and (iii) methods for efficient workload distribution. Apart from performances increase, the replication of kernels opens the possibility of supporting the multi-camera approach proposed in [Gut15], in which each camera’s input is processed by a dedicated kernel instantiation. • Multi-core and many-core platforms are undoubtedly the emerging trend in high performance computing (HPC). The appearance of hybrid solutions that use these platforms in tandem with FPGAs [PCC+14, WHK+14, MJK12], suggests new ways of conducting design space exploration for SR mappings onto these platforms. It would be of great interest to implement a SR algorithm onto several of these platforms, compare results of the mapping, analyze the arising implementation challenges, and the algorithmicand architecture-level modifications their efficient mapping would require. • In this line, we could take advantage of the ANSI C and SystemC implementations of NUGP algorithm developed for this thesis, as this code provides straight forward route for implementation using C-based languages like OpenCL or CUDA. OpenCL would be the preferred choice as it is an open and well defined standard supported on many platforms (CPU, GPUs, MP-SoCs, etc.) including FPGAs [Alt15]. OpenCL is expected to allow re-using most of the code in mappings on various platforms (NVIDIA Tegra 1000, Kalray, Intel Xeon Phi, Movidius Myriad 2, Amd APP, Altera Cyclone V, Xilinx EX-850 etc). The methodology established in this work is expected to facilitate this process as ESL flows using C-based languages (e.g. SystemC) have been reported to be the preferred starting point for OpenCL implementations in FPGAs [Eco12]. • The established methodology has been successfully used for mapping of SVC deblocking filter onto FPGA [CEN+13a, CEN+13b]. These studies extend the established methodology by describing steps that lead to a successful HDL migration to ASIC technology. It would be of interest to provide a manual (HDL) mapping and compare the results of the ESL-centered and manual HDLs mapping, both onto FPGAs and ASICs. This would allow quantitative evaluation of the difference in quality of the resulting HDL descriptions produced by these flows. • Finally, it would be interesting to use the established methodology for implemen252 7.3 FUTURE WORK tation of other high complexity image processing algorithms and/or using different toolchains. This would (i) further reinforce the credibility of the methodology’s effectiveness, robustness and utility in other applications, and (ii) allow revising the methodology’s strengths and limitations. Additionally, it would be of use to compare the SystemC-based flow with flows centered on using other high-level of abstraction languages. 253 Appendix A Publications International conferences [C1] Tomasz Szydzik, G.M. Callico, and Antonio Nunez, Quantitative Modelling of Image Processing Algorithms for Hardware Implementation, in Design of Circuits and Integrated Circuits (DCIS), Conference on, November 2015, Accepted for presentation. [C2] Tomasz Szydzik, G.M. Callico, Antonio Nunez, F. Tobajas, R. Sarmiento, and E. Quevedo, Optimization of non-uniform grid projection image super-resolution algorithms by reduced granularity and modified addressing, in Design of Circuits and Integrated Circuits (DCIS), Conference on, vol., no., pp.1–6, 26–28 Nov. 2014, doi: 10.1109/DCIS.2014.7035527. [C3] Pedro P. Carballo, Omar Espino, Romen Neris, Pedro Hernandez-Fernandez, Tomasz M. Szydzik, and Antonio Nunez, Scalable Video Coding Deblocking Filter FPGA and ASIC Implementation Using High-Level Synthesis Methodology, in Digital System Design (DSD), Euromicro Conference on, vol., no., pp.415–422, 4–6 Sept. 2013, doi: 10.1109/DSD.2013.52. [C4] Pedro P. Carballo, Omar Espino, Romen Neris, Pedro Hernandez-Fernandez, Tomasz M. Szydzik, and Antonio Nunez, Implementation of scalable video coding deblocking filter from high-level SystemC description, VLSI Circuits and Systems VI, SPIE Conference on, vol. 8764, pp. 876408, 28 May 2013, doi: 10.1117/12.2016885. [C5] Tomasz Szydzik, G.M. Callico, and Antonio Nunez, Real-time Implementation in FPGA of a Super-Resolution Algorithm using Macro-Block Execution Flow, in De255 Appendix A. Publications sign of Circuits and Integrated Circuits (DCIS), Conference on, Poster abstract, Avignon France, Nov. 2012. [C6] T. Szydzik, G.M. Callico, and A. Nunez, RealTime Video Restoration using FPGA Devices, in ACACES 2012 Poster Abstracts, Belgium: Academia press, Ghent, pp. 181, Fiuggi, Italy, 8–14 July 2012, http://hdl.handle.net/1854/LU-3234245 . [C7] Tomasz Szydzik, G.M. Callico, and Antonio Nunez, A Low Memory Requirements Execution Flow for the Non-Uniform Grid Projection Super-Resolution Algorithm, in Multimedia (ISM), 2011 IEEE International Symposium on , vol., no., pp.85–90, 5–7 Dec. 2011, doi: 10.1109/ISM.2011.22. [C8] Tomasz Szydzik, G.M. Callico, and Antonio Nunez, Closing the gap between software and hardware super-resolution image reconstruction: provision of high-quality output, VLSI Circuits and Systems V, SPIE Conference on, vol. 8067, pp. 80670M, 3 May 2011, doi: 10.1117/12.888253. [C9] M. Thadani, T. Szydzik, P. P. Carballo, P. Hernandez, Gustavo M. Callico, and A. Nunez, ESL quantitative assessment of the optimization steps in an ESL flow for a hardware implementation of a H.264/AVC decoder, Digital System Design, Architectures, Methods and Tools, 12th Euromicro Conference, 2009. [C10] M. Thadani, P. P. Carballo, P. Hernandez, Gustavo M. Callico, T. Szydzik*, and A. Nunez, ESL flow for a hardware H.264/AVC decoder using TLM-2.0 and high level synthesis: a quantitative study, VLSI Circuits and Systems IV, SPIE Conference on, vol. 7363, pp. 73630K, 28 May 2009, doi: 10.1117/12.821647. Journals [J1] T. Szydzik, G. M. Callico, and A. Nunez, Efficient FPGA implementation of a highquality super-resolution algorithm with real-time performance, in Consumer Electronics, IEEE Transactions on, vol.57, no.2, pp.664–672, May 2011, JCR=1.045, Q2, doi: 10.1109/TCE.2011.5955206. 256 Appendix A. Publications Other publications [O1] Tomasz Szydzik, and D Moloney, Precision Refinement for Media-Processor SoCs: fp32 -> fp64 on Myriad, Hot Chips: A Symposium on High Performance Chips, IEEE Technical Committee on Microprocessors and Microcomputers, in cooperation with ACM SIGARCH, poster, Flint Center, Cupertino, CA, August 10–12, 2014, poster, acceptance rate < 4.9%. [O2] Tomasz Szydzik, Marius Farcas, Valeriu Ohan, and David Moloney, Level-3 BLAS on Myriad Multi-Core Media-Processor, Hot Chips: A Symposium on High Performance Chips, IEEE Technical Committee on Microprocessors and Microcomputers, in cooperation with ACM SIGARCH, poster, Flint Center, Cupertino, CA, August 10–12, 2014, acceptance rate < 4.9%. [O3] T. Szydzik, A. Nunez, L. De Paepe, and L. Albani, Contributions to visualization algorithm enabling GPU-accelerated image displaying for dual panel high dynamic range LCD display, in Image and Signal Processing and Analysis (ISPA), 8th International Symposium on, vol., no., pp.483–488, 4–6 Sept. 2013, doi: 10.1109/ISPA.2013.- 6703789. [O4] T. Szydzik, G.M. Callico, and A. Nunez, Real-time Implementation in FPGA of a Super-Resolution Algorithm using Macro-Block Execution Flow, in Design of Circuits and Integrated Circuits (DCIS), Conference on, poster, Avignon France, Nov. 2012. [O5] T. Szydzik, G.M. Callico, and A. Nunez, RealTime Video Restoration using FPGA Devices, ACACES 2012, poster, Fiuggi, Italy, 8–14 July 2012. [O6] M. Marrero-Martin, T. Szydzik, J. Garcia, B. Gonzalez, and A. Hernandez, Effect of separation and depth of N+ diffusions in the quality factor and tuning range of PN varactors, VLSI Circuits and Systems V, SPIE Conference on, vol. 8067, pp. 80670T, 3 May 2011, doi: 10.1117/12.888663. 257 Appendix B. Resumen en castellano por SR se puede llevar a cabo tal y como se muestra en la Fig. B.1. A menos que el movimiento sea conocido a priori, antes de combinar las imágenes éstas tienen que ser registradas para determinar las variaciones entre ellas. Desde el momento en el que las imágenes de entrada representen la misma escena y puedan ser registradas de forma satisfactoria, el flujo de SR permite que cualquier método de adquisición sea usado. El pos-procesado depende del núcleo de SR usado y del proceso de aplicación de la SR. La mayoría de los algoritmos del estado del arte llevan a cabo una regularización post-fusión y restauración/reconstrucción. B.2.3.2 Técnicas directas La base teórica de las técnicas directas es la teoría de muestreo no-uniforme, la cual permite la reconstrucción de funciones de muestras tomadas de una posición distribuida no uniformemente. Esos métodos siguen la estrategia SR clásica presentada en la Fig. B.1. Dado un conjunto de observaciones de LR, una de las imágenes de LR es elegida como imagen objetivo y las otras son registradas frente a ésta. A continuación, la imagen objetivo es escalada usando un factor de escala específico, y las otras imágenes de LR son mapeadas (es decir: escaladas y desplazadas) dentro de la cuadrícula objetivo, usando la información obtenida en el registro. Los datos perdidos son obtenidos por medio de una interpolación no uniforme. Finalmente, se podría aplicar un núcleo opcional de enfocado/refinamiento a la imagen resultante. El proceso de fusión es normalmente implementado como una suma de valores multiplicados por pesos, y es por lo que esta estrategia normalmente se referencia como desplazamiento y suma (shift and add) [FREM03, FREM04b, FREM04a, FM06, AIAOM13]. Muchos de los métodos directos recientes llevan a cabo el registro y la fusión descomponiendo la imagen en zonas menores o bloques. Eso significa que las imágenes de LR son primero dividas en bloques. Luego cada bloque de las imágenes objetivo es comparado con el conjunto de bloques, incluyendo el correspondiente bloque y algunos bloques de su alrededor, en las imágenes LR de referencia. Basándose en la similitud de esos bloques y en sus similitudes con respecto al bloque actual, se asociado un peso a cada bloque, indicando la contribución de ese bloque en el proceso de producir un bloque de salida. Esta estrategia ayuda a manejar mejor una oclusión, local y movimientos complejos (por ejemplo, cambios faciales de la expresión). La familia de métodos que hace uso de este tipo de procesamiento es normalmente referenciada como SR directa no-paramétrica [PE09, PETM09, TMPE09, BBV10, CCL10, HS12, ZLRG12]. 264 Appendix B. Resumen en castellano B.2.3.3 Súper-resolución de secuencias de vídeo Los métodos existentes para la SR de secuencias de vídeo pueden ser clasificados en las siguientes 4 categorías: (i) métodos basado ventana deslizante de fotogramas (Sliding Frame Window, SFW) [SKR07, NHBS07, NSLZ07, PJ07], (ii) método secuencial [EF99c, EF99b, FEM06, CB07], (iii) método simultáneo [BS99, ZM07, AMK03], and (iv) método basado en el aprendizaje [BBM03, DKA04, KHX+06]. Las 3 primeras estrategias requieren un contexto multi-imágenes/fotograma, mientras que la SR basada en aprendizaje es capaz de ejecutarse en el contexto de un solo imagen/fotograma. Considerando asegurada la disponibilidad de las secuencias de vídeo, la técnica de SR multi-imagen constituye la estrategia preferente, especialmente en el contexto de una implem entación hardware. Las estrategias multi-imagen requieren de un cierto número de fotogramas para ejecutarse. Este conjunto de fotogramas, llamado ventana de fotogramas (o ventana de trabajo), forma el contexto de la ejecución. La ventana de fotograma comprende dos tipos de imágenes: la imagen objetivo que está siendo súper-resuelta y un número de imágenes de referencia (o fotogramas de referencia (RF)). Una consecuencia de usar un contexto multi-imagen para la ejecución del algoritmo, es necesario registrar los fotogramas pertenecientes al contexto actual. En el contexto de SR de secuencias de vídeo hay dos tipos principales de estrategia que están siendo utilizadas para el registro de imagen y actualización del contexto: el método de fijación y el método progresivo. En el método de fijación, uno de los fotogramas del contexto es elegido como el fotograma de referencia y los otros fotogramas desalineados se registran en relación a este fotograma. En el caso del registro progresivo, el fotograma actual es registrado en relación a sus vecinos temporales inmediatos a partir de los fotogramas anteriores. El algoritmo NUGP implementa la estrategia de SR basada en la ventana deslizante ilustrada en la Fig. B.2. Esta estrategia utiliza un registro fijo sobre un conjunto de fotogramas de baja resolución consecutivos, los cuales son combinados más tarde produciendo un fotograma de alta resolución. El mayor inconveniente de esta estrategia es que la correlación temporal entre las sucesivas imágenes de alta resolución reconstruidas queda prácticamente sin explotar. Además, la ocupación de memoria asociada con esta estrategia es alta, debido a que todos los datos requeridos del contexto tienen que estar disponibles para su lectura de la memoria. Normalmente se requiere que el procesamiento SR de secuencias de vídeo sea dinámico, es decir, que produzca una secuencia de salida que (al menos) mantenga el número y la tasa de fotogramas de la secuencia de entrada. Debe notarse que la SR dinámica de secuencias de vídeo tiene una demanda significante mayor en términos de rendimiento, velocidad de 265 Appendix B. Resumen en castellano High resolution sequence Low resolution sequence Fig. B.2 Estrategia de ventana deslizante de fotogramas para la SR de secuencias de vídeo ejecutándose en un contexto multi-fotograma. Basado en [TM11]. ejecución y memoria ocupada comparada con la SR de imágenes estáticas [Goh13]. En el contexto de la SR dinámica de secuencias de vídeo, a causa de su enorme carga computacional, se usan algoritmos de menor complejidad para poder mejorar la velocidad de ejecución y por tanto facilitar la posibilidad de poder llevar a cabo el procesamiento en el tiempo destinado para ello. Las implementaciones de los algoritmos basados en interpolación, los métodos directos y los métodos IBP (Iterative Back Projection), resultan ser más satisfactorios con respecto a las limitaciones de ejecución en tiempo real. B.3 El algoritmo de proyección sobre la cuadrícula no-uniforme En este trabajo se aborda el desafío de proporcionar la implementación de un algoritmo de SR capaz de procesar vídeo y trabajando en tiempo real, usando un algoritmo de proyección sobre una cuadrícula no-uniforme , tal y como el que ha sido propuesto en [MC03]. B.3.1 Flujo de ejecución del algoritmo NUGP NUGPA es un algoritmo de fusión que utiliza los datos no redundantes encontrados en un conjunto de imágenes afectadas por deformación asociada con el movimiento y por el aliasing causado por la limitación en banda de los sensores. Usando la clasificación mencionada, NUGPA encaja con la fusión basada en la interpolación de la familia de métodos directos que aplican la transformación en el dominio espacial, usando para ello múltiples 266 Appendix B. Resumen en castellano imágenes. La principal ventaja de esta clase de algoritmos de SR (la baja carga computacional y que hacen posible crear aplicaciones en tiempo real) tiene como contrapartida un modelo de degradación limitado. NUGPA lleva a cabo la reconstrucción de la imagen de súper-resolución en tres etapas: (1) la estimación de movimiento, (2) la reconstrucción basada en la fusión y (3) la interpolación no-uniforme. B.3.1.1 Estimación de movimiento En la primera etapa se estima la función de deformación y las métricas que la caracterizan. Con el objetivo de registrar las imágenes y encontrar las regiones con mayor probabilidad de contener información adicional, el algoritmo de súper-resolución NUGP utiliza una versión de la estimación de movimiento basada en coincidencias de bloques. Durante el registro de la imagen un conjunto de píxeles que forma un bloque (en adelante, un macro-bloque (MB)) del fotograma que está siendo procesado se compara con los bloques de otros fotogramas de la SFW. Para cada MB sólo los píxeles que pertenecen a un conjunto limitado (llamado área de búsqueda, SA) y delimitado a cierta proximidad espacial (llamado radio del área de búsqueda, SAR; expresado en número de píxeles) participa en el proceso de creación de un conjunto de candidatos que serán usados en el proceso de registro de la imagen. El proceso de estimación de movimiento calcula los llamados híper-datos que, en nuestro caso constan del vector de movimiento (MV) y de los índices de similitud asociados basado en la suma de diferencias absolutas SAD (Sum of Absolute Differences). Estos vectores de movimiento son vectores que identifican al MB cuyos índices de similitud tienen el valor más próximo a los que se buscan. En la B.3 se muestra un ejemplo de coincidencia de bloques para un SFW que contiene dos fotogramas de referencia, donde el radio del área de búsqueda es igual al ancho del MB (MBwidth), y contiene 2×MBwidth +1)2candidatos (para mayor claridad sólo se muestran nueve candidatos). B.3.1.2 Reconstrucción basada en la fusión Los híper-datos generados durante la etapa de estimación de movimiento se introducen en el núcleo de SR. EL proceso de reconstrucción se inicia con la creación de la cuadrícula HR. Las dimensiones de la cuadrícula se determinan en base a la precisión usada en el registro de imágenes (en lo sucesivo precisionme). En primer lugar, la cuadrícula de HR se rellena con los píxeles del fotograma de LR que está siendo procesado. Las nuevas coordenadas espaciales de HR se calculan multiplicando las coordenadas LR por la precisión de la estimación de movimiento. 267 Appendix B. Resumen en castellano SAD SAD Sliding Frame Window Search Area Motion Vector Patch being processed Best Match Candidate Fig. B.3 Estimación de movimiento basada en la coincidencia de bloques utilizado por NUGPA usando el criterio de similitud basado en la suma de diferencias absolutas SAD (Sum of Absolute Differences). Tras colocar todos los píxeles LR en la cuadrícula de HR, se consideran los píxeles que provienen de un conjunto de los fotogramas contenidos en el SFW. En este punto, la mayoría de las coordenadas de la cuadrícula de HR no contienen datos válidos. Esas coordenadas se conocen como huecos (holes). Los valores de los huecos están creados en base de la fusión de los datos extraídos de los fotogramas de la ventana de trabajo. Durante la fusión, los híper-datos recibidos (es decir, MVs y SADs) determinan qué datos contribuirán a la estimación del valor de los huecos y con qué peso. Se permite más de un valor para tomar parte en el proceso de la formación de un nuevo valor súper-resuelto del hueco. B.3.1.3 Interpolación no-uniforme En la mayoría de los casos, no todos los huecos de la cuadrícula de HR se rellenan durante la fusión. A los huecos que no han obtenido un valor durante el proceso de fusión se les 268 Appendix B. Resumen en castellano Decimation SRiuma Benchmark Image Quality Assessment Reference Input Output CIF YUV 4:2:0 Decimated YUV 4:2:0 Original Fig. B.4 Flujo de datos del experimento utilizado para evaluar la calidad de la salida de SR. asigna un valor estimado mediante interpolación. Por último, la cuadrícula de HR postinterpolación se ajusta a las dimensiones esperadas especificadas por el factor de escala (en lo sucesivo scale). La escala determina la relación entre la dimensión de la imagen de salida súper-resuelta y la dimensión de la imagen de entrada de baja resolución. Para valores de escala más pequeños que el valor de la precisión (ME) debe reducirse el muestreo de la cuadrícula de HR. Cuando los valores de los parámetros anteriormente mencionados son iguales, el ajuste no es necesario y la cuadrícula post-interpolación se convierte directamente en la imagen súper-resuelta producida por NUGPA. B.3.2 Evaluación cuantitativa de la calidad de la imagen súper-resuelta del algoritmo NUGP Para evaluar el rendimiento en términos de la calidad de la salida percibida, se aplica el NUGPA a los componentes de luma (componente Y del formato YUV) y se utiliza el flujo que se muestra en la Fig. B.4. La secuencia de referencia (y por lo tanto también la de salida) será una secuencia en formato YUV 4:2:0p de 8 bits y tamaño CIF (288x352 píxeles). El uso de este formato supone la entrada de LR correspondiente al formato QCIF (144x176 píxeles) o un formato 72x88 píxeles para los valores del factor de escala de la súper-resolución (scale) de 2 y 4, respectivamente. El caso común es utilizar imágenes de todavía menor resolución (32x32, 64x64, 128x128) [ESHEK12]. Las secuencias utilizadas en los experimentos se describen brevemente en la Tabla B.1. Estas secuencias forman parte de la ‘Xiph.org video Test Media - derf’s collection’ disponible on-line [Xip15]. En la Tabla B.2 se muestra la media de valores de PSNR de 90 fotogramas iniciales 269 Appendix B. Resumen en castellano TABLE B.1 Secuencias de test utilizadas en los experimentos. Sequence Visualization Description Foreman Presenta un obrero realizando rápidos movimientos de la cabeza al azar mientras habla a la cámara. La secuencia cambia a presentar el sitio de trabajo. Esta secuencia presenta movimiento local y global rápido y un cambio de contexto acentuado. Mobile Presenta un tren eléctrico en movimiento en la mesa, un sistema planetario en forma de móvil y un calendario colgado en la pared. El tren se mueve horizontalmente, mientras que el calendario se mueve de arriba abajo presentando un movimiento vertical. Una bola que gira presenta una mezcla de ambos tipos de movimiento. Existe un movimiento global ya que la cámara se mueve hacia la izquierda. Esta secuencia contiene gran cantidad de textura con muchas detalles presentes. Paris Presenta una pareja hablando sentada junto a una mesa. Mientras hablan, la mujer está tirando una pequeña pelota al aire, y el hombre está haciendo mover un lápiz entre sus manos. La secuencia presenta movimientos locales rápidos. El movimiento global está ausente. TABLE B.2 PSNR observado de la secuencia de salida (media sobre los 90 fotogramas iniciales; CIF; componente luma; MBwidth = 4, RF=4 (2+2), SAR = 16, y scale=2). Implementation Media de PSNR [dB] SRiuma 22.38 29.42 22.79 IBP con pesos (20 iteraciones) [BB08] 20.46 27.80 22.22 Nearest-neighbor interpolation 18.78 26.68 19.89 Bilinear interpolation 20.46 28.07 21.23 Bi-cubic interpolation 20.29 28.10 20.92 mobile foreman paris 270 Appendix B. Resumen en castellano de las secuencias paris, foreman y mobile (componentes de luma) mejoradas usando varios métodos de SR y de interpolación. La Tabla B.2 proporciona una comparación con el estado del arte de las implementaciones hardware [BB08], demostrando que NUGPA es capaz de proporcionar imágenes de muy alta calidad. En nuestros experimentos solo se ha comparado la calidad del componente luma (Y) el cual se ha llevado a cabo el proceso de súper-resolución. Los resultados presentados en [WBS05] indican que la incorporación de los componentes de color no altera el rendimiento del modelo de evaluación de una forma significativa. B.3.3 Flujo de ejecución con granularidad fina Los principales factores que limitan la implementación hardware son la naturaleza iterativa de los procesos y el consumo de recursos. La selección del algoritmo NUGP soluciona la primera limitación, permitiendo una implementación hardware eficiente. Desafortunadamente, la ocupación de memoria de la versión no iterativa del algoritmo mostrado en el capítulo anterior hace prohibitiva su implementación en una FPGA que no emplee memoria externa. La referencia software usada originalmente (software SRiuma) trabajaba a nivel de fotograma y por este motivo sus requerimientos del almacenamiento en memoria son muy elevados, resultando en que la implementación hardware usando solo bloques de memoria RAM (Block RAM, BRAM) sea inviable. Una rápida evaluación de la referencia software y del flujo de datos del algoritmo identifica los buffers de fotograma completo usados por el núcleo de fusión como los principales consumidores de memoria. El flujo de ejecución implementado requiere que al menos dos de los fotogramas que están siendo procesados estén disponibles para ser leídos en la memoria durante el tiempo en que se lleva a cabo el proceso de fusión. Por lo tanto, para eliminar los buffers de fotogramas completos el flujo de ejecución necesita ser modificado. En concreto, la granularidad del flujo de datos debe ser ‘refinada’ para reducir los requerimientos de almacenamiento en memoria del algoritmo. Para afrontar este problema hemos propuesto un flujo que internamente aplique un proceso completo de SR a un pequeño conjunto de píxeles (llamado macro-bloque, MB) y que por lo tanto represente un punto medio entre los flujos de proceso con granularidad “gruesa” y “fina”. En comparación con el procesamiento a nivel de fotograma, el cual requiere que la representación de HR de todos los fotogramas de la SFW estén disponibles en la memoria de salida durante todo el proceso de súper-resolución, la estrategia a nivel de macro-bloque requiere que solo las áreas de búsqueda (SA), en lugar de todo el fotograma, 271 Appendix B. Resumen en castellano estén disponible en la memoria. Como resultado, se reducen los requerimientos de almacenamiento de memoria con el inconveniente de incrementar el número de accesos a memoria. El procesado de un bloque de píxeles en lugar de un solo píxel permite optimizar la lógica de síntesis, consiguiendo así que ocupe un menor espacio y permitiendo la implementación del algoritmo completo en el dispositivo FPGA seleccionado. El flujo de ejecución propuesto para el NUGPA es similar al flujo a nivel de fotograma. La principal diferencia entre los flujos son: (i) los tamaños de los elementos básicos guardados en las memorias locales, (ii) la sustitución de los fotogramas de referencia por áreas de búsqueda, (iii) el hecho de que la transición de un paso del algoritmo al siguiente se produce cuando un macro-bloque, y no un fotograma, es procesado, y, (iv) la necesidad de tener en cuenta dependencias entre macro-bloques, lo cual no ocurre cuando el flujo de ejecución se lleva a cabo a nivel de fotograma, y que complica significativamente el proceso de la localización y rellenado de los huecos (holes filling). La elección del flujo de ejecución ha supuesto un impacto importante en las características de implementación, como la ocupación de la memoria y la lógica, el rendimiento y la escalabilidad del sistema, entre las principales. Estos cambios resultan en diferentes estrategias de almacenamiento en memorias intermedia, incrementan la complejidad del manejo de la memoria y mitigan el incremento en el tráfico de la memoria. B.3.4 Evaluación del flujo de ejecución a nivel de macro bloque Los factores de reducción en ocupación de la memoria y del incremento en el tráfico de las comunicaciones asociado con el cambio de flujo realizado desde el nivel de fotograma al nivel de macro-bloque, han sido cuantitativamente evaluados. También han sido cuantitativamente evaluados el factor de reducción de ocupación de la memoria (Total Memory storage Requirements, TMR) y del incremento del tráfico de las comunicaciones (Total Memory Access Count, TMAC), asociados con el cambio del flujo de nivel de fotograma a nivel de macro bloque. B.3.4.1 Estudio de la metodología empleada El estudio de la reducción en la ocupación de memoria debido al cambio de flujo desde el nivel de fotograma al flujo a nivel de macro-bloque se ha llevado a cabo tal y como se explica a continuación. En primer lugar, se ha identificado la memoria empleada en cada etapa del algoritmo, y se ha clasificado dentro de cada uno de los tipos de memoria definidos. A continuación, para cada uno de los tipos de memoria definidos, se ha evaluado el tamaño de 272 Appendix B. Resumen en castellano la memoria y las dependencias con los parámetros del algoritmo. Para evaluar el impacto del cambio de flujo desde nivel de fotograma al flujo a nivel de macro-bloque, en el tráfico de la memoria se han introducido, en primer lugar, figuras para modelar los accesos condicionales a la memoria. A continuación, se han identificado los patrones de acceso a memoria de cada etapa del algoritmo, así como sus dependencias con respecto a los parámetros del algoritmo. Esto ha llevado a la creación de ecuaciones generales para modelar el número de accesos a memoria realizado en cada etapa del algoritmo. Basándonos en las ecuaciones definidas, la ocupación de memoria y el número de accesos a memoria de cada etapa del algoritmo han sido cuantificados, para las dos estrategias de flujo de ejecución (flujo a nivel de fotograma y flujo a nivel de macro-bloque) para las 96 combinaciones de los parámetros del algoritmo NUGP. El conjunto de las 96 combinaciones empleadas para llevar a cabo las pruebas ha sido construido combinando los distintos valores de tamaño de macro-bloque (MBwidth = 4, 8 y 16), radio del área de búsqueda (SAR = 2, 4, 8 y 16), número de fotogramas de referencia soportados (RF = 2, 4, 8 y 16), y el aumento del factor de escalado espacial (scale = 2 y 4). B.3.4.2 Evaluación cuantitativa de la reducción de la memoria y del aumento del tráfico Los resultados de este estudio han demostrado que, para el formato de fotograma QCIF, el cambio desde el nivel de fotograma al esquema de ejecución a nivel de MB, puede conducir a una reducción de la memoria en un factor de entre 6, 8 y 40, y un incremento en los accesos totales a memoria en un factor de entre 1, 22 y 117, dependiendo de los valores de los parámetros del algoritmo NUGP. Los requisitos de memoria totales mínimos y máximos calculados para el flujo de MBL son de 122 KB y de 1051 KB. Dado que los conjuntos de parámetros del algoritmo que se han analizado y que se han utilizado en ambos estudios han sido los mismos, es posible confrontar los valores de ambos factores y analizar la relación de la reducción de la memoria frente al aumento del tráfico (un factor adecuado de mérito, en adelante TMR/TMACratio) como una función de los valores de los parámetros del algoritmo. El estudio realizado ha demostrado que para 20 de las 96 combinaciones probadas, el factor observado de aumento en el tráfico de memoria ha sido mayor que el factor de reducción observado en la ocupación de la memoria. Esto era esperado, dado que el rango de valores calculados anteriormente (de 1,22 a 117) es mayor que el rango de valores finalmente calculados (de 6,8 a 40). Un análisis de estas 20 combinaciones ha demostrado que en trece de ellas se ha observado un tamaño de MB 4, lo que significa que las ventajas aportadas por la granularidad más fina del flujo del SRIR 273 Appendix B. Resumen en castellano Decimation Compare (luma only) Reference Input Output Sequence SR Kernel iuma (Frame-level) (MB-level) Design under test Decimation SR iuma SR Config SR Kernel Fig. B.7 Flujo de datos del experimento utilizado para verificar el modelo funcional. responsables de la carga/almacenamiento de las imágenes, el diezmado de la imagen y la estimación de movimiento. Este hecho condujo a la configuración utilizada en la etapa de validación que se presenta en la Fig. B.7. El formato de la imagen de referencia (y por lo tanto de salida) es una secuencia YUV 4:2:0p con codificación de 8 bits y tamaño CIF (288x352 píxeles). La resolución de entrada de LR depende del valor del factor de escala de súper-resolución. Una vez el procesamiento se llevó a cabo, las imágenes de salida producidas por la implementación se compararon con las imágenes del software de referencia. Para considerar que la verificación era correcta, se buscaba que los resultados observados fueran idénticos para las luminancias de todos los fotogramas. Sólo después de que las salidas generadas por todas las simulaciones hubieran sido idénticas, el modelo del sistema fue considerado como verificado. B.4.1.3 Validación del modelo ESL Una vez validado el modelo funcional del sistema con el núcleo SR propuesto frente a la referencia, se pasó a la creación de la descripción SystemC del nuevo sistema. Para validar la descripción SystemC se usaron los datos de referencia obtenidos durante la validación a nivel algorítmico. La validación se llevó a cabo siguiendo una configuración de validación genérica del HDL. De acuerdo a esta configuración, a nivel más alto de la organización, el sistema comprende dos módulos: el módulo objeto de la verificación DUT (Desing Under Test) y el módulo del banco de pruebas TB (Test Bench), que encapsula el resto del código utilizado en la validación. Es importante destacar que el banco de pruebas no tiene por qué ser sintetizable en ningún punto del desarrollo. La depuración del diseño está basada en el almacenamiento y revisión de los resultados intermedios de procesamiento. Para ello, se han utilizado interfaces adicionales de salida, a través de las que se trasmiten al TB los datos recogidos en los llamados “puntos de va280 Appendix B. Resumen en castellano lidación”. En la mayoría de los casos, los puntos de validación se han situado en lugares intermedios entre los distintos módulos, recolectando y traspasando datos de unos a otros. Estos puntos se han implementado a través de la replicación de canales de transferencia a las interfaces externas. En algunos casos, con el fin de seguir el progreso del procesamiento dentro de los módulos de computación y/o poder investigar los valores de variables y señales internas, los puntos de validación se instanciaron como una parte de estos módulos. El receptor de los datos recogidos en los puntos de verificación fue el hilo de monitorización implementado por el módulo del banco de pruebas. Una vez recibida, la información se comprobaba para determinar la correcta ejecución, proporcionando la información por consola y/o almacenando dicha información. La mayor parte del código del testbench ha sido trasladada desde el código C usado para la validación del modelo de algoritmo. Con el fin de disminuir los tiempos de simulación, el proceso de estimación de movimiento no se ha implementado como código, sino que se ha emulado mediante la carga de los híper-datos (métricas ME) a partir de los archivos obtenidos en las simulaciones del algoritmo a nivel funcional. La salida observada para el modelo de algoritmo de la implementación a nivel de MB se considera como la salida esperada y ha sido usada en la verificación de la implementación a nivel de MB de los modelos ESL y RTL. El resto de cambios han sido realizados con el fin de adaptarse al protocolo de comunicaciones a través de las interfaces a nivel de pin o TLM. B.4.1.4 RTL Una vez verificada, la descripción ESL se convirtió en sintetizable y se usó como entrada en la síntesis de alto nivel. La descripción del sistema generada es la descripción RTL en VHDL o EDIF. Con el fin de validar el HDL sintetizado, se llevaron a cabo distintas simulaciones. La configuración de la validación ESL presentada anteriormente fue reutilizada en las simulaciones RTL. Esto fue posible a través de una simulación de lenguaje mixto en la cual las unidades de diseño VHDL se instanciaron en el diseño en SystemC. La instanciación de los módulos HDL en el SystemC se llevó a cabo mediante unos “encapsulados” que comunican las interfaces RTL con las de SystemC. Una vez que la descripción HDL fue compilada y los encapsulados fueron introducidos en el sistema, se usó el entorno de verificación presentado en la Fig. B.8 para validar la exactitud de la descripción HDL incrustada en el código SystemC. 281 Appendix B. Resumen en castellano Output receiver Input loader Execution monitor Design Under Test module Testbench module Top SystemC module Motion estimation metrics Input pels Design configuration Debug dumps Output pels Verification points VHDL VHDL VHDL SystemC SystemC SystemC Fig. B.8 Entorno de verificación RTL basado en encapsulados SystemC de VHDL. B.5 Implementación en FPGA En este capítulo se explica en detalle la implementación hardware. Esta implementación se ha llevado a cabo a través de la metodología establecida en la sección anterior. De acuerdo a esta metodología, se han modelado distintos niveles de abstracción para alcanzar el nivel de descripción usado en la síntesis lógica. El resultado de la síntesis ha permitido determinar el cuello de botella del sistema y ha servido de guía para el proceso de refinamiento de la arquitectura y la organización de los flujos de procesamiento. El proceso para alcanzar prestaciones de tiempo real implica asumir distintos retos de implementación y múltiples iteraciones de refinamiento. B.5.1 Visión general del sistema de súper-resolución El sistema de súper-resolución basado en el algoritmo NUGP, en el nivel más alto de organización, comprende tres bloques: (i) la estimación de movimiento, (ii) el núcleo de súper-resolución SRK (Super Resolution Kernel), y (iii) el código que gestiona los interfaces de entra/salida. Como ya se ha mencionado, la implementación de la estimación de movimiento está fuera del objetivo de este trabajo. Para mejorar el rendimiento de la simulación, la estimación de movimiento se ha emulado mediante la introducción de código adicional que gestiona la carga de los resultados de estimación de movimiento a partir de archivos con los resultados de estimaciones de movimiento. Como ya se ha indicado, la implementación hardware del algoritmo de BM se ha desarrollado en paralelo [HN10]. El 282 Appendix B. Resumen en castellano Ingress (BM) Block matching Adapter (frame-to-MB) Super-resolution kernel (MB-level) Ingress (SRK) Adapter (MB-to-frame) Egress (SRK) Block matching Adapter (frame-to-MB) Super-resolution kernel (MB-level) Adapter (MB-to-frame) Super-resolution coreSuper-resolution systemI/O Fig. B.9 Visión general del sistema de súper-resolución con el núcleo de procesamiento propuesto. módulo SRK incluye la lógica de gestión del NUGPA, esto es, dado el macro-bloque a procesar y los datos asociados (las métricas del ME, las áreas de búsqueda, etc.), genera la representación súper-resuelta del MB. El código original del SRK se transformó con el fin de operar a nivel de MBs, dado que el algoritmo de estimación de movimiento original operaba a nivel de fotograma, es decir, computaba las métricas fotograma tras fotograma. Tras las modificaciones, el código del SRK se volvió incompatible con el resto del sistema. Para hacer frente a esta situación, se desarrollaron dos módulos de adaptación: un primer módulo para recibir y reordenar las salidas de la estimación de movimiento, poniéndolas en el orden adecuado para el procesamiento MB a MB, y un segundo módulo para reconstruir el fotograma súper-resuelto de los MBs generados por el SRK. La implementación de las modificaciones indicadas estructura el código de referencia de súper-resolución como se indica en la Fig. B.9. En referencia a la nomenclatura establecida, el núcleo de súper-resolución y el código de adaptación forman el denominado dispositivo bajo test, mientras que la gestión de las I/O con la emulación del código BM se ha encapsulado formando un prototipo del módulo de banco de pruebas. In reference the established nomenclature, the super resolution kernel and the adapters code formed the design under test, and the I/O management code with BM emulation code was encapsulated forming a prototype of the testbench module. 283 Appendix B. Resumen en castellano B.5.2 Desafíos en la implementación El objetivo de prestaciones para la ejecución ha sido establecido para que el sistema sea capaz de súper-resolver 24 fotogramas por segundo (fps), de tamaño QCIF y con formato YUV 4:2:0 progresivo, que se corresponde con el formato 24p [Ski05]. A la frecuencia objetivo ftargeted se le asignó el valor de 109 MHz. Este valor fue presentado en las publicaciones [TSC+09, TCH+09] para la implementación de un decodificador H.264 baseline, que siguió el flujo de diseño ESL establecido y que iba dirigido a la misma familia de dispositivos FPGA. El número máximo de ciclos limitcycles disponibles para un procesamiento de un MB 4x4 se estimó dividiendo la estimación de la frecuencia objetivo, por el número de MBs (MBsnr) que vayan a procesarse en un segundo, como se muestra en la ecuación (B.1). Basándose en las prestaciones requeridas de 1584 MBs 4x4 por segundo y una estimación de la frecuencia de operación de 109 MHz, el límite de ciclos se estimó en limitcycles ciclos. La idea de procesar un MB abarca los siguientes elementos: el procesamiento interno del módulo, la sincronización y la comunicación entre módulos. limitcycles =ftargeted MBsnr ∗f ps=109∗106 1584∗24=2867 (B.1) Con el fin de cumplir con los objetivos de ejecución establecidos, se han abordado una serie de desafíos que se presentan brevemente a continuación: •Implicaciones a nivel de sistema de los valores de los parámetros del SRK. El valor de los parámetros del SRK tiene un impacto directo en el rendimiento y en la eficiencia de ejecución de la implementación del núcleo de súper-resolución, NUGPA. A fin de tener ese impacto en cuenta, el código del diseño se ha parametrizado usando las macros que representan valores determinables en tiempo de síntesis para la sincronización, la definición de la memoria, la definición de los iteradores de bucle y la replicación del módulo condicional. •Desafíos relacionados con la memoria. El rendimiento alcanzado por una implementación en dispositivos FPGA de un algoritmo de procesamiento de imágenes por ordenador, es más probable que esté limitado por los recursos de memoria disponibles que por los recursos computacionales. –Implementación de memoria usando recursos eficientes. El uso de look-up tables y de memorias distribuidas para la implementación de las matrices es altamente ineficiente. Con el fin de forzar al sintetizador a que infiera correctamente 284 Appendix B. Resumen en castellano el uso de los dispositivos de bloques de RAM (BRAM) disponibles se han utilizado macros específicas del compilador para la replicación de memorias de acceso aleatorio. –Comunicación entre módulos que permita la ejecución de módulos en paralelo. Con el fin de aprovechar el paralelismo ofrecido por el hardware, el sistema tiene que descomponerse en módulos que puedan ejecutarse en paralelo, haciendo uso de un esquema de comunicación basado en memoria compartida con espacios de direcciones separados de escritura y lectura, intercambiados después de los eventos de sincronización. –Garantía del rendimiento de memoria. Con el fin de aumentar el rendimiento de la memoria, la implementación propuesta utiliza una mezcla de replicación de memoria y de accesos a memoria coalescentes. Se logran ganancias adicionales mediante el uso de memorias compartidas con ampliación de doble mapa de memoria e interfaces de comunicación y replicación de memoria. •Dependencia de datos en el rellenado de los huecos. El proceso de interpolación requiere que las coordenadas de una vecindad de píxeles en forma de cuadrado, cuyo valor está siendo calculado, estén disponibles de forma inmediata para su acceso. En nuestra implementación, la carga de trabajo actual se determina en función de la distancia (en coordenadas de alta resolución) de los bordes de los MB a partir de los bordes del fotograma utilizando un conjunto de indicadores. Los píxels que no puedan ser procesados en el momento de la recepción debido a dependencias de datos se almacenan temporalmente hasta que las dependencias de datos se puedan resolver. •Tamaño variable de la zona de búsqueda. Después de cambiar a procesamiento a nivel de MB, el buffer de fotograma se eliminó. Esto dio lugar a la necesidad de extracción de los píxeles del área de búsqueda a partir de fuentes diferentes. En nuestra implementación, los píxeles del área de búsqueda cuyas coordenadas son coordenadas válidas de fotograma se extraen fuera del núcleo súper-resolución y se realimentaron de nuevo como entradas. El SA se reconstruye internamente en base a un conjunto definido de indicadores que señalan la posición del SA dentro del fotograma. •Construcción y gestión de las cuadrículas. Las cuadrículas de alta resolución del sistema se construyen a partir de una imagen de LR por medio de interpolación con el factor de escala igual a la precisión. En nuestra implementación, el problema con las cuadrículas de alta resolución del SA fue resuelto usando un esquema de direc285 Appendix B. Resumen en castellano cionamiento de memoria intermedia que, junto con las coordenadas de HR de la cuadrícula, fue capaz de proporcionar datos correctos operando en la representación de LR. •Determinación del tamaño de la ventana de fotogramas. El concepto de ventana deslizable de fotogramas se utilizó con el fin de proporcionar capacidades de súperresolución dinámica. El número exacto de fotogramas que se incluyen en la SFW varía a lo largo de la ejecución. En nuestra implementación, un número variable de fotogramas en la SFW se gestiona por medio de contadores internos y de señalizadores externos que señalan transiciones de estado de la SFW. •Operación de división. La operación de división no es una construcción sintetizable por las herramientas de síntesis de alto nivel utilizadas. En nuestra implementación, esta operación fue emulada tanto por medio de otras operaciones (es decir, desplazamientos de bits (+ sumas) o multiplicaciones) como implementándola utilizando un módulo hardware personalizado. B.5.3 Arquitectura del núcleo de súper-resolución: arquitectura implementada y organización La arquitectura funcional del núcleo de súper-resolución que alcanza las prestaciones y la funcionalidad perseguidas se presenta en la Fig. B.10. Internamente, el núcleo mantiene la división en tres partes principales: (i) la gestión del buffer de reordenación, (ii) el núcleo de súper-resolución, y (iii) la reconstrucción del fotograma. Sin embargo, con el fin de cumplir con el número de ciclos de ejecución establecido, muchos de los bloques funcionales tuvieron que ser descompuestos en módulos más pequeños. La tarea de gestión de reordenación del búfer (Reorder Buffer management) tuvo que ser dividida en dos etapas, encapsuladas como dos módulos: el ReorderBufferReap y el ReorderBufferSew. El primero de los módulos opera a nivel de fotograma, almacenando el resultado del ME en la memoria siguiendo un orden adecuado para procesar a nivel de MB. Por lo tanto, el módulo ReorderBufferReap controla la recepción de las métricas de estimación de movimiento, el reordenamiento de datos y el acceso de almacenamiento al buffer de reordenamiento. Sólo una vez que se han recibido todos los datos de los fotogramas de referencia del SFW actual, el ReorderBufferSew puede procesarlos. El ReorderBufferSew gestiona las solicitudes de datos que surgen dentro del núcleo de súper-resolución, lleva a cabo la carga de datos del buffer de reordenamiento y proporciona datos leídos al ShadPrep del SRK. 286 Appendix B. Resumen en castellano BM Output Reorder Buffer Reap Reorder buffer Reorder Buffer Sew Reorder Buffer Reap Reorder buffer Reorder Buffer Sew Shad Prep UpGrid MBHR UpHoles SALR Shad Step MBHR UpHoles SALR Shad Step MBHR UpHoles SALR Shad Step Input (SRK) MBHR Reap Steps MBp HR pel sums MBw HR weights sums Normalize MB Divider /8 / 13* MBHR Holes Filling Prep Intbuffer Pel Map Pel Map Pel Map Pel Map Mean Nearest Mean Nearest Divider Divider MB∗ HR MB∗ HR MB Collector Frame buffer MB Collector Frame buffer MB∗ HR MB∗ HR / 8 / 8 / 32 / 32 / 32 / 32 / 32 / 32 / 32 / 32 / 8 / 8 / 8 / 8 /8 /8 / 13* / 13* Output (SRK) / 32 / 8 / 8 / 8 / 8 / 8 / 32 / 8 / 8 / / 8 / 8 / 8 / 8 / 8 / 8 / 8 / 8 / 8 / 8 / 8 / 8 / 8 / 8 /8 /8 / 13* /8 / 13* / 8 / 8 / 32 /32 Grids construction and fusion Reorder Buffer management Frame reconstruction Holes filling Super-resolution core 1 . . . RF I/O Fig. B.10 Visión general de la organización del modelo pin-accurate. El núcleo de súper-resolución lleva a cabo dos tareas principales: (i) construcción de la cuadrícula y fusión, y (ii) interpolación. De este modo, los módulos que componen el SRK pueden verse como pertenecientes a uno de los dos grupos, dependiendo de sus tareas. Los módulos UpGrid,ShadPrep,UpHoles yShadStep, llevan a cabo la construcción de las cuadrículas. La fusión de las cuadrículas es realizada por los módulos ReapSteps yNormalizeMB. El proceso de rellenado de huecos se lleva a cabo por los módulos HolesFillingPrep yMeanNearest. La entrada al SRK comprende los píxeles LR del macro-bloque y de sus áreas de búsqueda (uno por cada fotograma en el SFW), y un conjunto de indicadores (flags). Estos indicadores describen el macro-bloque y el área de búsqueda dentro del fotograma e identifican las dependencias de datos entre macro-bloques. El módulo ShadPrep recibe los indicadores y solicita resultados de ME (MVs y SADs) del módulo ReorderBufferSew. Una vez que estos parámetros están disponibles, se calcula la localización actual del macrobloque dentro de la zona de búsqueda y los pesos para cada área de búsqueda de los actuales SFW. 287 Appendix B. Resumen en castellano Las coordenadas de localización calculadas y los pesos, junto con los MVs y los indicadores, se transmiten a los módulos ShadStep. Los píxeles del macro-bloque de LR se proporcionan desde la entrada principal de súper-resolución al módulo UpGrid. Cada píxel recibido se somete a la transformación zeroes2ones y se almacena en la memoria de entrada del módulo ReapSteps. Como resultado, la memoria de entrada de ReapSteps mantiene la representación de HR (interpolados 4x) de los MBs de LR. Después de haber terminado su tarea, el módulo UpGrid envía una señal de sincronización al módulo ReapSteps. Los píxeles del área de búsqueda son recibidos por los módulos UpHoles. La tarea de cada uno de los módulos de UpHoles es almacenar los píxeles recibidos en la memoria de entrada del módulo ShadStep correspondiente. Los módulos ShadSteps examinan los indicadores, identifican (en base a los MVs y a las coordenadas recibidas del ShadPrep) y extraen los píxeles que se utilizan para la creación de la cuadrícula del SA. La cuadrícula del SA construida se almacena en la memoria de entrada del módulo ReapSteps. Después de haber terminado su tarea, cada ShadStep envía una señal de sincronización al módulo ReapSteps Habiendo recibido las señales de sincronización de los ShadSteps y el módulo UpGrid, el ReapSteps inicia su funcionamiento. La tarea del módulo ReapSteps consiste en la carga de datos de las cuadrículas, la fusión de las mismas, y el almacenamiento del resultado en la memoria de entrada del módulo NormalizeMB. En primer lugar, recibe los indicadores y los pesos de uno de los módulos ShadSteps. Luego, para cada coordenada de HR se cargan píxeles de las memorias de entrada (la cuadrícula de HR de MB y la cuadrícula del SA) y las fusiona. Los resultados de la operación de fusión para cada coordenada de HR son los valores de la suma ponderada de los píxeles con valor no-nulo cargados desde las memorias de entrada, y la suma de los pesos de los píxeles no-nulos que contribuyeron a la operación de fusión. Dichos valores se almacenan en la memoria de entrada del módulo NormalizeMB. Para cada coordenada de HR, el módulo NormalizeMB carga los dos valores producidos por los módulos ReapSteps. Entonces, se divide la suma ponderada por la suma de los pesos y se almacena el resultado en la memoria de entrada del módulo HolesFillingPrep. Los módulos MeanNearest yHolesFillingPrep llevan a cabo el proceso de rellenado de huecos. Estos módulos calculan los valores de píxel para las coordenadas de HR de la cuadrícula, cuyos valores no fueron estimados durante el proceso de fusión. El resultado producido por las partes de construcción de la cuadrícula y fusión del SRK se transmite al módulo HolesFillingPrep. Los datos se utilizan para construir una estructura llamada el mapa de píxeles. Esta estructura se compone de píxeles del macro-bloque que se está 288 Appendix B. Resumen en castellano procesando actualmente y de los MBs procesados previamente. Los elementos de la imagen que se utilizan para ver el mapa de píxeles de la construcción se determinan en función de la ubicación del macro-bloque dentro del fotograma. Las dependencias de datos entre macro-bloques se gestionan mediante la exclusión del procesamiento de algunas partes del mapa de píxeles. Los píxeles que se van a procesar por el módulo MeanNearest forman una región continua llamada la región de trabajo. La región de trabajo se identifica por sus coordenadas ubicadas en la esquina superior izquierda e inferior derecha. Sólo los píxeles pertenecientes a la región de trabajo son procesados por los módulos MeanNearest. Los píxeles con dependencias entre macro-bloque no resueltas son excluidos de la región de trabajo. Una copia de estos píxeles se almacena en las memorias locales hasta que se resuelven las dependencias de datos. Una vez creado, el mapa de píxeles se almacena en la memoria de entrada del módulo MeanNearest. Sobre la base de las coordenadas recibidas, cada módulo MeanNearest determina su región de trabajo y la procesa píxel a píxel. Para cada píxel cargado se calculan las coordenadas en las que se almacenará en la memoria de entrada del siguiente módulo. Si el valor del píxel es diferente de cero, se almacena en la memoria de entrada del módulo MBcollector, y se carga el siguiente píxel. En caso contrario, antes de almacenar un nuevo píxel, el valor de píxel se calcula usando una interpolación del vecino más próximo (mean nearest neighbors interpolation). El módulo MBcollector recibe los píxeles súper-resueltos y lleva a cabo la reconstrucción del fotograma. Debido a las dependencias entre macro-bloques, el número de píxeles súper-resueltos en las memorias de entrada del módulo MBcollector es variable. El número exacto de píxeles súper-resueltos que se tienen que extraer, su ubicación en la memoria de entrada y la memoria buffer del fotograma son determinados en base a los indicadores recibidos y del valor del parámetro del factor de escalado. Habiendo recibido y procesado todos los MBs que pertenecen a un fotograma, el módulo envía el contenido de las memorias de fotograma a través del puerto de píxeles súper-resueltos. B.5.4 Resultados La Tabla B.3 resume los resultados de la implementación hardware, y muestra una comparación con la implementación del estado de la técnica presentada en [BB08]. El sistema ha sido implementado en un dispositivo FPGA avanzado[Xil09]. La implementación resultante (etiquetada como Tech1en la Tabla B.3). para tamaño de MB de 4x4 píxeles, factor de escala igual a 2, y SFW con 3 fotogramas (RF = 2) ocupó 10291 LUTs, 16 bloques DSP, 289 References [BAA05] D. Barreto, L. D. Alvarez, and J. Abad. Motion estimation techniques in super-resolution image reconstruction. A performance evaluation. In Virtual Observatory: Plate Content Digitization, Archive Mining Image Sequence Processing, pages 254–268. Heron Press, 2005. [BAHH92] James R. Bergen, P. Anandan, Keith J. Hanna, and Rajesh Hingorani. Hierarchical model-based motion estimation. In The Second European Conference on Computer Vision, Proc. of, European Conference on Computer Vision (ECCV) ’92, pages 237–252, 1992. [BB95] S. S. Beauchemin and J. L. Barron. The computation of optical flow. ACM Comput. Surv., 27(3):433–466, September 1995. [BB08] O. Bowen and C.-S. Bouganis. Real-time image super resolution using an FPGA. In Field Programmable Logic and Applications, 2008. FPL 2008. International Conference on, pages 89–94, September 2008. [BB09] Wilhelm Burger and Mark J. Burge. Principles of Digital Image Processing: Core Algorithms (Undergraduate Topics in Computer Science). Springer, 2009 edition, March 2009. [BBM03] Christopher M. Bishop, Andrew Blake, and Bhaskara Marthi. Superresolution enhancement of video. In Artificial Intelligence and Statistics, In Proc. of, 2003. [BBV10] S.V. Basavaraja, A.S. Bopardikar, and S. Velusamy. Detail warping based video super-resolution using image guides. In 17th IEEE International Conference on Image Processing (ICIP), pages 2009–2012, September 2010. [BEZN04] M. Ben-Ezra, A. Zomet, and S.K. Nayar. Jitter camera: high resolution video from a low resolution detector. In IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR), Proc. of, volume 2, pages II–135–II–142, June 2004. [BEZN05] M. Ben-Ezra, A. Zomet, and S.K. Nayar. Video super-resolution using controlled subpixel detector shifts. Pattern Analysis and Machine Intelligence, IEEE Transactions on, 27(6):977–987, June 2005. [BK00a] S. Baker and T. Kanade. Hallucinating faces. In Automatic Face and Gesture Recognition, 2000. Proceedings. Fourth IEEE International Conference on, pages 83–88, 2000. [BK00b] S. Baker and T. Kanade. Limits on super-resolution and how to break them. In Computer Vision and Pattern Recognition, 2000. Proceedings. IEEE Conference on, volume 2, pages 372–379, 2000. [BM10] Brian Bailey and Grant Martin. ESL Models and their Application: Electronic System Level Design and Verification in Practice (Embedded Systems). Springer, December 2010. [BM11] R. Sudheer Babu and K. E. Sreenivasa Murthy. A survey on the methods of super-resolution image reconstruction. International Journal of Computer Applications, 15(2):1–6, February 2011. [Boc09] A. Bock. Video Compression Systems: From First Principles to Concatenated Codecs (Iet Telecommuncations). The Institution of Engineering and Technology, July 2009. 296 References [Bro92] Lisa Gottesfeld Brown. A survey of image registration techniques. ACM Comput. Surv., 24(4):325–376, December 1992. [BS72] Daniel I. Barnea and H.F. Silverman. A class of algorithms for fast digital image registration. Computers, IEEE Transactions on, C-21(2):179–186, February 1972. [BS99] S. Borman and R.L. Stevenson. Simultaneous multi-frame MAP superresolution video enhancement using spatio-temporal priors. In IEEE International Conference on Image Processing (ICIP), Proceedings of, volume 3, pages 469–473, 1999. [BS00] S. Bhattacharjee and Malur K. Sundareshan. Modeling and extrapolation of prior scene information for set theoretic restoration and super-resolution of diffraction-limited images. In Image Processing, 2000. Proceedings. 2000 International Conference on, volume 2, pages 347–350, September 2000. [Cad10] Cadence. Cadence C-to-Silicon Compiler User Guide. Cadence, September 2010. [CB06] G.H. Costa and J.C.M. Bermudez. Statistical analysis of the lms algorithm applied to super-resolution video reconstruction. In Acoustics, Speech and Signal Processing, 2006. ICASSP 2006 Proceedings. 2006 IEEE International Conference on, volume 3, page III, May 2006. [CB07] G.H. Costa and J.C.M. Bermudez. Statistical analysis of the lms algorithm applied to super-resolution image reconstruction. Signal Processing, IEEE Transactions on, 55(5):2084–2095, May 2007. [CB08] G.H. Costa and J.C.M. Bermudez. Informed choice of the lms parameters in super-resolution video reconstruction applications. Signal Processing, IEEE Transactions on, 56(2):555–564, February 2008. [CCL10] Ming-Hui Cheng, Hsuan-Ying Chen, and Jin-Jang Leou. Video superresolution reconstruction using a mobile search strategy and adaptive patch size. In Signal Processing and Multimedia Applications (SIGMAP), Proceedings of the 2010 International Conference on, pages 108–111, July 2010. [CD00] B. Cohen and I. Dinstein. Polyphase back-projection filtering for image resolution enhancement. Vision, Image and Signal Processing, IEE Proceedings of, 147(4):318–322, August 2000. [Cel06] Celoxica. Agility Compiler manual, 2006. [CEN+13a] Pedro P. Carballo, Omar Espino, Romén Neris, Pedro Hernández-Fernández, Tomasz M. Szydzik, and Antonio Núñez. Implementation of scalable video coding deblocking filter from high-level SystemC description. In Proc. SPIE, volume 8764, pages 876408–1–876408–10, 2013. [CEN+13b] P.P. Carballo, O. Espino, R. Neris, P. Hernandez-Fernandez, T.M. Szydzik, and A. Nunez. Scalable video coding deblocking filter FPGA and asic implementation using high-level synthesis methodology. In Digital System Design (DSD), 2013 Euromicro Conference on, pages 415–422, September 2013. 297 References [CG03] Lukai Cai and Daniel Gajski. Transaction level modeling: an overview. In The 1st IEEE/ACM/IFIP international conference on Hardware/software codesign and system synthesis (CODES+ISSS), Proc. of, pages 19–24, 2003. [CJ06] Woong Il Choi and Byeungwoo Jeon. Hierarchical motion search for h.264 variable-block-size motion compensation. Optical Engineering, 45(1):017002–017002–9, 2006. [CKK+96] Peter Cheeseman, Bob Kanefsky, Richard Kraft, John Stutz, and Robin Hanson. Super-resolved surface reconstruction from multiple images. In Maximum Entropy and Bayesian Methods, volume 62 of Fundamental Theories of Physics, pages 293–308. Springer Netherlands, 1996. [CLL+06] Gustavo M. Callico, Rafael Peset Llopis, Sebastian Lopez, Jose Fco. Lopez, Antonio Nunez, Ramanathan Sethuraman, and Roberto Sarmiento. Lowcost super-resolution algorithms implementation over a hw/sw video compression platform, user’s guide and reference manual. EURASIP Journal on Applied Signal Processing, 2006:1–29, 2006. [CLS+08] G. Callico, S. Lopez, O. Sosa, J.F. Lopez, and R. Sarmiento. Analysis of fast block matching motion estimation algorithms for video super-resolution systems. Consumer Electronics, IEEE Transactions on, 54(3):1430–1438, August 2008. [CLT+08] G.M. Callico, S. Lopez, K. Tarajano, J. Lopez, and R. Sarmiento. Impact of fast motion estimation algorithms on super-resolved video sequences. Consumer Electronics, 2008. ICCE 2008. Digest of Technical Papers. International Conference on, pages P3–19, January 2008. [CMMC08] Jérôme Cornet, Florence Maraninchi, and Laurent Maillet-Contoz. A method for the efficient development of timed and untimed transaction-level models of systems-on-chip. In The conference on Design, Automation and Test in Europe (DATE), Proc. of, pages 9–14, 2008. [CN07] Gustavo M. Callico and Antonio Nunez. Mobile receiver design increases the resolution of compressed video. In SPIE newsroom Electronic Imaging & Signal Processing. SPIE, 2007. [CNYA12] Jin Chen, J. Nunez-Yanez, and A. Achim. Video super-resolution using generalized gaussian markov random fields. Signal Processing Letters, IEEE, 19(2):63–66, February 2012. [CTH03] G. Caner, A.M. Tekalp, and W. Heinzelman. Super resolution recovery for multi-camera surveillance imaging. In Multimedia and Expo, 2003. ICME ’03. Proceedings. 2003 International Conference on, volume 1, pages I– 109–12, July 2003. [CY13] Yuzhang Chen and Kecheng Yang. MAP-regularized robust reconstruction for underwater imaging detection. Optik - International Journal for Light and Electron Optics, 124(20):4514–4518, 2013. [CYX+11] Yuzhang Chen, Bofei Yang, Min Xia, Wei Li, Kecheng Yang, and Hiaohui Zhang. Model-based super-resolution reconstruction techniques for underwater imaging. In Photonics and Optoelectronics Meetings (POEM) 2011: Optoelectronic Sensing and Imaging, Proc. of SPIE, volume 8332, pages 83320G–83320G–10, November 2011. 298 References [CZ01] D. Capel and A. Zisserman. Super-resolution from multiple views using learnt image models. In IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR), Proc. of, volume 2, pages II–627– II–634, 2001. [DBS06] Jean-Pierre Deschamps, Gery Jean Antoine Bioul, and Gustavo D. Sutter. Synthesis of arithmetic circuits : FPGA, ASIC and embedded systems. John Wiley & Sons, Inc., Hoboken, New Jersey, March 2006. [DH97] Frank H. Durgin and Alexander C. Huk. Texture density aftereffects in the perception of artificial and natural textures. Vision Research, 37(23):3273 – 3282, 1997. [DHWG07] Shengyang Dai, Mei Han, Ying Wu, and Yihong Gong. Bilateral backprojection for single image super resolution. In Multimedia and Expo, 2007 IEEE International Conference on, pages 1039–1042, July 2007. [DKA04] G. Dedeoglu, T. Kanade, and J. August. High-zoom video hallucination by exploiting spatio-temporal regularities. In IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR), Proc. ofn, volume 2, pages II–151–II–158, June 2004. [Dur96] FrankH. Durgin. Visual aftereffect of texture density contingent on color of frame. Perception & Psychophysics, 58(2):207–223, 1996. [Eco12] G. Economakos. ESL as a gateway from OpenCL to FPGAs: Basic ideas and methodology evaluation. In Informatics (PCI), 2012 16th Panhellenic Conference on, pages 80–85, October 2012. [EF95] A.M. Eskicioglu and P.S. Fisher. Image quality measures and their performance. Communications, IEEE Transactions on, 43(12):2959–2965, December 1995. [EF97] M. Elad and A. Feuer. Restoration of a single superresolution image from several blurred, noisy, and undersampled measured images. Image Processing, IEEE Transactions on, 6(12):1646–1658, December 1997. [EF99a] M. Elad and A. Feuer. Super-resolution reconstruction of continuous image sequences. In IEEE International Conference on Image Processing (ICIP), Proceedings of, volume 3, pages 459–463, 1999. [EF99b] M. Elad and A. Feuer. Super-resolution reconstruction of image sequences. Pattern Analysis and Machine Intelligence, IEEE Transactions on, 21(9):817–834, September 1999. [EF99c] M. Elad and A. Feuer. Superresolution restoration of an image sequence: adaptive filtering approach. Image Processing, IEEE Transactions on, 8(3):387–395, March 1999. [EHO01] M. Elad and Y. Hel-Or. A fast super-resolution reconstruction algorithm for pure translational motion and common space-invariant blur. Image Processing, IEEE Transactions on, 10(8):1187–1193, August 2001. [ESHEK12] Fathi E. Abd El-Samie, Mohiy M. Hadhoud, and Said E. El-Khamy. Image Super-Resolution and Applications. CRC Press, December 2012. 299 References [Fat07] Raanan Fattal. Image upsampling via imposed edge statistics. ACM Trans. Graph., 26(3), July 2007. [FD87] DavidJ. Finlay and PeterC. Dodwell. Speed of apparent motion and the wagon-wheel effect. Perception & Psychophysics, 41(1):29–34, 1987. [FDC84] David Finlay, Peter Dodwell, and Terry Caelli. The waggon-wheel effect. Perception, 13(3):237, 1984. [FEM06] Sina Farsiu, Michael Elad, and Peyman Milanfar. Video-to-video dynamic super-resolution for grayscale and color sequences. EURASIP J. Appl. Signal Process., 2006:232–232, January 2006. [FM06] Michael Farsiu, S.and Elad and Peyman Milanfar. A practical approach to super-resolution. In Visual Communications and Image Processing, Proc. of SPIE, pages 24–38, 2006. [Fos10] Luca Fossati. TLM 2.0 standard into action: Designing efficient processor simulators, 2010. [Fos11] Luca Fossati. Electronic system level design at esa a hardware perspective. Online, September 2011. [FREM03] S. Farsiu, D. Robinson, M. Elad, and P Milanfar. Robust shift and add approach to super-resolution. In Conference on Applications of Digital Signal and Image Processing, Proceedings of SPIE, pages 121–130, 2003. [FREM04a] S. Farsiu, D. Robinson, M. Elad, and P Milanfar. Fast and robust multi-frame super-resolution. In Proceedings SPIE Conference on Image Reconstruction from Incomplete Data, number 13(10), pages 1327–1344, 2004. [FREM04b] Sina Farsiu, Dirk M. Robinson, Michael Elad, and Peyman Milanfar. Dynamic demosaicing and color superresolution of video sequences. In Conference on Image Reconstruction from Incomplete Data III, Proc. of SPIE, volume 5562, pages 169–178, October 2004. [Gau07] Matthew Gaubatz. Metrix mux visual quality assessment package. Online, 2007. http://foulard.ece.cornell.edu/gaubatz/metrix_mux. [GBI09] D. Glasner, S. Bagon, and M. Irani. Super-resolution from a single image. In Computer Vision, 2009 IEEE 12th International Conference on, pages 349–356, September 2009. [GDAG09] Giaime Ginesu, Tiziana Dessi, Luigi Atzori, and Daniele D. Giusto. Superresolution reconstruction of video sequences based on back-projection and motion estimation. In Proceedings of the 5th International ICST Mobile Multimedia Communications Conference, Mobimedia ’09, pages 25:1–25:7, ICST, Brussels, Belgium, Belgium, 2009. ICST (Institute for Computer Sciences, Social-Informatics and Telecommunications Engineering). [Ger74] R.W. Gerchberg. Super-resolution through error energy reduction. J. Mod. Opt., (21(9)):709–720, 1974. [Get11] Pascal Getreuer. Linear Methods for Image Interpolation. Image Processing On Line, 1, 2011. Online. 300 References [GH97] J.J. Green and B.R. Hunt. Super-resolution in a synthetic aperture imaging system. In International Conference on Image Processing, Proceedings of, volume 1, pages 865–868, October 1997. [Gha90] M. Ghanbari. The cross-search algorithm for motion estimation [image coding]. Communications, IEEE Transactions on, 38(7):950–953, July 1990. [Ghe06] Frank Ghenassia. Transaction-Level Modeling with SystemC: TLM Concepts and Applications for Embedded Systems. Springer-Verlag New York, Inc., Secaucus, NJ, USA, 2006. [Gib50] James J. Gibson. The Perception of the Visual World. Houghton Mifflin, 1950. [GJZGSZ00] D. D. Gajski, Domer Jianwen Zhu, R. Gerstlauer, and A. Shuqing Zhao. SPECC: Specification Language and Methodology. Springer, first edition, 2000. [Goh13] S. Gohshi. Limitation of super resolution image reconstruction for video. In Computational Intelligence, Communication Systems and Networks (CICSyN), 2013 Fifth International Conference on, pages 217–221, June 2013. [Goh14] Seiichi Gohshi. 4k-to-8k tv up-converter with super resolution. In Annual Technical Conference Exhibition, SMPTE 2014, pages 1–6, October 2014. [Gro02] Thorsten Grotker. System Design with SystemC. Kluwer Academic Publishers, Norwell, MA, USA, 2002. [Gro09] OSCI TLM Working Group. OSCI TLM-2.0 language reference manual. Open SystemC Initiative (OSCI), July 2009. [GSH+11] T. Goto, S. Suzuki, S. Hirano, M. Sakurai, and T.Q. Nguyen. Fast and high quality learning-based super-resolution utilizing tv regularization method. In 18th IEEE International Conference on Image Processing (ICIP), pages 1185–1188, September 2011. [GSP86] A. Goshtasby, George C. Stockman, and Carl V. Page. A region-based approach to digital image registration with subpixel accuracy. Geoscience and Remote Sensing, IEEE Transactions on, GE-24(3):390–399, May 1986. [Gut15] Eduardo Quevedo Gutierrez. Contribuciones al proceso de SuperResolutcion mediante tecnicas de filtros selectivos, topologia de MacroBloques y sistemas Multi-Camara. PhD thesis, University of Las Palmas de Gran Canaria, April 2015. [GW08] R.C. Gonzalez and R.E. Woods. Digital Image Processing, pages 65–68. Prentice Hall, Inc., 3rd edition, 2008. [Has14a] Hasselblad. Hasselblad h5d 200c. Online, 2014. http://static.hasselblad. com/2014/11/H5D-200c-MS_Datasheet_EN_v3.pdf. [Has14b] Hasselblad. Hasselblad multi-shot cameras. Online, 2014. http://www. hasselblad.com/medium-format/h5d-multi-shot. [HBA97] R.C. Hardie, K.J. Barnard, and E.E. Armstrong. Joint MAP registration and high-resolution image estimation using a sequence of undersampled images. Image Processing, IEEE Transactions on, 6(12):1621–1633, December 1997. 301 [Document text truncated for crawler view.]