IEEE TRANSACTIONS ON COMPUTATIONAL IMAGING, VOL. X, NO. X, XXX 20XX

1

Optimal Threshold Design for Quanta Image Sensor

arXiv:1704.03886v2 [cs.CV] 30 Oct 2017

Omar A. Elgendy, Student Member, IEEE and Stanley H. Chan, Member, IEEE

Abstract—Quanta Image Sensor (QIS) is a binary imaging
device envisioned to be the next generation image sensor after
CCD and CMOS. Equipped with a massive number of single
photon detectors, the sensor has a threshold q above which the
number of arriving photons will trigger a binary response “1”, or
“0” otherwise. Existing methods in the device literature typically
assume that q = 1 uniformly. We argue that a spatially varying
threshold can significantly improve the signal-to-noise ratio of
the reconstructed image. In this paper, we present an optimal
threshold design framework. We make two contributions. First,
we derive a set of oracle results to theoretically inform the maximally achievable performance. We show that the oracle threshold
should match exactly with the underlying pixel intensity. Second,
we show that around the oracle threshold there exists a set of
thresholds that give asymptotically unbiased reconstructions. The
asymptotic unbiasedness has a phase transition behavior which
allows us to develop a practical threshold update scheme using
a bisection method. Experimentally, the new threshold design
method achieves better rate of convergence than existing methods.

(a) Observed, q = 3

(b) Observed, q = q ∗ (c)

(c) Observed, q = 12

(d) Reconstruction, q = 3

(e) Reconstruction, q = 12

(d) Reconstruction, q = q ∗ (c)

(e) Ground Truth

Index Terms—Quanta image sensor, single-photon imaging,
high dynamic range, binary quantization, maximum likelihood.

I. I NTRODUCTION
A. Threshold Design for Quanta Image Sensor
Quanta Image Sensor (QIS) is a class of solid-state image
sensors envisioned to be the next generation imaging device
after CCD and CMOS. Originally proposed by Eric Fossum
in 2005 [1], the sensor has gained significant momentum in
the past decade, both in terms of hardware design [2]–[4]
and image processing [5]–[9]. The advantage of QIS over the
mainstream CCD and CMOS is attributed to its high spatial
resolution (e.g., 109 pixels per sensor with 200nm pitch per
pixel [10]) and high speed (e.g., 100k fps as reported in
[11]). However, in order to simplify circuit, minimize power
and reduce data transfer, QIS is operated in a binary mode:
When the number of photons arriving at the sensor exceeds a
threshold q, the sensor generates a binary bit “1”. When the
number of photons is less than q, the sensor generates a “0”.
The goal of this paper is to address the question of how to
optimally choose q to maximize the signal-to-noise ratio of
the reconstructed image.
Optimal threshold design for QIS is important as it directly
affects the dynamic range of an image. Figure 1 illustrates an
example where we simulate the raw binary data acquired by
The authors are with the School of Electrical and Computer Engineering,
Purdue University, West Lafayette, IN 47907, USA. Email: { oelgendy,
stanchan}@purdue.edu. The work was supported, in part, by the U.S.
National Science Foundation under Grant CCF-1718007. A preliminary version of this paper was presented at the 2016 IEEE International Conference
on Image Processing (ICIP), Phoenix, AZ.
This paper follows the concept of reproducible research. All the results and
examples presented in the paper are reproducible using the code and images
available online at http://engineering.purdue.edu/ChanGroup/.

Fig. 1. Simulated QIS data and the reconstructed gray-scale images using
different thresholds. Top row: The binary measurements obtained using
thresholds q = 3, q = q ∗ (c), and q = 12. Bottom figures: The maximum
likelihood estimates obtained from the binary measurements, with comparison
to the ground truth. The results show that our spatially varying threshold q ∗ (c)
offers the best reconstruction. In this experiment, we spatially oversample each
pixel by K = 2 × 2 binary bits and we use T = 25 independent temporal
measurements.

a QIS using a uniform threshold q. When q is low, most of
the bits in the raw input are “1”. The reconstructed image is
therefore an over-exposed image. On the other hand, when
q is high, most of the bits in the raw input are “0”. The
reconstructed image is then under-exposed. In both cases, it
is evident from the simulation that a uniform threshold has
limited performance. A better way is to allow q to vary spatially so that a pixel (or a group of pixels) has its own threshold
value. The result in Figure 1(d) shows the reconstruction result
using a spatially varying threshold obtained from our proposed
technique, which is clearly better than the uniform thresholds.
B. Scope and Contributions
The goal of this paper is to present an optimal threshold
design methodology and provide theoretical justifications. The
two major contributions are summarized as follows.
First, we provide a rigorous theoretical analysis of the
performance limit of the image reconstruction as a function of
the threshold. These results form the basis of our subsequent
discussions of the threshold update scheme. Some results are

2

IEEE TRANSACTIONS ON COMPUTATIONAL IMAGING, VOL. X, NO. X, XXX 20XX

TABLE I
L IST OF QIS P ROTOTYPES AND PARAMETERS
Camera

Canon 5D CMOS

EMCCD [12]

GMAPD [13]

SPC SPAD [14]

SwissSPAD [11]

Fossum QIS [15]

Price

$5, 000

$20, 000

Prototype

Prototype

Prototype

Prototype

Resolution

4096 × 2160

1024 × 1024

256 × 256

320 × 240

512 × 128

1376 × 768

Full-well capacity

69 ke- (@ISO100)

180 ke-

-

-

Frames per second (fps)

6
88.6 Mbps

8 × 103

1 − 250 e-

Sensor data rate

26 − 92

56 − 125 e-

Pixel Pitch (µm)

2.3

13

0.48 Gbps

25

0.52 Gbps

known, e.g., the signal-to-noise ratio is a function of the Fisher
Information [16], [17], but a number of new results are shown.
In particular, we show that (i) the maximum likelihood estimate has a closed-form expression in terms of the incomplete
Gamma function (Section III.B); (ii) the oracle threshold can
be derived in closed-form by maximizing the signal-to-noise
ratio (Section III.C); (iii) the image reconstruction has a phase
transition behavior (Section IV.A - Section IV.D).
Second, we propose an efficient threshold update scheme
based on our theoretical results. The new scheme is a bisection
method which iteratively updates the threshold without the
need of reconstructing the image. By checking whether the
proportion of one’s and zero’s approaches 0.5 in a spatialtemporal block, the threshold is guaranteed to be near optimal.
Compared to other existing threshold update schemes such as
[18] and [19]–[21], the new scheme offers significantly faster
rate of convergence (Section IV.E). We also demonstrate how
the dynamic range can be extended for high dynamic range
(HDR) imaging (Section IV.F).
A preliminary version of this paper was presented in ICIP
2016 [8]. This journal version contains significantly more details including complete proofs of major results, more comprehensive comparisons with existing methods, and discussions of
HDR imaging.
II. BACKGROUND
A. Current State of QIS
Quanta Image Sensor (QIS) belongs to the family of
photon-counting devices. These photon-counting devices have
been known for a long time. Some better-known examples
are the electron-multiplying charge-coupled device (EMCCD)
[22], [23], single-photon avalanche diode (SPAD) [11], [14],
[24], Geiger-mode avalanche photodiode (GMAPD) [13], etc.
The common feature of these devices is their single photon
sensitivity, which makes them useful in medical imaging
[25]–[27], astronomy [28], defense [29], nuclear engineering
[30], depth and reflectivity reconstruction [31], ultra-fast lowlight tracking [32], and recently in quantum random number
generation used in cryptography [33], [34].
The concept of QIS was first proposed by Fossum in 2005
as a solution for sub-diffraction limit pixels. The sensor was
called the digital film sensor, and later the quanta image sensor
[15], [35], [36]. After the introduction of QIS, researchers
in EPFL developed a similar concept called the Gigavision
camera [6], [37], [38]. Recently, teams at the University of
Edingburgh [14], [24], [39] and EPFL [33], [40] have made
new progresses in QIS using binary single photon detectors.

8

2 × 104

1.54 Gbps

24

1.56 × 105

10.24 Gbps

3.6

1 × 103
1 Gbps

In the industry, Rambus Inc. (Sunnyvale, CA) has developed
binary image sensors for high dynamic range imaging [19]–
[21]. Table I lists several recent QIS prototypes that are
available or are currently being developed. As a comparison
we also show a Canon 5D Mark III CMOS camera. Among
many different features, the most noticeable is the frame
rate. For example, SPS SPAD can be operated at 20k fps.
SwissSPAD can even achieve 156k fps. Both are significantly
faster than a standard CMOS camera.
B. Related Work on Threshold Design
Existing work on QIS threshold design study can be summarized into three classes of methods.
• Markov Chain [18]. The Markov Chain method developed by Hu and Lu [18] is a time-sequential update
scheme. A Markov Chain probability is used to control
how easy the threshold should be increased or decreased.
While the method has provable convergence, the threshold of each single photon detector of the QIS has to be
updated sequentially in time. In contrast, our proposed
method allows a group of single photon detectors to share
the same threshold. As a result, our proposed method has
significantly faster rate of convergence.
• Conditional Reset [19]–[21]. The conditional reset
method is a hardware solution proposed by Vogelsang
and colleagues. The idea is to take a sequence of images
with ascending (or descending) thresholds, and digitally
integrate the sequence to form an image. The drawback
of the method, besides the additional hardware cost of
the per-pixel reset transistors, is the limited quality of
the reconstructed image. For the same number of frames,
our proposed method produces better images.
• Checkerboard Threshold [16]. This method constructs a
checkerboard of thresholds by alternating two threshold
values q1 and q2 . The optimality criterion of q1 and q2
is based on minimizing the Cramér-Rao lower bound
(CRLB) integrated over a range of light intensities, which
is essentially an average case result. Our proposed method
obtains the optimal threshold for each pixel. This perpixel optimization has higher reconstruction performance
compared to checkerboard threshold.
C. QIS Imaging Model
In this subsection we provide an overview of the QIS
imaging model. The model has been previously discussed in
several papers, e.g., [6]–[9]. Readers interested in details can
refer to these papers for further explanations.

ELGENDY-CHAN: OPTIMAL THRESHOLD DESIGN FOR QUANTA IMAGE SENSOR

Fig. 2.

×

K

gk

θm
Poisson

Ym

•

( K

n

θ = αGc

α
cn

•

•

3

bm,T

T

T

M

}|

(z

{

bm,2

bm,1

Reconstruction

−−−→

Block diagram illustrating the image formation process of QIS.

1) Spatial Oversampling: We denote the discrete version
of the light intensity as a vector c = [c0 , . . . , cN −1 ]T , where
n = 0, . . . , N − 1 specify the spatial coordinates. We assume
that cn is normalized to the range [0, 1] for all n so that there
is no scaling ambiguity. To model the actual light intensity,
we multiply cn by a constant α to yield αcn , where α > 0 is
a fixed scalar constant.
Given the N -dimensional vector c, QIS uses M  N tiny
def
pixels called jots to sample c. The ratio K = M/N is known
as the spatial oversampling factor. The oversampling process
is illustrated in Figure 2, where it first upsamples the vector
c by a factor of K, and then filters the output by a lowpass
filter {gk }. Mathematically, the process can be expressed as
θ = αGc,

(1)

where θ = [θ0 , . . . , θM −1 ]T denotes the light intensity sampled at the M jots, and the matrix G is defined as
G=

1
I N ×N ⊗ 1K×1 ,
K

(2)

where 1K×1 is a vector of all ones and ⊗ denotes the
Kronecker product. Note that the choice of G in (2) is the
result of simplifying the model by assuming that the lowpass
filter is gk = 1/K for all k. This assumption is typically
reasonable, because on each QIS jot there is a micro-lens to
focus the incident light. Although previous papers, e.g., [6],
[7], do not make such assumption, in this paper we decide to
use a simplified G, for otherwise the theoretical analysis will
become very complicated. Nevertheless, in the Supplementary
Material we show comparison between a general G and the
simplified G. The gap is usually insignificant.
2) Truncated Poisson Process: We assume that the operating speed of QIS is significantly faster than the scene motion.
Therefore, for a given scene c (and also θ), we are able to
acquire a set of T independent measurements. We illustrate
this using the T channels in Figure 2.
The oversampled signal θ generates a sequence of Poisson
random variables according to the distribution
y

P(Ym,t = ym,t ) =

θmm,t e−θm
,
ym,t !

(3)

where m = 0, 1, . . . , M −1 denotes the m-th jot of the QIS and
t = 0, 1, . . . , T − 1 denotes the t-th independent measurement
in time. Denoting q ∈ N as the quantization threshold, the final
observed binary measurement Bm,t is a truncation of Ym,t :
(
0, if Ym,t < q.
Bm,t =
1, if Ym,t ≥ q

z

N

}|

{

Fig. 3. Image reconstruction of QIS data. Given the binary bit planes, the
reconstruction algorithm estimates the gray-scale image shown on the right.

The probability mass function of Bm,t is given by
q−1
k −θm
P θm

e

, if bm,t = 0,

k!
P(Bm,t = bm,t ) = k=0
∞
k
−θ
P θm e m


, if bm,t = 1.

k!

(4)

k=q

The goal of image reconstruction is to recover the underlying image c from the binary measurements B = {Bm,t | m =
0, . . . , M − 1, and t = 0, . . . , T − 1}. A pictorial illustration
of the reconstruction is shown in Figure 3.
3) Properties of Truncated Poisson Processes: The probability mass function of Bm,t in (4) is Bernoulli. However, the
right hand side of (4) involves infinite sums which are difficult
to interpret. To simplify the equations, we consider the upper
incomplete Gamma function Ψq : R+ → [0, 1] defined in [41]
as:
Z ∞
1
def
tq−1 e−t dt, for θ > 0, q ∈ N.
Ψq (θ) =
Γ(q) θ

where Γ(q) = (q − 1)! is the standard Gamma function. The
incomplete Gamma function allows us to rewrite the infinite
sums in (4) using the following identity [41]:
Ψq (θ) =

q−1
X
θk −θ
e .
k!

(5)

k=0

Consequently, the probabilities in (4) become
P(Bm,t = 0) = Ψq (θm ),
P(Bm,t = 1) = 1 − Ψq (θm ).

(6)

Example 1. In the special case of q = 1, we obtain:
Z ∞
1
P(Bm,t = 0) =
t0 e−t dt = e−θm ,
Γ(1) θm

which coincides with the results shown in [6] and [7].
The incomplete Gamma function Ψq (θ) is a decreasing
function of θ because the first order derivative of Ψq (θ) with
respect to θ is negative:
d
−θq−1 e−θ
Ψq (θ) =
< 0,
dθ
Γ(q)

∀q ∈ N, and θ > 0.

(7)

The limiting behavior of Ψq (θ) is important. For a fixed q,
the function Ψq (θ) → 1 as θ → 0 and Ψq (θ) → 0 as
θ → ∞. While Ψ−1
q still exists in these situations because Ψq
is monotonically decreasing, for a given z the value Ψ−1
q (z)
could be numerically very difficult to evaluate. To characterize
the sets of θ and q that Ψq is (numerically) invertible, we
define the θ-admissible set and the q-admissible set.

4

IEEE TRANSACTIONS ON COMPUTATIONAL IMAGING, VOL. X, NO. X, XXX 20XX

Definition 1. The θ-admissible set and q-admissible set of the
incomplete Gamma function are
def

Θq = {θ | ε ≤ Ψq (θ) ≤ 1 − ε},
def

Qθ = {q | ε ≤ Ψq (θ) ≤ 1 − ε},
respectively, where 0 < ε <

1
2

(8)

is a constant.

More discussions of the incomplete Gamma function can
be found in the Supplementary Material.
Remark 1. In this paper, we assume that QIS is noise-free, i.e.,
the only source of randomness is the truncated Poisson random
variable. In real sensors, there will be readout noise, photoresponse non-uniformity caused by conversion gain variation,
dark count rate (a.k.a. dark current), optical crosstalk and
electronic crosstalk. See [42] for details.
III. O PTIMAL T HRESHOLD : T HEORY
A. Image Reconstruction by MLE
We begin the optimal threshold design by discussing image reconstruction because the optimality of the threshold is
measured with respect to the reconstructed image. However,
since QIS is a new device, the number of reconstruction
methods is limited. A few examples that can be found in the
literature are the gradient descent [6], dynamic programming
[43], ADMM [7], and Transform-Denoise method [9], and
neural network [44]. In this paper, we shall focus on the
maximum likelihood estimation (MLE) approach as it provides
closed-form expressions.
Given B, MLE solves the following optimization problem:
(a)

cb = argmax
c

(b)

= argmax
c

−1
TY
−1 M
Y

P[Bm,t = 1 ; θm ]bm,t

t=0 m=0

−1 n
T
−1 M
X
X
t=0 m=0

× P[Bm,t = 0 ; θm ]1−bm,t

bm,t log(1 − Ψq (θm ))

o
+ (1 − bm,t ) log Ψq (θm ) ,

(9)

subject to the constraint that θ = αGc. Here, the right
hand side of (a) is the likelihood function of a Bernoulli
random variable, and (b) follows from taking the logarithm.
With the G defined in (2), we can partition B into N blocks
{B1 , . . . , BN } where each block is
def

Bn = {BKn+k,t | k = 0, . . . , K − 1, t = 0, . . . , T − 1}.

Then, the pixel b
cn can be estimated as follows.

Proposition 1 (Closed-form ML Estimate). The solution of
the MLE in (9) is


Sn
K −1
1−
,
(10)
b
cn = Ψq
α
KT
def PT −1 PK−1
where Sn = t=0
k=0 BKn+k,t is the sum of bits in the
n-th block Bn .
Proof: See [9].

B. Signal-to-Noise Ratio of ML Estimate
In order to determine the optimal threshold, we need to
quantify the performance of the ML estimate. The performance
metric we use is the signal-to-noise ratio of the ML estimate
at every pixel b
cn . Considering each b
cn individually is allowed
here because they are independently determined according to
(10). For notation simplicity we drop the subscript n in the
subsequent discussions.
Definition 2. The signal-to-noise ratio (SNR) of the ML
estimate b
c is defined as
def

SNRq (c) = 10 log10

c2
,
E[(b
c − c)2 ]

(11)

where the expectation is taken over the probability mass
function of the binary measurements in (6).
The difficulty of working with SNRq (c) is that it does not
have a simple closed-form expression. In view of this, Lu [17]
showed that the SNR is asymptotically linear to the log of the
Fisher Information.
Proposition 2. As KT → ∞,


SNRq (c) ≈ 10 log10 c2 Iq (c) + 10 log10 KT,
(12)
i
h
2
def
where Iq (c) = EB −∂
∂c2 log P(B = b; θ) is the Fisher Information measuring the amount of information that the random
variable B carries about the unknown value c.
Proof: See [17].
While the asymptotic result shown in Proposition 2 has
significantly simplified the SNR, we still need to determine
the Fisher Information. The following proposition gives a new
result of the Fisher Information with arbitrary q.
Proposition 3. The Fisher Information Iq (c) of the probability
mass function in (6) under a threshold q is:
2q−2
αc
 α 2
e−2( K ) αc
K

 .
Iq (c) =
(13)
K
Γ2 (q)Ψq αc
1 − Ψq αc
K
K

Proof: See Appendix A-A.
Substituting (13) into (12), we observe that the SNR can be
approximated as
2q
αc
KT e−2( K ) αc
K

 , (14)
SNRq (c) ≈ 10 log10
1 − Ψq αc
Γ(q)2 Ψq αc
K
K

which is characterized by the unknown pixel value c, the
threshold q, the spatial oversampling ratio K and the number
of temporal measurements T . To understand the behavior of
(14), we show in Figure 4 SNRq (c) as a function of c for
different thresholds q ∈ {1, . . . , 16}. For a fixed q, SNRq (c)
is a convex function with a unique maximum. The goal of
optimal threshold design is to determine a q which maximizes
SNRq (c) for a fixed c.
Remark 2. The SNRq (c) in (14) can also be derived from
a concept in the device literature called the exposure-referred
SNR [45]. See Supplementary Material for discussions.

ELGENDY-CHAN: OPTIMAL THRESHOLD DESIGN FOR QUANTA IMAGE SENSOR
40

IV. O PTIMAL T HRESHOLD : P RACTICE

q increases

20

SNRq (c)

0
-20
-40
-60
-80
-100
10−6

10−4

10−2

100

c

Fig. 4.
SNRq (c) for different thresholds q ∈ {1, . . . , 16}. In this
experiment, we set α = 400, K = 4, and T = 30. For fixed q, SNRq (c) is
always a convex function.

C. Oracle Threshold
We now discuss the optimal threshold design in the oracle
setting. We call the result oracle because the optimal threshold
depends on the unknown pixel intensity c. The practical
threshold design scheme will be discussed in Section IV.
Using the definition of the signal-to-noise ratio, the optimal
threshold is determined by maximizing SNRq (c) with respect
to q:
q ∗ = argmax SNRq (c) = argmax log(c2 Iq (c)).
q∈N

5

(15)

q∈N

The second equality follows from Proposition 2. Substituting
(13) yields an expression of the right hand side of (15). To
further simplify the expression we derive the following lower
bound.
Proposition 4. The function log(c2 Iq (c)) is lower bounded as
follows.


αc
αc
log(c2 Iq (c)) ≥ 2 log 2 −
+ q log
− log Γ(q) .
K
|
{z K
}
def

=Lq (c)

Proof: See Appendix A-B.
Using this lower bound, we can derive the optimal threshold
q as follows 1 .
Proposition 5. The optimal threshold q ∗ (c) is
j αc k
q ∗ (c) = argmax Lq (c) =
+ 1,
K
q∈N

(16)

where b·c denotes the flooring operator that returns the largest
integer smaller than or equal to the argument.
Proof: See Appendix A-C.
The result of Proposition 5 is important as it states that
the oracle threshold is exactly the same as the light intensity
αc/K. The flooring operation and the addition of a constant
1 are not crucial here because they are only used to ensure
that q is an integer. In [18], a special where α = 1 was
demonstrated experimentally. Proposition 5 now provides a
theoretical justification.
1 Straightly speaking, the result shown in Proposition 5 is a “near-optimal”
result because we are minimizing the lower bound. From our experience,
the gap between the near-optimality and the exact optimality is typically
insignificant.

The oracle threshold derived in the previous section provides
a theoretical foundation but is practically infeasible as it
requires knowledge of the ground truth c. In this section,
we present an alternative solution by relaxing the optimality
criteria. Our strategy is to consider a set of thresholds which
are close to the oracle threshold q ∗ (c), and show that they
are asymptotically unbiased when the number of observed bits
approaches infinity (Section IV.A). This result will allow us to
characterize the estimate b
c (Section IV.B). We will then show
that there exists a phase transition region where the asymptotic
unbiasedness is maintained as q stays within a certain range
around q ∗ (c), and is lost rapidly as q falls outside this range
(Section IV.C - IV.D). Based on these observations, we will
present a practical threshold update scheme (Section IV.E).
A. Asymptotic Unbiasedness
In order to derive an alternative threshold that does not
require the ground truth, we start by reconsidering the ML
estimate b
c in Proposition 1. For a spatial-temporal block
B = {Bk,t | 0 ≤ k < K − 1, 0 ≤ t < T − 1}, the ML
estimate b
c satisfies the condition
 
S
αb
c
=1−
,
(17)
Ψq
K
KT
P
where S = k,t Bk,t is the sum of bits in B. The right hand
side of this equation is an important quantity. We denote it as
S
.
(18)
KT
In the device literature (e.g., [45]), the term 1−γq (c) is known
as the bit-density as it is the proportion of ones in B. Note
that γq (c) is a random variable because S is the sum of KT
i.i.d. random binary bits. Therefore, if we want to understand
(17), we must first derive the the mean and variance of γq (c).
def

γq (c) = 1 −

Proposition 6. The mean and variance of γq (c) are
 αc 
E[γq (c)] = Ψq
, and
K

 αc i
1
αc  h
Var[γq (c)] =
Ψq
1 − Ψq
,
KT
K
K
respectively.

(19)

Proof: See Appendix A-D.
We can now look at the asymptotic behavior of γq (c) to see
if it offers any insight about the optimal threshold. Applying
the strong law of large number to S/KT , we can show that
as KT → ∞,
a.s.

γq (c) = 1 − S/KT → 1 − E[Bk,t ] = Ψq (αc/K).

(20)

Going back to (17)-(18), the ML estimate b
c should have the
expectation:

(a) K  −1
E[ b
c] =
E Ψq (γq (c))
α
 αc  (c)
(b) K −1
→ Ψq Ψq
= c.
(21)
α
K
where (a) follows from the definition of b
c, (b) follows from
(20), and (c) holds because Ψq and Ψ−1
cancels
each other.
q

IEEE TRANSACTIONS ON COMPUTATIONAL IMAGING, VOL. X, NO. X, XXX 20XX

What is the implication of (21)? It shows that the ML
estimate b
c is asymptotically unbiased. That is, as the number of
independent measurements grows, the estimate b
c approaches
to the ground truth c. In other words, as long as KT is large
enough, the random variable b
c would be an accurate estimate
of the ground truth. How can this be used to determine the
threshold q? Let us look at Qθ .
B. Set of Admissible Thresholds Qθ

The result in (17)-(21) shows that for a given S (or
equivalently γq (c)), the ML estimate can be found by
K
(γq (c)) .
(22)
b
c = Ψ−1
α q
When this happens, the b
c given by (22) is asymptotically
unbiased. However, the inversion Ψ−1
is not always allowed.
q
There is a set of q’s that can make Ψq invertible, which
is defined as Qθ in Definition 1. The following proposition
relates Qθ to γq (c).
Proposition 7. Let 0 < δ < 1 be a constant. Then, for any
(
1
1 )
  KT
  KT
δ
δ
def
≤ Ψq (θ) ≤
, (23)
q ∈ Qθ = q 1 −
2
2

the random variable γq (c) will not attain 0 or 1 with probability at least 1 − δ, i.e.,
P[0 < γq (c) < 1] > 1 − δ.
In this case, the ML estimate b
c is uniquely defined by (22).

Proof: See Appendix A-E.
Before we proceed, let us look at some rough magnitude of
the parameters in the following example.

Example 2. Let the ground truth pixel value be c = 0.5. The
sensor parameters are set as T = 50, K = 4, α = 300.
For a constant δ = 2 × 10−4 , the tolerance level is ε =
1−(δ/2)1/KT = 0.045. Therefore, as long as q ∈ {q | 0.045 ≤
Ψq (θ) ≤ 1 − 0.045}, which is the set {q | 28 ≤ q ≤ 48}, the
probability that γq (c) equals to 0 or 1 is upper bounded by
δ = 2 × 10−4 .
C. Gap between Qθ and q ∗

The result in the previous subsection shows that as long as
q ∈ Qθ , the ML estimate is asymptotic unbiased. However,
how is a q ∈ Qθ compared to the oracle threshold q ∗ ? We
answer this question in three parts.
First, does an asymptotically unbiased estimate maximize
the SNR? The answer is no, because Proposition 5 states that
if q ∗ is the optimal threshold, then SNRq∗ (c) ≥ SNRq (c) for
any q 6= q ∗ . Therefore, moving from the exact optimal q ∗
to an asymptotically unbiased threshold is a relaxation of the
optimality criteria.
If asymptotic unbiasedness is a relaxed optimality criteria,
how much SNR drop will there be if we choose a q ∈ Qθ
but not necessarily q = q ∗ ? We show in Figure 5 the plot of
a typical experiment with setup discussed in Example 2. As
shown in the figure, the green zone is the set Qθ = {q | 28 ≤

q ∈ Qθ
q∈
/ Qθ
Exact E[ b
c ]/c
Monte-Carlo E[ b
c ]/c
Exact 1 − E[γq (c)]
Monte-Carlo 1 − E[γq (c)]

1
0.8
0.6
0.4
0.2
0

10

20

30

40

50

60

70

80

90

100

70

80

90

100

Threshold q
40

SNR in dB

6

30
20
10
0

10

20

30

40

50

60

Threshold q

Fig. 5. Phase transition of the ML estimate and its relationship to the average
bit density 1 − E[γq (c)]. The red region is where it is impossible to recover
c, whereas the green region is where we can have perfect recovery.

q ≤ 48}, or equivalently Qθ = {q | 0.045 ≤ Ψq (θ) ≤ 0.9955}.
For any q in this Qθ , the reconstruction has a SNR at least
30dB. If we further tighten Qθ so that Qθ = {q | 35 ≤ q ≤
42}, or equivalently Qθ = {q | 0.25 ≤ Ψq (θ) ≤ 0.6}, the SNR
stays in the range 36.15dB ≤ SNRq (c) ≤ 36.65dB, which is
reasonably narrow.
How tight should Qθ be? Ideally we want Qθ to be as tight
as possible. But knowing the fact that the incomplete Gamma
function has a rapid transition (See the black line in Figure 5),
Qθ can be much wider. In fact, we can choose Qθ such that
1 − γq (c) stays close to 0.5, so that we are guaranteed to
obtain a near optimal threshold. From an information theoretic
point of view, 1 − γq (c) ≈ 0.5 is where the bit density attains
the maximum information — if q is too high then most bits
become 0 whereas if q is too low then most bits become 1. It
is maximum when q leads to 50% zeros and 50% ones. 2
D. Phase Transition Phenomenon
We can now point out a very interesting phenomenon in
Figure 5. In the upper plot of Figure 5 we show two sets of
curves: blue curves (solid and dotted), and black curves (solid
and dotted). The black curves represent the ratio E[ b
c ]/c, and
the black curves represent the average bit density 1−E[γq (c)].
For both sets of curves, we use dotted lines to illustrate the
Monte-Carlo simulation using 10,000 random samples, where
each sample refers to a spatial-temporal block Bn containing
KT = 200 binary bits. Notice that these dotted lines overlap
exactly with their expectations, and hence (17)-(21) are valid.
Let us take a closer look at the blue curve E[ b
c ]/c. Let
Qθ = {q | qL ≤ q ≤ qH }, where qL and qH are the smallest
and the largest integers in Qθ respectively. There are three
distinct phases:
• When q < qL , the threshold is low and so most bits
become 1. Therefore, γc (q) → 0 and hence b
c → ∞. Thus,
E[ b
c ]/c → ∞ as q decreases.

2 The exact optimal value of 1 − γ (c) at q ∗ is slightly lower than 0.5 due
q
to the nonlinearity of the Gamma function. See Supplementary Material for
additional discussion.

ELGENDY-CHAN: OPTIMAL THRESHOLD DESIGN FOR QUANTA IMAGE SENSOR

Algorithm 1 Bisection Threshold Update Scheme
Initial thresholds qA and qB such that 1 − γqA > 0.5 and
1 − γqB < 0.5.
Compute qM = d(qA + qB )/2e, where d·e denotes the
ceiling operator.
while |γqM − 0.5| < tol do
If γqM < 0.5, then set qA = qM . Else, set qB = qA .
Compute qM = d(qA + qB )/2e.
end while
return qM

1
0.9
iter 2

0.8

Bit density 1 − γq (c)

0.7
iter 4
optimal

0.6
0.5
0.4
iter 3

0.3

7

iter 1

0.2
0.1
0

0

2

4

6

8

10

12

14

16

Threshold q

iter 1, 27.4 dB

iter 2, 37.1 dB

iter 3, 38.8 dB

iter 4, 39.1 dB

Fig. 6. The proposed bisection update scheme adjusts the threshold q such
that the bit density 1 − γq (c) approaches 0.5. The upper graph illustrates
the bisection steps. Bottom row shows cropped patches from reconstructed
images using threshold maps at different iterations and the PSNRs.

• When q > qH , the threshold high and so most bits become
0. Therefore, γc (q) → 1 and hence b
c → 0. Thus, E[ b
c ]/c →
0 as q increases.
• When qL ≤ q ≤ qH , the ML estimate b
c is asymptotically
unbiased. Therefore, E[ b
c ]/c = 1.
Essentially, Figure 5 demonstrates a phase transition behavior
of the threshold. Such phase transition exists because Ψq is
only invertible when q ∈ Qθ .
E. Bisection Threshold Update Scheme
Now we present a practical threshold update scheme. As
we discussed in Section IV.C, the oracle threshold q ∗ can be
obtained when bit density γq (c) is close to 0.5. Therefore, a
practical procedure to determine q is to sweep through a range
of q until the bit density reaches 0.5. To achieve this objective,
we propose a bisection method illustrated in Figure 6 and
Algorithm 1. Starting with initial thresholds qA and qB , we
check whether the bit density satisfies 1 − γqA > 0.5 and
1 − γqB < 0.5. If this is the case, then we find a mid point
qM = (qA +qB )/2 and check whether 1−γqM is greater or less
than 0.5. If 1 − γqM > 0.5, we replace qA by qM , otherwise
we replace qB by qM . The process repeats until 1 − γqM is
sufficiently close to 0.5.
In our proposed threshold update scheme, we assume that
the image has been partitioned into N blocks {Bn | n =
0, . . . , N − 1}. Each Bn contains KT binary bits and is
used to estimate one pixel value cn . This setting results in
N different thresholds, one for every pixel. To generalize the
setting, it is also possible to allow multiple pixels to share a
common threshold. Figure 7 shows an example. The advantage
of sharing a threshold for multiple pixels is that circuits
associated with the sensor can be simplified. In terms of
performance, since neighboring pixels are typically correlated,
sharing the threshold causes little drop in the resulting SNR.

Fig. 7. Concept of shared thresholds. (Left) binary measurements, spatial
oversampling K = 3 × 3, Temporal oversampling T = 5 . (Right) Threshold
map, one threshold value is shared by 6 × 6 jots.

The price that the proposed bisection algorithm has to pay
is the number of frames it requires to determine a good q.
For every evaluation of γqM , the sensor has to physically
acquire one frame and compute the bit density in each of the
N blocks. Therefore, the more bisection steps we need, the
more frames that the sensor has to physically acquire. The rate
of convergence of the proposed method and existing methods
will be compared in Section V.
F. Extension to High Dynamic Range
While QIS is a photon counting device, it is designed to
count a few photons to keep the full-well capacity small, e.g.
20 photoelectrons as reported in [46]. Therefore, for practical
imaging tasks, we need to extend the dynamic range for QIS.
There are two ways to enable dynamic range extension:
• Bright Scenes: Reduce Duty Cycle. In the signal processing block diagram shown in Figure 2, we can replace
the constant α by a fraction as ατ , where 0 ≤ τ ≤ 1
determines the ratio between the actual integration time
and the readout scan time. It can also be referred to the
shutter duty cycle because the shutter is opened to collect
photons during this proportion of time [47]. For very
bright scenes, a low duty cycle will prevent QIS from
saturating early.
• Dark Scenes: Multiple Measurements. For dark scenes,
multiple measurements can be taken to ensure enough
photons over the measurement period. This, however, is
different from conventional HDR imaging. In conventional HDR imaging, the multiple shots are taken at different shutter speeds, e.g., 1/8192, 1/2048, 1/512, 1/128,
1/32, 1/8, 1/2 seconds [48], which is redundant. QIS’s
multiple shot functions more similar to burst photography

8

IEEE TRANSACTIONS ON COMPUTATIONAL IMAGING, VOL. X, NO. X, XXX 20XX
60
50

q
q
q
q
q

= 1, τ = 1
= 1, τ = 0.2
= 1, τ = 0.04
= 1, τ = 0.008
= 1, overall SNR

q
q
q
q
q

= 25, τ = 1
= 25, τ = 0.2
= 25, τ = 0.04
= 25, τ = 0.008
= 25, overall SNR

q∗ (c), τ = 1
q∗ (c), τ = 0.2
q∗ (c), τ = 0.04
q∗ (c), τ = 0.008
q∗ (c), overall SNR

40

SNRq (c)

30
20
16 dB

54 dB

10
0
-10
-20
-30
10−2

10−1

100

101

102

103

V. E XPERIMENTAL R ESULTS
In this section we evaluate the proposed threshold update
scheme by comparing it with existing methods. We consider
two evaluation metrics: (1) convergence rate of the threshold
update methods; (2) quality of the reconstructed images. For
reconstruction evaluation, we create our own Purdue dataset
comprising 77 images captured by a Canon EOS Rebel T6i
camera. For HDR imaging, we use the HDR-Eye dataset by
Nemoto et al. [51], [52]. In all experiments, we fix the spatial
over-sampling factor as K = 4 × 4 = 16, and number of
temporal frames as T = 13. The maximum threshold level is
set as qmax = 16 to ensure that it is realistic for today’s QIS.

104

θ = αc/K

Fig. 8. SNR in dB vs. exposure θ for HDR imaging mode obtained by
fusion of frames with shutter duty cycles τ ∈ {1, 0.2, 0.04, 0.008}. Three
scenarios are shown: constant threshold with q = 1 (black), q = 25 (red) and
an optimal spatially varying threshold (blue).

[49]. The amount of acquisition time is significantly less
than the conventional HDR imaging.
These two methods can be used for any threshold scheme,
including ours and others. The benefit of using our proposed
threshold scheme is that it supports a much wider dynamic
range extension. In Figure 8, we illustrate the total dynamic
range that can be covered using 4 multiple measurements at
duty cycles τ = 1, τ = 0.2, τ = 0.04, and τ = 0.008. The
maximum threshold level is qmax = 25, and the minimum
threshold level is qmin = 1. It can be seen from the figure
that with the optimal threshold q ∗ , the dynamic range is
significantly more than the non-optimal ones. In particular,
we observe a 16dB and a 54dB improvement compared to
qmin = 1 and qmax = 25, respectively. Experimental results
will be shown in Section V.C.
G. Hardware Consideration
Concerning the hardware implementation, we anticipate that
future QIS will be equipped with per-pixel FPGAs to perform
the proposed threshold update scheme. On-sensor FPGA is
an actively developing technology. For example, MIT Lincoln
Lab’s digital focal plane array can achieve on-sensor image
stabilization and edge detection [50] . For QIS threshold
update, the complexity is low because we are only counting
the number of ones in the bisection. More specifically, in order
to perform the bisection, we only need
K additions to compute
P
PK−1
K−1
b
;
one
comparison
k=0 Kn+k,t
k=0 bKn+k,t ≥ 0.5; one
addition and one multiplication (with a constant 0.5) to update
the threshold qM = d(qA +qB )/2e. The dominating factor here
is the K additions, which can be implemented efficiently by
shifting bits in a buffer.
We should also point out that the proposed bisection method
can be flexibly adjusted spatially and temporally for different
hardware configurations. For example, we can use a spatialtemporal window 4 × 4 × 1 for low-resolution high-speed
imaging, or 1 × 1 × 16 for high-resolution low-speed imaging. This flexibility offers additional advantages of QIS over
conventional CCD and CMOS cameras.

A. Convergence
We compare the proposed threshold update scheme with
the Markov Chain (MC) adaptation proposed by Hu and Lu
[18]. The Markov Chain adaptation models the threshold as
a variable with 2L states. These 2L states can be regarded
as 2L steps before reaching to the next threshold level. The
probability of changing from one state to another is controlled
by a parameter 1 − β with 0 < β < 1. When a bit arrives,
the state will be updated (increased or decreased) or will
remain unchanged. Once the state is increased by 2L times,
the threshold will be increased by one.
When comparing Markov Chain adaptation with the proposed bisection algorithm, one should be aware of the difference between the two methods. Markov Chain adaptation
is a per-jot update scheme whereas the proposed bisection
algorithm is a per-pixel update scheme. For a pixel with
K × K jots, Markov Chain adaptation needs K 2 iterations
to update the threshold sequentially. In contrast, the proposed
bisection algorithm updates a common threshold for all K 2
jots simultaneously. Thus in practice our bisection algorithm
is significantly less complex to implement in hardware than
the Markov Chain. In order to take the different forms of
updates into account, we treat the K 2 iterations of Markov
Chain adaptation as one “major iteration” and compare it with
the one bisection step of the proposed algorithm.
The first comparison we make is to check the threshold
at different jots. Figure 9 shows the results of three typical
runs with underlying optimal thresholds q ∗ = 1, 8, 16. In this
experiment, we generate 100 random binary blocks of size K×
K and estimate the threshold at each major iteration. We report
the average of these 100 estimates to minimize the randomness
of the data. The results show that one iteration of the proposed
bisection algorithm works as good as the K 2 iterations of the
Markov Chain adaptation. In some cases, Markov Chain tends
to oscillate whereas the bisection result is stable.
The second comparison we make is to check how close
the estimated threshold is compared to the optimal threshold.
The optimal threshold q ∗ is obtained using the oracle scheme.
In Figure 10, we plot the mean squared error between the
estimated threshold and the oracle threshold. For fairness, we
show the results of the MSE averaged over the 77 images of
our dataset, and 50 random samples per image. One threshold
is shared by K × K jots, and each K × K jots correspond

ELGENDY-CHAN: OPTIMAL THRESHOLD DESIGN FOR QUANTA IMAGE SENSOR
[18], L = 4

Threshold

16

[18], L = 3

[18], L = 2

Proposed

q∗ (c)

14

9

TABLE II
AVERAGE PSNR AND S TANDARD DEVIATION OF 77 RECOVERED IMAGES
USING DIFFERENT Q- MAPS AND 50 RANDOM SAMPLES .

12

Configuration

Average
PSNR

Std

Uniform Threshold

q=1
q=5
q = 10
q = 16

10.30
28.80
23.22
12.95

0.01
0.04
0.02
0.01

Conditional Reset [21]

Ascending q sequence
Descending q sequence

23.77
24.95

0.52
0.53

Proposed Method

2K 2 × 2K 2
K2 × K2
K×K

30.14
31.18
32.78

0.06
0.06
0.02

10
8

Threshold

9

1

2

3

4

5

6

7

8

9

10

11

8
7
6
1

2

3

4

5

6

7

8

9

10

11

Threshold

8
6
4
2
1

2

3

4

5

6

7

8

9

10

11

Major Iteration

Fig. 9. Convergence of the threshold at 3 jots. Each curve is averaged
over 100 random samples. The red curve indicates the proposed bisection
method. The blue curves are the Markov chain adaptation [18] with β =
0.25. Note that one major iteration of Markov Chain adaptation corresponds
to K 2 sequential updates, and one major iteration of the bisection method
corresponds to a single update to K 2 jots simultaneously.

Mean Square Error

15

[18], L = 4
[18], L = 3
[18], L = 2
Proposed

10

5

0

1

2

3

4

5

6

7

8

9

10

11

Major Iteration

Fig. 10. Mean square error between the estimated threshold and the ideal
oracle threshold. Each curve is averaged over 50 random samples and 77
images. The red curve indicates the proposed bisection method. The blue
curves are the Markov chain adaptation [18] with β = 0.25.

to one pixel. The result is consistent with the ones shown in
Figure 9.
B. Image Reconstruction Quality
The convergence comparison in the previous subsection is
only useful to compare threshold update methods that actually
return a threshold. In the QIS literature, there are methods
that implicitly update the threshold, e.g., the conditional reset
method [21]. For comparison with these methods, we have
to compare the quality of the image reconstructed from the
binary raw data. The image reconstruction is done using the
closed-form ML estimate in Section III-A.
We consider three classes of methods:
• Uniform Threshold. Uniform threshold is commonly used
in the device literature [5]–[7]. A uniform threshold is a
single threshold applied to all pixels in the image. In this
experiment, we consider the following choices of uniform
thresholds: q = 1, q = 5, q = 10 and q = 16.
• Conditional Reset [21]. Conditional reset counts the number of photons and is reset when it is above the threshold.

The threshold in conditional reset is sequentially increasing or decreasing. The reconstructed image is obtained
by digitally integrating the raw binary frames.
• Proposed Method. As we discussed in Section IV-E, the
proposed method can be implemented to let multiple
pixels share a common threshold. Thus, in this experiment we consider three sharing strategies: (1) Share a
threshold between a neighborhood of K×K jots (i.e., one
threshold for one pixel); (2) Share a threshold between
a neighborhood of K 2 × K 2 jots (i.e., one threshold
for K × K pixel); (3) Share a threshold between a
neighborhood of 2K 2 × 2K 2 jots (i.e., one threshold for
2K × 2K pixels).
The result of the experiment is shown in Table II. The PSNR
values reported are averaged over 77 images in our dataset.
Each image generates 50 random realizations, and the PSNR
of an image is averaged over these 50 random realizations
to minimize the randomness. As shown in the table, while
conditional reset generally performs better than a uniform
threshold, it performs significantly worse than the proposed
threshold update scheme.
C. Influence of QIS Threshold on HDR Imaging
Since QIS does not have sufficient full well capacity to
accumulate photons for HDR imaging, we apply the dynamic
range extension method discussed in Section IV-F. When
different threshold schemes are used, the reconstructed HDR
images will be affected. The objective of this experiment is to
evaluate the influence of the threshold in HDR imaging.
In this experiment, we consider the HDR-Eye image dataset
[51], [52]. Each HDR image in this dataset contains 9 images
acquired at different exposure settings (−2.7,−2, −1.3, −0.7,
0, 0.7, 1.3, 2, and 2.7 EV). A snapshot of these images
are shown in Figure 11. From each exposure, we simulate
the photon counts resulting from the luminance channel. The
sensor gain is set as α = K 2 (qmax − 1) to ensure proper
number of photons, where K = 4×4 = 16 and qmax = 16. On
the reconstruction side, we reconstruct the 9 images using the
MLE discussed in Section III-A. Tone mapping and exposure
fusion [53] are applied to the 9 imags to generate an HDR
image. As a reference, we apply the same tone mapping and

10

IEEE TRANSACTIONS ON COMPUTATIONAL IMAGING, VOL. X, NO. X, XXX 20XX

Fig. 11.

Bracketed images with different exposure settings. From Left to Right: −2.7,−2, −1.3, −0.7, 0, 0.7, 1.3, 2, and 2.7 EV.

Fig. 12.

The reconstructed HDR images using different thresholds. See supplementary material for additional results.

q = 1, PSNR = 17.94 dB

q = 16, PSNR = 20.77 dB

fusion algorithm to the 9 ground truth images. PSNR between
the reference and the estimated is then recorded.
The result of this experiment is shown in Figure 12.
With the proposed threshold update scheme, the reconstructed
images achieve the highest PSNR value and visual quality.
When q = 1, which is too low, the image appears underexposed. When q = 16, which is too high, the image appears
over-exposed. The spatially varying property of the proposed
method mitigates the issue by allowing multiple thresholds.
In practice, one would typically add image denoisers to
handle the randomness in the ML estimate and potentially
other types of noise. This can be done using methods such
as [9]. In HDR literature, there are also optical approaches
that reduce the number of exposures, e.g., [54], [55]. These
techniques are complementary to QIS, because QIS is a sensor
of similar functionality of a CMOS. Thus optical techniques
can always be added.
VI. C ONCLUSION
Quanta Image Sensor is a new image sensor for high speed,
high resolution and high dynamic range imaging. The sensor
has a threshold which needs to be carefully adjusted so that the
dynamic range can be maximized. We studied the threshold
design problem by establishing several theoretical results.
First, we showed that an oracle threshold can be obtained
assuming that we know the underlying pixel value. Our result
showed that the oracle threshold must match with the pixel
value in order to maximize the signal-to-noise ratio. Second,
we showed that around the oracle threshold, there exists a
set of thresholds that can produce asymptotically unbiased
estimates of the pixel value. Within this set of threshold, the
signal-to-noise ratio stays very close to the oracle case. Third,
we developed a bisection method to update the threshold
scheme. We also discussed how QIS can be used in HDR
imaging, and its advantages compared to standard sensors.
Experimental results showed the effectiveness of our proposed
approach compared to the standard approach that uses uniform
threshold for all pixels.

Proposed, PSNR = 31.46 dB

ACKNOWLEDGMENT
The authors thank Professor Eric Fossum, Jiaju Ma and
Saleh Masoodian at Dartmouth College for many insightful
discussions about the physics and circuits of QIS.
A PPENDIX A
A. Proof of Proposition 3
The Fisher Information metric is defined as:

 2
−∂
def
log
P(B
=
b;
θ,
q)
,
Iq (c) = EB
∂c2

(24)

where θ = αc/K. Using the chain rule, we can derive the
Fisher Information as follows

 2
 α 2
−∂
Iq (c) =
EB
log
P(B
=
b;
θ,
q)
.
(25)
K
∂θ2

The expectation can be calculated as follows

 α 2  −∂ 2
Iq (c) =
log P(B = 1; θ, q) P(B = 1; θ, q)
K
∂θ2


 α 2 −∂ 2
+
log
P(B
=
0;
θ,
q)
P(B = 0; θ, q)
K
∂θ2
(26)
Using (7) to differentiate the 1st term, we get:
∂2
∂2
log
P(B
=
1;
θ,
q)
=
log (1 − Ψq (θ))
∂θ2
∂θ2
0
2
R (1 − Ψq (θ)) − R /Γ(q)
,
=
2
Γ(q) (1 − Ψq (θ))

(27)

where R = e−θ θq−1 and R0 = ∂R/∂θ. Similarly, the second
term is
∂2
∂2
log P(B = 0; θ, q) = 2 log Ψq (θ)
2
∂θ
∂θ
R0 Ψq (θ) + R2 /Γ(q)
=−
.
2
Γ(q) (Ψq (θ))

(28)

ELGENDY-CHAN: OPTIMAL THRESHOLD DESIGN FOR QUANTA IMAGE SENSOR

E. Proof of Proposition 7

Substitute (27) and (28) in (26) yields
 α 2 h R0 Γ(q)(1 − Ψ (θ)) − R2
q
−
Iq (θ) =
K
Γ2 (q) (1 − Ψq (θ))
R0 Γ(q)Ψq (θ) + R2 i
+
Γ2 (q)Ψq (θ)
 α 2
e−2θ θ2q−2
=
.
2
K
Γ (q)Ψq (θ) (1 − Ψq (θ))

The probability P[0 < γq (c) < 1] can be evaluated by
checking the complement when γq (c) = 0 or γq (c) = 1:
P[0 < γq (c) < 1] = 1 − P[γq (c) = 0] − P[γq (c) = 1]
= 1 − P[S = 0] − P[S = KT ]
(a)

= 1 − Ψq (θ)KT − [1 − Ψq (θ)]KT ,

B. Proof of Proposition 4
The lower bound is obtained by observing that the
product Ψq (θ) (1 − Ψq (θ)) attains its maximum value
when Ψq (θ) = 1/2. Substituting with the upper bound
Ψq (θ) (1 − Ψq (θ)) ≤ 1/4, we get:

 
e−2θ θ2q−2
αc 2
log(c2 Iq (c)) = log
K
Γ2 (q)Ψq (θ) (1 − Ψq (θ))
e−2θ θ2q
= log 2
Γ (q)Ψq (θ) (1 − Ψq (θ))
4e−2θ θ2q
≥ log
Γ2 (q)
= 2 log 2 − 2θ + 2q log θ − 2 log Γ(q)


αc
αc
= 2 log 2 −
+ q log
− log Γ(q) .
K
K
C. Proof of Proposition 5
Using the definition of Gamma function Γ(q) = (q − 1)!
and θ = αc
K , we can rewrite the lower bound in Proposition 4
as follows.
Lq (c) = 2 (log 2 − θ + q log θ − log(q − 1)!)
= 2 log 2 − θ + (q − 1) log θ + log θ − log
= 2 log 2 − θ +

q−1
X

k=1

log(θ/k) + log θ

!

q−1
Y

k=1

11

k

!

The only dependence on q is in the second term, so we take a
closer look at it. When q−1 < bθc, all summands log(θ/k) are
positive because k < bθc. Hence, the total sum increases by
increasing q. On the other hand, when q − 1 > bθc, we start to
add negative summands log(θ/k) because k > θ. Therefore,
the total sum decreases on increasing q − 1 over bθc. Thus,
maximum is obtained at q = bθc + 1 = b αc
K c + 1.
D. Proof of Proposition 6
def PT −1 PK−1
By definition, S = t=0
k=0 Bk,t is the summation of
KT independent i.i.d. Bernoulli random variables. Therefore,
def
S is a binomial random variable with parameters n = KT
def
and p = 1 − Ψ(αc/K). The mean and variance of a binomial
random variable is E[S] = np, and Var[S] = np(1 − p).
Therefore, we have
 αc 
E[S]
E [γq (c)] = 1 −
= Ψq
, and
KT
K
 αc  
 αc 
Var [S]
1
Var [γq (c)] =
=
Ψq
1 − Ψq
.
2
2
K T
KT
K
K

where (a) follows from the fact that S, which is a sum of i.i.d.
Bernoulli random variables, is a binomial random variable.
Let 0 < δ < 1. If
1
1
  KT
  KT
δ
δ
1−
≤ Ψq (θ) ≤
,
2
2

then we have

Ψq (θ)KT <

δ
δ
and [1 − Ψq (θ)]KT < .
2
2

Thus, it holds that
1 − Ψq (θ)KT − [1 − Ψq (θ)]KT > 1 − δ.
R EFERENCES
[1] E. R. Fossum, “What to do with sub-diffraction-limit (SDL) pixels?–
A proposal for a gigapixel digital film sensor (DFS),” in Proc. IEEE
Workshop CCDs Adv. Image Sensors, Sep. 2005, pp. 214–217.
[2] J. Ma, D. Hondongwa, and E. R. Fossum, “Jot devices and the quanta
image sensor,” in Proc. IEEE Int. Electron Devices Meeting, Dec. 2014,
pp. 10.1.1–10.1.4.
[3] J. Ma and E. R. Fossum, “A pump-gate jot device with high conversion
gain for a quanta image sensor,” IEEE J. Electron Devices Soc., vol. 3,
no. 2, pp. 73–77, Mar. 2015.
[4] J. Ma, L. Anzagira, and E. R. Fossum, “A 1 µm-pitch quanta image
sensor jot device with shared readout,” IEEE J. Electron Devices Soc.,
vol. 4, no. 2, pp. 83–89, Mar. 2016.
[5] F. Yang, Y. M. Lu, L. Sbaiz, and M. Vetterli, “An optimal algorithm for
reconstructing images from binary measurements,” in Proc. IS&T/SPIE
Electronic Imaging, Computational Imaging VIII, Jan. 2010, vol. 7533,
pp. 75330K–75330K–12.
[6] F. Yang, Y. M. Lu, L. Sbaiz, and M. Vetterli, “Bits from photons:
Oversampled image acquisition using binary poisson statistics,” IEEE
Trans. Image Process., vol. 21, no. 4, pp. 1421–1436, Apr. 2012.
[7] S. H. Chan and Y. M. Lu, “Efficient image reconstruction for gigapixel
quantum image sensors,” in Proc. IEEE Global Conf. Signal and
Information Processing (GlobalSIP’14), Dec. 2014, pp. 312–316.
[8] O. A. Elgendy and S. H. Chan, “Image reconstruction and threshold
design for quanta image sensors,” in Proc. IEEE Int. Conf. Image
Process. (ICIP’16), Sep. 2016, pp. 978–982.
[9] S. H. Chan, O. A. Elgendy, and X. Wang, “Images from bits: Noniterative image reconstruction for quanta image sensors,” MDPI Sensors,
vol. 16, no. 11, pp. 1961, 2016.
[10] E. R. Fossum, J. Ma, and S. Masoodian, “Quanta image sensor:
concepts and progress,” in Proc. SPIE 9858, Advanced Photon Counting
Techniques X, 2016, vol. 9858, pp. 985804–985804–14.
[11] I. M. Antolovic, S. Burri, C. Bruschini, R. Hoebe, and E. Charbon,
“Nonuniformity analysis of a 65k pixel CMOS SPAD imager,” IEEE
Trans. Electron Devices, vol. 63, no. 1, pp. 57–64, Jan. 2016.
[12] “Andor ixon ultra 888 specifications,” http://www.andor.com/pdfs/
specifications/iXon Ultra 888 EMCCD Specifications.pdf, Accessed:
2016-08-03.
[13] B. F. Aull, D. R. Schuette, D. J. Young, D. M. Craig, B. J. Felton, and
K. Warner, “A study of crosstalk in a 256×256 photon counting imager
based on silicon Geiger-mode avalanche photodiodes,” IEEE Sens. J.,
vol. 15, no. 4, pp. 2123–2132, Apr. 2015.
[14] N. A. W. Dutton, I. Gyongy, L. Parmesan, and R. K. Henderson, “Single
photon counting performance and noise analysis of CMOS SPAD-based
image sensors,” MDPI Sensors, vol. 16, no. 7, pp. 1122, Jul. 2016.
[15] S. Masoodian, A. Rao, J. Ma, K. Odame, and E. R. Fossum, “A 2.5 pj/b
binary image sensor as a pathfinder for quanta image sensors,” IEEE
Trans. Electron Devices, vol. 63, no. 1, pp. 100–105, Jan. 2016.

12

IEEE TRANSACTIONS ON COMPUTATIONAL IMAGING, VOL. X, NO. X, XXX 20XX

[16] F. Yang, Bits from Photons: Oversampled Binary Image Acquisition,
Ph.D. thesis, École Polytechnique Fédérale De Lausanne, 2012.
[17] Y. M. Lu, “Adaptive sensing and inference for single-photon imaging,”
in Proc. 47th Annual Conf. Inf. Sci. and Sys. (CISS’13), Mar. 2013, pp.
1–6.
[18] C. Hu and Y. M. Lu, “Adaptive time-sequential binary sensing for high
dynamic range imaging,” in Proc. SPIE Conf. Adv. Photon Counting,
2012, vol. 8375, pp. 83750A–1.
[19] T. Vogelsang and D. G. Stork, “High-dynamic-range binary pixel
processing using non-destructive reads and variable oversampling and
thresholds,” in Proc. IEEE Sens, Oct. 2012, pp. 1–4.
[20] T. Vogelsang, M. Guidash, and S. Xue, “Overcoming the full well
capacity limit: high dynamic range imaging using multi-bit temporal
oversampling and conditional reset,” in Proc. Int. Image Sens. Workshop
(IISW), Snowbird, UT, 2013.
[21] T. Vogelsang, D. G. Stork, and M. Guidash, “Hardware validated unified
model of multibit temporally and spatially oversampled image sensors
with conditional reset,” J. Electronic Imaging, vol. 23, no. 1, pp. 013021,
2014.
[22] J. Hynecek, “Impactron-a new solid state image intensifier,” IEEE Trans.
Electron Devices, vol. 48, no. 10, pp. 2238–2241, Oct. 2001.
[23] M. S. Robbins and B. J. Hadwen, “The noise performance of electron
multiplying charge-coupled devices,” IEEE Trans. Electron Devices, vol.
50, no. 5, pp. 1227–1232, May 2003.
[24] N. A. W. Dutton, I. Gyongy, L. Parmesan, S. Gnecchi, N. Calder, B. R.
Rae, S. Pellegrini, L. A. Grant, and R. K. Henderson, “A SPAD-based
QVGA image sensor for single-photon counting and quanta imaging,”
IEEE Trans. Electron Devices, vol. 63, no. 1, pp. 189–196, Jan. 2016.
[25] L. Liang, H. Shen, P. De Camilli, and J. S. Duncan, “A novel
multiple hypothesis based particle tracking method for clathrin mediated
endocytosis analysis using fluorescence microscopy,” IEEE Trans. Image
Process., vol. 23, no. 4, pp. 1844–1857, Apr. 2014.
[26] L. H. C. Braga, L. Gasparini, L. Grant, R. K. Henderson, N. Massari,
M. Perenzoni, D. Stoppa, and R. Walker, “A fully digital 8 × 16 SiPM
array for PET applications with per-pixel TDCs and real-time energy
output,” IEEE J. Solid-State Circuits, vol. 49, no. 1, pp. 301–314, Jan.
2014.
[27] S. P. Poland, N. Krstajić, J. Monypenny, S. Coelho, D. Tyndall, R. J.
Walker, V. Devauges, J. Richardson, N. Dutton, P. Barber, D. D. Li,
K. Suhling, T. Ng, R. K. Henderson, and S. M. Ameer-Beg, “A high
speed multifocal multiphoton fluorescence lifetime imaging microscope
for live-cell fret imaging,” Biomed. Opt. Express, vol. 6, no. 2, pp.
277–296, Feb. 2015.
[28] G. Grubbs, R. Michell, M. Samara, D. Hampton, and J.-M. Jahn, “A
synthesis of star calibration techniques for ground-based narrowband
electron-multiplying charge-coupled device imagers used in auroral
photometry,” J. Geophysical Research: Space Physics, vol. 121, no.
6, pp. 5991–6002, 2016, 2015JA022186.
[29] P. Seitz and A. J. Theuwissen, Single-photon imaging, vol. 160, Springer
Science & Business Media, 2011.
[30] L. J. Meng, “An intensified EMCCD camera for low energy Gamma
ray imaging applications,” IEEE Trans. Nuclear Science, vol. 53, no. 4,
pp. 2376–2384, Aug. 2006.
[31] D. Shin, F. Xu, D. Venkatraman, R. Lussana, F. Villa, F. Zappa, V. K.
Goyal, F. N. C. Wong, and J. H. Shapiro, “Photon-efficient imaging
with a single-photon camera,” Nature Communications, vol. 7, pp. –,
2016.
[32] I. Gyongy, T. Al Abbas, N. A.W. Dutton, and R. K. Henderson, “Object
tracking and reconstruction with a quanta image sensor,” in Proc. Int.
Image Sens. Workshop (IISW), 2017.
[33] S. Burri, Y. Maruyama, X. Michalet, F. Regazzoni, C. Bruschini, and
E. Charbon, “Architecture and applications of a high resolution gated
spad image sensor,” Opt. Express, vol. 22, no. 14, pp. 17573–17589,
July 2014.
[34] E Amri, Y. Felk, D. Stucki, J. Ma, and E. R. Fossum, “Quantum random
number generation using a quanta image sensor,” MDPI Sensors, vol.
16, no. 7, pp. 1002, Jun. 2016.
[35] E. R. Fossum, “The Quanta Image Sensor (QIS): Concepts and
Challenges,” in Proc. OSA Topical Mtg Computational Optical Sensing
and Imaging, Jul. 2011, paper JTuE1.
[36] J. Ma and E. R. Fossum, “Quanta image sensor jot with sub 0.3er.m.s. read noise and photon counting capability,” IEEE Electron Device
Letters, vol. 36, no. 9, pp. 926–928, Sep. 2015.
[37] L. Sbaiz, F. Yang, E. Charbon, S. Susstrunk, and M. Vetterli, “The
gigavision camera,” in Proc. 2009 IEEE Int. Conf. Acoustics, Speech
and Signal Processing (ICASSP’09), Apr. 2009, pp. 1093–1096.

[38] F. Yang, L. Sbaiz, E. Charbon, S. Süsstrunk, and M. Vetterli, “On
pixel detection threshold in the gigavision camera,” in Proc. IS&T/SPIE
Electronic Imaging, Digital Photography VI, 2010, vol. 7537, pp.
75370G–75370G–8.
[39] N. A. W. Dutton, L. Parmesan, A. J. Holmes, L. A. Grant, and R. K.
Henderson, “320 × 240 oversampled digital single photon counting
image sensor,” in Proc. Symp VLSI Circuits, Jun. 2014, pp. 1–2.
[40] I. M. Antolovic, S. Burri, R. A. Hoebe, Y. Maruyama, C. Bruschini,
and E. Charbon, “Photon-counting arrays for time-resolved imaging,”
MDPI Sensors, vol. 16, no. 7, pp. 1005, 2016.
[41] M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions:
with Formulas, Graphs, and Mathematical Tables, Number 55. Courier
Corporation, 1964.
[42] E. R. Fossum, “Photon counting error rates in single-bit and multi-bit
quanta image sensors,” IEEE J. Electron Devices Soc., vol. 4, no. 3, pp.
136–143, May 2016.
[43] F. Yang, L. Sbaiz, E. Charbon, S. Ssstrunk, and M. Vetterli, “Image
reconstruction in the gigavision camera,” in Proc. IEEE 12th Int. Conf.
on Computer Vision Workshops (ICCV Workshops), 2009, Sep. 2009,
pp. 2212–2219.
[44] R. A. Rojas, W. Luo, V. Murray, and Y. M. Lu, “Learning optimal
parameters for binary sensing image reconstruction algorithms,” in Proc.
IEEE Int. Conf. Image Process. (ICIP’17), Sep. 2017, pp. 2791–2795.
[45] E. R. Fossum, “Modeling the performance of single-bit and multi-bit
quanta image sensors,” IEEE J. Electron Devices Soc., vol. 1, no. 9, pp.
166–174, Sep. 2013.
[46] J. Ma, D. Starkey, A. Rao, K. Odame, and E. R. Fossum, “Characterization of quanta image sensor pump-gate jots with deep sub-electron
read noise,” IEEE J. Electron Devices Soc., vol. 3, no. 6, pp. 472–480,
Nov. 2015.
[47] E. R. Fossum, “Multi-bit quanta image sensors,” in Proc. Int. Image
Sens. Workshop (IISW), Jun. 2015, pp. 292–295.
[48] I. Sprow, D. Kuepper, Z. Barańczuk, and P. Zolliker, “Image quality
assessment using a high dynamic range display,” in Proc. 12th Congress
of the International Colour Association, Newcastle upon Tyne, England,
Jul. 2013.
[49] S. Hasinoff, D. Sharlet, R. Geiss, A. Adams, J. T. Barron, F. Kainz,
J. Chen, and M. Levoy, “Burst photography for high dynamic range and
low-light imaging on mobile cameras,” ACM Transactions on Graphics
(Proc. SIGGRAPH Asia), 2016.
[50] K. I. Schultz, M. W. Kelly, J. J. Baker, M. H. Blackwell, M. G. Brown,
C. B. Colonero, C. L. David, B. M. Tyrrell, and J. R. Wey, “Digitalpixel focal plane array technology,” Lincoln Lab J, vol. 20, no. 2, pp.
36–51, 2014.
[51] “HDR-eye dataset,” http://mmspg.epfl.ch/hdr-eye, Accessed: 2016-1130.
[52] H. Nemoto, P. Korshunov, P. Hanhart, and T. Ebrahimi, “Visual attention
in LDR and HDR images,” in 9th Int. Workshop on Video Process. and
Quality Metrics for Consumer Electronics (VPQM), 2015.
[53] T. Mertens, J. Kautz, and F. V. Reeth, “Exposure fusion,” in Proc. 15th
Pacific Conf. Comp. Graphics and App., 2007. PG ’07., Oct 2007, pp.
382–390.
[54] S. K. Nayar and T. Mitsunaga, “High dynamic range imaging: spatially
varying pixel exposures,” in Proc. IEEE Conf Computer Vis Pattern
Recog (CVPR), 2000, vol. 1, pp. 472–479 vol.1.
[55] C. Aguerrebere, A. Almansa, J. Delon, Y. Gousseau, and P. Muse, “A
Bayesian hyperprior approach for joint image denoising and interpolation, with an application to HDR imaging,” IEEE Trans Comput Imag,
vol. PP, no. 99, pp. 1–1, 2017.

1

Optimal Threshold Design for Quanta Image Sensor
(Supplementary Material)
Omar A. Elgendy, Student Member, IEEE and Stanley H. Chan, Member, IEEE

Abstract

arXiv:1704.03886v2 [cs.CV] 30 Oct 2017

•

•
•
•
•

•

This supplementary report provides the following additional information of the main article.
Derivation of SNRq (c) from exposure-referred SNR,
Properties of the incomplete Gamma function,
Comparison with the threshold design scheme by Yang [1],
Phase transition under different configurations,
Influence of Non-Boxcar Kernel G, and
Additional results for HDR image reconstruction.

I. D ERIVATION OF SNRq (c) FROM EXPOSURE - REFERRED SNR
In the literature of QIS devices, one metric to quantify the image quality is the exposure-referred signal-to-noise
[2]. In image processing, however, exposure-referred SNR is not commonly used. The goal of this section is to
show that the SNR we showed in the main article is equivalent to the exposure-referred SNR.

•
θ
Input

Fig. 1.

Poisson

•

•

BK−1,T −1
B0,1

Y

+

B0,0

S
Output
F (.)

Block diagram illustrating a QIS with input-output relation output = F (input)

To understand the exposure-referred SNR, we have to first understand two common ways of defining a signal
to noise ratio. Consider the truncated Poisson part of the QIS model shown in Figure 1. The input to this model
is the over-sampled measurement θ. The truncated Poisson process can be considered as a black box function F
which takes an input θ and generates an output S , defined as
S=

T
−1 K−1
X
X

Bk,t ,

(1)

t=0 k=0

where Bn = {Bk,t | k = 0, 1, . . . , K − 1, t = 0, 1, . . . , T − 1} is the spatial-temporal block containing all binary
bits corresponding to θ. As shown in the main article, the mean and variance of S are
E[S] = KT (1 − Ψq (θ)),

Var[S] = KT Ψq (θ)(1 − Ψq (θ)),

(2)

respectively.
The first notion of signal-to-noise, which is the one used in CCD and CMOS, is called the output-referred SNR.
SNROR is defined as the ratio between the output signal and the photon shot noise. Referring to Figure 1, this is
s
1 − Ψq (θ)
output signal
E[S]
SNROR =
=p
= KT
.
(3)
noise
Ψq (θ)
Var[S]

The authors are with the School of Electrical and Computer Engineering, Purdue University, West Lafayette, IN 47907, USA. Email: {
oelgendy, stanchan}@purdue.edu.

2
40

100

q increases

20

q increases

0

SNRER (c)

SNROR (c)

50

0

-20
-40
-60

-50

-80
-100
10−6

10−4

100

10−2

-100
10−6

10−4

(a) SNROR
Fig. 2.

10−2

100

c

c

(b) SNRER

Comparison of the SNRs for q ∈ {1, . . . , 16}. In this experiment, we fix α = 400, K = 4, and T = 30.

However, SNROR fails to work for QIS because the shot noise is arbitrarily small if all bits are 1 or 0. In [2], Fossum
called it squeezing of the noise. If we plot SNROR as a function of θ, then we observe that SNROR approaches to
infinity as θ grows.
The second notion of signal-to-noise, which is a modification of SNROR , is the exposure-referred SNR. SNRER
is the ratio between the exposure signal θ and the exposure-referred noise. This noise is defined as [2]:
dθ p
Var[S]
Exposure-referred noise =
dE[S]

The factor dEdθ[S] can be considered as the “inverse” transfer function from the output to the input.
determined by taking derivative of the expectation in (2) with respect to E[S]

dθ
dE[S]

can be

dKT (1 − Ψq (θ))
dE[S]
=
dE[S]
dE[S]

Using chain rule, we observe that
1 = −KT

Since

d
dθ Ψq (θ)

=

−e−θ θq−1
Γ(q) ,

dθ
d
Ψq (θ)
dθ
dE[S]

it holds that
1 = −KT

Hence,



−e−θ θq−1
Γ(q)



dθ
dE[S]

dθ
Γ(q)
=
dE[S]
KT e−θ θq−1

The exposure-referred SNR is defined as
exposure signal
exposure-referred noise
θ
=p
Var[S] dEdθ[S]
s
e−θ θq
KT
.
=
Γ(q) Ψq (θ) 1 − Ψq (θ)

SNRER =

Taking logarithm shows that SNRER is identical to the SNR derived from the Fisher Information shown in the
main article.

3
Ψq (θ)

1

1
θ ∈ Θq , q = 8
θ∈
/ Θq , q = 8
Ψq (θ), q = 1
Ψq (θ), q = 8
Ψq (θ), q = 16

0.9
0.8

0.9
0.8

0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0

2

4

6

8

10

12

14

16

18

20

Ψq (θ)

0

q ∈ Qθ , θ = 8
q∈
/ Qθ , θ = 8
Ψq (θ), θ = 1
Ψq (θ), θ = 8
Ψq (θ), θ = 16

2

4

6

8

Intensity θ

(a) Ψq (θ) vs. θ
Fig. 3.

10

12

14

16

18

20

Threshold q

(b) Ψq (θ) vs. q

Ψq (θ) as a function of θ and q. In defining, Qθ and Θq , we set  = 0.01.

II. P ROPERTIES OF THE INCOMPLETE G AMMA FUNCTION
In the main article, we used the incomplete Gamma function for QIS analysis. In this section, we provide more
details about the properties of the incomplete Gamma function.
First, we recall that the normalized upper incomplete Gamma function Ψq : R+ → [0, 1] is defined as
Z ∞
1
def
Ψq (θ) =
tq−1 e−t dt, for θ > 0, q ∈ N.
(4)
Γ(q) θ

where Γ(q) = (q − 1)! is the standard Gamma function.
In this equation, we note that Ψq (θ) depends on two variables: q and θ.
• As a function of θ . As we showed in the main article, Ψq (θ) is a monotonically decreasing function of θ
because the derivative is negative:
d
−θq−1 e−θ
Ψq (θ) =
< 0.
dθ
Γ(q)

However, Ψq (θ) is very close to 1 when θ is small, and is very close to 0 when θ is large. Therefore, there
exists a range of θ in which Ψq (θ) can attain a reasonably good inverse. We define this set as the θ-admissible
set
def
Θq = {θ | ε ≤ Ψq (θ) ≤ 1 − ε},
(5)
•

for any fixed q and a tolerance level ε. An illustration of Θq is shown in Figure 3.
As a function of q . The incomplete Gamma function Ψq (θ) can also be considered as a function of q . In this
case, Ψq (θ) is only defined for integer values of q . We illustrate the behavior of Ψq (θ) as a function of q in
Figure 3. The set of q in which Ψq (θ) is sufficiently away from 0 and 1 is defined as the q -admissible set.
def

Qθ = {q | ε ≤ Ψq (θ) ≤ 1 − ε}.

(6)

III. C OMPARISON WITH THE THRESHOLD DESIGN SCHEME BY YANG [1]
In this section, we compare our threshold scheme with the one in [1].
First, we recall that the optimality of our method is based on a lower-bound Lq (c) for the per-pixel SNR:
j αc k
q ∗ (c) = argmax SNRq (c) ≈ argmax Lq (c) =
K
q∈N
q∈N

(7)

Therefore, the optimal threshold is a function of c, which changes in space and in time.
In contrast, [1] uses a checkerboard pattern by alternating two thresholds (q1∗ , q2∗ ). These two thresholds are
obtained by maximizing the Cramér-Rao lower bound (CRLB) over a range of light intensity values [cmin , cmax ]:
Z cmax
∗ ∗
CRLB(q1 , q2 , c) dc.
(8)
(q1 , q2 ) = argmin
1≤q1 ,q2 ≤qmax

cmin

4

(a) q = 1, 15.97dB

(b) q = 5, 27.80dB

(c) q = 8, 25.84dB

(d) q = 10, 22.36dB

(e) q = 15, 15.35dB

(f) [1], 27.32dB

(g) Proposed, 28.15dB

(h) Ground Truth

Fig. 4. Spatial oversampling K = 4. Temporal oversampling T = 20. Quadratic B-spline kernel is used in synthesis and reconstruction
models. Gradient descent is used to obtain the ML estimate. For bisection threshold map, 8 frames are used for adapting the map, and 12
frames are used for reconstruction. For all other maps, the whole 20 frames are used for reconstruction.

As a result, the threshold is optimal in the average sense. To compare the two approaches, we followed the same
steps in [1] to obtain CRLB(q1 , q2 , c) for a checkerboard pattern in terms of Ψq (c) as follows
CRLB(q1 , q2 , c) =

2
X
e−2θ θ(2(qi −1))
α2
2K Γ(qi )2 Ψqi (θ) [1 − Ψqi (θ)]

(9)

i=1

where θ = αc/K . Using the parameters α = K(qmax −1 ), qmax = 16, K = 4, and using trapezoidal technique for
numerical integration over c, we obtain that q1∗ = 4 and q2∗ = 12. Figure 4 shows the reconstructed images using
uniform threshold maps with thresholds q ∈ {1, 5, 8, 10, 15}, the checkerboard threshold map in [1] with q1∗ = 4
and q2∗ = 12, and the oracle threshold map obtained by (7). In this experiment, our proposed method achieves 28.15
dB, which is 0.83 dB higher than the checkerboard pattern.
IV. P HASE TRANSITION UNDER DIFFERENT CONFIGURATIONS
In the main article, we showed the phase transition behavior of the ML estimate using K = 4, T = 50, and
δ = 2 × 10−4 . In this section, we study the effect of changing K , T , and δ on the phase transition region width.
As a function of T . Figure 5-Figure 6 illustrate the phase transition behavior when T = 10, 25, 50, and 100. As
T increases, the width of the green region increases. However, if we fix the range of the bit density 1 − E[γq (c)],
we observe that the SNR does not vary significantly even as T changes.
As a function of K . The spatial oversampling K affects both the threshold q ∗ (c) = bαc/Kc + 1 and the phase
transition width. Figure 7(a) illustrates the behavior of the threshold q ∗ as a function of K . As K increases, q ∗
decreases. However, the optimal q ∗ still stays within the set Qθ .
As a function of δ . The constant δ is used to define the set Qθ :
(
  1
  1 )
δ KT
δ KT
def
Qθ = q 1 −
≤ Ψq (θ) ≤
.
(10)
2
2
The constant δ is the tolerance level. When δ increases, the size of the set Qθ should also increase. This result is
shown in Figure 7(b).
Using the closed form expression of the average bit density 1 − Ψq (θ), we can calculate the average bit density
at the optimal threshold q ∗ = bθc + 1, which is shown in Figure 8. We notice that as long as θ ≥ 1, the average
bit density is between 0.264 and 0.630. Within this range, we observe from Figure 5-Figure 6 that the SNR does
not vary significantly if the estimated threshold is deviated from the optimal threshold. This observation relaxes
the requirement of the bisection method from obtaining the exact optimal threshold to obtaining a threshold that
make the bit density equal to 0.5. Since 0.5 ∈ [0.264, 0.630], we guarantee to achieve an SNR which is sufficiently
close to the optimal SNR.

5

Controlling θ ≥ 1 can be achieved by tuning the constant α. Tuning α can be hardware-implemented by increasing
the exposure period. Intuitively what θ ≥ 1 requires is that the average number of impinging photons per jot must
be at least one. If θ is less than one, then most bits will become zeros. Increasing exposure period (i.e., increasing
α) will ensure sufficient number of photons.
q ∈ Qθ
q∈
/ Qθ
Exact E[ b
c ]/c
Monte-Carlo E[ b
c ]/c
Exact 1 − E[γq (c)]
Monte-Carlo 1 − E[γq (c)]

1
0.8
0.6

0.8
0.6

0.4

0.4

0.2

0.2
10

20

30

40
q ∗ (c)

50

60

70

80

90

0

100

40

40

30

30

20
10
0

10

20

30

40

10

20

30

Threshold q

SNR in dB

SNR in dB

0

q ∈ Qθ
q∈
/ Qθ
Exact E[ b
c ]/c
Monte-Carlo E[ b
c ]/c
Exact 1 − E[γq (c)]
Monte-Carlo 1 − E[γq (c)]

1

50

60

70

80

90

50

60

70

80

90

100

70

80

90

100

Threshold q

20
10
0

100

40
q ∗ (c)

10

20

30

40

Threshold q

50

60

Threshold q

(a) T = 10, SNR ∈ [28.93, 29.57]

(b) T = 25, SNR ∈ [32.98, 33.74]

Fig. 5. Phase transition for T = 10 and T = 25. SNR range is shown for average bit density 1 − E[γq (c)] in the range [0.264, 0.630]. For
all cases, we set δ = 2 × 10−4 , and K = 4.

q ∈ Qθ
q∈
/ Qθ
Exact E[ b
c ]/c
Monte-Carlo E[ b
c ]/c
Exact 1 − E[γq (c)]
Monte-Carlo 1 − E[γq (c)]

1
0.8
0.6

0.8
0.6

0.4

0.4

0.2

0.2

0

10

20

30

40

50

60

70

80

90

q ∈ Qθ
q∈
/ Qθ
Exact E[ b
c ]/c
Monte-Carlo E[ b
c ]/c
Exact 1 − E[γq (c)]
Monte-Carlo 1 − E[γq (c)]

1

0

100

10

20

30

40

40

30

30

SNR in dB

SNR in dB

Threshold q

20
10
0

10

20

30

40

50

60

70

80

Threshold q

(a) T = 50, SNR ∈ [36.15, 36.80]

90

100

40
q ∗ (c)

50

60

70

80

90

100

70

80

90

100

Threshold q

20
10
0

10

20

30

40

50

60

Threshold q

(b) T = 100, SNR ∈ [39.35, 39.82]

Fig. 6. Phase transition for T = 50 and T = 100. SNR range is shown for average bit density 1 − E[γq (c)] in the range [0.264, 0.630].
For all cases, we set δ = 2 × 10−4 , and K = 4.

V. I NFLUENCE OF N ON -B OXCAR K ERNEL G
In this section, we discuss the boxcar kernel assumption in QIS model, i.e., G = K1 I N ×N ⊗ 1K×1 . We also
study the effect of assuming a general kernel G on our results.
On QIS, we typically assume that there are micro-lenses on top of each jot or a group of jots. These micro-lenses
ensure that the incident light converges onto the sensing site with no (or very minor) interference with adjacent
jots or groups. As a result, we can model the incoming light using the boxcar kernel. This assumption is perhaps
strong in some perspective, but it allows us to significantly simplify the theory and offer efficient implementations.
What if there is a mismatch between the physical model (e.g., using B-spline or Gaussian kernel G) and the
reconstruction (e.g., using boxcar)? To see the effect of this mismatch on the reconstruction quality, we conduct
two sets of experiments.

6

50

35

q ∗ (c) = ⌊αc/K⌋ + 1
q ∈ Qθ , T = 10
q ∈ Qθ , T = 100

40

30

30

Width of Qθ

Threshold q

25

20

20
15
10
δ
δ
δ
δ

10
5

0

4

9

16

25

36

49

64

81

100

0

0

Spatial Oversampling K

200

400

800

1000

Spatio-Temporal Oversampling KT

(a)

(b)

(a) The threshold q and Qθ as K increases. (b) The width of Qθ as KT and δ changes.
Average Bit Density 1 − E[γq (c)] at q ∗ (c)

Fig. 7.

600

= 10−3
= 10−4
= 10−5
= 10−6

0.7
0.630
0.6
0.5
0.4
0.3
0.264

0.2
0.1
0

0

10

20

30

40

50

Exposure θ = αc/K

Fig. 8.

Average bit density 1 − E[γq (c)] calculated at optimal threshold q ∗ = bθc + 1.

1D Signal: We consider a 1D signal with 10 coefficients. These 10 coefficients are modulated with boxcar
kernels and B-spline kernels to generate two sets of incident light. On the QIS simulator, we set the spatial
and temporal oversampling factors as K = 9 and T = 30, respectively. Then we use the oracle threshold map
for quantization. To reconstruct the images, we use boxcar kernel for both cases so that we have one matching
case and one mismatching case. Figure 9 shows the reconstructed signals. As expected, when the forward
model matches with the reconstruction model, the reconstructed image has the highest PSNR. However, the
gap between the cases are not significant.
• 2D Signal: Figure 10 shows a 2D example. Similar to the 1D case, boxcar kernel leads to the best reconstruction
but its gap with the other cases are not significant.
The reader might think why we do not use B-spline on the reconstruction so that it will match with the forward
model? In principle this is doable, but we need an iterative algorithm to compute the ML estimate such as gradient
descent as reported in [1]. In contrast, the boxcar assumption allows us to use a closed-form ML estimate, which
is practically much more affordable.
•

7

0.9

0.7

0.6
0.5
0.4
0.3

0.6
0.5
0.4
0.3

0.2

0.2

0.1

0.1

0

0

0.2

0.4

0.6

Ground Truth cn
Ground Truth θm
Reconstructed b
cn
Reconstruction θbm
Linear Kernels

0.8

Light Field Intensity

0.7

Light Field Intensity

0.9

Ground Truth cn
Ground Truth θm
Reconstructed b
cn
Reconstruction θbm
Box-Car Kernels

0.8

0.8

0

1

0

0.2

0.4

t

(a) Boxcar kernel, PSNR= 33.74 dB
0.9

0.7

0.6
0.5
0.4
0.3

0.6
0.5
0.4
0.3

0.2

0.2

0.1

0.1

0

0

0.2

0.4

1

0.6

Ground Truth cn
Ground Truth θm
Reconstructed b
cn
Reconstruction θbm
Cubic Kernels

0.8

Light Field Intensity

Light Field Intensity

0.7

0.8

(b) Linear B-spline, PSNR= 31.62 dB
0.9

Ground Truth cn
Ground Truth θm
Reconstructed b
cn
Reconstruction θbm
Quadratic Kernels

0.8

0.6

t

0.8

1

0

0

0.2

t

0.4

0.6

0.8

1

t

(c) Quadratic B-spline, PSNR= 31.66 dB

(d) Cubic B-spline, PSNR= 32.15 dB

Fig. 9. Spatial oversampling K = 9. Temporal oversampling T = 30. Oracle threshold map is used for quantization. Different kernels are
used in synthesis and boxcar kernel is used in reconstruction. ML closed-form is used for reconstruction

(a)

Boxcar
Ground Truth

(b)

Linear B-spline
Ground Truth

(c)

Quadratic B-spline
Ground Truth

(d)

Cubic B-spline
Ground Truth

Reconstruction
Reconstruction
Reconstruction
Reconstruction
(f)
(g)
(h)
PSNR= 33.33
PSNR= 30.67
PSNR= 32.54
PSNR= 33.74
Fig. 10. Ground truth and reconstructed images using simulated binary measurements synthesized by (a)(e) Boxcar, (b)(f) Linear B-spline,
(c)(g) Quadratic B-spline, and (d)(h) Cubic B-spline kernels. In this experiment, we spatially oversample each pixel by K = 4 × 4 binary
bits and we use T = 15 independent temporal measurements. We use 8 frames for learning the threshold map using bisection method, and
the remaining 7 frames are used for image reconstruction using the ML closed-form by the boxcar kernel assumption.
(e)

8

VI. S UPPLEMENTARY HDR

RESULTS

In this section, we show more results for HDR image reconstruction using our method compared to the fixed
threshold approach. Figure 11 show reconstructed HDR images using adapted Q-map by the bisection algorithm, and
fixed Q-maps with low threshold (q = 1) and high threshold (qmax = 16). The spatial and temporal oversampling
factors are K = 4, and T = 13, respectively. Sensor gain is α = K 2 /(qmax − 1).

Fig. 11.

q = 1, 15.94 dB

q = 16, 20.77 dB

Proposed, 29.97 dB

Ground Truth

q = 1, 17.94 dB

q = 16, 20.77 dB

Proposed, 31.46 dB

Ground Truth

q = 1, 15.74 dB

q = 16, 20.01 dB

Proposed, 31.65 dB

Ground Truth

Reconstructed HDR images using different threshold maps

R EFERENCES
[1] F. Yang, Bits from Photons: Oversampled Binary Image Acquisition, Ph.D. thesis, École Polytechnique Fédérale De Lausanne, 2012.
[2] E. R. Fossum, “Modeling the performance of single-bit and multi-bit quanta image sensors,” IEEE J. Electron Devices Soc., vol. 1, no.
9, pp. 166–174, Sep. 2013.

