Seeing into Darkness:
Scotopic Visual Recognition

Our computational framework starts from a model of image capture. Each pixel in an image reports the brightness estimate of a cone of visual space by counting photons coming from that direction. The estimate improves over time. Starting from a probabilistic assumption of the imaging process and of the target classification application, we derive an approach that allows for the best tradeoff between exposure time and classification accuracy.

We make three assumptions and relax them later:

• •

  The world is stationary during the imaging process. This may be justified as many photon-counting sensors sample the world at $>1kHz$ [8, 9]. Later we test the model under different camera movements and show robust performance.

• •

  Photon arrival times follow a homogeneous Poisson process. This assumption is only used to develop the model. We will evaluate the model in Sec. 4.4 using observations from realistic noise sources.

• •

  A probabilistic classifier based on photon counts is available. We discuss how to obtain such a model in Sec. 3.4.

Formally, the input ${\bf X}_{1:t}$ is a stream of photons incident on the sensors during time $[0,t\Delta)$, where time has been discretized into bins of length $\Delta$. $X_{t,i}$ is the number of photons arrived at pixel $i$ in the $t$th time interval, i.e. $[(t-1)\Delta,t\Delta)$. The task of a scotopic visual recognition system is two fold: 1) computing the category $C\in\{0,1,\ldots,K\}$ of the underlying object, and 2) crucially, determining and minimizing the exposure time $t$ at which the observations are deemed sufficient.

Our decision strategy for trading off accuracy and speed is based on SPRT, for its simplicity and attractive optimality guarantees. Assume that a probabilistic model is available to predict the class label $C$ given a sensory input ${\bf X}_{1:t}$ of any duration $t$ – either provided by the application or learned from labeled data using techniques described in Sec. 3.4 – SPRT conducts a simple accumulation-to-threshold procedure to estimate the category $\hat{C}$:

Let $S_{c}({\bf X}_{1:t})\buildrel\triangle\over{=}\log\frac{P(C=c|{\bf X}_{1:t})}{P(C\neq c|{\bf X}_{1:t})}$ denote the class posterior probability ratio of the visual category $C$ for photon count input ${\bf X}_{1:t}$, $\forall c\in\{1,\ldots,K\}$, and let $\tau$ be an appropriately chosen threshold. SPRT conducts a simple accumulation-to-threshold procedure to estimate the category $\hat{C}$:

When a decision is made, the declared class $\hat{C}$ has a probability that is at least $\exp(\tau)$ times bigger thanthe probability of all the other classes combined, therefore the error rate of SPRT is at most $1-Sigm(\tau)$, where $Sigm$ is the sigmoid function: $Sigm(x)\buildrel\triangle\over{=}\frac{1}{1+\exp(x)}$.

For simple binary classification problems, SPRT is optimal in trading off speed versus accuracy in that no other algorithm can respond faster while achieving the same accuracy [18]. In the more realistic case where the categories are rich in intra-class variations, SPRT is shown to be asymptotically optimal, i.e. it gives optimal error rates as the exposure time becomes large [19]. Empirical studies suggest that even for short exposure times SPRT is near-optimal [45].

In essence, SPRT decides when to respond dynamically, based on the stream of observations accumulated so far. Therefore, the response time is different for each example. This regime is called “free-response” (FR), in contrast to the “interrogation” (INT) regime, typical of photopic vision, where a fixed-length observation is collected for each trial [46]. The observation length may be chosen according to a training set and fixed a priori. In both regimes, the length of observation should take into account the cost of errors, the cost of time, and the difficulty of the classification task.

Despite the striking similarity between the two regimes, SPRT (the FR regime) is provably optimal in the asymptotical tradeoff between response time and error rate, while such proofs do not exist for the INT regime. We will empirically evaluate both regimes in Sec. 4.3.3.

The challenge of applying SPRT is to compute $S_{c}({\bf X}_{1:t})$ for class $c$ and the input stream ${\bf X}_{1:t}$ of variable exposure time $t$, or in a more information-relevant unit, variable PPP levels. Thanks to the Poisson noise model (Eq. 1), the sufficient statistics for observation ${\bf X}_{1:t}$ is the cumulative count $\boldsymbol{N}_{t}=\sum_{t^{\prime}=1}^{t}{\bf X}_{t^{\prime}}$ (visualized in  Fig. 1), and the observation duration length $t$, therefore we may rewrite $S_{c}({\bf X}_{1:t})$ as $S_{c}(\boldsymbol{N}_{t},t)$. We further shorthand the notation to $S_{c}(\boldsymbol{N}_{t})$ since the exposure time is evident from the subscript. Since counts at different PPPs (and, equivalently, exposure times) have different statistics, it would appear that a specialized system is required for each PPP. This leads to the naive ensemble approach. Instead, we also propose a network called WaldNet that can process images at all PPPs and has the size of only a single specialized system. We describe the two approaches below.

Our goal is to train WaldNet to optimally trade off the expected exposure time (or PPP) and error rate in the FR regime. Optimality is defined by the Bayes risk $R$ [18]:

where $\mathbb{E}$[PPP] is the expected photon count required for classification, $\mathbb{E}[C\neq\hat{C}]$ is the error rate, and $\eta$ describes the user’s cost of photons per pixel (PPP) versus error rate. WaldNet asymptotically optimizes the Bayes risk provided that it can faithfully capture the class log posterior ratio $S_{c}(\boldsymbol{N}_{t})$, and selects the correct threshold $\tau$ (Eq. 3) based on the tradeoff parameter $\eta$. Sweeping $\eta$ traverses the optimal time versus error tradeoff (Fig. 2c).

Since picking the optimal threshold according to $\eta$ is independent from training a ConvNet to approximately compute the log posterior ratio $S_{c}(\boldsymbol{N}_{t})$, the same ConvNet is shared across multiple $\eta$’s. This suggests the following two-step learning algorithm.

Both scotopic algorithms (ensemble and WaldNet) assume knowledge of the light-level PPP in order to choose the right model parameters. This knowledge is easy to acquire when the illuminance is constant over time, in which case PPP is linear in the exposure time $t$ (Eq. 2), which may be measured by an internal clock.

However, in situations where the illuminance is dynamic and unknown, the linear relationship between PPP and exposure time is lost. In this case we propose to estimate PPP directly from the photon stream itself, as follows. Given a cumulative photon count image $\boldsymbol{N}$ ($t$, the time it takes to accumulate the photons, is no longer relevant as the illuminance is unknown), we examine local neighbors that receive high photon counts, and pool the photon counts as a proxy for PPP. In detail, we (1) convolve $\boldsymbol{N}$ using an $s\times s$ box filter, (2) compute the median of the top $k$ responses, and (3) fit a second order polynomial to regress the median response towards the true PPP. Here $s$ and $k$ are parameters, which are learned from a training set consisting of $(\boldsymbol{N},$PPP$)$ pairs. Despite its simplicity, this estimation procedure works well in practice, as we will see in  Sec. 4.

One major challenge of scotopic systems is to compute log posterior ratio computations as quickly as photons stream in. Photon-counting sensors [8, 9] sample the world at $1k-10kHz$, while the fastest reported throughput of ConvNet [48] is $2kHz$ for $32\times 32$ color images, and $800Hz$ for $100\times 100$ color images. Fortunately, the reported throughputs are based on independent images, whereas the scotopic systems operate on photon streams, where temporal coherence may be leveraged for accelerated processing. Since the photon arrival events within any time bin is sparse, changes to the input and the internal states of a scotopic system are small. An efficient implementation thus could model the changes and propagate only those that are above a certain magnitude.

One such implementation relies on spiking recurrent hidden units. A spiking recurrent hidden unit is characterized by two aspects: 1) computation at time $t$ reuses the unit’s state at time $t-1$ and 2) only changes above a certain level will be propagated to layers above.

Specifically, the first hidden layer of WaldNet may be exactly implemented using the following recurrent network where the features $S^{H}(\boldsymbol{N}_{t})$ (Eq. 4) are represented by membrane voltages $\boldsymbol{V}(t)\in\mathbb{R}^{n_{H}}$. The dynamics of $\boldsymbol{V}(t)$ is:

where $r(t)\buildrel\triangle\over{=}\frac{\alpha(t)}{\alpha(t-1)}$ is a damping factor, $\boldsymbol{l}(t)\buildrel\triangle\over{=}\boldsymbol{\beta}(t)-r(t)\boldsymbol{\beta}(t-1)$ is a leakage term (derivations in Sec. 0.A.3). The photon counts ${\bf X}_{t}$ in $[(t-1)\Delta,t\Delta)$ is sparse, thus the computation $\boldsymbol{W}{\bf X}_{t}$ is more efficient than computing $S^{H}(\boldsymbol{N}_{t})$ from scratch.

To propagate only large changes to the layers above, we use a similar thresholding mechanism as (Eq. 3). For each hidden unit $j$, we associate a ‘positive’ and a ’negative’ neuron that communicate with the layer above. For each time bin $t$:

where $\tau_{dis}>0$ is a discretization threshold. By taking the difference between the spike counts from the ‘positive’ and the ‘negative’ neuron, the layers above can reconstruct a discretized version of $V_{j}(t)$. Hidden units from higher layers in WaldNet may be approximated using spiking recurrent units (Eq. 9) in a similar fashion.

The discretization threshold affects not only the number of communication spikes, but also the quality of the discretization, and in turn the classification accuracy. For spike-based hardwares [49], the number of spikes is an indirect measure of the energy consumption (Fig. 4(B) of [49]). For non-spiking hardwares, the number of spikes also translate to the number of floating point multiplications required for the layers above. Therefore, the $\tau_{dis}$ controls the tradeoff between accuracy and power / computational cost. We will empirically evaluate this tradeoff in Sec. 4.

We compare WaldNet against the following baselines, under both the INT regime and the FR regime:

Ensemble. We construct an ensemble of $4$ specialists with PPPs from $\{.22,2.2,22,220\}$ respectively. The performance of the specialists at their respective PPPs gives a lower bound on the optimal performance by ConvNets of the same architecture.

Photopic classifier. An intuitive idea is to take a network trained in normal lighting conditions, and apply it to properly rescaled lowlight images. We choose the specialist with PPP$=220$ as the photopic classifier as it achieves the same accuracy as a network trained with $8$-bit images.

Rate classifier. A ConvNet on the time-normalized image (rate) without weight adaptation. The first hidden layer is computed as $S^{H}_{j}(\boldsymbol{N}_{t})\approx\boldsymbol{W}\boldsymbol{N}_{t}/t+\boldsymbol{b}^{H}$. Note the similarity with the WaldNet approximation used in Eq. 4.

For all models above we assume that an internal clock measures the exposure time $t$, and the illuminance $\lambda^{\phi}$ is known and constant over time. We remove this assumption for the model below for unknown and varying illuminance:

WaldNet with estimated light-levels. A WaldNet that is trained on constant illuminance $\lambda^{\phi}$, but tested in environments with unknown and dynamic illuminance. In this case the linear relationship between exposure time $t$ and PPP (Sec. 2) is lost. Instead, the light-level is first estimated according to Sec. 3.5 directly from the photon count image $\boldsymbol{N}$. The estimate $\hat{PPP}$ is then converted to an ‘equivalent’ exposure time $\hat{t}$ using $\hat{t}=\frac{\hat{PPP}}{\lambda^{\phi}\Delta}$ (by inverting Eq. 2), which is used to adapt the first hidden layer of WaldNet in Sec. 4, i.e. $S^{H}(\boldsymbol{N})\approx\alpha(\hat{t})\boldsymbol{W}\boldsymbol{N}+\boldsymbol{\beta}(\hat{t})$.

We consider two standard datasets: MNIST [17] and CIFAR10 [16]. We simulate lowlight image sequences using Eq. 1.

MNIST contains gray-scaled $28\times 28$ images of $10$ hand-written digits. It has $60k$ training and $10k$ test images. We treat the pixel values as the ground truth intensity⁷⁷7The brightest image we synthesize has about $2^{8}$ photons, which corresponds to a pixel-wise maximum signal-to-noise ratio of $16$ ($4$-bit accuracy), whereas the original MNIST images has ($7$ to $8$-bit accuracy) that corresponds to $2^{14}$ to $2^{16}$ photons. . Dark current $\epsilon_{dc}=3\%$. We use the default LeNet architecture from the MatConvNet package [50] with batch normalization [51] after each convolution layer. The architecture is $784$-$20$-$50$-$500$-$10$⁸⁸8The first and last number represent the input and output dimension, each number in between represents the number of feature maps used for that layer. The number of units is the product of the number of features maps with the size of the input. with $5\times 5$ receptive fields and $2\times 2$ pooling.

CIFAR10 contains $32\times 32$ color images of $10$ visual categories. It has $50k$ training and $10k$ test images. We use the same sythensis procedure above to each color channel⁹⁹9For simplicity we do not model the Bayer filter mosaic.. We again use the default $1024$-$32$-$32$-$64$-$10$ LeNet architecture [22] with batch normalization. We use the same setting prescribed in [22] to achieve $18\%$ test error on normal lighting conditions. [22] uses local contrast normalization and ZCA whitening as preprocessing steps. We estimate the local contrast and ZCA from normal lighting images and transforming them according to the lowlight model to preprocess scotopic images.

The speed versus accuracy tradeoff curves in the INT regime are shown in Fig. 3a,b. Median PPP versus accuracy tradeoffs for all models in the FR regime are shown in Fig. 3c,d. All models use constant thresholds for producing the tradeoff curves. In Fig. 4a are average PPP versus accuracy curves when the models use optimized dynamic thresholds described in Sec. 3.4, step-two.

How robust is the speedup observed in Sec. 4.3 affected by sensor noise? For MNIST and CIFAR10, we take WaldNet and vary independently the dark current, the read noise and the fixed pattern noise. We also introduce a rotational jitter noise to investigate the model’s robustness to motion. The jitter applies a random rotation to the camera (or equivalently to the object being imaged by a stationary camera) where the rotation at PPP follows a normal distribution: $\Delta\theta\sim\mathcal{N}(0,\left(\frac{\sigma_{\theta}PPP}{220}\right)^{2})$, where $\sigma_{\theta}$ controls the level of jitter. e.g. $\sigma_{\theta}=22^{\circ}$ means that at $PPP=220$, the total amount of rotation applied to the image has an std of $22^{\circ}$. The result is shown in Fig. 5a,b.

First, the effect of dark current and fixed pattern noise is minimal. Even an $11\%$ dark current (i.e. photon emission rate of the darkest pixel is $10\%$ of that of the brightest pixel) merely doubles the exposure time with little loss in accuracy. The multiplicative fixed pattern noise does not affect performance because WaldNet in general makes use of very few photons. Second, current industry standard of read noise ($\sigma_{r}=22\%$ [9]) guarantees no performance loss for MNIST and minor loss for CIFAR10, suggesting the need for improvement in both the algorithm and the photon-counting sensors. The fact that $\sigma_{r}=50\%$ hurts performance also suggests that single-photon resolution is vital for scotopic vision. Lastly, while WaldNet provides certain tolerance to rotational jitter, drastic movement ($22^{\circ}$ at $220$ PPP) could cause significant drop in performance, suggesting that future scotopic recognition systems and photon-counting sensors should not ignore camera / object motion.

Finally, we inspect the power efficiency of the spiking network implementation (Eq. 9, 10) on the MNIST dataset. Our baseline implementation (“Continuous”) runs a ConvNet from end-to-end every time the input is refreshed. As a proxy for power efficiency we use the number of multiplications [49], normalized by the total number in the baseline. For simplicity we vary the discretization threshold $\tau_{dis}$ for inducing spiking events (Eq. 10), and the threshold is common across all layers.

The power, speed and accuracy results shown in Fig. 6a,b suggest that for $\tau_{dis}\leq 0.2$ discretization not only faithfully preserves the SAT of WaldNet (Fig. 6a), but also could be optimized to consume only $35\%$ of the total multiplications, i.e. the spiking implementation provides a $3\times$ power reducation. The amount of spiking events starts high and tappers off gradually (Fig. 6c) as the noise in the hidden unit estimates (Eq. 9) improves over time. Thus most of the computational savings comes at the later stage ($PPP\geq 22$). Further savings may reside in optimizing the discretization thresholds per layer or over time, which we reserve for future investigations.