Megapixel Photon-Counting Color Imaging using Quanta Image Sensor

Quanta Image Sensor (QIS) is a single-photon detector designed for extremely low light imaging conditions. Majority of the existing QIS prototypes are monochrome based on single-photon avalanche diodes (SPAD). Color imaging has not been demonstrated with single-photon detectors due to the intrinsic difficulty of shrinking the pixel size and increasing the spatial resolution while maintaining acceptable intra-pixel cross-talk. In this paper, we present image reconstruction of the first color QIS with a resolution of $1024 \times 1024$ pixels, supporting both single-bit and multi-bit photon counting capability. Our color image reconstruction is enabled by a customized joint demosaicing-denoising algorithm, leveraging truncated Poisson statistics and variance stabilizing transforms. Experimental results of the new sensor and algorithm demonstrate superior color imaging performance for very low-light conditions with a mean exposure of as low as a few photons per pixel.


Introduction
Quanta image sensor (QIS) is a concept proposed by E. Fossum in 2005 aiming to overcome the limited sensor performance with the continuously shrinking pixel size [1,2].In a CMOS image sensor [3], images are acquired by accumulating hundreds to thousands of photoelectrons in each pixel and converting the charges to voltages and then to digital signals.As pixel size shrinks, the amount of photons that fall on individual pixels drops.Under low-light (or high-speed imaging) conditions, photon shot noise dominates the acquisition process.As the photon becomes scarce, eventually the signal-to-noise ratio will drop to a point that no meaningful images can be generated.Small pixel also reduces the full-well capacity which limits the dynamic range.QIS overcomes the problem by spatial-temporally oversampling the scene.In a QIS, each "image pixel" is partitioned into a group of smaller pixels called jots, and each jot is sensitive enough for photon counting.Over the past decade, numerous studies have demonstrated the feasibility of QIS, with a few prototypes developed [4][5][6][7][8].On the image reconstruction side, various methods have also been proposed [9][10][11][12].
Despite the advances in QIS hardware and image reconstruction algorithms, one important open problem is the color imaging capability of QIS.Color imaging using QIS is difficult for two reasons.First, the resolution of a typical single-photon detector array is not high, e.g., 320 × 240 pixels [13], or most recently 512 × 512 pixels [8].If a color filter array such as the standard Bayer pattern is implemented, the effective resolution will be halved to 160 × 120 or 256 × 256.For detection/tracking applications [14] this might be sufficient; but for photography the resolution could be low compared to the mainstream CMOS image sensors.Second, the color reconstruction of QIS data is different from CMOS image sensors because the QIS data is acquired under sparse-photon conditions with much lower signal-to-noise ratio.The demosaicing algorithm of QIS data needs to overcome the photon shot noise and accurately reproduce the color information.As we will discuss, QIS requires a joint reconstruction and demosaicing on top of the truncated Poisson statistics, a problem that has not been studied.
In this paper, we report the first demonstration of QIS color imaging.The QIS Pathfinder camera module developed at Gigajot Technology was used in this demonstration.It has a spatial resolution of 1024 × 1024, single-bit and multi-bit photon-counting outputs, and up to 1040 fps speed in single-bit mode.Figure 1 shows a snapshot of the work presented in this paper.The key innovation, besides the sensor, is a new customized image reconstruction algorithm which enables joint reconstruction and demosaicing of raw QIS data.Our method leverages a few recently developed techniques in the QIS literature, including the transform-denoise method [10] and the Plug-and-Play ADMM [15].We demonstrate both our new sensor and the new algorithm on real data, and make quantitative evaluation based on synthetic data.

Megapixel Quanta Image Sensor
When the concept of QIS was first proposed, the two mainstream design methodologies to realize single-photon detection were the electron-multiplying charge-coupled device (EMCCD) [6] and the single-photon avalanche diode (SPAD) [16,17].Both devices amplify the signal from a single photoelectron with a physical phenomenon called electron avalanche multiplication, where high electrical voltage (typically higher than 20V) is used to accelerate the photoelectron to produce more free electrons through collision.SPAD is more widely adopted in consumer applications because of its excellent low-noise performance and temporal resolution.The SwissSPAD [8,18,19] and the SPC-SPAD [7,13] are two better known examples.However, due to the relatively large pixel size (typically more than 5µm pitch), high dark count, and low quantum efficiency, SPADs are not ideal candidates for implementing the QIS devices.The QIS we use in this paper is based on a new type of single-photon detector proposed by Ma et al. [2,4,5,20,21].The approach is to enhance the voltage signal generated by a single photoelectron by reducing the capacitance of the pixel output node so that the single-photon signal can overcome the background thermal noise.
Table 1 shows a comparison between a conventional CMOS image sensor, SPAD, and the prototype QIS used in this work.Among the listed different sensor parameters, the most noticeable difference between the QIS and SPADs is the small pixel pitch, where the QIS realized 1.1 µm pixel size, but SPADs typically have a pitch size around 5-20 µm.The smaller pixel size makes high spatial resolution (e.g.10s or 100s of megapixel) possible with a small optical format (e.g.1/4").Another advantage of QIS over SPADs is the high quantum efficiency.Quantum efficiency measures the percentage of photons that are detectable by the sensor, which is critical for low-light imaging.SPAD cameras generally have a low quantum efficiency (e.g.<50%) because it requires read-out electronics to initiate and quench the electron avalanche amplification.This reduces the photosensitive areal percentage in the SPAD pixels and limits the quantum Fig. 2. Quanta Image Sensor (QIS) vs. CMOS Image Sensor (CIS).The top row shows a simulated CIS data at a photon level same as the QIS.The second row shows the real analog sensor data obtained from our prototype camera.For each photon flux level, we show both the raw input data and the denoised data.The CIS is assumed to have a read noise of 1.2e − r.m.s., and we do not considered dark count.With dark count, the performance of CIS will deteriorate even more.efficiency.In contrast, since electron avalanche gain is not required in the QIS device, it has a quantum efficiency over 80% within the visible wavelengths.

QIS Color Imaging Model
The imaging model of a color QIS is illustrated in Figure 3. Relative to the previous QIS imaging models, e.g., [10][11][12], the new model considers a color filter array which allows us to select wavelengths and multiple color channels.There are a few essential components in the model.When the scene image arrives at the sensor, the color filter array first selects the wavelength according to the colors.Each color pixel is then sensed using a photon-detector.In single-bit mode the detector reports a binary value, and in multi-bit mode the detector reports an integer up to the saturation limit.The measured data contains three subsampled sequences, each representing a measurement in the color channel.channels, respectively.QIS was originally designed as a spatial-temporally oversampling device which samples each pixel using K jots spatially.This means every pixel has K measurements, and totally there are M def = NK measurements.Thus, the model of this oversampling process is

Spatial
where G ∈ R 3M×3N represents the K-fold upsampling and interpolation filtering.The constant α is a gain factor.In our current prototype, the oversampling factor is K = 1 as we want to maximize the spatial resolution.In this case, we have G = I , and M = N.
Color Filter Array.To obtain color information, a color filter is placed on the top of each jot to collect light for only one of the RGB colors, and filter out the remaining two colors.This process is represented by three mutually exclusive M × M diagonal matrices or masks: S r , S g and S B for the red, green and blue channels, respectively.These matrices contain either zeros or ones on the diagonal, and S r + S g + S b = I .Incorporating the color filter array, the light exposure on the M jots of QIS is represented by We call θ the mosaiced image.If we let S def = [S r , S g , S b ] ∈ R M×3M , then θ and c are related by Photon Arrival.The photon arrival of QIS is modeled as a Poisson process.Let Y ∈ N M be a non-negative random integer denoting the number of photons arrived at each jot during one integration period.Then, the probability of observing For single-bit QIS, we put a voltage comparator to threshold the arrived photon count into a binary decision.Thus, the read out of each jot is a binary random variable where q > 0 is a threshold.
For multi-bit QIS, we use an L-bit quantizer to quantize the measurement.In this case, In other words, we count the photons until the number reaches a threshold.When the threshold is reached, we report the maximum number before saturation.Assumptions.While the above model is theoretically tractable, we made three assumptions.The first is regarding the read noise associated with Y .If read noise is present, the measured photon count will become y + η where η is typically i.i.d.Gaussian.We assume that the read noise is negligible compared to the photon count, and can be eliminated by the thresholding process.We find that this is often the case with the QIS device used in this paper, which provides ultra-low read noise (0.24e-rms at room temperature) and accurate photon-counting capability.See Figure 5 for an actual measurement of the photon count.
The second assumption we make is that the multi-bit saturation does not happen.That is, during the multi-bit sensing, we will never reach the upper threshold q and so B m = Y m , which is a standard Poisson random variable.In terms of hardware, this can be achieved by implementing an automatic exposure control scheme so that the threshold q is increased when the scene is too bright.We will leave this to a follow up paper as it is outside the scope of this work.Fig. 6.The image reconstruction pipeline consists of (i) a temporal binning step to sum the input frames, (ii) a variance stabilizing transform T to transform the measurement so that the variance is stabilized, (iii) a joint reconstruction and demosaicing algorithm to recovery the color, (iv) an inverse transform to compensate the forward transform, and (v) a tone mapping operation to correct the contrast.
The third assumption we make is that the camera is fast enough to catch up the scene movement, and the scene remains static during the multiple temporal samplings.This allows us to utilize multiple independent measurements over time to improve the statistics.In this case, the probability distribution of Y becomes where y m,t is the Poisson count at pixel m in time t.

QIS Color Image Reconstruction
The task of image reconstruction is to recover color scene c from the measurements B. In the gray-scale setting, we can formulate the problem as maximum-likelihood and solve it using convex optimization tools [9,11,12,23].We can also use learning-based methods, e.g., [24][25][26][27] to reconstruct the signal.The method we present here is based on the transform-denoise approach by Chan et al. [10].We prefer transform-denoise because it is a physics-based approach and is robust to different sensor configurations.For example, in learning-based approaches, if we change the number of frames to sum, we need to train a different model or neural network.

Reconstruction Pipeline
The pipeline of the proposed reconstruction algorithm is shown in Figure 6.Given the random variable B = {B m,t } with realizations b = {b m,t }, we first sum the frames to generate a single image Z = {Z m }: In the future, this step can be integrated into the hardware of the camera, so that the output of the camera will be the sum of multiple frames.Depending on whether we are using the single-bit mode or the multi-bit mode, the statistics of {Z m } follows either a Binomial distribution or a Poisson distribution.In single-bit mode, the distribution of {Z m } is where Ψ q (θ) θ j e −θ / j! is the probability of obtaining a one in the binary variable B. As shown in [11], Ψ q is the incomplete Gamma function.In multi-bit mode, since a sum of Poisson remains a Poisson, the distribution of {Z m } is If we assume that there is no processing of the bits besides temporal summation, then a simple maximum-likelihood of the probability will return us the solution: The operation M can be interpreted as a tone-mapping.

Variance Stabilizing Transform
Now, we want to process the data besides just summing the frames.However, in either single-bit or multi-bit mode, the Binomial or the Poisson statistics are not easy as the variance changes with the mean.This prohibits the use any off-the-shelf algorithms that are based on i.i.d.Gaussian assumptions, e.g., denoising or inpainting.The variance stabilizing transform (VST) is a statistical technique that tries to alleviate the changing variance problem [28][29][30].Among the different types of VSTs, we are particularly interested in the two transforms originally proposed by F. Anscombe in 1948 [31]: Here, the first transform is called the Anscombe Binomial, and the second is called the Anscombe Poisson.Both transforms are pixel-wise, meaning that we can apply the transform to each pixel independently.As proved in [10] and [31], the variance of T (z) is 1/4 for both the single-bit and the multi-bit case.The inverse T −1 is simply the algebraic inverse (or the unbiased inverse if we follow [32]).The important question for color imaging is now that the photon count z is the result of applying a color filter array matrix S = [S r , S g , S b ].That is, instead of observing the full color we only observe a subsampled version z = S z.However, since T is point-wise and S is a concatenation of three non-overlapping binary matrices, it holds that Consequently, if we define T (z) as the observed measurement (which is available from the sensor) and treat x def = T ( z) as the unknown quantity to be determined, then the problem can be formulated as where g is some regularization function controlling the smoothness of x.Note that ( 12) is a standard demosaicing-denoising problem assuming i.i.d.Gaussian noise.

Joint Reconstruction and Demosaicing
The optimization problem in ( 12) is a standard least squares with regularization function g.Thus, most convex optimization algorithm can be used as long as g is convex.In this paper, we adopt a variation of the alternating direction method of multiplier (ADMM) by replacing g with an off-the-shelf image denoiser.Such algorithm is coined the name Plug-and-Play (PnP) [33] (and different versions thereafter [15]).
For the particular problem in (12), the PnP ADMM algorithm iteratively updates the following two steps: Demosaicing Module: Denoising Module: and updates the Lagrange multiplier by u (k+1) = u (k) − (x (k+1) − v (k+1) ).Readers interested in the detailed derivation can consult, e.g., [15].Here, ρ is an internal parameter that controls the convergence.The operator D is an off-the-shelf image denoiser, e.g., BM3D or deep neural network denoisers.The subscript ρ/λ denotes the denoising strength, i.e., the hypothesized "noise variance".Since S T S is a diagonal matrix, the inversion is pointwise.

Non-Iterative Algorithm
The 1.1µm pixel pitch of the proposed QIS can potentially lead to a spatial resolution as high as or even higher than a conventional CMOS sensor.When this happens, in certain applications we can trade-off the color reconstruction efficiency and the resolution.For example, instead of using one jot for one pixel, we can use four jots for one pixel as shown in Figure 7.In the future when QIS can achieve higher spatial resolution, we can use four color jots to reconstruct one pixel.In this case, we can bypass the iterative ADMM algorithm and use a one-shot denoising method.
Using four jots for one pixel allows us to bypass the iterative ADMM steps because there is no more missing pixel problem.In this case, the matrix S ∈ R M×3M will become S = diag{ 1  4 I, 1 2 I, 1 4 I } ∈ R 3M×3M , and hence ( 12) is simplified to a denoising problem with different noise levels for the three channels.In particular, the green channel has half of the variance of the red and the blue.For implementation, we can modify a denoiser, e.g., BM3D to accommodate the different noise variances.Since there is no more ADMM iteration, the algorithm is significantly faster.

Comparisons
We compare the proposed method with several existing methods on a synthetic dataset shown in Figure 8.We simulate the raw color QIS data by assuming a Bayer pattern and using the  [34]; (e) Hirakawa's PSDD method [36], with a built-in wavelet shrinkage denoiser; (f) Proposed method with BM3D.
image formation pipeline described in the previous section.We demosaic the images using: (a) a baseline method using MATLAB's demosaic preceded by gray-scale BM3D denoising of R, G 1 , G 2 and B channels and followed by color BM3D denoising; (b) LSLCD method [34], which has a built-in BM3D denoiser; (c) Hirakawa's PSDD method [35,36], which does joint denoising and demosaicing for Poisson noise; and (d) the proposed method using BM3D with (λ, ρ) = (0.001, 5).We apply variance stabilizing transform, except for PSDD which is designed for Poisson noise.
The results show that the proposed method has a significantly better performance, both in terms of PSNR and visual quality.

Synthesized QIS Data
We conduct a synthetic experiment to provide a quantitative evaluation of the performance of the proposed algorithm.To this end, we simulate the image formation pipeline by passing through color images to generate the QIS raw input data, with different number of bits. Figure 9 shows one example, and in the Supplementary Report we have additional examples.
In our simulation, we assume that the number of QIS frames is T = 4, and the average number of photon per pixel is 0.28, 0.85, 1.98 and 4.23 photons / frame for 1-bit, 2-bit, 3-bit and 4-bit QIS, respectively.On the measurement side, we generate single-bit and multi-bit data by thresholding the raw sensor output.To reconstruct the image, we use the proposed method with PnP and BM3D.The parameters are set as ρ = 1 and λ = 0.007, 0.003, 0.002 and 0.0007 for 1-bit, 2-bit, 3-bit and 4-bit QIS, respectively.With as low as 1-bit, the reconstructed image in Figure 9 is already capturing most of the features.As the number of bits increases, the visual quality improves.

Real QIS Data
For real images, we first show the reconstruction of an image of the Digital SG ColorChecker chart.We generate two sets of measurements: (a) a set of 50 one-bit frames, quantized with a threshold q = 4, and (b) a set of 10 five-bit frames.We use PnP and Bm3D with (λ, ρ) = (0.15, 10) and (0.01, 10) for the 1-bit and 5-bit data, respectively.After reconstruction, the results are multiplied by a 3 × 3 color correction matrix to mitigate any color cross-talk.This matrix is generated by linear least square regression of the 24 Macbeth color patches that lies inside the SG ColorChecker chart.The results in Figure 1 suggest that while 1-bit mode has color discrepancy with the ground truth, the 5-bit mode is producing a reasonably-high color accuracy.Next, we show the result of imaging real scenes.See Figure 10.In this experiment, the exposure time for each frame is set to 50µs.The average number of photons per pixel is approximately 4.2 photons for the "QIS" sign image, 3 photons for the Pathfinder image, 1.9 photons for the duck image, and 2.9 photons for the mushroom image.For all images, we collect the data using a 5-bit QIS.
We demonstrate two algorithms: (i) the proposed transform-denoise framework using PnP + BM3D, and (ii) assuming four-jots to one-pixel scenario by trading half of the spatial resolution.The typical runtime on an unoptimized MATLAB code is approximately 4 minutes for PnP, and 10 seconds for the four-jot to one-pixel method.An interesting observation is that even in the lower-resolution case, the details are not significantly deteriorated unless we zoom-in.However, the speed up we get is substantial.

Conclusion and Discussion
We demonstrated the first color imaging for a mega-pixel Quanta Image Sensor (QIS).We designed a joint reconstruction-demosaicing algorithm for the quantized Poisson (and Poisson) statistics of the sensor data.Our solution involves a variance stabilizing transform and an iterative (and a non-iterative) algorithm.By integrating the sensor and the new algorithm, the sensor is able to acquire color images at a photon level of 1e − or less, a level that would be difficult for standard CMOS image sensors.The average number of photons per frame is 4.2, 3.0, 1.9, and 2.9 for each image respectively.Both methods use 4 frames for reconstruction.The raw data has a resolution of 1024 × 1024 pixels.The ADMM method retains the resolution, whereas the non-iterative method reduces the resolution to 512 × 512.Reconstruction using both the methods are shown at the same size for easier visual comparison.Notice that the non-iterative algorithm is able to achieve a visual quality almost similar to the ADMM method.

Fig. 1 .
Fig. 1.(a) One 1-bit frame.(b) Reconstructed color image using 50 frames of 1-bit input with threshold q = 4. (c) One 5-bit frame.(d) Reconstructed color image using 10 frames of 5-bit input.The average number of photons per frame is 5.
Oversampling.We represent the discretized scene image as a column vector c = [c r ; c g ; c b ] ∈ [0, 1] 3N , where c r , c g , and c b are the light intensities of the red, green and blue

Fig. 3 .
Fig.3.QIS Imaging Model.When the scene image arrives at the sensor, the color filter array first selects the wavelength according to the colors.Each color pixel is then sensed using a photon-detector.In single-bit mode the detector reports a binary value, and in multi-bit mode the detector reports an integer up to the saturation limit.The measured data contains three subsampled sequences, each representing a measurement in the color channel.

Fig. 4 .Fig. 5 .
Fig. 4. Gigajot prototype QIS camera module and the QIS sensor chip.Note that there is no additional optics required besides a simple focusing lens.The camera is powered by standard 5V DC input, and has a USB3 data interface to transmit data to external storage.

Fig. 7 .
Fig.7.In the future when QIS can achieve higher spatial resolution, we can use four color jots to reconstruct one pixel.In this case, we can bypass the iterative ADMM algorithm and use a one-shot denoising method.

Fig. 8 .
Fig.8.Simulated QIS experiment.The goal of this experiment is to compare the proposed iterative algorithm with existing methods.We assume the observed Bayer RGB image is from a 3-bit QIS sensor.(a) Ground Truth; (b) One 3-bit QIS frame; (c) MATLAB demosaic preceded and followed by BM3D; (d) LSLCD[34]; (e) Hirakawa's PSDD method[36], with a built-in wavelet shrinkage denoiser; (f) Proposed method with BM3D.

Fig. 9 .
Fig. 9. Synthetic experiment for quantitative evaluation.[Top row]: One frame of the QIS measurements using different number of bits.[Bottom row]: Reconstructed images using the proposed method with 20 frames of QIS data.The average photon counts per pixel are 0.25, 0.75, 1.75 and 3.75 for 1-bit, 2-bit, 3-bit and 4-bit QIS, respectively.
(a) Raw QIS 5-bit Data. 4 frames with exposure time of 50 µs were obtained for each scene.(b) Reconstruction using the proposed method (iterative) 1024 × 1024 (c) Reconstruction using the proposed method (non-iterative) 512 × 512

Fig. 10 .
Fig. 10.Real QIS Image Reconstruction.The exposure time for each frame is 50 µs.The average number of photons per frame is 4.2, 3.0, 1.9, and 2.9 for each image respectively.Both methods use 4 frames for reconstruction.The raw data has a resolution of 1024 × 1024 pixels.The ADMM method retains the resolution, whereas the non-iterative method reduces the resolution to 512 × 512.Reconstruction using both the methods are shown at the same size for easier visual comparison.Notice that the non-iterative algorithm is able to achieve a visual quality almost similar to the ADMM method.

Table 1 .
Comparison Mean Dark Count Rate> 10 e − /s 47 e − /pix/s 75 e − /pix/s 0.068 e − /s of the available image sensor technologies.Comparison between a standard CMOS sensor, single-photon avalanche diode (SPAD) cameras, and QIS.The QIS we present in this work has significantly smaller pixel pitch compared to the two existing SPAD cameras, and is half of a CMOS pixel.In addition, we support color imaging which are not currently available on SPAD cameras.