{'fraction_non_alphanumeric': 0.0492135865473103, 'fraction_numerical': 0.014724336986530579, 'mean_word_length': 4.271554575523704, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bits from Photons: Oversampled Image Acquisition Using Binary Poisson Statistics

We study a new image sensor that is reminiscent of traditional photographic film. Each pixel in the sensor has a binary response, giving only a one-bit quantized measurement of the local light intensity. To analyze its performance, we formulate the oversampled binary sensing scheme as a parameter estimation problem based on quantized Poisson statistics. We show that, with a single-photon quantization threshold and large oversampling factors, the Cram\'er-Rao lower bound (CRLB) of the estimation variance approaches that of an ideal unquantized sensor, that is, as if there were no quantization in the sensor measurements. Furthermore, the CRLB is shown to be asymptotically achievable by the maximum likelihood estimator (MLE). By showing that the log-likelihood function of our problem is concave, we guarantee the global optimality of iterative algorithms in finding the MLE. Numerical results on both synthetic data and images taken by a prototype sensor verify our theoretical analysis and demonstrate the effectiveness of our image reconstruction algorithm. They also suggest the potential application of the oversampled binary sensing scheme in high dynamic range photography.


I. INTRODUCTION
Before the advent of digital image sensors, photography, for the most part of its history, used film to record light information. At the heart of every photographic film are a large number of light-sensitive grains of silver-halide crystals [1]. During exposure, each micron-sized grain has a binary fate: Either it is struck by some incident photons and becomes "exposed", or it is missed by the photon bombardment and remains "unexposed". In the subsequent film development process, exposed grains, due to their altered chemical properties, are converted to silver metal, contributing to opaque spots on the film; unexposed grains are washed away in a chemical bath, leaving behind them transparent regions on the film. Thus, in essence, photographic film is a binary imaging medium, using local densities of opaque silver grains to encode the original light intensity information. Thanks to the small size and large number of these grains, one hardly notices this quantized nature of film when viewing it at a distance, observing only a continuous gray tone.
In this work, we study a new digital image sensor that is reminiscent of photographic film. Each pixel in the sensor has a binary response, giving only a one-bit quantized measurement of the local light intensity.
At the start of the exposure period, all pixels are set to 0. A pixel is then set to 1 if the number of photons reaching it during the exposure is at least equal to a given threshold q. One way to build such binary sensors is to modify standard memory chip technology, where each memory bit cell is designed to be sensitive to visible light [2]. With current CMOS technology, the level of integration of such systems can exceed 10 9 ∼ 10 10 (i.e., 1 giga to 10 giga) pixels per chip. In this case, the corresponding pixel sizes (around 50 nm [3]) are far below the diffraction limit of light (see Section II for more details), and thus the image sensor is oversampling the optical resolution of the light field. Intuitively, one can exploit this spatial redundancy to compensate for the information loss due to one-bit quantizations, as is classic in oversampled analog-to-digital (A/D) conversions [4]- [7].
Building a binary sensor that emulates the photographic film process was first envisioned by Fossum [8], who coined the name "digital film sensor". The original motivation was mainly out of technical necessity.
The miniaturization of camera systems calls for the continuous shrinking of pixel sizes. At a certain point, however, the limited full-well capacity (i.e., the maximum photon-electrons a pixel can hold) of small pixels becomes a bottleneck, yielding very low signal-to-noise ratios (SNRs) and poor dynamic ranges. In contrast, a binary sensor whose pixels only need to detect a few photon-electrons around a small threshold q has much less requirement for full-well capacities, allowing pixel sizes to shrink further.
In this paper, we present a theoretical analysis of the performance of the binary image sensor, and propose an efficient and optimal algorithm to reconstruct images from the binary sensor measurements.
Our analysis and numerical simulations demonstrate that the dynamic ranges of the binary sensors can be orders of magnitude higher than those of conventional image sensors, thus providing one more motivation for considering this binary sensing scheme.
Since photon arrivals at each pixel can be well-approximated by a Poisson random process whose rate is determined by the local light intensity, we formulate the binary sensing and subsequent image reconstruction as a parameter estimation problem based on quantized Poisson statistics. Image estimation from Poisson statistics has been extensively studied in the past, with applications in biomedical and astrophysical imaging.
Previous work in the literature has used linear models [9], multiscale models [10], [11], and nonlinear piecewise smooth models [12], [13] to describe the underlying images, leading to different (penalized) maximum likelihood and/or Bayesian reconstruction algorithms. The main difference between our work and previous works is that we have only access to one-bit quantized Poisson statistics. The binary quantization and spatial oversampling in the sensing scheme add interesting dimensions to the original problem. As we will show in Section III, the performance of the binary sensor depends on the intricate interplay of three parameters: the average light intensity, the quantization threshold q, and the oversampling factor.
The binary sensing scheme studied in this paper also bears resemblance to oversampled A/D conversion schemes with quantizations (see, e.g., [4]- [7]). Previous work on one-bit A/D conversions considers bandlimited signals or, in general, signals living in the range space of some overcomplete representations. The effect of quantization is often approximated by additive noise, which is then mitigated through noise shaping [4], [6], or dithering [7], followed by linear reconstruction. In our work, the binary sensor measurements are modeled as one-bit quantized versions of correlated Poisson random variables (instead of deterministic signals), and we directly solve the statistical inverse problem by using maximum likelihood estimation, without any additive noise approximation.
The rest of the paper is organized as follows. After a precise description of the binary sensing model in Section II, we present three main contributions in this paper: 1. Estimation performance: In Section III, we analyze the performance of the proposed binary sensor in estimating a piecewise constant light intensity function. In what might be viewed as a surprising result, we show that, when the quantization threshold q = 1 and with large oversampling factors, the Cramér-Rao lower bound (CRLB) [14] of the estimation variance approaches that of unquantized Poisson intensity estimation, that is, as if there were no quantization in the sensor measurements. Furthermore, the CRLB can be asymptotically achieved by a maximum likelihood estimator (MLE), for large oversampling factors.
Combined, these two results establish the feasibility of trading spatial resolutions for higher quantization bit depth.

Advantage over traditional sensors:
We compare the oversampled binary sensing scheme with traditional image sensors in Section III-C. Our analysis shows that, with sufficiently large oversampling factors, the new binary sensor can have higher dynamic ranges, making it particularly attractive in acquiring scenes containing both bright and dark regions.
3. Image reconstruction: Section IV presents an MLE-based algorithm to reconstruct the light intensity field from the binary sensor measurements. As an important result in this work, we show that the loglikelihood function in our problem is always concave for arbitrary linear field models, thus ensuring the achievement of global optimal solutions by iterative algorithms. For numerically solving the MLE, we present a gradient method, and derive efficient implementations based on fast signal processing algorithms in the polyphase domain [15], [16]. This attention to computational efficiency is important in practice, due to extremely large spatial resolutions of the binary sensors.
Section V presents numerical results on both synthetic data and images taken by a prototype device [17].
These results verify our theoretical analysis on the binary sensing scheme, demonstrate the effectiveness of our image reconstruction algorithm, and showcase the benefit of using the new binary sensor in acquiring scenes with high dynamic ranges.
To simplify the presentation we base our discussions on a one-dimensional (1-D) sensor array, but all the results can be easily extended to the 2-D case.

A. Diffraction Limit and Linear Light Field Models
In this section, we describe the binary imaging scheme studied in this work. Consider a simplified camera model shown in Fig. 1(a). We denote by λ 0 (x) the incoming light intensity field (i.e., the radiance map).
By assuming that light intensities remain constant within a short exposure period, we model the field as only a function of the spatial variable x. Without loss of generality, we assume that the dimension of the sensor array is of one spatial unit, i.e., 0 ≤ x ≤ 1.
After passing through the optical system, the original light field λ 0 (x) gets filtered by the lens, which acts like a linear system with a given impulse response. Due to imperfections (e.g., aberrations) in the lens, the impulse response, a.k.a. the point spread function (PSF) of the optical system, cannot be a Dirac delta, thus, imposing a limit on the resolution of the observable light field. However, a more fundamental physical limit is due to light diffraction [18]. As a result, even if the lens is ideal, the PSF is still unavoidably a small blurry spot [see, for example, Fig. 1(b)]. In optics, such diffraction-limited spot is often called the Airy disk [18], whose radius R a can be computed as where w is the wavelength of the light and f is the F-number of the optical system. Example 1: At wavelength w = 420 nm (i.e., for blue visible light) and f = 2.8, the radius of the Airy disk is 1.43 µm. Two objects with distance smaller than R a cannot be clearly separated by the imaging system as their Airy disks on the image sensor start blurring together. Current CMOS technology can already make standard pixels smaller than R a , reaching sizes ranging from 0.5 µm to 0.7 µm [19]. In the case of binary sensors, the simplicity of each pixel allows the feature size to be further reduced. For example, based on standard memory technology, each memory bit-cell (i.e., pixel) can have sizes around 50 nm [3], making it possible to substantially oversample the light field.
In what follows, we denote by λ(x) the diffraction-limited (i.e., "observable") light intensity field, which is the outcome of passing the original light field λ 0 (x) through the lens. Due to the lowpass (smoothing) nature of the PSF, the resulting λ(x) has a finite spatial-resolution, i.e., it has a finite number of degrees of freedom per unit space.

Definition 1 (Linear field model):
In this work, we model the diffraction-limited light intensity field as where ϕ(x) is a nonnegative interpolation kernel, N is a given integer, τ is the exposure time, and {c n : c n ≥ 0} is a set of free variables.

Remark 1:
The constant N/τ in front of the summation is not essential, but its inclusion here leads to simpler expressions in our later analysis.
The function λ(x) as defined in (1) has N degrees of freedom. To guarantee that the resulting light fields are physically meaningful, we require both the interpolation kernel ϕ(x) and the expansion coefficients {c n } to be nonnegative. Some examples of the interpolation kernels ϕ(x) include the box function, cardinal B-splines [20], and the squared sinc function,

B. Sampling the Light Intensity Field
The image sensor in Fig. 1 where β(x) is the box function defined in (2) and ·, · represents the standard L 2 -inner product. Substitute the light field model (1) into the above equality, where (5) is obtained through a change of variables (N x − n) → x.

Definition 2:
The spatial oversampling factor, denoted by K, is the ratio between the number of pixels per unit space and the number of degrees of freedom needed to specify the light field λ(x) in (1), i.e., In this work, we are interested in the "oversampled" case where K > 1. Furthermore, we assume that K is an integer for simplicity of notation. Using (6), and by introducing a discrete filter

Photon Counting Quantization
Binary sensing Light exposure Fig. 2. The signal processing block diagram of the imaging model studied in this paper. In the first step, the light exposure value sm at the mth pixel is related to the expansion coefficients cn through a concatenation of upsampling and filtering operations.
Subsequently, the image sensor converts {sm} into quantized measurements {bm} (see Fig. 3 and the discussions in Section II-C for details of this second step).
we can simplify (5) as The above equality specifies a simple linear mapping from the expansion coefficients {c n } of the light field to the light exposure values {s m } accumulated by the image sensor. Readers familiar with multirate signal processing [15], [16] will immediately recognize that the relation in (8) can be implemented via a concatenation of upsampling and filtering, as shown in the left part of Fig. 2. This observation can also be verified by expressing (8) in the z-transform domain and using the fact that C(z K ) is the z-transform of the K-fold upsampled version of c n . In Section IV, we will further study the signal processing block diagram in Fig. 2 to derive efficient implementations of the proposed image reconstruction algorithm.

Example 2:
The discrete filter g m is completely specified by the interpolation kernel ϕ(x) and the oversampling factor K. As a simple case, when the kernel ϕ(x) = β(x), we can compute from (7) that of a Poisson random variable Y m , whose intensity parameter is equal to s m , i.e.,

C. Binary Sensing and One-Bit Poisson Statistics
It is a well-known property of the Poisson process that E[Y m ] = s m . Thus, the average number of photons captured by a given pixel is equal to the local light exposure s m .
As a photosensitive device, each pixel in the image sensor converts photons to electrical signals, whose amplitude is proportional to the number of photons impinging on that pixel. 1 In a conventional sensor design, the analog electrical signals are then quantized by an A/D converter into 8 to 14 bits (usually the more bits the better). In this work, we study a new sensor design using the following binary (i.e., one-bit) quantization scheme.

Definition 3 (Binary Quantization):
Let q ≥ 1 be an integer threshold. A binary quantizer is a mapping In we can write Remark 2: The noise model considered in this paper is that of Poisson noise. In practice, the performance of image sensors is also influenced by thermal noise, which in our case can be modeled as random bitflipping in the binary sensor measurements. Due to space constraints, we leave further discussions on this additional noise source and its impact on reconstruction performance to a follow-up work.

D. Multiple Exposures and Temporal Oversampling
Our previous discussions focus on the case of acquiring a single frame of quantized measurements during the exposure time [0, τ ]. As an extension, we can consider multiple exposures and acquire J consecutive and independent frames. The exposure time for each frame is set to τ /J, so that the total acquisition time remains the same as the single exposure case. In what follows, we call J the temporal oversampling factor.
As before, we assume that τ ≪ 1 and thus light intensities λ(x) stay constant within the entire acquisition For the jth frame (0 ≤ j < J), we denote by s j,m the light exposure at the mth pixel. Following the same derivations as in Section II-B, we can show that where {c n } are the expansion coefficients of the light field λ(x), and g m is the discrete filter defined in (7). The only difference between (14) and (8) is the extra factor of 1/J, due to the change of exposure time from τ to τ /J. In the z-domain, similar to (9), In what follows, we establish an equivalence between temporal oversampling and spatial oversampling.
More precisely, we will show that an M -pixel sensor taking J independent exposures (i.e., with J-times In multirate signal processing, the above construction is called the polyphase representation [15], [16], and its alternating subsequences {s j,m } J−1 j=0 the polyphase components. Proposition 1: Let g m be a filter whose z-transform where G(z) is the z-transform of the filter g m defined in (7). Then, Proof: See Appendix A.
Remark 3: Proposition 1 formally establishes the equivalence between spatial and temporal oversampling.
We note that (18)

Remark 4:
Taking multiple exposures is a very effective way to increase the total oversampling factor of the binary sensing scheme. The key assumption in our analysis is that, during the J consecutive exposures, the light field remains constant over time. To make sure this assumption holds for arbitrary values of J, we set the exposure time for each frame to τ /J, for a fixed and small τ . Consequently, the maximum temporal oversampling factor we can achieve in practice will be limited by the readout speed of the binary sensor.
Thanks to the equivalence between spatial and temporal oversampling, we only need to focus on the single exposure case in our following discussions on the performance of the binary sensor and image reconstruction algorithms. All the results we obtain extend directly to the multiple exposure case.

III. PERFORMANCE ANALYSIS
In this section, we study the performance of the binary image sensor in estimating light intensity information, analyze the influence of the quantization threshold and oversampling factors, and demonstrate the new sensor's advantage over traditional sensors in terms of higher dynamic ranges. In our analysis, we assume that the light field is piecewise constant, i.e., the interpolation kernel ϕ(x) in (1) is the box function β(x). This simplifying assumption allows us to derive closed-form expressions for several important performance measures of interest. Numerical results in Section V suggest that the results and conclusions we obtain in this section applies to the general linear field model in (1) with different interpolation kernels.

A. The Cramér-Rao Lower Bound (CRLB) of Estimation Variances
From Definition 1, reconstructing the light intensity field λ(x) boils down to estimating the unknown deterministic parameters {c n }. Input to our estimation problem is a sequence of binary sensor measurements Assume that the light field λ(x) is piecewise constant. We have computed in Example 2 that, under this case, the discrete filter g m used in (8) We In what follows, and without loss of generality, we focus on estimating c 0 from the block of binary For notational simplicity, we will drop the subscript in c 0 and use c instead. To analyze the performance of the binary sensing scheme, we first compute the CRLB [14], which provides a theoretical lower bound on the variance of any unbiased estimator.
Denote by L b (c) the likelihood function of observing K binary sensor measurement b. Then, where (20) is due to the independence of the photon counting processes at different pixel locations, and (21) follows from (13) and (19). Defining K 1 (0 ≤ K 1 < K) to be the number of "1"s in the binary sequence b, we can simplify (21) as Proposition 2: The CRLB of estimating the light intensity c from K binary sensor measurements with threshold q ≥ 1 is Proof: See Appendix B.
It will be interesting to compare the performance of our binary image sensor with that of an ideal sensor which does not use quantization at all. To that end, consider the same situation as before, where we use K pixels to observe a constant light intensity value c. The light exposure s m at each pixel is equal to c/K, as in (19). Now, unlike the binary sensor which only takes one-bit measurements, consider an ideal sensor that can perfectly record the number of photon arrivals at each pixel. By referring to Fig. 3, the sensor measurements in this case will be {y m }, whose probability distributions are given in (11).
In Appendix C, we compute the CRLB of this unquantized sensing scheme as which is natural and reflects the fact that the variance of a Poisson random variable is equal to its mean (i.e., c, in our case).
To be sure, we always have CRLB bin (K, q) > CRLB ideal (K), for arbitrary oversampling factor K and quantization threshold q. This is not surprising, as we lose information by one-bit quantizations. In practice, the ratio between the two CRLBs provides a measure of performance degradations incurred by the binary sensors. What is surprising is that the two quantities can be made arbitrarily close, when q = 1 and K is large, as shown by the following proposition.
Proof: Specializing the expression (23) for q = 1, we get and thus (25). The statements for cases when q ≥ 2 are shown in Appendix D.
Proposition 3 indicates that it is feasible to use oversampling to compensate for information loss due to binary quantizations. It follows from (25) that, with large oversampling factors, the binary sensor operates as if there were no quantization in its measurements. It is also important to note that this desirable tradeoff between spatial resolution and estimation variance only works for a single-photon threshold (i.e., q = 1). For other choices of the quantization threshold, the "gap" between CRLB bin (K, q) and CRLB ideal (K), measured in terms of their ratio, cannot be made arbitrarily small, as shown in (26). In fact, it quickly tends to infinity as the oversampling factor K increases.
The results of Proposition 3 can be intuitively understood as follows: The expected number of photons collected by each pixel during light exposure is equal to s m = c/K. As the oversampling factor K goes to infinity, the mean value of the Poisson distribution tends to zero. Consequently, most pixels on the sensor will only get zero or one photon, with the probability of receiving two or more photons at a pixel close to zero. In this case, with high probability, a binary quantization scheme with threshold q = 1 does not lose information. In contrast, if q ≥ 2, the binary sensor measurements will be almost uniformly zero, making it nearly impossible to differentiate between different light intensities.

B. Asymptotic Achievability of the CRLB
In what follows, we show that, when q = 1, the CRLB derived in (23) can be asymptotically achieved by a simple maximum likelihood estimator (MLE). Given a sequence of K binary measurements b, the MLE we seek is the parameter that maximizes the likelihood function L b (c) in (22). More specifically, where we substitute p 1 (c/K) in (22) where p is the inverse function of p 0 (x). Remark 5: From the definition in (12), we can easily verify that d dx p 0 (x) < 0 for all x > 0. It follows that the function p 0 (x) is strictly decreasing for x > 0 and that the inverse p . In this particular case, and for It follows that we can use the sum of the K binary measurements as a first-order approximation of the light intensity estimation.
Proof: At the two extreme cases, when K 1 = 0 or K 1 = K, it is easy to see that (28) is indeed the solution to (27). Next, we assume that 0 < K 1 < K.
Computing the derivative of L b (c) and setting it to zero, we can verify that the equation d dc L b (c) = 0 has a single solution at Since L b (c) ≥ 0 and L b (0) = lim c→∞ L b (c) = 0, we conclude that the likelihood function L b (c) achieves its maximum value at c max . Finally, the MLE solution c ML = min { c max , S}, and thus, we have (28).
Theorem 1: When q = 1, we have where |ε 1 | ≤ 2c e 1−c ec Meanwhile, the mean squared error (MSE) of the estimator approaches CRLB ideal , i.e., where |ε 2 | ≤ 2c(c + 1)e −c ec Remark 6: It is easy to verify that, for fixed c, the two terms ε 1 and ε 2 converge (very quickly) to 0 as S tends to infinity. It then follows from (29) and (30) that the MLE is asymptotically unbiased and efficient, in the sense that We leave the formal proof of this theorem to Appendix E. Its main idea can be summarized as follows. As K goes to infinity, the area of each pixel tends to zero, so does the average number of photons arriving at that pixel. As a result, most pixels on the sensor will get only zero or one photon during exposure. A single-photon binary quantization scheme can record perfectly the pattens of "0"s and "1"s on the sensor.
It loses information only when a pixel receives two or more photons, but the probability of such events tends to zero as K increases.
Suppose, now, that we use a quantization threshold q ≥ 2. In this case, as K tends to infinity, the binary responses of different pixels will almost always be "0", essentially obfuscating the actual light intensity values. This problem leads to poor performance in the MLE. As stated in the following proposition, the asymptotic MSE for q ≥ 2 becomes c 2 instead of c.
Proposition 4: When q ≥ 2, the MLE is asymptotically biased, that is, for any fixed c and S, Meanwhile, the MSE becomes Proof: See Appendix F.

C. Advantages over Traditional Sensors
In what follows, we demonstrate the advantage of the oversampled binary sensing scheme, denoted by "BIN", in achieving higher dynamic ranges. We focus on the case where the quantization threshold is set to q = 1. For comparisons, we also consider the following two alternative sensing schemes: The first, denoted by "IDEAL", uses a single pixel to estimate the light exposure parameter (i.e., nonoversampled), but that pixel can record perfectly the number of photon arrivals during exposure; The second scheme, denoted by "SAT", is very similar to the first, with the addition of a saturation point C max , beyond which the pixel can hold no more photons. Note that in our discussions, the "SAT" scheme serves as an idealized model of conventional image sensors, for which the saturation is caused by the limited full-well capacity of the semiconductor device. The general trend of conventional image sensor design has been to pack more pixels per chip by reducing pixel sizes, leading to lower full-well capacities and thus lower saturation values.
where c is the estimation of the light exposure value we obtain from each of the sensing schemes.
We observe that the "IDEAL" scheme (the red dash-dot line in the figure) represents an upper-bound of the estimation performance. To see this, denote by y the number of photons that arrive at the pixel during exposure. Then y is a realization of a Poisson random variable Y with intensity equal to the light exposure value c, i.e., P(Y = y; c) = c y e −c y! .
Maximizing this function over c, we can compute the MLE for the "IDEAL" scheme as c IDEAL (y) = y.
It is easy to verify that this estimator is unbiased, i. which is the full well capacity of the image sensor reported in [21]. The sensor measurement in this case is y SAT def = min{y, C max }, and the estimator we use is We can see that the "SAT" scheme initially has the same performance as "IDEAL". It remains this way until the light exposure value c approaches the saturation point C max , after which there is a drastic drop 2 in SNR. Denoting by SNR min the minimum acceptable SNR in a given application, we can then define the dynamic range of a sensor as the range of c for which the sensor achieves at least SNR min . For example, if we choose SNR min = 20 dB, then, as shown in the figure, the SAT scheme has a dynamic range from c = 10 2 to c ≈ 10 4 , or, if measured in terms of ratios, 100 : 1.
Finally, the three dashed lines represent the "BIN" scheme with q = 1 and increasing oversampling factors (from left to right: K = 2 13 , 2 14 , 2 15 and 2 16 , respectively). We use the MLE given in (28) and plot the corresponding estimation SNRs. We see that, within a large range of c, the performance of the "BIN" scheme is very close to that of the "IDEAL" scheme that does not use quantization. This verifies our analysis in Theorem 1, which states that the "BIN" scheme with a single-photon threshold can approach the ideal unquantized CRLB when the oversampling factor is large enough. Furthermore, when compared with the "SAT" scheme, the "BIN" scheme has a more gradual decrease in SNR when the light exposure values increase, and has a higher dynamic range. For example, when K = 2 16 , the dynamic range of the "BIN" scheme spans from c = 10 2 to c = 10 5.8 , about two orders of magnitude higher than that of "SAT".
In Section V, we will present a numerical experiment that points to a potential application of the binary sensor in high dynamic range photography. 2 The estimator in (33) is biased around c = Cmax. For a very narrow range of light intensity values centered around Cmax, the MSE of this biased estimator is lower than the ideal CRLB. Thus, there is actually a short "spike" in SNR right before the drop.

Remark 7:
Note that K is the product of the spatial oversampling factor and the temporal oversampling factor. For example, the pixel pitch of the image sensor reported in [21] is 1.65µm. If the binary sensor is built on memory chip technology, with a pitch size of 50 nm [3], then the maximum spatial oversampling factor is about 1089. To achieve K = 2 13 , 2 14 , 2 15 and 2 16 , respectively, as required in Fig. 4, we then need to have temporal oversampling factors ranging from 8 to 60. Unlike traditional sensors which require multi-bit quantizers, the binary sensors only need one-bit comparators. This simplicity in hardware can potentially lead to faster readout speeds, making it practical to apply temporal oversampling.

IV. OPTIMAL IMAGE RECONSTRUCTION AND EFFICIENT IMPLEMENTATIONS
In the previous section, we studied the performance of the binary image sensor, and derived the MLE for where e m is the mth standard Euclidean basis vector. 3 Remark 8: In using the above notations, we do not distinguish between single exposure and multiple exposures, whose equivalence has been established by Proposition 1 in Section II-D. In the case of multiple exposures, the essential structure of G-upsampling followed by filtering-remains the same. All we need 3 Here we use zero-based indexing. Thus, e0 to do is to replace s by the interlaced sequence { s m } constructed in (16), the oversampling factor K by KJ, and the filter g m by g m in (17). (20) and (21), the likelihood function given M binary measurements

Similar to our derivations in
where (36) follows from (12) and (35). In our subsequent discussions, it is more convenient to work with the log-likelihood function, defined as For any given observation b, the MLE we seek is the parameter that maximizes L b (c), or equivalently, The constraint c ∈ [0, S] N means that every parameter c n should satisfy 0 ≤ c n ≤ S, for some preset maximum value S.
Example 3: As discussed in Section III, when the light field is piecewise-constant, different light field parameters {c n } can be estimated independently. In that case, the likelihood function has only one variable [see (22)] and can be easily visualized. In Fig. 5, we plot L b (c) in (22)  show that this result is general, namely, the log-likelihood functions in the form of (37) are always concave.
Proof: It follows from the definition in (12) that, for any b m ∈ {0, 1}, the function log p bm (s) is either We can apply Lemma 2 in both cases, and show that {log p bm (s)} are concave functions for s ≥ 0. Since the sum of concave functions is still concave and the composition of a concave function with a linear mapping (s m = e T m Gc) is still concave, we conclude that the log-likelihood function defined in (37) is concave. In general, there is no closed-form solution to the maximization problem in (38). An MLE solution has to be found through numerical algorithms. Theorem 2 guarantees the global convergence of these iterative numerical methods.

B. Iterative Algorithm and Efficient Implementations
We compute the numerical solution of the MLE by using a standard gradient ascent method. Denote by c (k) the estimation of the unknown parameter c at the kth step. The estimation c (k+1) at the next step is obtained by where ▽l b (c (k) ) is the gradient of the log-likelihood function evaluated at c (k) , γ k is the step-size at the current iteration, and P D is the projection onto the search domain D Taking the derivative of the log-likelihood function l b (c) in (37), we can compute the gradient as The choice of the step size γ k has significant influence over the speed of convergence of the above iterative algorithm. We follow [9] by choosing, at each step, a γ k so that the gradient vectors at the current and the next iterations are approximately orthogonal to each other. By assuming that the estimates s (k+1) and s (k) at consecutive iterations are close to each other, we can use the following first-order approximation

It follows that
Assuming that the gradient update c (k) + γ k ▽l b (c (k) ) is inside of the constraint set D, we can neglect the projection operator P D in (40), and write Substituting the above equality into (42), we get Finally, by requiring that ▽l b (c (k+1) ) be orthogonal to ▽l b (c (k) ), we compute the optimal step size as Remark 9: By definition, H b (s) (for b = 0, 1) are the second-order derivatives of concave functions (see Lemma 2), and are thus nonpositive. Consequently, the terms −H b (s) in the denominator of (43) are At every iteration of the gradient algorithm, we need to update the gradient and the step size γ k . We see from (41) and (43) that the computations always involve matrix-vector products in the form of Ga and G T b, for some vectors a, b. The matrix G is of size M × N , where M is the total number of pixels. In practice, M will be in the range of 10 9 ∼ 10 10 (i.e., gigapixel per chip), making it impossible to directly implement the matrix operations. Fortunately, the matrix G used in both formulae is highly structured, and can be implemented as upsampling followed by filtering (see our discussions in Section II-B and the expression (8) for details). Similarly, the transpose G T can be implemented by filtering (by g −m ) followed by downsampling, essentially "flipping" all the operations in G. Fig. 6(a) and Fig. 6(b) summarizes these operations.
We note that the implementations illustrated in Fig. 6(a) and Fig. 6(b) are not yet optimized: For example, the input to the filter g m in Fig. 6(a) is an upsampled sequence, containing mostly zero elements; In Fig. 6(b), we compute a full filtering operation (by g −m ), only to discard most of the filtering results in the subsequent downsampling step. All these deficiencies can be eliminated by using the tool of polyphase representations from multirate signal processing [15], [16], as follows.
First, we split the filter g m into K non-overlapping polyphase components g 0,m , g 1,m , . . . , g K−1,m , defined as g k,m = g Km+k , for 0 ≤ k < K.
Intuitively, the polyphase components specified in (44) are simply downsampled versions of the original filter g m , with the sampling locations of all these polyphase components forming a complete partition. The mapping between the filter g m and its polyphase components is one-to-one. To reconstruct g m , we can easily verify that, in the z-domain, Following the same steps as above, we can also split the sequences u m and b m in Fig. 6 and  Fig. 6(b). Recall that K is the oversampling factor of our image sensor. As we operate in highly oversampled regimes (e.g., K = 1024) to compensate for information loss due to one-bit quantizations, the above improvements make our algorithms orders of magnitude faster.

V. NUMERICAL RESULTS
We present several numerical results in this section to verify our theoretical analysis and the effectiveness of the proposed image reconstruction algorithm.

A. 1-D Synthetic Signals
Consider a 1-D light field λ(x) shown in Fig. 7(a). The interpolation filter ϕ(x) we use is the cubic B-spline function β 3 (x) defined in (3). We can see that λ(x) is a linear combination of the shifted kernels, with the expansion coefficients {c n } shown as blue dots in the figure.
We simulate a binary sensor with threshold q = 1 and oversampling factor K = 256. Applying the proposed MLE-based algorithm in Section IV, we obtain a reconstructed light field (the red dashed curve) shown in Fig. 7(b), together with the original "ground truth" (the blue solid curve). We observe that the low-light regions are well-reconstructed but there exist large "overshoots" in the high-light regions.
We can substantially improve the reconstruction quality by increasing the oversampling factor of the sensor. Fig. 7(c) shows the result obtained by increasing the spatial oversampling factor to K = 2048.
Alternatively, we show in Fig. 7(d)  This observation agrees with our previous theoretical analysis in Section II-D on the equivalence between spatial and temporal oversampling schemes.

B. Acquiring Scenes with High Dynamic Ranges
A well-known difficulty in photography is the limited dynamic ranges of the image sensors. Capturing both very bright and very dark regions faithfully in a single image is difficult. For example, Fig. 8(a) shows several images taken inside of a church with different exposure times [22]. The scene contains both sun-lit areas and shadow regions, with the former over a thousand times brighter than the latter. Such high dynamic ranges are well-beyond the capabilities of conventional image sensors. As a result, these images are either overexposed or underexposed, with no single image rendering details in both areas. In light of this problem, an active area of research in computational photography is to reconstruct a high dynamic range radiance map by combining multiple images with different exposure settings (see, e.g., [22], [23]). While producing successful results, such multi-exposure approaches can be time-consuming.
In Section III-C, we have shown that the binary sensor studied in this work can achieve higher dynamic ranges than conventional image sensors. To demonstrate this advantage, we use the high dynamic range radiance map obtained in [22] as the ground truth data [i.e., the light field λ(x) as defined in (1)], and simulate the acquisition of this scene by using a binary sensor with a single photon threshold. The spatial oversampling factor of the binary sensor is set to 32 × 32, and the temporal oversampling factor is 256 (i.e., 256 independent frames). Similar to our previous experiment on 1-D signals, we use a cubic B-spline kernel [i.e., ϕ(x) = β 3 (x)] along each of the spatial dimensions. Fig. 8(b) shows the reconstructed radiance map using our algorithm described in Section IV. Since the radiance map has a dynamic range of 3.3 × 10 5 : 1, the image is shown in logarithmic scale. To have a visually more pleasing result, we also shown in Fig. 8(c) a tone-mapped [23] version of the reconstruction. We can see from Fig. 8(b) and Fig. 8(c) that details in both light and shadow regions have been faithfully preserved in the reconstructed radiance map, suggesting the potential application of the binary sensor in high dynamic range photography.

C. Results on Real Sensor Data
We have also applied our reconstruction algorithm to images taken by an experimental sensor based on single photon avalanche diodes (SPADs) [17]. The sensor has binary-valued pixels with single-photon sensitivities, i.e., the quantization threshold is q = 1. Due to its experimental nature, the sensor has limited spatial resolution, containing an array of only 32×32 detectors. To emulate the effect of spatial oversampling, we apply temporal oversampling and acquire 4096 independent binary frames of a static scene. In this case, we can estimate the light intensity at each pixel independently by using the closed-form MLE solution in (28). Fig. 9 shows 50 such binary images, together with the final reconstruction result (at the lower-right corner). The quality of reconstruction verifies our theoretical model and analysis.

VI. CONCLUSIONS
We presented a theoretical study of a new image sensor that acquires light information using one-bit pixels-a scheme reminiscent of traditional photographic film. By formulating the binary sensing scheme By substituting (15) into (48) and using the definition (17), we can simplify (48) as Finally, since C(z KJ ) is the z-transform of the sequence n c n δ m−KJn , it follows from (49) that s m = ( n c n δ m−KJn ) * g m , and thus (18).

B. The CRLB of Binary Sensors
We first compute the Fisher information, defined as Using (22), we get where p ′ 0 (x) = d dx p 0 (x) and p ′′ 0 (x) = d 2 dx 2 p 0 (x) are the first and second order derivative of p 0 (x), respectively. In reaching (50), we have also used the fact that p 1 (x) = 1 − p 0 (x) and thus p ′ 0 (x) = − d dx p 1 (x) and p ′′ 0 (x) = − d 2 dx 2 p 1 (x). Note that K 1 = 0≤m<K b m is a binomial random variable, and thus its mean can be computed as On substituting the above expression into (50), the Fisher information can be simplified as Using the definition of p 0 (x) in (12), the derivative in the numerator of (51) can be computed as Finally, since CRLB bin (K, q) = 1/I(c), we reach (23) by substituting (12) and (52) into (51), and after some straightforward manipulations.

C. The CRLB of Ideal Unquantized Sensors
Without quantization, the sensor measurements are y def = [y 0 , y 1 , . . . , y K−1 ] T , i.e., the number of photons collected at each pixel. The likelihood function in this case is where (53) follows from the independence of {Y m } and the expressions (11) and (19).
Computing the Fisher information I(c) = E[− ∂ 2 ∂c 2 log L y (c)] in this case, we get where in reaching (56) we have used the Cauchy-Schwarz inequality.
It is easy to verify (through a change of variables t → x 2 t ) that Consequently, the term on the right-hand side of (56) is equal to xe 2x 1 . Through a change of variables t → t 2 , we have, where (57) is obtained through another change of variables t → u 2 + x + u 2 4 . We can easily verify that (57) is a monotonically increasing function with respect to x > 0. So, for For 0 ≤ x < 1, we can obtain the following inequalities by keeping the first two terms in p 0 (x) and the first term in p 1 (x): It is easy to see that 2(q−1)! q is a monotonically increasing function with respect to q ≥ 2. Therefore, Finally, for q = 2, we keep the first two terms in p 0 (x) and p 1 (x) and get Combining (58), (59) and (60), we reach the desired result in (26).
Finally, we show that, for q ≥ 2, the "gap" between CRLB bin (K, q) and CRLB ideal (K) will only get bigger as the oversampling factor K grows. To that end, we notice that when K → ∞, the variable x = c/K → 0.
Keeping the first terms in p 0 (x) and p 1 (x), we have For q ≥ 2, the above quantity goes to infinity as x → 0. Therefore, lim K→∞ CRLB bin (K, q)/CRLB ideal (K) = ∞.

E. Proof of Theorem 1
When q = 1, we have p 0 (x) = e −x and thus p 0 (x) = − log(x). In this case, the MLE solution in (28) can be rewritten as We note that −K log . . and that lim K→∞ K(1 − e −S/K ) = S. Thus, for sufficiently large K, the above MLE solution can be further simplified as In what follows, we derive bounds for the quantities on the right-hand side of the above inequality. First, consider the probability P(K 1 = n). Since K 1 is a binomial random variable, we have For every n < S, it is easy to verify that K(K−1)...(K−n+1) Next, consider the second term on the right-hand side of (62).
for all c < S, where (66) follows from (64) and the inequality (67) is due to the Chernoff bound on the tail of Poisson distributions [24]. Similarly, the third term on the right-hand side of (62) can be rewritten as where the inequality is again an application of the Chernoff bound. Finally, on substituting (65), (67) and (68) into (62), and after some simple manipulations, we reach (29).
The proof for the mean-squared error formula (30) is similar. Using (61), we have where in reaching (69), we have used the estimation (64) of the Binomial probabilities. We note that the variance of a Poisson random variable is equal to its mean. Thus, c = ∞ n=0 (n − c) 2 c n e −c n! . On combining this identity with (69), Applying the Chernoff bound to the above inequality, we get (30).
For q ≥ 2 and any fixed constant S, the above quantity converges to 0 as K tends to infinity. As a result, for sufficiently large K, the MLE solution in (28) can be simplified as where we have also used the fact that p 0 (0) = 1 and thus p We have K 1 = 0 when all the pixel responses are uniformly 0. The probability of seeing such an event is P(K 1 = 0) = p 0 (c/K) K = e −c 1 + c K + . . . + c q−1 K q−1 (q − 1)! K , which converges to 1 as K tends to infinity, i.e., lim K→∞ P(K 1 = 0) = 1.

G. Proof of Lemma 2
The function h(x) def = log j k=i x k e −x k! is continuously differentiable on the interval (0, ∞). Therefore, to establish its concavity, we just need to show that its second derivative is nonpositive. To that end, we first introduce a sequence of functions {r k (x)} k∈Z∪{∞} , defined as It is straightforward to verify that d dx r k (x) = r k−1 (x) for all k ∈ Z ∪ {∞}. Now, rewriting h(x) as log j k=i r k (x) − x and computing its second derivative, we get where we have omitted the function argument x in r k (x), r k−1 (x) and r k−2 (x) for notational simplicity.
Recall that our goal is to show that d 2 dx 2 h(x) ≤ 0, for x > 0. Since the denominator of (74) is always positive, we just need to focus on its numerator. Using the identities i≤k≤j r k = i≤k≤j r k−1 + r j − r i and i≤k≤j r k−1 = i≤k≤j r k−2 + r j−1 − r i−1 , we can simplify the numerator of (74) as follows: In what follows, we show that for arbitrary choices of x ≥ 0 and i ≤ k ≤ j, where 0 ≤ i ≤ j < ∞ or 0 ≤ i < j = ∞. Note that, when k < 2 or j = ∞, the left-hand side of (76) becomes −r k−1 (x)r j−1 (x) and thus (76) automatically holds. Now, assume that k ≥ 2 and j < ∞. From the definition in (73), the left-hand side of (76) is Using similar arguments, we can also show that On substituting the inequalities (76) and (77) into (75), we verify that the numerator of (74) is nonpositive, and therefore d 2 dx 2 h(x) ≤ 0, for all x > 0.

H. Proof of Proposition 5
Expressing the signal processing operations in Fig. 6 in the z-domain, we have where A(z K ) in the first equality is the z-transform of the K-times upsampled version of a m , and (78) follows from (45). Similar to (45), we can expand U (z) in terms of the z-transforms of its polyphase components, as U (z) = U 0 (z K ) + z −1 U 1 (z K ) + . . . + z −(K−1) U K−1 (z K ).
Comparing (78) against (79) and using the uniqueness of the polyphase decomposition, we conclude that U k (z) = A(z)G k (z), for all 0 ≤ k < K.
Now, consider Fig. 6(b). We note that the z-transform of g −m is G(z −1 ). Denote by d m the output of the filtering operation. Then, its z-transform can be computed as When downsampling d m by K, only the first term on the right-hand side of (80) will be retained; the second term is "invisible" to the sampling operation due to mismatched supports. Therefore, we have, after downsampling, V (z) = K−1 k=0 B k (z)G k (z −1 ).