MATRIX STRUCTURES AND PARALLEL ALGORITHMS FOR IMAGE SUPERRESOLUTION RECONSTRUCTION QIANG ZHANG * , RICHARD T. GUY † , AND ROBERT J. PLEMMONS ‡ Abstract. Computational resolution enhancement (superresolution) is generally regarded as a memory intensive process due to the large matrix-vector calculations involved. In this paper, a detailed study of the structure of the n 2 × n 2 superresolution matrix is used to decompose the matrix into 9 matrices of size l 2 × l 2 where l is the upsampling factor. As a result, previously large martix vector products can be broken into many small, parallelizable products. An algorithm is presented that utilizes the structural results to perform superresolution on compact, highly parallel architectures such as Field-Programmable Gate Arrays. Key words. image superresolution, FPGA, parallel computation, structured matrices AMS subject classifications. 65R32, 65F10, 65F50, 94A08 1. Introduction. Computational methods for resolution improvement (super- resolution) have attracted much attention lately due in part to their ability to over- come the optical limitations of inexpensive, lower resolution sensors. See, for instance, [6, 14, 16]. Superresolution (SR) is based on the idea that slight variations in the in- formation encoded in a series of low resolution (LR) images can be used to recover a high resolution (HR) image. The basic superresolution problem can be posed as an inverse problem [1, 6], min f ||DH i S i f − g i || 2 2 ,i =1,...,l 2 , (1.1) where f is the vectorized true high resolution image, g i is a vectorized lower resolution image, D is the decimation matrix, H i is a blurring matrix, S i is a shift matrix and l is the upsampling factor. In the models that follow, the decimation matrix D is a local averaging matrix that aggregates values of non-intersecting small neighborhoods of HR pixels to produce LR pixel values. The shift matrix S i , also called the interpolation matrix, assigns weights according to a bilinear interpolation of HR pixel values to perform a rigid translation of the original image. The blurring matrix H i is generated from a point spread function (PSF) and represents distortion from atmospheric and other sources. As it will be better explained in Section 2, usually the l 2 matrices DH i S i are stacked to create one large least squares problem min f ||Af − g|| 2 2 , (1.2) where, using the MATLAB notation, A =[DH 1 S 1 ; ... ; DH l 2 S l 2 ], g =[g 1 ; ... ; g l 2 ], and A ∈ R n 2 ×n 2 , being n × n the dimension of the true high resolution image f . The dimensionality of the problem is usually quite large. Given a moderate HR image size of 256 × 256 with l = 4, the na¨ ıve way to construct A would require 2l 2 matrices H i and S i , for i =1,...,l 2 , each 65536 × 65536, plus one smaller matrix D that is 4096 × 65536, assuming l = 4. The system matrix A is sparse but is of size 65536 × 65536. * Department of Biostatistical Sciences, Wake Forest University Health Sciences, Medical Center Boulevard, Winston-Salem, NC 27157 (qizhang@wfubmc.edu). † Department of Mathematics, Wake Forest University, Winston-Salem, NC 27106 ‡ Departments of Mathematics and Computer Science, Wake Forest University, Winston-Salem, NC 27106. 1