Compressive Imaging Using RIP-Compliant CMOS Imager Architecture and Landweber Reconstruction

IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS FOR VIDEO TECHNOLOGY, VOL. 30, NO. 2, FEBRUARY 2020 387 Compressive Imaging Using RIP-Compliant CMOS Imager Architecture and Landweber Reconstruction Marco Trevisi , Ali Akbari , Maria Trocan, Senior Member, IEEE, Ángel Rodríguez-Vázquez, Fellow, IEEE, and Ricardo Carmona-Galán ,Senior Member, IEEE Abstract—In this paper, we present a new image sensor architecture for fast and accurate compressive sensing (CS) of natural images. Measurement matrices usually employed in CS CMOS image sensors are recursive pseudo-random binary matrices. We have proved that the restricted isometry property of these matrices is limited by a low sparsity constant. The quality of these matrices is also affected by the non-idealities of pseudo-random number generators (PRNG). To overcome these limitations, we propose a hardware-friendly pseudo-random ternary measurement matrix generated on-chip by means of class III elementary cellular automata (ECA). These ECA present a chaotic behavior that emulates random CS measurement matrices better than other PRNG. We have combined this new architecture with a block-based CS smoothed-projected Landweber reconstruction algorithm. By means of single value decomposition, we have adapted this algorithm to perform fast and precise reconstruction while operating with binary and ternary matrices. Simulations are provided to qualify the approach. Index Terms—CMOS image sensor architecture, compressive sensing, Landweber reconstruction, power spectral density, random binary matrix RIP proof, single value decomposition, ternary measurement matrix. I. INTRODUCTION COMPRESSIVE Sensing (CS) is a data acquisition technique that can be used to represent the content of an image with fewer samples than required by Shannon-Nyquist theorem. These samples are called compressed samples. Manuscript received June 18, 2018; revised September 21, 2018 and December 4, 2018; accepted December 25, 2018. Date of publication January 11, 2019; date of current version February 5, 2020. This work was supported in part by the Spanish Government (European Region Development Fund, ERDF/FEDER) under Grant TEC2015-66878-C3-1-R MINECO, in part by the Junta de Andalucía under Grant TIC 2338-2013 CEICE, in part by the Office of Naval Research, USA, under Grant N000141410355, and in part by the EU H2020 MSCA ACHIEVE-ITN under Grant 765866. This paper was recommended by Associate Editor X. Li. (Marco Trevisi and Ali Akbari contributed equally to this work.) (Corresponding author: Ricardo Carmona-Galan.) M. Trevisi, Á. Rodríguez-Vázquez, and R. Carmona-Galán are with the Instituto de Microelectrónica de Sevilla (IMSE-CNM), CSIC-Universidad de Sevilla, Seville 41092, Spain (e-mail: [email protected]; angel@ imse-cnm.csic.es; [email protected]). A. Akbari is with the Centre for Vision, Speech and Signal Processing, University of Surrey, Guildford GU2 7XH, U.K. (e-mail: [email protected]). M. Trocan is with the Department of Signal, Images and Telecommunications, Institut Supérieur d’Électronique de Paris, Paris 92130, France (e-mail: [email protected]). Color versions of one or more of the figures in this article are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TCSVT.2019.2892178 Reconstruction algorithms exploit the sparsity of the original image in order to recover it by finding the unique solution to an underdetermined linear system [1]. This system can be expressed in matrix notation as: y=x(1) where y∈RMis the set of Mcompressed samples, x∈RN represents the Npixels values folded in a column vector, i.e. the unknowns of the system of equations, and ∈ RM×N(MN)is referred to as measurement matrix. It has been shown in [2] that, if holds the Restricted Isometry Property (RIP), then the solution to this inverse linear problem, by means of convex optimization, is unique and error free [3]. CS is a field of research that is located in-between the areas of analog microelectronics and of signal processing. The link between these areas is given by the feasible yet RIP-compliant design of measurement matrices. Measurement matrices affect both the performance of reconstruction algorithms as well as the architecture of CS CMOS image sensors (CS-CIS). CS differs from standard acquisition-plus-compression techniques in that the image is not captured and then digitized before compression. To benefit from the characteristics of CS, the image needs to be sensed and compressed simultaneously. For that reason, in CS-CIS, compressed samples are formed in the analog domain and digital conversion occurs just before their delivery. Since each pixel has different coefficients for each compressed sample, if the coefficients of a measurement matrix ()were real numbers, a CS-CIS would require the implementation of in-pixel analog multipliers. This solution would be impractical in terms of pixel sensitivity and spatial resolution. Furthermore, each compressed sample is a linear combination of the weighted readings of all pixels (N), which means it is the sum of Nproducts. The upper bound on the amount of bits required to represent theses linear combinations is: BCS =log2N+(B+BI)(2) where ·denotes the smallest integer greater than the argument and BIand Bare the number of bits used to describe the pixel values and the coefficients of the measurement matrix, respectively. It is virtually impossible to design analog to digital converters (ADCs) with such resolutions in standard technologies. In fact, one of the main drawbacks of most CS-CIS 1051-8215 © 2019 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See https://www.ieee.org/publications/rights/index.html for more information. Authorized licensed use limited to: Universidad de Sevilla. Downloaded on March 09,2020 at 14:01:04 UTC from IEEE Xplore. Restrictions apply. 388 IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS FOR VIDEO TECHNOLOGY, VOL. 30, NO. 2, FEBRUARY 2020 implementations [4]–[6] is the lack of dynamic range to properly represent compressed samples. The use of binary measurement matrices and block based compressed sampling (BCS) [7], which divides a pixel array into smaller sub-arrays that are digitized independently, is essential for the implementation of practical measurement matrices in CS-CISs. In signal processing literature, there are many works that study deterministic binary measurement matrices; i.e. cyclic matrices [8], sparse bipartite graphs [9] or other forms of deterministic matrices that try to match the performance of Gaussian random matrices [10]. These matrices are usually coined hardware-friendly by their authors, but, even if they are binary, given the technology at our disposal, they cannot be implemented in application specific integrated circuits (ASICs). As compressive sampling an acquisition technique, the fact that we cannot apply these theories to actively design sensors is, in itself, the limitation of these works. This happens because, even if these matrices are binary in nature, they still disregards the fact that, on-chip, they need hardware resources to be generated, transmitted and/or stored. To store a binary measurement matrix in its integrity a large memory would be required. This exceeds nowadays CS-CIS capabilities in term of area occupation. If, on the contrary, we were to introduce a connection between sensor and reconstruction system to continuously send the rows of the matrix, the amount of exchanged data would void CS benefits as a compression technique. For this reason, measurement matrices designed in analog microelectronics [5], [6] solely focus on the use of pseudo random binary numbers generators (PRNG) to recursively create the rows of the matrix sample after sample. This solution employs the minimum possible amount of on-chip resources and does not need data feedback between sensor and reconstruction system. These matrices though do not ascertain the quality of the resulting samples because their primary focus is placed almost entirely on electrical parameters such as power consumption and compactness of the implementation. PRNG are non-linear spatially discrete and temporally discrete dynamic systems made of binary logic elements that show high sensitivity to initial conditions and evolve in time according to a divergent and fractal behavior [11]. A new state in their time evolution is derived from their actual state using a feedback mechanism that promotes this instability. Each state corresponds to a row of the binary measurement matrix. Each row of a measurement matrix is generated by the CS-CIS from the previous row and there is no need to store the whole binary matrix on chip or to receive it from outside the chip. In this paper, the performance of binary measurement matrices obtained using binary PRNG is evaluated in terms of RIP. We introduce a new method for the generation of a hardware-friendly measurement matrix that takes into account both technology limits and quality of the resulting matrix. A differential pixel readout system and an elementary cellular automaton (ECA) are proposed to recursively create a ternary measurement matrix in a row by row fashion. The performance of the generated ternary measurement matrix is analyzed and compared with other binary measurement matrices typically found in analog microelectronics. To do so we introduce a new reconstruction algorithm that uses this matrix to deliver reconstruction results close to those of the more refined matrices typically found in signal processing. II. RIP OF BINARY MATRICES PRODUCED BY PSEUDO RANDOM BINARY GENERATORS PRNG commonly used to implement measurement matrices into CS-CIS are Linear Feedback Shift Registers (LFSRs) [12], [13]. LFSRs consist of a series of sequentially connected flip-flops. Some of the outputs of these flip-flops, besides being connected to the input of the flip-flop that they precede, are also connected, by means of XOR logic gates, to the input of the very first one. When considered individually, the output of each flip-flop in this sequence evolves in time with a behavior that resembles a symmetric probability distribution. Using one flip-flop for each pixel it is possible to sequentially create the pseudo-random binary coefficients of the measurement matrix one row after the other. The resulting matrix is an approximation of a very efficient type of measurement matrix popular in CS: a matrix in which each element is derived from the normalized Gaussian distribution. Early implementations of CS-CIS binary measurement matrices were carried out incorporating a flip-flop inside each pixel. These flip-flops were then connected sequentially to form LFSRs. This solution though posed a heavy burden on the design specifications of image sensors due to the large area required at pixel level. Recent solutions have overcome this problem placing LFSRs around the pixel array as row and column drivers (Fig. 1) [6], [14], [15]. Inside of each pixel, the flip-flops are replaced by simple and more compact selectors. These selectors are driven by LFSRs surrounding the array. Selectors can be implemented with few transistors, using for instance switches or gated inverters. This simplification was possible because, at pixel level, we can consider random row and column selection as the multiplication of two independent binary random variables. The multiplication of random variables still delivers random outputs. Using PRNG, such as LFSRs, the recursive row and column drivers also reduce the amount of on-chip memory needed to store their elements to Ng=2√Nbits. Since this recursive pseudo-random sequence can be reproduced from an initial Ngbits seed, it allows the independent generation of the same measurement matrix in the system in charge of reconstruction, without continuous transmission of the measurement matrix from the sensor. It is important to notice that LFSRs produce binary matrices that are neither normalized nor orthogonalized. To understand the quality of the compressed samples that they produce, it is necessary to analyze them in terms of RIP. Suppose that an area of size √N×√Npixels is sampled by a measurement matrix generated row by row as described above. Let us consider the simplest possible case in which an element of the pseudo random generator is assigned to each pixel of the imager so that each row of the measurement matrix corresponds to a time step in the discrete evolution of the generator itself. Mathematically, this sampling process Authorized licensed use limited to: Universidad de Sevilla. Downloaded on March 09,2020 at 14:01:04 UTC from IEEE Xplore. Restrictions apply. TREVISI et al.: COMPRESSIVE IMAGING USING RIP-COMPLIANT CMOS IMAGER ARCHITECTURE AND LANDWEBER RECONSTRUCTION 389 Fig. 1. Conceptual floor plan of a CS-CIS. is formulated as Eq.(1), where x∈RNis the √N×√N array of pixels folded in vector form. The reconstruction of the original signal xfrom the compressed samples yis an ill-posed problem. By making use of sparsity, a well-known characteristic existing in natural signals, a k-sparse1signal, like images and videos, can be reconstructed from a few samples obtained using an appropriate measurement matrix. A sufficient condition for the unique and exact recovery of the signal is the RIP of the measurement matrix [2]. Amatrix∈RM×Nsatisfies the RIP of order kif there is constant δk(0<δ k<1)such that, for all vectors x∈RN with x0≤k(i.e. k-sparse signals2), it holds: (1−δk)x2 2≤x2 2≤(1+δk)x2 2(3) The smallest non-negative value for δkis called restricted isometry constant (RIC) of order k. The construction of a measurement matrix which satisfies the RIP is a central problem in CS. If sparse signals with maximum possible sparsity level kcan be recovered exactly and stably, it is said the measurement matrix has sparsity order k. An upper bound for the sparsity level is k≤CM/log(N/M)where C is a constant [16]. The signals with kabove this bound can only be approximately reconstructed. In practice, to recover a signal xwith a large k, is desirable to have a measurement matrix with a δkas small as possible. The limit imposed on the RIP of random binary measurement matrices has already been studied in [17] Such binary matrices satisfy the RIP with: δk=(3k−2)m (k−2)m+4dif 3 ≤d≤M/2and2≤m≤d(4) where dis the estimated number of non-zero elements in each column of and mis the maximum inner product 1xis called k-sparse or a signal with sparsity level k,ifithasnomore than knon-zero components with kN. 2x0norm denotes the number of non-zero components in x. between two distinct columns. The distribution of 1’s among the elements of random binary matrices, generated by pseudo random generators, can be associated with a Bernoulli probability distribution in which each element has a probability P of being 1 and ¯ P=1−P of being 0. Since these pseudo random measurement matrices try to emulate the normalized Gaussian distribution, which is a symmetric probability distribution, usually, circuit designers choose LFSR because their distribution is akin to P =¯ P=0.5. Let us consider the ideal case of an ideal pseudo random generator capable of recursively producing independent binary elements with a symmetric Bernoulli probability distribution having P =0.5. If we consider a large enough pseudo random binary generator (a large enough number of pixels N),since each of its elements, by construction, is independent from the others and the probability distribution is symmetrical, we will obtain a measurement matrix with an equal amount of 0s and of 1s row-wise. If the acquisition process is long enough (a large enough number of acquired compressed samples M), we will have an equal amount of 0s and of 1s columnwise as well. Such measurement matrix, whose entries are randomly drawn from pseudo random generators, satisfies the RIP with Eq.(4). Following [17], let us define the mutual coherence of as: μ()=md(5) We can derive mand das functions of the number of pixels in the image N, the sampling subrate of the sensor S=MN, and the probability of each element of being different from 0, P. From the definition of subrate (S), it is straight forward that M=S×N. Moreover, using the probability that each element of the matrix has of being different from 0, we can define the number of non-zero elements in a column of as: d=S·N·P(6) Mutual coherence μ()of a matrix represents the maximum absolute value of the cross-correlations among its normalized columns and is defined as: μ()=maxj=kϕH jϕk(7) where ϕkis the k-th normalized column of matrix and ϕH j is the conjugate transpose of the j-th normalized column of the same matrix. For a matrix to be column-wise normalized it means that for each column: ϕH jϕj=1. Given the definition of mutual coherence and the randomness of ϕH jand ϕk, to approximate m, we will make use of the most probable outcome of a dot product between normalized columns of : m=EϕH jϕkd(8) Remember that ϕH jϕkis the sum of Melement by element products. Each one of these Mproducts will be different from zero if and only if the two elements that are being multiplied differ from 0 as well. By construction, we are multiplying independent random variables; the probability that their product be different from 0 is equal to the joint Authorized licensed use limited to: Universidad de Sevilla. Downloaded on March 09,2020 at 14:01:04 UTC from IEEE Xplore. Restrictions apply. 390 IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS FOR VIDEO TECHNOLOGY, VOL. 30, NO. 2, FEBRUARY 2020 probability of the single elements. For this reason, the expected amount of non-zeros among the Mproducts is: ETr ϕH j·ϕk=P2·S·N(9) where Tr (·)is the trace of the matrix obtained by the vector multiplication of ϕH jand ϕk. Furthermore we can derive the value of each product in terms of P, Nand Sby using the definition of normalization and remembering that is binary. Normalizing a vector involves dividing each one of its non-zero elements by the Euclidean norm of the vector itself. Since each non-zero element in a binary vector has to take the value of 1, all elements of a normalized binary vector will be equal to the inverse of the vector Euclidean norm itself. For this reason, if the i-th elements in column ϕjand column ϕk are non-zeros, then their normalized product will be: ϕH jiϕki=1 P·N·S(10) Combining Eq.(9) and Eq.(10), we can deduce the most probable cross-correlation between two columns of a binary matrix This correlation only depends on P and it is independent from the number of pixels and the sampling subrate: EϕH jϕk=P2·S·N P·N·S=P (11) Joning Eq.(4), Eq.(5) and Eq. (11) we obtain: δk=(3k−2) (k−2)+4/Pif 3/(S·N)≤P≤1/2and2/(S·N)≤P2≤P (12) For P =0.5andN≥6, it is possible to see that all of the conditions hold and that random binary matrices that follow a symmetric probability distribution will have RIP (δk≤1) only if k≤4. A recovery algorithm can deliver error free reconstructions only if the measurement matrix holds RIP. Eq. (12) poses a harsh limitation on matrices that need to be implemented in CS-CISs. Block based CS-CIS, in many occasions, need small blocks in order to maintain acceptable errors during ADC. ADC dynamic range is determined by noise and the maximum signal range. It is very difficult to allocate large dynamic range. Small blocks reduce the requirement for signal range and consequently making it feasible for implementation. This in turn introduces asymmetries in the pixel array. These asymmetries create artefacts in the captured images in correspondence of the borders of the blocks. Reconstruction of block based compressive sampled images usually entails smoothing of those regions. When using a CS-CIS the smoothing process needs to take care of physical discontinuities as well. Because these two effects combine themselves they can potentially deteriorate the quality of the recovered images. It would be desirable to implement a measurement matrix that did not add errors during reconstruction to facilitate the smoothing process. But, since the RIP of random binary matrices holds only if k≤4, errors during reconstruction are to be expected. III. TERNARY MATRICES GENERATED USING CLASS III ELEMENTARY CELLULAR AUTOMATA We have bound the RIP of binary measurement matrices usually implemented in a CS-CIS with a limit (k≤4)that, in practice, is too small to grant the error free reconstruction of natural images. As explained in the introduction, PRNG generated measurement matrices are necessary to design feasible CS-CIS architectures. However, to overcome the limitation imposed by Eq.(12), we propose to use them to generate ternary measurement matrices instead of the commonly used binary ones. The resulting ternary matrices present smaller coherence than the aforementioned binary ones of same size thus improving RIP and diminishing reconstruction errors. Similar results have been obtained using binary matrices with binary correlations between columns applied to CS in [18]. Ternary measurement matrices can be generated by delivering two independent binary coefficients to each pixel: one to define the number either 0 or 1 and the other to define the sign of the contribution. As mentioned in section II, if one set of row/column selectors (or driving signal) can be used to produce a binary matrix having one bit per element we could use two sets (or driving signals) to produce a ternary matrix that needs two bits per element. If each pixel received two driving signals, it would be possible to use one of them to determine if the pixel took part in a compressed sample and, in case it did, use the other to select the sign of its contribution. Since recent CS-CIS examples [6], [14], [15] have shown that as few as Ng=2√Nbits of information are needed to generate binary measurement matrices on-chip, the amount of resources (number of transistors) needed to incorporate a ternary matrix would still be implementable in a CS-CIS design. Each coefficient of the resulting matrix would have a probability distribution of: ⎧ ⎪ ⎨ ⎪ ⎩ P+1=0.25 P0=0.5 P−1=0.25 (13) where P−1represents the probability that a pixel has of contributing negatively in a compress sample, P+1represents the probability that a pixel has of contributing positively and P0represents the probability that a pixel has of not contributing at all. Furthermore, a ternary measurement matrix would have a positive effect on the limit, imposed by Eq. (2), on the amount of bits required to represent a compressed sample. In hardware, the sign associated to a pixel contribution can be implemented by means of a differential readout system that routes the output of the pixels through one of two output lines [5]. These lines, outside of the pixel array can then be digitized separately or used as input of differential circuits, such as analog subtractors or transimpedance amplifiers. Using a differential readout system to divide the pixel contribution to different output lines to be treated separately would in fact reduce the amount of bits needed for ADC thus relaxing the converter design parameters. Another important factor to study in order to improve the quality of CS-CISs measurement matrices are LFSRs Authorized licensed use limited to: Universidad de Sevilla. Downloaded on March 09,2020 at 14:01:04 UTC from IEEE Xplore. Restrictions apply. TREVISI et al.: COMPRESSIVE IMAGING USING RIP-COMPLIANT CMOS IMAGER ARCHITECTURE AND LANDWEBER RECONSTRUCTION 391 performances. Even if LFSRs are good candidates to emulate random probability distributions thanks to their low impact on on-chip area consumption, ultimately the matrices that they produce still follow a row-by-row deterministic pattern. For this reason it is necessary to analyze this pattern to understand how well the generated measurement matrix approximates the functionality of one that is extracted directly from a real binary probability distribution. Moreover, LFSRs are not the only example of circuits that can serve as PRNG in a CS-CIS. A valid alternative are elementary cellular automata (ECA) [19]. An ECA is a spatially and temporally discrete dynamic system made of identical interconnected cells. Each ECA cell is a triple-input single-output binary unit. Based on this setup there are 256 possible ECA configurations. Each configuration, also known as rule, corresponds to a particular logic implementation. We will focus our attention on rules that belong to class III, in particular rule 30 [20]. We will study rule 30 hardware implementation and the matrices that it produces in order to establish if it is better suited than LFSRs to be used as PRNG in the design of CS-CIS. Let us consider the simple case of binary measurement matrices. This simplification poses no loss of generality because the elements of the aforementioned ternary matrices are generated following the same pattern used for binary measurement matrices and the bits are grouped in sets of two only at pixel level. Let us define G∈{0,1}M×Nas a binary measurement matrix obtained using a PRNG and R∈ {0,1}M×Nas a binary measurement matrix obtained selecting each element at random from a binary probability distribution with P =¯ P=0.5. Gwill be a good approximation of Rif it holds three characteristics: the number of non-zero elements in a row of Gmust approximate P ≈0.5; the temporal evolution of the elements of the PRNG do not present repeating patterns; and the temporal evolution of the elements of the PRNG present no correlation with one another. The first and second requirements bind the number of non-zero elements in Gwith the average of the probability distribution used to create R. The second and third requirements can be used to bind the mutual coherence of Gas: E[μ(G)]≈E[μ(R)]=P (14) We will analyze how binary measurement matrices generated using LFSRs and ECA fulfill these criteria in order to compare the performances of the two types of PRNGs in CS imagers. A. Density of the PRNG Output To study the number of non-zero elements in a row of G let us define the density of a PRNG output (DO(t))as the average of all the states of its binary elements at a given discrete time step t(following a notation similar to [21]): DO (t)=1 NN h=1ϕh(t)where t=1,2,···,M(15) being ϕh(t)the state of the h-th element of the PRNG at time tas well as the t-th element of the h-th column of G, Nthe number of elements of the PRNG as well as the number of columns in Gand Mthe number of discrete time steps of the PRNG as well as the number of rows in G.Lett=1 Fig. 2. Density of the output over a varying random initial seed for rule-30 ECA and LFSR. be the initial condition when the first row of Gcoincide with the seed that has been loaded in the memory on-chip. Since the elements of Rare extracted from a symmetric probability distribution, an optimal seed should have half of its elements set to 0 and half set to 1 and as such DO (1)=0.5. To analyze the performance of LFSRs and ECA it is necessary to evaluate how fast their states can reach DO (t)=0.5in the eventuality of having a suboptimal initial seed in which the number of elements set to 0 differs from the number of those set to 1 (Fig. 2). To do so we devise a MATLAB experiment in which we evolve a 64-cells rule-30 ECA and a 64-flip-flops LFSR using 65 different initial seeds having an density that varies from DO (1)=0 all the way to DO (1)=1. In the result of the experiment we represent DO(t)using greyscale elements in which DO(t)=0 is colored black and DO(t)=1 is colored white. In Fig. 2 the top line represents these initial seeds with increasing density, from the leftmost DO(1)=0 to the rightmost DO (1)=1. Each step underneath shows the temporal evolution of the PRNG output density evolving from each initial seed configuration. Fig. 2 shows that LFSRs density changes in time at a much slower rate than class III ECA. This implies that in the presence of sub-optimal initial seed, the DO of LFSRs would take several time steps longer than ECAs to achieve a stable situation centered around DO (t)=0.5. As reported in [22], the behavior of LFSRs is more similar to class IV ECA, where randomness in the evolution is linked to the randomness of the initial configuration, rather than class III ECA where the chaotic evolution is introduced by the rule itself. For this reason, in order to approximate Rusing LFSRs much care should be placed in the selection of an appropriate initial seed. ECA appears to be a safer choice to generate a Grather than LFSRs at least in terms of number of distribution of non-zero elements starting from seeds that are not ideal. Fig. 2 also shows an interesting behavior for ECA when its seed is either DO (1)=0orDO (1)=1andfor LFSR when its seed is DO (1)=0. These particular seeds represent the forbidden states of these PRNGs. All PRNGs have forbidden states. A forbidden state is a state in which the PRNG continuously delivers a fixed output. These states though do not pose any issues in real-world implementations because it is impossible for a PRNG to stumble upon Authorized licensed use limited to: Universidad de Sevilla. Downloaded on March 09,2020 at 14:01:04 UTC from IEEE Xplore. Restrictions apply. 392 IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS FOR VIDEO TECHNOLOGY, VOL. 30, NO. 2, FEBRUARY 2020 Fig. 3. Power spectral density of rule-30 ECA and LFSR. them accidentally. Even though these states do not pose any problem during standard operation, to use PRNGs in CS-CIS, it is important to be aware of their presence in order to avoid sending them accidentally. The only particular concern regarding forbidden states happens when power is first applied to a circuit. Since each register can randomly start up containing either a zero or a one, a PRNG could power up containing its forbidden state, but this can be quickly taken care of initializing the PRNG with an opportune seed value. B. Power Spectral Density To understand if the rows in Gpresent repeating patterns and to compare dynamical behaviors of LFSRs and rule-30 ECAs we will use Power Spectral Density (PSD) analysis. This technique has been extensively applied to analyze discrete dynamical systems such as ECA [21] or LFSRs [13]. For a spatially discrete and temporally discrete dynamic system, the DFT can be expressed as: φh(f)=1 M N  i=1 ϕh(t)e−i2πtf/Mwhere t=1,2,···,M (16) where, following a notation similar to [21], φh(f)is the DFT value of the h-th element of the PRNG at frequency f,Monce again is the number of discrete time steps of the PRNG as well as the number of rows in Gand ϕh(t)is the state of the h-th element of the PRNG at time tas well as the t-th element of the h-th column of G. PSD expresses the distribution of the energy of a waveform among its different frequency components. Any peak in a graphic of PSD(f)over fat a given frequency fpwould represent a strong repeating pattern of period 1fpamong the rows of GGiven all φh(f)PSD can be computed as: PSD (f)=1 NN e=1|φh(f)|2(17) where Nis the number of elements of the PRNG as well as the number of columns in G. The PSD profile of a suitable and efficient PRNG for CS should closely resemble white noise since its energy should equally spread throughout the entire spectrum. Once again we devise a MATLAB experiment in which we evolve a 64-cells rule-30 ECA and a 64-flip-flop LFSR. Fig. 3 shows the PSD profiles of rule-30 ECA and LFSR. These graphics have been produced for t=1000 time steps starting from an initial random binary seed with DO(1)=0.5. The first thing that we can notice is that both PSD profiles are similar to white noise. That was to be expected because otherwise neither of them would be suitable for CS but, the amplitude of the pattern of the LFSR presents a greater variation than that of rule-30 ECA as well as some unwanted fluctuation. Rule-30 ECA appears to a better choice over LFSRs to extract quasi-independent coefficients from a random like distribution. Furthermore, as we said in the introduction, a necessary requirement to design a measurement matrix for CS is that said matrix respects the RIP which means that it must be nearly orthonormal, at least when operating on sparse vectors [2]. Orthogonal matrices have uncorrelated columns and rows. For this reason, even though RIP and mutual coherence are two different concepts, a matrix G, compliant with RIP, and its transpose T Gwill both have low mutual coherence by definition [23]. Class III ECA PSD profile offers better guaranties over LFSRs that the mutual coherence of the generated matrix will be low. We opted to use a single class III rule elementary cellular automaton over other solutions like multi-rule hybrid cellular automata for two reasons: to facilitate the circuit design and because we are interested in the patterns that the individual cells produce and not the generated number in its decimal representation. C. Hardware Implementation To analyze the correlations among temporal evolutions of the elements of a PRNG is equivalent to analyze mutual coherence of T G. To do so it is necessary to study how these patterns are generated by the feedback mechanism that promotes the PRNG’s instability. Fig. 4 is an example that represents the feedback needed to evolve an 8-bit LFSR. This feedback is imposed on the first flip-flop of the sequence through XOR gates tapping flip-flops 4, 5, 6 and 8. The flip-flops that must be tapped to allow an N-bit LFSR to cycle through all 2N−1 combinations are tied to primitive polynomials with binary coefficients and degree equal to N[24]. On the contrary, the feedbacks needed for the evolution of a cell in a cellular automaton are spatially local. The block diagram of Fig. 5 is an example of a cell of a rule-30 ECA, in it we have defined as S the actual state of the cell. L and R are the actual states of its two closest neighbors, to the left and to the right respectively. NS corresponds to the cell’s next state. This is in sharp contrast with the feedback of LFSRs which is unique for all its elements. The only drawback that ECA present with respect to LFSRs is that while the basic element of a LFSR is a flip-flop, to implement a cell of a rule-30 ECA it takes two extra logic gates. Even though ECA require more resources (transistors) to be implemented than LFSRs it has been shown in [22] that, due to the shifting of data, inside large LFSRs, the outputs of the individual flip-flops have strong correlation between each other. A consequence of this fact is that the mutual coherence of G Authorized licensed use limited to: Universidad de Sevilla. Downloaded on March 09,2020 at 14:01:04 UTC from IEEE Xplore. Restrictions apply. TREVISI et al.: COMPRESSIVE IMAGING USING RIP-COMPLIANT CMOS IMAGER ARCHITECTURE AND LANDWEBER RECONSTRUCTION 393 Fig. 4. 8-bit LFSR. Fig. 5. ECA cell implementation of rule 30. could inadvertently increase, thus reducing k, the maximum amount of non-zero coefficients, in an appropriate domain, that the signal to be sampled can have before incurring in reconstruction errors. IV. RECOVERY ALGORITHM FOR PRNG GENERATED MEASUREMENT MATRICES When the signal to be acquired is a large image, the receiver needs a large memory to store the generated measurement matrix. Furthermore, the reconstruction algorithm undergoes a huge computational burden due to this large size. BCS [7], [25], [26] solves these issues by dividing the scene to be acquired into Hnon-overlapping blocks of size √Nb×√Nb. With BCS, the acquisition is achieved independently, block by block, using the same measurement matrix ∈RMb×Nb. For each block we can reformulate Eq. (1) as: yi=xi(18) where xi∈RNbis the vector representation of the i-th block of the image and yi∈RMbare the corresponding samples. At the receiver side, each block can be reconstructed separately. In order to avoid creating unpleasant blocking artefacts, recently, several BCS reconstruction algorithms have been designed [25], [27]–[29]. BCS with smoothed-projected Landweber reconstruction (BCS-SPL) is a fast and efficient reconstruction algorithm [25]. It is based on successive projections and thresholding in the transform domain and provides a good trade-off between computational complexity and reconstruction quality. While this method is simple to derive, convergence is only guaranteed for a limited class of matrices that hold the strong RIP condition, e.g. Gaussian random measurement matrices with real coefficients. Recently, Akbari and Trocan [26] proposed a way to improve the performance of the BCS-SPL algorithm to a wide range of random measurement matrices, including Bernoulli and Gaussian measurement matrices. Binary PRNG such as those used in hardware implementations deliver matrices that, a priori, would destabilize the BCS-SPL algorithm. In this paper, we propose one way to adapt this algorithm to work with these generators. To this purpose, we combine the BCS-SPL algorithm with singular value decomposition in order to design a more robust and stable reconstruction algorithm, called BCS-SVD. Let ˆx[k] ibe the estimated image at the k-th iteration, partitioned into H blocks. Furthermore, suppose that the SVD decomposition of the measurement matrix is given by: =UDVT(19) where Uand Vare orthogonal matrices of size Mb×Mb and Nb×Nbrespectively, and their column are known as the singular vectors of the matrix .Dis a diagonal matrix of size Mb×Nband its element along its diagonal are known as the singular values [30]. The proposed algorithm updates the i-th block ˆx[k] iiteratively by means of the following approximating operator: ˆx[k+1] i=ˆx[k] i+TUSUTyi−ˆx[k] i(20) where Sis a diagonal matrix of size Mb×Mbwhose entries are obtained by taking the reciprocal of the corresponding squared non-zero elements of the matrix D. The most useful feature of the approximating operator in BCS-SVD is that it partially generalizes matrix inversion of the measurement matrix . This operator is based on the error between the original samples of each block and the corresponding approximated samples, obtained by the estimated block at the previous iteration. Using this residual, it forms a correction equation for which this process is repeated. In the next step, all the reconstructed blocks at iteration k are put back into the full-image plane to update the reconstructed image ˆx[k+1]using the same smoothing and thresholding operations, proposed in the BCS-SPL algorithm [25]. The smoothing operation is done by Wiener filtering which reduces the blocking artefacts. The aim of thresholding process is controlling the local sparsity that is well-known in the natural image and video signals. To this purpose, soft thresholding is done in the transform domain by: ˜x[k+1]=˜x[k+1]˜x[k+1]≥τ 0else (21) where τis a predetermined threshold and ˜x[k+1]is the transformed image. A wide range of sparse transforms can be used, such as discrete cosine and wavelet transforms. The most Authorized licensed use limited to: Universidad de Sevilla. Downloaded on March 09,2020 at 14:01:04 UTC from IEEE Xplore. Restrictions apply. 394 IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS FOR VIDEO TECHNOLOGY, VOL. 30, NO. 2, FEBRUARY 2020 TABLE I BCS-SVD ALGORITHM effective one, as proposed in [25], is the dual-tree discrete wavelet transform (DDWT). Such a procedure is iterated until the difference of successive approximations of the image is a very small value, i.e.  ˆx[k+1]−ˆx[k] ≤ε. The details of the proposed BCS-SVD algorithm are summarized in Table 1.The operator R(·)splits the image into Hnon-overlapping blocks and the operator R−1(·)puts back the reconstructed blocks into the corresponding positions in the reconstructed image, padded with zeros elsewhere. Eq. (20) is derived to solve the following minimization problem argminxi{yi−xi2+λxi1(22) where λis the regularization parameter. Eq. (20) is the enhanced version of the projected Landweber step proposed in BCS-SPL [25] which, with the thresholding operator of Eq. (21), iteratively improves the reconstruction quality. This iterative threshing process has been analyzed extensively in the literature. Different theoretical analyses of the BCS-SPL algorithm applied to the CS recovery problem have been given in [31]–[33]. BCS-SPL is guaranteed to converge to a local minimum of Eq. (22). The argument for BCS-SVD follows the same line as that in [31]. V. EXPERIMENTAL RESULTS In this section, the performance of the CS measurement matrix, described in sections III, and the BCS-SVD algorithm, proposed in the section IV, will be evaluated via a suite of simulations carried out on 8-bit grayscale standard images of size 512×512 pixels. Here, as shown in Fig. 6, we report the results of four classic images often used in image processing: Boat, Lena, Livingroom, and Baboon.3The Peak Signal-to-Noise Ratio (PSNR) is chosen as objective quality measure for our experimental framework. PSNR is an approximation to human perception of reconstruction quality from lossy compression. When used for this purpose the original data is considered the 3The MATLAB implementation of the ternary measurement matrix obtained using ECA and the BCS-SVD recovery algorithm can be found, along with more test at the address: www2.imse-cnm.csic.es/icaveats/tcsvt Fig. 6. Test images: (a)Boat, (b)Lena, (c) Livingroom, (d)Baboon. signal and the error introduced by compression is the noise: PSNR =10log10 max (x)2 1 MN  x−ˆx 2 (23) being max(x)the highest possible pixel value. The simulations were carried out in MATLAB. In each set of simulations, 30 trials are performed due to the random nature of generated measurement matrices and the average PSNR is computed. For the BCS sampling, we use the block dimension of size 32 ×32, i.e. Nb=1024. The value of threshold in the reconstruction algorithm is set to ε=10−4.Weremind the reader that dividing the pixel array into smaller subarrays that can be digitized independently is essential for two reasons. It helps the implementation of practical measurement matrices in CS-CISs and it reduces software resources needed during reconstruction. Usually fast and accurate reconstruction algorithms opt for small blocks, typically 8 ×8. On the other hands, blocks that are too small deteriorate the quality of the sensor introducing asymmetries within the pixel array. We chose a block size of 32×32 as a compromise between the two divergent necessities. We have set the value of threshold in the reconstruction algorithm to be ε=10−4because we saw that it was a good compromise between reconstruction time and accuracy of the reconstruction result. Experiments are designed with two goals: first, the performance of ECA-based ternary matrices for CS should be compared with other types of matrices; second, the stability and robustness of the proposed measurement matrices with different recovery algorithms should be evaluated. This test set up includes 480 trials. More examples can be found in our online repository. A. Validation of the Proposed Measurement Matrices In this subsection, the performance of the proposed ternary and binary pseudo random measurement matrices based on Authorized licensed use limited to: Universidad de Sevilla. Downloaded on March 09,2020 at 14:01:04 UTC from IEEE Xplore. Restrictions apply. TREVISI et al.: COMPRESSIVE IMAGING USING RIP-COMPLIANT CMOS IMAGER ARCHITECTURE AND LANDWEBER RECONSTRUCTION 395 Fig. 7. PSNR of the reconstructed images as a function of subrate (S). ECA is compared to performance of the pseudo random measurement matrices based on LFSR. In addition, we have included the Gaussian random matrices, whose entries are randomly selected from a Gaussian distribution, as a reference. We conducted a series of test during which we applied the CS using four matrices: Gaussian a Gaussian random matrix chosen i.i.d. from N(0,1)and then normalized; B LFSR a binary pseudo random measurement matrix based on LFSR; B ECA a Binary pseudo random measurement matrix based on ECA; and T ECA a ternary pseudo random measurement matrix basedonECA. First, the test images are divided into non-overlapping blocks of size 32 ×32 pixels and each block is compressed with the abovementioned measurement matrices. Then, the test images are recovered using the proposed recovery algorithm in section IV. This process is repeated several times per image and then the results are averaged. Fig. 7 demonstrates the quality of reconstructed images in terms of PSNR as a function of subrate (S), where S varies from 0.1 to 0.5. As shown in Fig. 7, ternary matrices can be used to compressive sample the images and their performances degrade when the subrate decreases. An interesting observation is that the ternary matrix T ECA outperforms binary matrices almost reaching the Gaussian random matrix Gaussian. Fig. 9 represent Lena reconstructed using a ternary measurement matrix generated using ECA with a subrate of 0.1. We will provide more images and results along with the MATLAB code used to obtain them in the online repository. Depending on the subrate, on average, the proposed ternary matrix improves the PSNR of the reconstructed images from 0.08 dB to 0.35 dB when compared with the binary measurement matrix B ECA and from 0.27 dB to 1.90 dB when compared with the binary measurement matrix B ECA. As can be seen, the performance the binary measurement matrix B ECA is lower than that of the ternary matrix T ECA at lower subrates. However, it shows a good performance when higher subrates are considered. For example, at S =0.5, the PSNR of reconstructed images using the samples obtained from the matrix B ECA matrix is almost 0.08 dB lower on average than that of the matrix T ECA. Finally, it should be particularly noted that the ternary matrix T ECA competes with the performance of the Gaussian random matrix, which is an optima theoretical measurement matrix. It should be reminded that the proposed ternary matrix is much easier to implement in hardware and that adds a valuable advantage to it. Authorized licensed use limited to: Universidad de Sevilla. Downloaded on March 09,2020 at 14:01:04 UTC from IEEE Xplore. Restrictions apply.