Spatio-temporal Compressed Sensing with Coded Apertures and Keyed Exposures

1.3 MB·English

by  Zachary T. Harmany

1
Spatio-temporal Compressed Sensing with Coded
Apertures and Keyed Exposures
Zachary T. Harmany, Student Member, IEEE, Roummel F. Marcia, Member, IEEE,
and Rebecca M. Willett, Senior Member, IEEE
Abstract—Opticalsystemswhichmeasureindependentrandom the compressive data acquisition. Our approach is based on
projections of a scene according to compressed sensing (CS) two imaging techniques called coded apertures and keyed
theory face a myriad of practical challenges related to the size
exposures, which we explain next.
of the physical platform, photon efficiency, the need for high
2
temporal resolution, and fast reconstruction in video settings. Coded apertures [14], [15] have a long history in low-
1
This paper describes a coded aperture and keyed exposure light astronomical applications. Coded apertures were first
0
approach to compressive measurement in optical systems. The
developed to increase the amount of light hitting a detector
2
proposed projections satisfy the Restricted Isometry Property
in an optical system without sacrificing resolution (by, say,
n for sufficiently sparse scenes, and hence are compatible with
a theoreticalguaranteesonthevideoreconstructionquality.These increasing the diameter of an opening in a pinhole camera).
J concepts can be implemented in both space and time via The basic idea is to use a mask, i.e., an opaque rectangu-
6 either amplitude modulation or phase shifting, and this paper lar plate with a specified pattern of openings, that allows
describes the relative merits of the two approaches in terms of
significantly brighter observations with higher signal-to-noise
theoreticalperformance,noiseandhardwareconsiderations,and
] ratiothanthosefromconventionalpinholecameras[14],[15].
P experimental results. Fast numerical algorithms which account
A forthenonnegativityoftheprojectionsandtemporalcorrelations These masks encode the image before detection, inducing a
inavideosequencearedevelopedandappliedtomicroscopyand more complicated point spread function than that associated
.
t short-wave infrared data. with a pinhole aperture. The original scene is then recovered
a
t Index Terms—Compressive sensing, coded apertures, convex from the encoded observations in post-processing using an
s
[ optimization, sparsity, coded exposure, Toeplitz matrices appropriate reconstruction algorithm which exploits the mask
pattern.Thesemultiplexingtechniquesareparticularlypopular
2
in astronomical [16], [17] and medical [18], [19], [20] appli-
v I. INTRODUCTION
cations because of their efficacy at wavelengths where lenses
7 THE theory of compressed sensing (CS) suggests that we
4 cannot be used, but recent work has also demonstrated their
can collect high-resolution imagery with relatively few
2 utility for collecting both high resolution images and object
photodetectors or a small focal plane array (FPA) when the
7 depth information simultaneously [21].
. sceneissparseorcompressibleinsomedictionaryorbasisand
1 the measurements are chosen appropriately [1], [2]. There has Keyed exposure (also called coded exposure [22], flutter
1
1 been significant recent interest in building imaging systems shutter [23], or coded strobing [24]) imaging was initially
1 in a variety of contexts to exploit CS (cf. [3], [4], [5], [6], developed to facilitate motion deblurring in video using a
: [7], [8], [9], [10], [11], [12]). By designing optical sensors to relatively low frame rate. In some cases motion has been
v
i collect measurements of a scene according to CS theory, we inferredfromasinglekeyedexposuresnapshot.Thebasicidea
X can use sophisticated computational methods to infer critical is that the camera sensor continuously collects light while the
r scene structure and content. One particularly famous example shutter is rapidly opened and closed; the shutter movement
a
in optical imaging is the Single Pixel Camera [13], which effectively modulates the motion blur point spread function,
collectsindividualprojectionsofthescenesequentially.While and with well-chosen shutter movement patterns it becomes
these measurements are supported by the CS literature, there possible to deblur moving objects. Similar effects can be
are several practical challenges associated with the tradeoff achieved using a strobe light instead of moving a shutter.
between temporal resolution physical system footprint. In this Despitetheutilityoftheabovemethodsinspecificsettings,
paper we describe an alternative approach to designing a low they both face some limitations. The design of conventional
frame-rate snapshot CS camera which naturally parallelizes coded apertures does not account for the inherent structure
and compressibility of natural scenes, nor the potential for
This research is supported by DARPA Contract No. HR0011-04-C-0111,
nonlinear reconstruction algorithms. Likewise, existing keyed
ONR Grant No. N00014-06-1-0610, DARPA Contract No. HR0011-06-C-
0109,andNSF-DMS-08-11062.Portionsofthisworkwerepresentedatthe exposure methods focus on direct (uncoded) measurements of
IEEE International Conference on Acoustic, Speech, and Signal Processing, thespatialcontentofthesceneandhavelimitedreconstruction
March2008,andattheSPIEElectronicImagingConference,January2009,
capabilities, as we detail below. The Coded Aperture Keyed
andSPIEPhotonicsEurope,April2010.
Z.T.HarmanyandR.M.WillettarewiththeDepartmentofElectricaland Exposure (CAKE) sensing paradigm we propose in this paper
Computer Engineering, Duke University, Durham, NC 27708 USA (e-mail: is designed to allow nonlinear high-resolution video recon-
[e-mail],[e-mail]).
struction from relatively few measurements in more general
R.F.MarciaiswiththeSchoolofNaturalSciences,UniversityofCalifornia,
Merced,CA95343USA(e-mail:[e-mail]). settings.
2
A. Problem Formulation includingamean-subtractionpre-processingstepwhichallows
We consider the problem of reconstructing an N-frame us to sidestep challenging theoretical aspects associated with
video sequence f(cid:63), where each frame is an n × n two- nonnegative sensing matrices (such as amplitude modulating
1 2
dimensional image denoted f(cid:63). Using standard vector repre- coded apertures). This paper builds substantially upon earlier
t
sentation, we have that f(cid:63) ∈ Rn for t = 1,...,N where preliminary studies by the authors [25], [26], [27] and related
t
n(cid:44)n n isthetotalnumberofpixels.Asaresult,thevector independent work by Romberg [28].
1 2
representation of the video sequence is f(cid:63) = (f(cid:63),...,f(cid:63)) ∈
1 N
RnN.
C. Organization of the Paper
The observations y of f(cid:63) are also acquired as a video
The paper is organized as follows. Section II-A describes
sequence.Wedonotassumethattheobservationsareacquired
at the same rate at which we will ultimately reconstruct f(cid:63). conventional coded aperture imaging techniques and Sec-
tion II-Bdescribes keyedexposure techniquescurrently inthe
In general, we assume y is an M-frame video sequence, with
each frame y of size m ×m . Similarly to f(cid:63), we have literature. We describe the compressive sensing problem and
k 1 2
y ∈ Rm for k = 1,...,M, where m (cid:44) m m , therefore formulate it mathematically in Section II-C. In Section III,
k 1 2
y =(y ,...,y )∈RmM. we show how CS theory can be used for constructing coded
1 M
We observe f(cid:63) via a spatio-temporal sensing matrix A ∈ aperture masks that can easily be implemented for improving
RmM×nN which linearly projects the spatio-temporal scene image reconstruction resolution in snapshot imaging; this
includes theoretical results, a discussion of implementation
onto an mM-dimensional set of observations:
details and tradeoffs in practical optical systems, and experi-
y =Af(cid:63)+w, (1) mental results. We consider applications to video compressed
where w ∈ RmM is noise associated with the physics of the sensing using the full CAKE paradigm in Section IV, includ-
ing theory and experimental results.
sensor.
CS optical imaging systems must be designed to meet
several competing objectives: II. BACKGROUND
• ThesensingmatrixAmustsatisfysomenecessarycriterion
Prior to detailing our main contributions, we first review
(suchastheRIP,definedbelow)whichprovidestheoretical
pertinent background material. This review touches upon the
guaranteesontheaccuracywithwhichwecanestimatef(cid:63)
development of coded aperture imaging, coded exposure pho-
from y.
tography,andabriefreviewofsalientconceptsincompressed
• The total number of measurements, mM, must be lower sensing theory.
than the total number of pixels to be reconstructed, nN.
This is achievable via compressive spatial acquisition
A. Coded Aperture Imaging
(m<n), frame rate reduction (M <N), or simultaneous
spatio-temporal compression. Seminal work in coded aperture imaging includes the de-
velopmentofmasksbasedonHadamardtransformoptics[29]
• The measurements y must be causally related to the
temporal scene f(cid:63), which restricts the structure of the and pseudorandom phase masks [30]. Modified Uniformly
Redundant Arrays (MURAs) [31] are generally accepted as
projections A.
optimal mask patterns for coded aperture imaging. These
• The optical measurements modeled by A must be imple- mask patterns (which we denote by hMURA) are binary,
mentable in a way that results in a smaller, cheaper, more
square patterns, whose grid size matches the spatial resolu-
robust, or lower power system.
tion of the photo-detector. Each mask pattern is specifically
• Thesensingmatrixstructuremustfacilitatefastreconstruc- designed to have a complementary pattern hMURA such that
tion algorithms. hMURA∗hMURA is a single peak with flat side-lobes (i.e., a
This paper demonstrates that compressive Coded Aperture Kronecker δ function).
Keyed Exposure systems achieve all these objectives. In practice, the resolution of a detector array dictates the
properties of the mask pattern and hence resolution at which
B. Contributions f(cid:63) can be reconstructed. We model this effect as f(cid:63) being
downsampled to the resolution of the detector array and then
The primary contribution of this paper is the design and
convolved with the mask pattern hMURA, which has the same
theoretical characterization of compressive Coded Aperture
resolution as the FPA and the downsampled f(cid:63), i.e.,
Keyed Exposure (CAKE) sensing. We explore amplitude
modulating and phase shifting masks and describe theoretical y =(D f(cid:63))∗hMURA+w, (2)
int
and implementation aspects of both. We further describe how
keyed exposure ideas can be used in conjunction with coded where ∗ denotes circular convolution, w corresponds to noise
apertures to increase both the spatial and temporal resolution associated with the physics of the sensor, and D f(cid:63) is the
int
ofvideofromrelativelyfewmeasurements.Weprovehitherto integration downsampling of the scene, which consists of
unknown theoretical properties of such systems and demon- partitioning f(cid:63) into uniformly sized d × d blocks, where
1 2
strate their efficacy in several simulations. In addition, we d (cid:44) n /m for i = 1,2, and measuring the total intensity in
i i i
discussseveralimportantalgorithmicaspectsofourapproach, each block.
3
BecauseoftheconstructionofhMURA andhMURA,D f(cid:63) The assumption of a periodic video makes it possible to
int
can be reconstructed using apply much more general reconstruction algorithms that do
not require background subtraction or separating different
f(cid:98)=y∗hMURA.
moving components. However, it is a very strong assumption,
However, the resulting resolution is often lower than what is which places some limits in its applicability to real-world
necessary to capture some of the desired details in the image. settings. The approach described in this paper has similar
Clearly, the estimates from MURA reconstruction are limited performance guarantees but operates on much more general
by the spatial resolution of the photo-detector. Thus, high video sequences.
resolution reconstructions cannot generally be obtained from
low-resolution MURA-coded observations. C. Compressed Sensing
It can be shown that this mask design and reconstruction
In this section we briefly define the Restricted Isometry
result in minimal reconstruction errors at the FPA resolution
Property (RIP) and explain its significance to reconstruction
and subject to the constraint that linear, convolution-based
performance. In subsequent sections, we demonstrate our
reconstruction methods would be used. However, when the
primary theoretical contribution, which is to prove the RIP
scene of interest is sparse or compressible, and nonlinear
for compressive coded aperture and keyed exposure systems.
sparse reconstruction methods may be employed, then CS
ideascanbeusedtodesigncodedaperturewhichyieldhigher Definition 1 (Restricted Isometry Property (RIP) [32]). A
resolution images. Before describing the details of this, we matrix A satisfies the RIP of order s if there exists a constant
briefly review two key relevant concepts from CS. δs ∈(0,1) for which
(1−δ )(cid:107)f(cid:107)2 ≤(cid:107)Af(cid:107)2 ≤(1+δ )(cid:107)f(cid:107)2. (3)
s 2 2 s 2
B. Coded (Keyed) Exposure Imaging
holds for all s-sparse f ∈ Rn. If this property holds, we say
Coded (or keyed) exposures were developed recently in the
A is RIP(s,δ ).
s
computational photography community. Initial work in this
area was focused on engineering the temporal component of Matrices which satisfy the RIP are called CS matrices;
a motion blur point spread function by rapidly opening and and when combined with sparse recovery algorithms, they
closingtheshutterduringasingleexposureorasmallnumber are guaranteed to yield accurate estimates of the underlying
of exposures at a low frame rate [22], [23]. That is, function f(cid:63):
y =AKEf(cid:63)+w =(cid:88)f(cid:63)+w, Theorem 1 (Sparse Recovery with RIP Matrices [33√], [34]).
i
Let A be a matrix satisfying RIP(2s,δ ) with δ < 2−1,
i∈S 2s 2s
and let y = Af + w be a vector of noisy observations of
where the keyed exposure (KE) measurement matrix AKE any signal f ∈ Rn, where the w is a noise or error term
selects the subset of frames during which the shutter is open.
with (cid:107)w(cid:107) ≤ (cid:15) for some (cid:15) > 0. Let f be the best s-sparse
2 s
We refer to this subset as the exposure code S ⊆{1,...,N}.
approximation of f; that is, f is the approximation obtained
s
Ifanobjectismovingduringimageacquisition,thenastatic
by keeping the s largest entries of f and setting the others to
shutterwouldinduceatypicalmotionblur,makingthemoving
zero. Then the estimate
object difficult to resolve with standard deblurring methods.
However, by “fluttering” the shutter during the exposure f(cid:98) = argmin (cid:107)f(cid:107)1
f∈Rn (4)
using carefully designed patterns, the induced motion blur
subject to (cid:107)y−Af(cid:107) ≤(cid:15),
can be made invertible and moving objects can be accurately 2
reconstructed. Instead of a moving shutter, more recent work obeys
useWshailsetrtohbisenliogvhetltaoppprrooadcuhcetoavsiidmeoilaarceqfufiescittio[2n4c]a.n produce (cid:107)f −f(cid:98)(cid:107)2 ≤C1,s(cid:15)+C1,s(cid:107)f −√sfs(cid:107)1,
veryaccuratedeblurredimagesofmovingobjects,thereissig-
where C and C are constants which depend on s but not
1,s 2,s
nificant overhead associated with the reconstruction process.
on n or m.
To see why, note that every object moving with a different
velocityortrajectorywillproduceadifferentmotionblur.This Note that the reconstruction (4) in Theorem 1 is equivalent to
means that (a) any stationary background must be removed
duringpreprocessingand(b)multiplemovingobjectsmustbe f(cid:98)=argmin12(cid:107)y−Af(cid:107)22+τ(cid:107)f(cid:107)1 (5)
f∈Rn
separated and processed individually.
where τ > 0, which depends on (cid:15), can be viewed as a
More recently, it was shown that these challenges could
regularization parameter.
be sidestepped when the video is temporally periodic (e.g.,
consider a video of an electronic toothbrush spinning) [24].
TheperiodicassumptionamountstoasparsetemporalFourier III. COMPRESSIVECODEDAPERTURESFORSNAPSHOT
transform of the video, and this approach, called coded strob- IMAGING
ing,isacompressiveacquisitioninthetemporaldomain.Asa This section first considers a snapshot acquisition model.
result,theauthorswereabletoleverageideasfromcompressed Our goal is to recover a static high-resolution scene from a
sensing to achieve high-quality video reconstruction. single image where all pixels are collected simultaneously. In
4
termsofthenotationinSec.I-A,wehaveN =M =1,sothat
Aisofsizem×n.Thesensingmatrixforcompressivecoded
aperture (CCA) systems can be modeled mathematically as
Af(cid:63) =D(f(cid:63)∗h); (6)
where h is a coding mask, and D is a subsampling operator
(detailed below). Here the coding mask, h, is at the size
and resolution at which f(cid:63) will be reconstructed; this is in
contradistinction to the MURA system, in which hMURA is at
the size and resolution of the FPA. Thus in (6), we model the
measurements as the scene being convolved with the coded
mask and then downsampled. Fig.1. Then×nmatrixF−1diag(Fh)Fisblock-circulantwithn2blocks
Using a similar model, Romberg [28] conducted related ineachrowandcolumn.Eachblockisn1×n1 andiscirculant.
work concurrent and independent of our initial investigations
[25]. We will summarize the key features of this model
and compare with our approach in Sec. III-B. While these R is consists of n ×n blocks,
2 2
modelssharecommonelements,wewillseeinSec.III-Cthat
 
R R ··· R R
there are important tradeoffs associated with the theory and 1 n2 3 2
implementation of each strategy. R2 R1 ··· R4 R3
RecentworkbyBajwaetal.[35],[36],showedthatrandom R= ... ... ... ... ... , (8)
circulant matrices (and Toeplitz matrices, in general) are R R ··· R R
sufficient to recover sparse f(cid:63) from y with high probability. n2 n2−1 2 1
In particular, they showed that a Toeplitz matrix whose first where each Rj ∈Rn1×n1 is circulant; i.e., Rj is of the form
row contains elements drawn independently from a Gaussian
 
distribution are RIP(s,δs) when m ≥ Cs2log(n) for some rj,1 rj,n1 ··· rj,3 rj,2
constant C. Here, we extend these results to pseudo-circulant rj,2 rj,1 ··· rj,4 rj,3
matrices and use them to motivate our mask design. Our Rj = ... ... ... ... ... ,
model differs from that in [28] in that we consider a different
r r ··· r r
generative model for the coded aperture mask, as well as a j,n1 j,n1−1 j,2 j,1
different subsampling strategy: (seeFig.1).Thisblock-circulantwithcirculant-block(BCCB)
• The elements of the coding mask h are generated iid structure of R is a direct result of the fact that the k-point
according to a particular generating distribution (e.g., an one-dimensional Fourier transform F diagonalizes any k×k
k
appropriately scaled Rademacher, uniform, or Gaussian circulant matrix (such as R with k =n ) and so F ≡F ⊗
j 1 n2
distribution). F diagonalizes block-circulant matrices (such as R). Here
n1
• Our analysis allows for deterministic downsampling in ⊗ denotes the Kronecker matrix product.
which we collect one sample per nonoverlapping block. We now examine the generation of a compressed sensing
In our notation, we distinguish mask patterns in this manner matrix from the convolution R by restricting the number of
using the abbreviations BS, US, and GS, where the first measurements collected of the vector Rf(cid:63). Here we simply
character denotes the distribution used to generate the mask, subsample the vector Rf(cid:63) by applying a pointwise subsam-
i.e.,binary,uniform,orGaussian,andthesecondisareminder pling matrix Dsub. The operation of applying Dsub consists of
that these masks are generated directly in the spatial domain. retaining only one measurement per uniformly sized d1×d2
block, i.e., we subsample by d in the first coordinate and d
1 2
inthesecondcoordinatesothattheresultisanm ×m image
1 2
A. CCAs Generated in the Spatial Domain
withm =n /d ,i=1,2.InmatrixformD canbethought
i i i sub
The two-dimensional convolution of h with an image f(cid:63) ofasretainingacertainnumberofrowsoftheidentitymatrix.
as in (6) can be represented as the application of the Fourier Because of the structure and deterministic nature of this type
transform to f(cid:63) and h, followed by element-wise multiplica- of downsampling, it is often more straightforward to realize
tionandapplicationoftheinverseFouriertransform.Inmatrix in practical imaging hardware (see Sec. III-C).
notation, this series of linear operations can be expressed as The resulting projection matrix A is then given by
vect(f(cid:63)∗h)=F−1diag(Fh)Ff(cid:63) =Rf(cid:63), (7) A=D R=D F−1diag(Fh)F.
sub sub
where vect(f) is a vectorized representation of an image
We now show that if the elements of h are drawn from an
f, F is the two-dimensional Fourier transform matrix, and
appropriate probability distribution, then A will be RIP(s,δ )
diag(Fh) is a diagonal matrix whose elements correspond s
with high probability.
to the transfer function, which is the Fourier transform of h.
The matrix product R = F−1diag(Fh)F ∈ Rn×n is block- Theorem 2 (Spatial-Domain CCA Sensing). Let hBS be a
circulantandeachblockisinturncirculant.Inmatrixnotation, mask with entries generated i.i.d. according to the scaled
5
Rademacher distribution B. CCAs Generated in the Fourier Domain
(cid:40) (cid:112) The random convolution in [28] follows the same structure
d/n with probability 1/2,
hBk1S,k2 = −(cid:112)d/n with probability 1/2, (9) aissgtehnaetrianteSderca.nIdIoI-mAl.yHinowtheevferre,qiunenthcaytdwoomrkainth.eMcoornevsopluectiiofin-
cally,theentriesofthetransferfunctioncorrespondtorandom
for k1 = 1,...,n1,k2 = 1,...,n2. Let ABS = DsubR be an phase shifts (with some constraints to keep the resulting
m×nmatrix,whereR∈Rn×n istheBCCBmatrixgenerated observation matrix real-valued). We denote the resulting con-
from hBS, and Dsub ∈ Rm×n is the pointwise subsampling volutionkernelashUP,where“UP”refersto“uniformphase”.
matrix. Then, there exists constants c1,c2 ≥ 0 depending on For simplicity of presentation, we describe the generating
δs such that for any distributionforaone-dimensionalconvolutionforeven-length
masks. In particular, let σ =FhUP, and generate σ such that
m≥c1s2log(n), (10) for k =1,...,n,

ABS is RIP(s,δs) with probability exceeding ±eiφ1 wwiitthh eφq∼ualUp(r0o,b2aπb)ility iiff k2≤=k1,≤n/2,
σ =
1−2n2e−c2m/s2. (11) k ±σ∗1 with equal probability iiff kn/=2+n/22≤+k1≤n.
n−k+2
Remark 1. This binary-valued distribution was selected to (13)
model coded apertures with two states – open and closed – Here U(0,2π) denotes the uniform distribution over [0,2π].
per mask element. We discuss the issue of implementing these The real-valued convolution kernel is then given by hUP =
masks in Section III-C1. F−1σ.Thisresultcanbeextendedeasilyfortwo-dimensional
convolutions.
Remark 2. It is straightforward to extend the result to other
Toformacompressedsensingmatrixfromthisrandomcon-
mask generating distributions, such as uniform and Gaussian,
volution,[28]considerstwodifferentdownsamplingstrategies:
which we denote hUS and hGS, with associated sensing
samplingatrandomlocations,andrandomdemodulation.The
matrices AUS and AGS.
firstmethodentailsselectingarandomsubsetΩ⊂{1,...,n}
We present the proof of Theorem 2 in Appendix A. of indices. We form a downsampling matrix DΩ by retaining
For successful recovery, the coded aperture masks are only the rows of the identity matrix indexed by Ω, hence the
designed to be satisfy RIP(2s,δ ) as described in (3) with resulting measurement matrix is given by
2s
high probability when m ≥ c1s2log(n). Note that this is AUP =D F−1diag(σ)F, (14)
Ω
somewhat weaker than what can be achieved by a purely
unstructured i.i.d. randomly generated sensing matrix which whereσ isgeneratedaccordingtothetwo-dimensionalanalog
satisfies (3) with high probability when m ≥ c˜ slog(n/s) of(13).Therandomdemodulationmethodmultipliestheresult
1
for some constant c˜ > 0. Intuition may lead one to believe of the random convolution by a random sign sequence s ∈
1
the extra factor of s is due to the fact that the m projections {−1,1}n, such that s = ±1 with equal probability for all
i
sensedusingtheamplitudemodulatedmaskframeworkexhibit i = 1,...,n, then performs an integration downsampling of
dependencies. However, in many settings this theoretical dis- the result. Therefore in this case the measurement matrix is
advantageisoffsetbyadvantageouspracticalimplementations.
AUP =D diag(s)F−1diag(σ)F. (15)
SparseGradientsScenes: Theabovetheoryisapplicableto int
reconstructingf(cid:63) whenf(cid:63) issparseinthepixelbasis;similar It can be shown that both of these strategies yield RIP-
results hold when the gradient of f(cid:63) is sparse. For simplicity satisfying matrices with high probability.
ofnotationwewillshowthenina1dsetting;theextensionto
Theorem 3 (Fourier-Domain CCA Sensing). Let AUP be an
2d is straightforward. Let ∇f denote the first-order gradient
m×n sensing matrix resulting from random convolution with
of f, so that
phase shifts followed by a downsampling strategy, and let W
 1 0 0 ··· 0  denote an arbitrary orthonormal basis. If the downsampling
−1 1 0 ... ...  ics r>an0dosmuchsutbhsaatmwpiltihngpraosbainbil(i1ty4)e,xtcheeenditnhger1e−exOist(sn−co1n),stfaonrt
  3
∇(cid:44) 0 −1 1 ... 0 . (12) m≥c δ−2min(slog6n,s2log2n),
  3
 ... ... ... ... 0  AUPW satisfies RIP(2s,δ2s) with δ2s ≤ δ. If the down-
0 ··· 0 −1 1 sampling is random demodulation as in (15), then there
exists constant c > 0 such that with probability exceeding
4
Asnotedin[35],f(cid:63) =∇−1θ,whereθisthesparsegradient 1−O(n−1), for
image. Thus if we sense y =A∇f(cid:63)+w =A∇∇−1θ+w = m≥c δ−2min(slog6n,s2logn),
4
Aθ+w,wecanexpecttorecoverθusingsparsereconstruction
methods. AUPW satisfies RIP(2s,δ ) with δ ≤δ.
2s 2s
6
These theoretical results are stronger than those in The- halftoning procedure, which is easiest to implement at the
orem 2, especially since they allow sparsity in arbitrary necessary resolution with chrome on quartz.
orthonormal bases. However, there are important differences Both Robucci et al. [37] and Majidzadeh et al. [38] have
between the observation models which have a significant proposed performing the analog, random convolution step in
impactonthefeasibilityandeaseofhardwareimplementation. complementary,metal-oxide-semiconductor(CMOS)electron-
WeelaborateonthisinSectionIII-C.Nevertheless,thetheory ics.Aclearadvantagetothisarchitectureisthattheadditional
developed in [28] lends important theoretical support to the optics required for spatial light modulation are removed in
general concept of compressive coded apertures. favor of additional circuitry, immediately reducting imager
size.
C. Hardware and Practical Implementation Considerations 2) Phase Shift Masks: Phase shifting masks for coded
aperture imaging have been implemented recently using a
In this section we describe how we shift from modeling
phase screen [39]. This approach allows one to account for
the coded aperture masks in a way that is compatible with
diffraction in the optical design. However, depending on
compressed sensing theory, to a model that describes their
the precise optical architecture, phase shift masks may be
actual implementation in an optical system. Our analysis does
much less photon-efficient than amplitude modulation masks.
not account for the bandwidth of the lenses; in particular,
Additionally, the mask generation distribution described in
we implicitly assume that the bandwidth of the lenses is
Eq. (13) will result in negative entries for the corresponding
high enough that band-limitation effects are negligible at the
PSF hUP = F−1σ. To compensate for this, the phase mask
resolutionofinterest.Inallofthehardwaresettingsdescribed
must be mean-shifted to make all entries nonnegative, so that
below, precise alignment of the optical components (e.g., the
we are actually implementing
mask and the focal plane array) is critical to the performance
oftheproposedsystem.Oftenahigh-resolutionFPAishelpful hUP =c(hUP−min(hUP)1), (16)
+
for alignment and calibration.
with the constant c selected so that the implementable PSF
In this paper we focus on incoherent light settings (con-
hUP is flux-preserving and 1 is a matrix of ones and of the
sistent with many applications in astronomy, microscopy, and +
same dimension as hUP.
infrared imaging). In this case, the coded aperture must be
3) ImplementableMasksforSceneswithSparseGradients:
real-valued and flux-preserving (i.e., the light intensity hitting
As described in Section III-A, theoretical results for scenes
the detection cannot exceed the light intensity of the source).
with sparse gradients hold for the sparse difference operator
In this section, we consider the following apertures:
• Binary Spatial Mask: h ∈ {0,1/n}n1×n2, drawn with ∇ defined in (12). The problem with this approach is that ∇
is invertible but not circulant – and hence the sensing matrix
equal probability,
A∇ cannot be implemented with a coded aperture system.
• Uniform Spatial Mask: h ∈ [0,1/n]n1×n2, where each We address this by noting that if the upper right element of
element is drawn independently from a uniform distribu-
∇ were −1 instead of 0, then the resulting sensing matrix,
tion, or
denoted ∇(cid:101), could be implemented physically and is a close
• Uniform Phase Mask: h = |Fp|2 for some p, where p approximation to the theoretically supported ∇. In short, the
correspondstoaphase-shiftingmaskinanincoherentlight theoretically supported solution is to set
setting.
y =A∇f(cid:63)+w
1) AmplitudeModulationMasks: Inaconventionallensless
coded aperture imaging setup, the point spread function asso- θ(cid:98)=argmin12(cid:107)y−Aθ(cid:107)22+τ(cid:107)θ(cid:107)1
ciated with the aperture is the mask pattern h itself. To shift a θ
RIP-satisfying aperture as in Theorem 2 to an implementable f(cid:98)=∇−1θ(cid:98)
aperture, one simply needs to apply an affine transform to
(cid:112) (cid:112) while an implementable close approximation is to set AG (cid:44)
h mapping [− d/n, d/n] to [0,1/n]. This transform en-
A∇(cid:101) and
sures that the resulting mask pattern is nonnegative and flux-
preserving. y =AGf(cid:63)+w
These amplitude modulating masks may be implemented
using a spatial light modulator (SLM) or placing chrome on f(cid:98)=argmin12(cid:107)y−AGf(cid:107)22+τ(cid:107)f(cid:107)TV,
f
quartz. The SLMs may be preferable in video settings where
the underlying scene contains motion and using a different where (cid:107)·(cid:107)TV is the total variation seminorm which causes f(cid:98)
mask pattern at each time step boosts performance. However, to have sparse gradients. In our numerical experiments we
in order for the proposed approach to work, the mask or compare the performance of TV regularization to sparsity
SLM used must be higher resolution than the FPA. Currently, penalization for the task of reconstructing static scenes.
very high resolution SLMs are still in development. Chrome 4) Downsampling Implementation: In developing RIP-
on quartz masks can be made with higher resolution than satisfying AMM coded apertures, Theorem 2 assumes the
many SLMs, but cannot be changed on the fly unless we subsampling operation selects one measurement per d ×d
1 2
mount a small number of fixed masks on a rotating wheel or block. From an implementation standpoint, this operation
translating stage. The uniform amplitude modulation masks effectivelydiscardsalargeportion((d−1)/d)oftheavailable
in particular could be constructed using a high-resolution light, which would result in much lower signal-to-noise ratios
7
at the detector. A more pragmatic approach is to use larger zero mean and, therefore, negatively impacts the performance
detector elements that essentially sum the intensity over each of these CS-based reconstruction algorithms. In this section,
d ×d block, making a better use of the available light. We we describe how we take full advantage of the structure of A
1 2
call this operation integration downsampling to distinguish it and address how we mitigate the negative effects of A having
from subsampling. The drawback to this approach is that we nonzero-mean.
lose many of the desirable features of the system in terms of 1) Sparsity-promotingMethods: Wenotethatmostpopular
the RIP. Integration downsampling causes a large coherence andeffectivemethodsforperformingtheabovereconstruction
between neighboring columns of the resulting sensing matrix areiterativealgorithmswithrepeatedapplicationsoftheoper-
A. An intermediate approach would randomly sum a fraction ators A and (A)T to scene estimates and residuals [42], [43].
of the elements in each size d block, which increases the InmostCSsettings,computingthesematrix-vectormultiplica-
signal-to-noise ratio versus subsampling, but yields smaller tionsisalargecomputationalburden.However,becauseofthe
expectedcoherence.Thisapproachismotivatedbytherandom circulantstructureofA,thiscomputationisveryefficientusing
demodulation proposed in [28] and described in Sec III-B, fast Fourier transform algorithms. In particular, we compared
wherebythesignalismultipliedbyarandomsequenceofsigns the computation time for reconstructing a 256 × 256 scene
{−1,+1}, then block-wise averaged. The pseudo-random using 128×128 CCA measurements with the time required
summation proposed here can be thought of as an optically to reconstruct the same scene using the same number of
realizable instantiation of the same idea where we multiply measurements computing using a dense random projection
by a random binary {0,1} sequence. We explore the effects matrix; our reconstruction was 250 times faster because of
of these choices in our numerical results section. our exploitation of the Toeplitz structure in A.
5) NoiseandQuantization: WhileCSisparticularlyuseful 2) Mean Subtraction: Generative models for random pro-
when the FPA needs to be kept compact, it should be noted jection matrices used in CS involve drawing elements inde-
thatCSismoresensitivetomeasurementerrorsandnoisethan pendently from a zero-mean probability distribution [1], [44],
more direct imaging techniques. The experiments conducted [32], [35], [45], and likewise a zero-mean distribution was
in this paper simulated very high signal-to-noise ratio (SNR) usedtoanalyzethecodedaperturemasksdescribedinSec.III.
settings and showed that CS methods can help resolve high However, as coded aperture masks with zero mean are not
resolution features in images. However, in low SNR settings physicallyrealizableinopticalsystems,wegenerateourphys-
CS reconstructions can exhibit significant artifacts that may icallyrealizablemasksfromanappropriatelyscaledBernoulli
even cause more distortion than the low-resolution effects distribution with values {0,1/n}. This shifting ensures that
associated with conventional coded aperture techniques such the coded aperture corresponds to an implementable (i.e.,
as MURA. nonnegativeandintensitypreserving)maskpatternwhichdoes
Similar observations are made in [40], which presents a satisfies the assumptions needed to model photon propagation
direct comparison of the noise robustness of CS in contrast to through an optical system.
conventional imaging techniques both in terms of bounds on Whilenecessarytoaccuratelymodelreal-worldopticalsys-
how reconstruction error decays with the number of measure- tems,thisshiftingnegativelyimpactstheperformanceofwell-
ments and in a simulation setup; the authors conclude that for established (cid:96)2-(cid:96)1 reconstruction algorithms for the following
most real-world images, CS yields the biggest gains in high reason. Several of these algorithms (e.g., [42], [43]) assume
signal-to-noise ratio (SNR) settings. Related theoretical work that the (cid:96)2 data fidelity term φ(f) = 21(cid:107)Af −y(cid:107)22 (propor-
in [41] show that in the presence of low SNR photon noise, tional to the negative log Gaussian likelihood) is such that
theoretical error bounds can be large, and thus the expected ∇2φ(x)=ATAcanbewell-approximatedbyαI,α>0.This
performanceofCSmaybelimitedunlessthenumberofavail- assumption is crucial to the performance of these algorithms.
ablephotonstosenseissufficientlyhigh.Theseconsiderations If we collect measurements using a realizable mask with
play an important role in choosing the type of downsampling elements in the range [0,1/n], the resulting matrix A does
to implement. not satisfy this condition due to the mean offset. Therefore in
the reconstruction we use a mean-shifted sensing matrix
Similar issues arise when considering the bit-depth of focal
plane arrays, which corresponds to measurement quantization A0 =A−µ 1 ,
A m×n
errors. Future efforts in designing optical CS systems must
where µ = (1/mn)1TA1 is the mean value of A. This
carefully consider the amount of noise anticipated in the A m n
new matrix is such that (A0)TA0 ≈αI is a valid assumption.
measurements to find the optimal tradeoff between the focal
However, to use this matrix in the reconstruction, we need to
plane array size and image quality.
adjust our estimate accordingly. To compensate for this, we
decompose our estimate f into its mean component µ and a
f
D. Reconstruction zero-mean deviation f0: f =f0+µ 1 . So then
f n
To solve the CS minimization problem (5), we use well- Af =(A0+µ 1 )(f0+µ 1 )
A m×n f n
established gradient-based optimization methods. They are
=A0f0+nµ µ 1 +µ A01 +µ (1Tf0)1
particularly suitable in our setting because the block-circulant A f m f n A n m
≈A0f0+nµ µ 1 .
structure of A described in the previous section allows for A f m
very fast matrix-vector products that are critical to the speed The data y can be decomposed in a similar fashion, so that
oftheirperformance.However,ourobservationmatixAisnot y =y0+µ 1 . Thus a straightforward estimate for the mean
y m
8
Ground Truth Conventional
of our solution is µ = µ /µ n, and we can use the zero-
f y A
mean quantities in our data fidelity term
φ0(f0)= 1(cid:107)A0f0−y0(cid:107)2 (17)
2 2
within the iterative algorithm, which now has the desired
property that ∇2φ0(f0)≈αI.
We have motivated this mean subtraction approach from an
estimationsettingwherewefirstestimatethemean,thendevi-
ations about the mean. However, it can be equally motivated
from a pure optimization perspective by viewing the mean
subtraction step as a preconditioning strategy. The gradient-
basedreconstructionmethodsconsideredareknowntoexhibit
poor convergence when the Hessian has large condition num-
Different mask distributions
ber, and the mean subtraction operation essentially improves
Binary Spatial Uniform Spatial Uniform Phase
the conditioning of the problem by eliminating the single
dominant eigenvalue that occurs due to the nonzero mean.
g
n
pli E. Experimental Comparison on Static Images
m
sa Here we present numerical experiments supporting the
n
w proposed architectures for snapshot compressive imaging. We
o
D compare the simulated performance of the described imaging
on architectures in the task of high-fidelity image reconstruction.
ati Weutilizea256×256pixelrealisticphantomimageofaneu-
r
g
e ron(groundtruthimageinFig.2)whichisdesignedtotestthe
nt
I resolution limits of our architectures with its high-resolution
features. This image is inspired by a real microscopic image
foundathttp://www.lsi.umich.edu/facultyresearch/labs/bingye/
research.Weexamineatraditionalimagingsystemandatotal
of six convolution-based compressive imaging systems. These
n
o systems are constructed by considering each combination of
mati the following design considerations:
m
u • the mask generation method, including binary spatial
S
(BS) and uniform phase (UP) models,
m
o • switching among integration downsamping, pseudo-
d
n random summation (randomly summing d/2 values in
a
R each d ×d block), and subsampling.
1 2
To analyze the effect of the number of measurements we
collect, we vary the scene-to-sensor downsampling ratio d to
be either 4, 16, or 64 (i.e., downsampling in both directions
by a factor of 2, 4, or 8). While each architecture has its own
unique considerations in terms of signal-to-noise efficiency
(see Section III-C), we choose to normalize the experiments
g
n in such a way as to keep a constant signal-to-noise ratio at
pli the detector. This is achieved by fixing the variance of the
m
a additive white Gaussian noise to be σ2 = var(Af)/16 for
s
b
u each architecture when we choose a downsampling ratio of
S
d = 4. This variance is then fixed for each architecture as
we vary the downsampling ratio, the rationale being that this
allows us to showcase the impact of the various subsampling
operations.
Fig. 2. Best performing reconstructions for the different downsampling
We consider a comparison of many penalization schemes
strategies (integration downsampling, random summation, and subsampling)
andforeachdownsamplingstrategy,thedifferentmaskdistributions(binary for image reconstruction, such as sparsity in an orthonormal
spatial, uniform spatial, and uniform phase). For each set of images, the wavelet (Haar) basis, isotropic and anisotropic total-variation
bottomshowstheregionindicatedbytheyellowsquareinthegroundtruth [46], as well as an (cid:96) sparsity penalty in the overcomplete
image. Clearly, significant gains in performance can be achieved by the 1
compressivearchitecturesoverconventionalimaging. curvelet basis [47]. Curvelets are similar to wavelets in that
they capture more spatially localized information than the
Fourier components, however they are designed to be more
9
well-adapted for capturing curvilinear singularities in images Theorem 4 (Coded Aperture Keyed Exposure Sensing). Let
(e.g., edges). To reconstruct the image, we use our own A=[A ···A ] be an m×nB sensing matrix for the CAKE
1 B
implementation of the SpaRSA algorithm [42] which utilizes systemwherethecodedaperturepatternforeachA isdrawn
t
themeansubtractionproceduredescribedinSectionIII-D2.A i.i.d.fromasuitableprobabilitydistribution.Thenthereexists
coarse initialization is found by solving an unpenalized least- constants c ,c ≥0 depending on δ such that for any
1 2 s
squares minimization via a conjugate gradient method using
m≥c s2log(nB), (20)
only the compressive data. We terminate the reconstruction 1
when the relative change in the iterates falls below a pre- A is RIP(s,δ ) with probability exceeding
s
specified tolerance of 1e-3. We choose the regularization 1−2n2B2e−c2m/s2. (21)
weighting to minimize the final reconstruction RMSE, mea-
sured by RMSE(f(cid:98)) = (cid:107)f(cid:98)− f(cid:63)(cid:107)2/(cid:107)f(cid:63)(cid:107)2. A table of the Note here that s denotes the sparsity of the first B frames
resulting RMSE values is presented in Table I. of the video sequence, rather than simply the sparsity of
As evidenced by both the RMSE values and the images, an individual frame as in Thm. 2. The proof of Thm. 4
significant gains in performance can be achieved by the is presented in Appendix B. Similar to Thm. 2, this proof
compressive architectures. The UP mask architectures slightly assumesthattheentriesofeachh ,t=1,...,B aregenerated
t
outperform the BS-based approaches in this simulation, but i.i.d. according to the scaled Rademacher distribution (9).
the type of subsampling has a greater effect on the quality Again, it is straightforward to extend the result to other mask
of the reconstruction. Focusing on the d = 4 case, we see generating distributions.
that subsampling performs best overall, followed by random
summation, then integration downsampling. This is readily A. Sparse Transformation Sensing
apparent from the RMSE values and the greater clarity by It is most often the case that the original frames f(cid:63) are
which we can resolve fine-scale features in the dendrites of not sparse, but can be sparsely represented by a temporal
the neuron. This is exactly predicted by the theory: larger transformoftheframes.Forsimplicity,considerthefirstcoded
coherence yields poorer performance. However, pushing the measurement.Withinthisacquisition,thecoefficientsequence
level of undersampling, we see that simple subsampling is
B
not always the solution, since the performance degrades quite θ(cid:63) =(cid:88)W f(cid:63), k =1,...,B,
rapidly, as the per-measurement signal to noise ratio does k k,t t
t=1
not naturally increase as it does with the other subsampling
may be more sparse for a well-chosen sparse temporal trans-
architectures. An optimal architecture must strike a careful
form W. Notationally this is equivalent to θ(cid:63) =(W ⊗I )f(cid:63),
balance between low coherence and robustness to noise. n
where ⊗ denotes the Kronecker matrix product. A preferred
sensing strategy would be to use our RIP-satisfying CAKE
IV. CODEDAPERTUREKEYEDEXPOSURESENSINGFOR acquisition to sense the coefficient sequence θ(cid:63):
DYNAMICSCENES
B
Here we detail our compressive video acquisition method y =(cid:88)A θ(cid:63)+w. (22)
k k
thatcombinescodedaperturesandkeyedexposurestoaddress
k=1
allthecompetingchallengesdetailedinSec.I-A.Specifically,
Surprisingly, this can be accomplished using an identical
in our CAKE imaging method, each observed frame y is
k architecture,withsomeslightadjustmenttothecodedaperture
given by an exposure of B high-rate coded observations:
mask patterns used during the keyed acquisition.
B Toseethis,weexaminetheresultingsensingmatrixinterms
(cid:88)
yk = A(k−1)B+tf((cid:63)k−1)B+t+wk, (18) of f(cid:63):
t=1 B B B B (cid:34) B (cid:35)
(cid:88) (cid:88) (cid:88) (cid:88) (cid:88)
where each At describes an AMM CCA sensing matrix. Note Akθk(cid:63) = Ak Wk,tft(cid:63) = Wk,tAk ft(cid:63)
that since in our theory, the downsampling operator D is a
sub k=1 k=1 t=1 t=1 k=1
structured nonrandom operator, we can rewrite the above as B
(cid:88)
(cid:34) B (cid:35) = AWt ft(cid:63),
(cid:88)
y =D R f(cid:63) +w . (19) t=1
k sub (k−1)B+t (k−1)B+t k
where we define
t=1
Hence all that is required to implement this sensing paradigm (cid:88)B
AW = W A . (23)
is modulating the coded aperture mask over the B-frame t k,t k
exposure time. Because of this, one can think of our system k=1
as performing a coded exposure acquisition for each coded Therefore (22) is also a CAKE acquisition using the sensing
aperture element. matrices AWt . Because of the linear dependence on the gen-
Since our sensing strategy is independent across each low- erating mask patterns, this amounts to using the an identical
architecture with adjusted mask patterns
resolution frame, and also for simplicity of presentation, we
onlyconsidertherecoveryofalengthB blockofframesfrom B
(cid:88)
a single snapshot image in our theoretical analysis. hWt = Wk,thk. (24)
k=1
10
SensingArchitecture ReconstructionRMSE(%)
Subsampler MaskModel d CGInit (cid:96)2-Haar (cid:96)2-Aniso.TV (cid:96)2-Iso.TV (cid:96)2-Curvelet
Integration BS 4 43.552652 31.908468 23.632185 23.712212 34.552177
downsampling 16 52.508761 47.604093 44.602132 43.339913 43.606292
64 60.435296 58.592028 57.537579 56.960604 56.747754
US 4 43.220564 31.861358 23.497272 23.643347 34.325474
16 52.126070 47.551547 44.400403 43.069598 43.610772
64 59.642746 58.145647 57.118346 56.684747 56.541953
UP 4 38.474720 30.797415 23.265628 23.214944 32.475710
16 50.284641 46.653984 44.638167 43.091206 44.231886
64 58.815932 58.168464 57.184195 56.785933 56.558017
Random BS 4 56.890247 34.101763 18.992251 20.535913 35.202771
Summation 16 62.274874 51.430633 45.573628 44.644012 45.537800
64 65.672366 60.160050 58.296549 57.789539 57.630494
US 4 56.544494 33.659520 17.484428 19.021174 34.322115
16 61.950795 50.977833 44.857052 43.641823 45.217990
64 65.350187 59.697311 57.573941 57.159371 57.096311
UP 4 54.473107 29.831022 16.033758 17.728624 33.651808
16 60.297045 48.902002 44.939914 43.733698 45.219903
64 64.317781 60.909768 57.360741 57.106131 56.933892
Pointwise BS 4 61.578677 39.742969 17.129559 18.611897 35.793504
Subsampling 16 67.787137 65.569753 51.912236 50.842530 52.625216
64 69.244807 69.202225 64.772033 64.693113 64.741188
US 4 61.578677 39.780250 17.200487 18.703763 35.813267
16 67.787137 65.570111 51.928327 50.883713 52.629357
64 69.244807 69.202225 64.791854 64.687147 64.741653
UP 4 61.005834 36.260055 14.108310 15.883293 34.282494
16 67.675516 64.606384 48.280444 47.379639 48.924798
64 69.187676 69.189077 63.591662 63.514672 63.350203
ConventionalImager 4 38.157537 36.654837 36.067380 32.946067 33.157436
16 50.137698 50.117667 49.836071 47.932575 45.466410
64 58.873484 58.873455 58.472819 57.806412 56.764936
TABLEI
SUMMARYOFTHERMSEVALUESFORTHEDIFFERENTMASKGENERATIONANDDOWNSAMPLINGSCHEMES,INCLUDINGTHECONJUGATEGRADIENT
INITIALIZATION,ANDTHEVARIOUSRECONSTRUCTIONMETHODSCONSIDERED.THEIMAGESASSOCIATEDWITHTHEBESTPERFORMINGMETHODS
INDICATEDBYTHEBOLDRMSEVALUESARESHOWNINFIG.2
If we denote hW = (hW,...,hW), then this can be written active users, where in our system, we are using it to infer a
1 B
moresimplyashW =(WT⊗I )h=(W⊗I )Th.Hencethe sequence of sparse frames, or a sequence of sparse difference
n n
CAKE system can adapt to any temporal sparse coding of the frames.
video sequence over the block of coded frames. In summary,
to incorporate the sparse transform W applied to the frames,
B. Implementation and Normalization Details
we simply apply the adjoint transform to the independently
generated mask patterns. Recall that these RIP-satisfying sensing matrices are gen-
A useful transformation that allows the exploitation of erated using a zero-mean probability distribution and cannot
dependencies between frames is to assume that the difference be directly implemented in an optical system. As before, we
between subsequent frames are sparse. In this case, we select apply an affine transformation to each of the A so that the
k
W =∇, as defined in (12). In particular, we sense using the resulting mask elements for each h are within the range
k
coded sequence of aperture patterns h∇ where [0,1/n]. If we are generating a sensing matrix that satisfies
(cid:40) RIP with respect to the difference frames, then even a binary
h −h k =1,...,B−1,
h∇ = k k+1 generating distribution will cause the resulting mask patterns
t h k =B,
B to have elements that are neither fully open nor fully closed.
and the generating aperture patterns h are drawn from a This can be accomplished using half-toning as discussed in
suitable distribution. Sec. III-C. If such an implementation is difficult, then one
It is of interest to note that since the CAKE sensing matrix strategy would be to set these intermediate values to 0 or 1/n
is a concatenation of downsampled Toeplitz matricies, this at random with equal probability.
sensing strategy has clear connections to [48] where they
consider concatenations of Toeplitz matricies as a sensing
C. Reconstruction for Video
matrix for performing multiuser detection in wireless sensor
networks. The important conceptual link is that their sensing Given measurements from the proposed CAKE imaging
matrix is used to determine a sparse set of simultaneously system,wherewehavedesignedthemaskpatternsforsparsity