Analysis and Optimization of Aperture Design in Computational Imaging

ANAL YSIS AND OPTIMIZA TION OF APER TURE DESIGN IN COMPUT A TIONAL IMA GING Adam Y edidia, Christos Thrampoulidis, Gr e gory W ornell Department of Electrical Engineering and Computer Science, MIT ABSTRA CT There is growing interest in the use of coded aperture imag- ing systems for a v ariety of applications. Using an analy- sis frame work based on mutual information, we e xamine the fundamental limits of such systems—and the associated op - timum aperture coding—under simple but meaningful propa- gation and sensor models. Among other results, we sho w that when thermal noise dominates, spectrally-flat masks, which hav e 50% transmissivity , are optimal, but that when shot noise dominates, randomly generated masks with lo wer transmis- sivity offer greater performance. W e also provide compar- isons to classical pinhole cameras. Index T erms — coded aperture cameras, computational photography , optical signal processing 1. INTR ODUCTION Digital signal processing plays an important role in modern imaging systems. Many modern imaging systems operating at optical and higher frequencies use coded apertures, whereby the traditional lens in the aperture is replaced with a spa- tial mask that selectively blocks portions of the light from reaching the sensor . Y et while this is an increasingly impor- tant imaging modality—and one with a long history dating back to the earliest pinhole cameras—typical mask designs are guided by heuristics and/or numerical procedures. As Figure 1 depicts, with an empty aperture, scene reco v- ery from measurements at the imaging plane is very poorly conditioned. Coded-aperture cameras seek to improv e the conditioning of the problem through the use of more compli- cated (and transmissive) masks than a pinhole in combination with suitably designed post-processing. In this paper , we dev elop a comparativ e analysis of these imaging systems, using mutual information as our perfor- mance measure. Moreover , we use far-field geometric optics to model propagation, and our sensor model at the imaging plane includes thermal and shot noise components. Among the earliest and simplest instances of coded- aperture imaging are those based on pinhole structure [ 1 , 2 ] and pinspeck (anti-pinhole) structure [ 3 ], though more com- plex structure is often used. Other methods in v olve cameras that uses a mask in addition to a lens to, e.g., facilitate depth estimation [ 4 ], deblur out-of-focus elements in an image [ 5 ], Fig. 1 : Three imaging systems (left, top-to-bottom): no aper- ture, a pinhole and a lens. Arrows indicate paths light from the scene takes to a particular point on the imaging plane. On the right is an arbitrary mask, an illustration of its discretiza- tion and the corresponding transfer matrix. enable motion deblurring [ 6 ], and/or recov er 4D lightfields [ 7 ]. Some forgo the lens altogether to decrease costs and/or meet physical constraints [ 8 ] [ 9 ]. Certain other systems, intended for non-line-of-sight ap- plications, rely on known structure in between the scene and the imaging plane to impro v e the conditioning of the problem [ 10 ], like windows [ 11 ] or corners of buildings [ 12 ]. These can be vie wed as instances of broader class of coded-aperture systems that we analyze, in which the mask is naturally oc- curring and not chosen. 2. MODEL Scene . Let I ( x ) [W/m] represent the intensity of the scene ov er space in one dimension: 0 ≤ x ≤ L . W e denote J = R I ( x )d x [W] the net power radiated. Assume a uniform dis- cretization of [0 , L ] into n bins of size ∆ = L/n each, and denote x 1 , x 2 , . . . , x n their centers. W e assume that the dis- cretization is fine enough that the intensity at each bin i ∈ [ n ] takes constant value I ( x i ) . Let f i = I ( x i ) · ∆ be the power radiated from each bin. W e model f = [ f 1 , . . . , f n ] as a mul- tiv ariate Gaussian distribution N ( µ 1 , Q ) with mean µ and cov ariance matrix Q . W e set µ = J /n to ensure that the av erage net power is E [ P i ∈ [ n ] f i ] = P i ∈ [ n ] µ = J . The Gaussian statistics model for images is frequently used, such as in [ 13 , 4 ]. In this paper, we consider the follo wing two cases: IID: W e assume that the f i ’ s are uncorrelated, i.e., Q = I . While natural scenes will exhibit correlations, studying the IID case is a means of performing a worst-case analysis. 1 /f -prior: W e follow a classical statistical model according to which the power spectrum of natural images depends as 1 /f over the spatial frequency [ 14 ], by taking Q = F ∗ n D ? F n , where F n is the normalized DFT matrix of size n and D ? is a diagonal matrix with the following entries: d ? i = d ? n/ 2+ i = 1 /i, for i = 1 , . . . b n/ 2 c . Imaging plane . The imaging plane consists of m adja- cent and equally-sized pixels. W e focus on the case where m = n . The power y j [W] measured at each pixel is y j = 1 n P n i =1 A j i · f i , where f i is the power radiated from the i th bin. The (1 /m ) –scaling is chosen to ensure preservation of energy: E [ P j y j ] = 1 m P j P i A j i · E [ x i ] ≤ 1 m · mn · J n = J. The measurement model is a reduction of a more complete forward model, which further accounts for distance atten- uation and cosine factors in light propagation [ 15 ]. This reduction corresponds to a scenario in which the scene is far enough from the imaging plane that the distance attenuation and cosine factors are well-approximated by constants. Apertur e . Denote by A the m × n transfer matrix whose en- tries A j i model the aperture. W e assume that a maximal inte- gration time is allowed, and normalize it so that the maximal value for each entry of A is 1 . W e let ρ denote the transmis- sivity of the aperture. For an on-off aperture, ρ measures the fraction of elements that transmit light (See Fig. 1 ). In gen- eral, we assume a circulant A ; that is equiv alent to assuming that the mask repeats a certain pattern (of length n ) twice: A j i = a ( i − j ) mod n where a T = ( a 0 , . . . , a n − 1 ) T is the first row of A . Noise . W e distinguish between two different types of noise. (Thermal noise): This includes noise sources that are inde- pendent of the contribution to the measurements due to the scene of interest. W e model it as additi v e Gaussian with v ari- ance W /m , i.e., constant net noise po wer W and each pixel absorbs power proportional to its size, giving rise to the 1 /m factor . (Shot noise): This includes measurement noise that depends on the contribution due to the scene of interest. This results in additiv e Gaussian noise of variance ρ · J m (proportional to the net power of light that goes through the aperture). Overall, the measurement at each pixel is modeled as y j = (1 /m ) P i ∈ [ n ] A j i f i + z j , where z j ∼ N (0 , ( W + ρ · J ) /m ) . Mutual information . The mutual information (MI) between the measurements y j , j ∈ [ m ] and the unknowns f i , i ∈ [ n ] of the imaging problem is gi ven as I = log det  1 W + ρ · J · 1 m · A QA T + I  . Recall that a circulant matrix is diagonalized by F n . Also, Q = F ∗ n DF n where D = I (IID scene) or D = D ? (1/f-prior). W ith these, I reduces to (recall m = n ) I = n X i =1 log  1 W + ρ · J · d i · | λ i ( A ) | 2 n + 1  , (1) where, λ i ( A ) denotes the eigenv alue of A corresponding to the i th frequency . W e often write λ i when clear from context. Apertur e T ypes . Here, we summarize sev eral types of aper- ture designs and their corresponding models. Pinhole: W e model a pinhole camera as an on-of f mask with only a single open element, i.e., A = I (or, any permutation of the identity). Also, for a pinhole: ρ = 1 /n . Spectrally-Flat patterns: The family includes pseudo-noise binary (0/1) patterns such as maximum length sequences (MLS) and uniform redundant array patterns such as URA and MURA. Onwards, we refer to patterns with the follo wing properties as spectrally-flat patterns: (i) ρ ≈ 1 / 2 (there is one more one than zero); (ii) they are spectrally flat with the exception of a DC term [ 16 , 17 , 18 , 19 ]. Random on-off patterns: W e study random patterns where each entry of a is generated IID Bern( p ) , for p ∈ (0 , 1] . For such random on-off patterns we use ρ = p , since for large n (which is our focus) the number of on-elements is ≈ np . Random uniform patterns: W e also study patterns consisting of elements that can partially absorb light, e.g., [ 20 , 7 ]. W e focus on random such patterns where each entry of a is IID Uniform ([0 , 1]) . For these patterns, the e xpected transmissi v- ity ρ = 1 / 2 . 3. RESUL TS 3.1. IID scene Throughout this section we study the IID scene model. It is con v enient to work with the normalized mutual information per pixel I := I /n. 3.1.1. Pinhole From ( 1 ) the (normalized) MI of a pinhole is giv en by I pinhole = log  1 n · W + J + 1  . By allo wing only a fraction of 1 /n of the light to go through, the formula justifies that the performance of a pinhole deteriorates drastically for large n (cf., MI goes to zero, unless W becomes negligible, e.g., un- less it scales inv ersely proportionally to 1 /n ). Note that this result applies only to a vanishingly small pinhole (decreas- ing in size as n increases); a pinhole of fixed size achie ves constant mutual information per pixel. 3.1.2. Spectrally-flat patterns The following proposition characterizes the MI of spectrally- flat patterns and shows that they maximize MI when thermal noise is dominant. See Appendix A for a proof sketch. Proposition 3.1. Consider the IID scene model. Let I ? be the mutual information of a spectrally-flat pattern for an odd n . 1 It holds that: lim n →∞ I ? = log  1 / 4 W + J / 2 + 1  . (2) Mor eover , if W  J , then given the mutual information I p of any on-off apertur e design with np “on” elements and p 6 = 1 2 , for lar ge enough n , it holds that I p < I ? . Remark 1 . For spectrally-flat occluders, λ 1 ≈ n 2 and | λ 2 | = . . . = | λ n | ≈ √ n 2 . The contribution of the first eigen value to the sum in ( 1 ) is O (log ( n ) /n ) , which captures the rate at which con vergence in ( 2 ) is true. Throughout, statements that in v olve n → ∞ are to be interpreted with the rest of parame- ters (such as W , J , ρ ) held constant (independent of n ). Remark 2 . The adv antage of spectrally-flat apertures when thermal noise is dominant is shown by concavity of log and applying Jensen’ s inequality to ( 1 ). More generally , Jensen’ s inequality is tight if f λ 2 d 2 = λ 3 d 3 = . . . = λ n d n . This leads to optimality of flat-spectrum patterns for d i = 1 (IID scenes), but the conclusion might be different for correlated scenes. Also, we require that W  J . In the next section, we show that if this is not the case then certain patterns with p < 1 / 2 can outperform the spectrally flat ones. 3.1.3. Random on-off patterns W e explicitly compute the asymptotic value of the MI for ran- dom on-off patterns. Our theoretical results use tools from random matrix theory (RMT) [ 21 , 22 ] and are thus asymptotic in nature. (Howe ver , numerical simulations suggest accuracy of the predictions for n on the order of a fe w hundreds.) A proof sketch is deferred to Appendix A . Proposition 3.2. Assume the IID scene model. Let X be a random variable with density f X ( x ) = | x | e − x 2 . The mu- tual information I p for a random on-off cir culant system with parameter 0 < p < 1 con ver ges in probability with n to: e I p = E X [log( p (1 − p ) W + pJ X 2 + 1)] . Remark 3 . Maximizing the formula of the proposition over p giv es the optimal choice of the transmissi vity parameter . Since log is increasing, it can be sho wn that the maximum occurs at p ? = ( W /J ) · ( p 1 + J /W − 1) . (3) In particular , when ambient noise is dominant ( W  J ), then using q 1 + J W ≈ 1 + J 2 W giv es p ? ≈ 1 2 . On the other hand, when shot noise is dominant ( J  W ), then p ? ≈ q 1 J ; thus, 1 Here, we implicitly assume that n is such that an MLS, or URA, or MURA pattern exists. For example, MURA patterns can be generated for any prime n that is of the form 4 d + 1 , d = 1 , 2 , . . . . Fig. 2 : A plot of the MI per pixel of a spectrally-flat occluder , a random on-off occluder with p = 0 . 5 , and a random on-off occluder with optimally chosen p ? . See Proposition 3.2 . fewer open holes in the aperture design are desirable. See Figure 2 for an illustration. For small v alues of 1 /W (relati ve to 1 /J ): e I p ? ≈ e I 1 2 , but e I p ? > e I 1 2 when 1 /J is small. Remark 4 . For Bern(1 / 2) patterns, an application of Jensen’ s inequality verifies that e I 1 2 < log( 1 / 4 W + J/ 2 E X [ X 2 ] + 1) = I ? , i.e., spectrally-flat patterns are superior . On the other hand, a random pattern Bern( p ? ) with optimal parameter given by ( 3 ) can outperform the spectrally-flat one. For example, this happens when shot noise is dominant, as illustrated in Fig- ure 2 . In the same figure, spectrally-flat patterns are superior when W  J as predicted by Proposition 3.1 . 3.1.4. Random uniform patterns Similar to Proposition 3.2 we le verage results of [ 21 ] to ev alu- ate the MI performance of random uniform patterns; we omit the details due to space limitations. Proposition 3.3. Consider the IID scene model. Let X be a random variable with density f X ( x ) = | x | e − x 2 . The nor- malized mutual information I uniform for a random uniform cir- culant system con ver ges in pr obability with n to: e I uniform = E X [log( 1 / 24 W + J/ 2 X 2 + 1)] . Comparing the formula of the proposition to Proposition 3.2 , reveals that e I uniform < e I p , for all p ∈ [ 1 2 − 1 √ 6 , 1 2 ] . Hence, random on-off masks in this range of p outperform random uniform masks. In short, if physical limitations pre- vent the use of apertures that can redirect light, but can only absorb it, then absorbing all (with appropriate p ) is better than partially (at least for random designs). 3.2. Correlated scene W e extend the “worst-case” analysis of the previous section regarding IID scenes to correlated ones. W e follow the 1/ f scene prior model. Due to space limitations, we restrict the exposition to spectrally-flat and random on-of f patterns. Spectrally-flat patterns: The MI of the spectrally-flat pat- terns for correlated scenes can be computed similar to ( 2 ). For large enough n , we find that I ? ≈ log ( n/ 4 W + J/ 2 + 1) + 2 P n − 1 2 k =2 log( 1 / 4 W + J/ 2 1 k + 1) ≈ log ( 1 / 4 W + J/ 2 · n ) + 1 / 2 W + J/ 2 (log( n/ 2) − 1) , where, for the first approximation: n − 1 n ≈ 1 , and, for the second one: log(1 + x ) ≈ x for | x |  1 and P n k =1 1 n ≈ log n . In contrast to the IID case where the MI scaled linearly with n , here it scales as O (log( n )) . Random on-off patterns: Contrary to the case of IID scenes where knowledge of the the limiting spectral density of A suf- fices to characterize the MI, for correlated scenes each eigen- value is weighted differently . Hence, the behavior of the MI depends on the statistics of each individual eigen v alue. Since A is circulant, the eigenv alues are exactly the Fourier coef- ficients of the entries of the generating vector a , i.e., λ 1 = P n − 1 ` =0 a ` , and, for k = 2 , . . . , n − 1 2 (assume n is odd for sim- plicity): λ 2 k = λ 2 n − k = g 2 k + h 2 k , where g k := P n − 1 ` =0 a ` · cos( `k 2 π n ) , h k := P n − 1 ` =0 a ` · sin( `k 2 π n ) . Next, observ e that if the a i ’ s were standard Gaussians then the following state- ments hold. (a) λ 1 is distributed N (0 , n ) . (b) g k ’ s and h k ’ s are IID N (0 , 1 / 2) ; therefore, λ 2 k iid ∼ 1 2 χ 2 2 where χ 2 2 denotes a chi-squared random variable with two degrees of freedom. This leads to the following conclusion: Lemma 3.1. Let the first r ow of a cir culant A have entries drawn IID standar d Gaussians and the MI be given as in ( 1 ) , for some γ := 1 W + ρ · J and for d i = d ? i . Then, E [ I ] equals E G ∼N (0 , 1) log  γ nG 2 + 1  + 2 P n − 1 2 k =2 E X ∼ χ 2 2 log  γ X 2 i + 1  . W e conjecture that the conclusion of Lemma 3.1 is uni- versal over the distribution of the entries of a T , i.e., it holds for entries that have zero mean, unit variance, and bounded third moment. Based on this assumption, we conjecture that the expected mutual information E [ I p ] for a random on-off circulant system with parameter 0 < p < 1 for the correlated scene model is giv en by: E G ∼N (0 , 1) log  ( p p (1 − p ) · G + p √ n ) 2 W + pJ + 1  + 2 n − 1 2 X k =2 E X ∼ χ 2 2 log  p (1 − p ) X W + pJ 1 2 k + 1  . (4) Figure 3 sho ws a comparison of the formula predicted by ( 4 ) against simulated data. It further re veals that ( 4 ) can be used to numerically ev aluate the optimal p = p ∗ . 4. DISCUSSION AND FUTURE W ORK Our framework allows to rigorously show that spectrally-flat patterns are optimal for IID scenes, and formalize the ar - guably unintuiti ve empirical claim (discussed, for instance, Fig. 3 : Analytical formula follows Eqn. ( 4 ). Simulated data are av erages of 1000 randomly generated apertures of size n = 250 for v arious different v alues p . W e set W = − 20 dB . in [ 4 ]) that the best masks tend to transmit half the light they receiv e. [ 7 ] raises the question of whether continuous-valued masks perform better than binary-valued ones; we plan to use our framew ork to find an answer in the future. In this work, we focused exclusi v ely on 1D masks, which are relev ant for example in de-blurring along one dimension [ 6 ]. W e leave extensions to 2D masks to future work. Howe v er , we men- tion in passing that that much of the analysis conducted here can be directly applied to study separable 2D apertures, i.e. ones that can be expressed as the outer product of two 1D apertures. A. PR OOF SKETCHES Pr oof sketch of Pr oposition 3.1 : For con venience set λ i := λ i ( A ) and γ = 1 W + pJ . W e treat the DC-term of the spectrum, i.e. λ 1 , separately from the rest. Note that A1 = ( np ) 1 ;s thus, λ 1 = np. Next, let us denote I ∼ 1 the MI in ( 1 ), ex- cluding the term that in volv es λ 1 . By concavity of log and Jensen’ s inequality , I ∼ 1 is upper bounded by n − 1 n log  γ ( n − 1) n n X i =2 | λ i ( A ) | 2 + 1  ≈ log  γ p (1 − p ) + 1  , (5) where the bound is tight iff | λ 2 | = | λ 3 | = . . . = | λ n | ; ( 5 ) uses the fact that P n i =2 | λ i ( A ) | 2 = k A k 2 F − λ 2 1 = n 2 p (1 − p ) and n ≈ n − 1 for large n . In particular , spectrally-flat patterns achiev e the upper bound, which gi ves I ? ≈ 1 n log( γ n 4 + 1) + log  γ 4 + 1  n →∞ − → log  γ 4 + 1  . Next assume W  J such that γ ≈ 1 W . The upper bound in ( 5 ) is then maximized for p = 1 / 2 . On the other hand, the contribution of λ 1 is at most 1 n log( γ n 2 + 1) , which goes to zero for large n . Pr oof sketch of Pr oposition 3.2 : The proof lev erages the fol- lowing result of [ 21 ]. Consider a rev erse circulant matrix 1 √ n B with entries B j i = b j + i − 2 mo d n and ( b 0 , b 1 , . . . , b n ) a sequence of IID random variables with mean zero, unit v ari- ance and bounded third moment. Then, the empirical spectral density (ESD) of B con ver ges to the limiting spectral distri- bution with density f X ( x ) . In our setting, we are interested on the ESD of AA T for A that has entries Bern( p ) . T o apply the result of [ 21 ], consider: e A = ( A − p 11 T ) / p p (1 − p ) . The entries of e A hav e no w zero mean and unit variance. More- ov er , λ j ( e A ) = λ j ( A ) / p p (1 − p ) for j = 2 , . . . , n . It can be shown that | λ j ( e A ) | 2 = λ 2 j ( B ) [ 21 , Lem. 1]. Applying these to ( 1 ) gi ves I = 1 n P n i =1 log  p (1 − p ) W + pJ · λ 2 i  1 √ n B  + 1  n →∞ → e I p , where the con vergence result follo ws from [ 21 ]. Refer ences [1] E. E. Fenimore and T . Cannon, “Coded aperture imag- ing with uniformly redundant arrays, ” Applied optics , vol. 17, no. 3, pp. 337–347, 1978. [2] M. Y oung, “Pinhole, ” Applied Optics , vol. 10, pp. 2763– 2767, 1971. [3] A. L. Cohen, “ Anti-pinhole imaging, ” Journal of Mod- ern Optics , vol. 29, no. 1, pp. 63–67, 1982. [4] A. Levin, R. Fergus, F . Durand, and W . T . Freeman, “Image and depth from a con ventional camera with a coded aperture, ” ACM transactions on graphics (TOG) , vol. 26, no. 3, p. 70, 2007. [5] C. Zhou, S. Lin, and S. Nayar, “Coded aperture pairs for depth from defocus and defocus deblurring, ” Inter- national Journal of Computer V ision , vol. 93, no. 1, pp. 53–72, 2011. [6] R. Raskar , A. Agraw al, and J. T umblin, “Coded ex- posure photography: motion deblurring using fluttered shutter , ” in ACM T ransactions on Graphics (TOG) , vol. 25, no. 3. A CM, 2006, pp. 795–804. [7] A. V eeraraghav an, R. Raskar, A. Agrawal, A. Mohan, and J. T umblin, “Dappled photography: mask enhanced cameras for heterodyned light fields and coded aperture refocusing, ” ACM T ransactions on Graphics (TOG) , vol. 26, no. 3, p. 69, 2007. [8] M. F . Duarte, M. A. Davenport, D. T akbar , J. N. Laska, T . Sun, K. F . Kelly , and R. G. Baraniuk, “Single-pixel imaging via compressi ve sampling, ” IEEE signal pro- cessing magazine , v ol. 25, no. 2, pp. 83–91, 2008. [9] M. S. Asif, A. A yremlou, A. Sankaranarayanan, A. V eeraraghavan, and R. Baraniuk, “Flatcam: Thin, bare-sensor cameras using coded aperture and compu- tation, ” arXiv pr eprint arXiv:1509.00116 , 2015. [10] C. Thrampoulidis, G. Shulkind, F . Xu, W . T . Freeman, J. H. Shapiro, A. T orralba, F . N. W ong, and G. W . W or- nell, “Exploiting occlusion in non-line-of-sight acti ve imaging, ” arXiv pr eprint arXiv:1711.06297 , 2017. [11] A. T orralba and W . T . Freeman, “ Accidental pinhole and pinspeck cameras: Re vealing the scene outside the picture. ” Computer V ision and P attern Recognition , pp. 374–381, 2012. [12] K. Bouman, V . Y e, A. Y edidia, F . Durand, G. W ornell, A. T orralba, and W . T . Freeman, “T urning corners into cameras: Principles and methods, ” International Con- fer ence on Computer V ision , 2017. [13] C. Bouman and K. Sauer , “ A generalized gaussian im- age model for edge-preserving map estimation, ” IEEE T ransactions on Image Pr ocessing , 1993. [14] A. V an Der Schaaf and J. H. van Hateren, “Modelling the power spectra of natural images: statistics and in- formation. ” V ision Resear c h , vol. 36, no. 17, pp. 2759– 2770, 1996. [15] K. Ikeuchi, Lambertian Reflectance , 2014. [16] F . J. MacW illiams and N. J. Sloane, “Pseudo-random sequences and arrays, ” Pr oceedings of the IEEE , vol. 64, no. 12, pp. 1715–1729, 1976. [17] S. R. Gottesman and E. Fenimore, “New family of bi- nary arrays for coded aperture imaging, ” Applied optics , vol. 28, no. 20, pp. 4344–4352, 1989. [18] M. J. Cie ´ slak, K. A. Gamage, and R. Glover , “Coded- aperture imaging systems: Past, present and fu- ture dev elopment–a revie w , ” Radiation Measur ements , vol. 92, pp. 59–71, 2016. [19] A. Busboom, H. Elders-Boll, and H. Schotten, “Uni- formly redundant arrays, ” Experimental Astr onomy , vol. 8, no. 2, pp. 97–123, 1998. [20] J. Ice, N. Narang, C. Whitelam, N. Kalka, L. Hornak, J. Dawson, and T . Bourlai, “Swir imaging for facial im- age capture through tinted materials. ” Pr oceedings of SPIE , vol. 8353, 2012. [21] A. Bose and J. Mitra, “Limiting spectral distribution of a special circulant, ” Statistics & pr obability letters , vol. 60, no. 1, pp. 111–120, 2002. [22] A. Bose, R. S. Hazra, K. Saha, et al. , “Limiting spectral distribution of circulant type matrices with dependent inputs, ” Electr on. J. Pr obab , vol. 14, no. 86, pp. 2463– 2491, 2009.