IEEE TRANSACTIONS ON COMPUTATIONAL IMAGING, VOL. X, NO. X, AUGUST XXXX

1

Measurement-Adaptive Sparse Image Sampling and
Recovery

arXiv:1706.03129v2 [] 23 Nov 2017

Ali Taimori, Farokh Marvasti, Senior Member, IEEE

Abstract—This paper presents an adaptive and intelligent
sparse model for digital image sampling and recovery. In the proposed sampler, we adaptively determine the number of required
samples for retrieving image based on space-frequency-gradient
information content of image patches. By leveraging texture
in space, sparsity locations in DCT domain, and directional
decomposition of gradients, the sampler structure consists of
a combination of uniform, random, and nonuniform sampling
strategies. For reconstruction, we model the recovery problem as
a two-state cellular automaton to iteratively restore image with
scalable windows from generation to generation. We demonstrate
the recovery algorithm quickly converges after a few generations
for an image with arbitrary degree of texture. For a given number
of measurements, extensive experiments on standard image-sets,
infra-red, and mega-pixel range imaging devices show that the
proposed measurement matrix considerably increases the overall
recovery performance, or equivalently decreases the number
of sampled pixels for a specific recovery quality compared
to random sampling matrix and Gaussian linear combinations
employed by the state-of-the-art compressive sensing methods. In
practice, the proposed measurement-adaptive sampling/recovery
framework includes various applications from intelligent compressive imaging-based acquisition devices to computer vision
and graphics, and image processing technology. Simulation codes
are available online for reproduction purposes.
Index Terms—Cellular automaton, compressive sensing, directional gradients, measurement-adaptive sampling, sparse recovery, sparsity location, texture.

I. I NTRODUCTION
URRENT digital imaging devices at first acquire images
and then separately compress them leading to big datarelated problems. Contrarily, the aim of emerging Compressive
Sensing (CS) theory is to merge sampling and compressing
into one step by introducing the concept of “compressing
during sensing” for reducing data-rate, acquisition time, power
consumption, and device manufacturing cost with demanding
applications to next-generation Infra-Red (IR), remote sensing
and Magnetic Resonance Imaging (MRI) systems [1–3]. The
theory of CS states that signals with a sparse representation
in a specific domain can be recovered at a rate less than
the traditional Shannon-Nyquist rate theorem [4, 5]. DCT,
wavelet, and gradient spaces represent paradigms in which
natural images are sparse or at least compressible. The later
is not fully invertible.

C

A. Considered Scenario and Related Arts
In the literature, significant attempts have been done to
compressively sample data at a very low sampling rate [6–
10]. Conventional CS-based sampling/recovery strategies only
A. Taimori and Farokh Marvasti are with the Electrical Engineering
Department, Sharif University of Technology, Tehran 14588-89694, Iran (emails: alitaimori@yahoo.com; marvasti@sharif.edu).

exploit the sparsity-prior of signals in an appropriate domain
[11–13]. On the other hand, theoretical compressive sampling
analyses rely on using Gaussian or random sampling matrices
to measure signals [14]. However, to design sensing matrix,
recent scientific studies discover considering signal modelprior in sampling phase can significantly improve recovery
quality [7, 8, 10]. In this regard, the present paper concentrates
on finding a way to adaptively and sparsely sense image scene
in order to simply and efficiently retrieve the image.
Ji et al. in [6] introduced an adaptive Bayesian CS (BCS)
with abilities such as determining sufficient number of measurements and confidence of recovery. The researchers in [7]
proposed a nonuniform CS-based image sampling/recovery
method in which a Hidden Markov Tree (HMT) is utilized
to consider the correlation among sparse wavelet coefficients
in both sampling and reconstructing phases, so-called uHMTNCS. In [8], a variable density sampling strategy was designed
in the frequency domain which exploits priori statistical distributions of images in the wavelet space. Malloy and Nowak
in [9] suggested an adaptive compressed sensing algorithm
that in comparison to standard random nonadaptive design
matrices requires nearly the same number of measurements
but succeeds at lower SNR levels. Yang et al. in [10] designed
an efficient Mixed Adaptive-Random (MAR) sensing matrix
which employs image edge information with a sensor array
16 times less than the ultimate length of recovered image.
The Iterative Method with Adaptive Thresholding (IMAT) [13]
and its modified version equipped with interpolation, IMATI
[15], are extensions of the Iterative Hard Thresholding (IHT)
family for sparse recovery [12]. Instead of a fixed thresholding,
they adaptively threshold coefficients in a transform domain
such as Discrete Cosine Transform (DCT) during iterations.
To recover sparse signals from Gaussian measurements, the
authors in [16] presented a general iterative framework based
on proximal methods called Iterative Sparsification-Projection
(ISP). ISP family contains the well-known SL0 algorithm [17]
as a special case. In [18], a learning-based approach by jointly
optimizing the sensing and overcomplete dictionary matrices
was introduced which performs better than random matrices
or the matrices that are optimized without learning dictionary.
In [19], a variational BCS in complex wavelet domain was
suggested. This model, called TSCW-GDP-HMT, considers
sparsity and structural information of images.
B. Motivations and Contributions
Although aforementioned efforts result in performance improvement somewhat, it seems that the role of Artificial
Intelligence (AI) in the context of CS-based sampling/recovery

IEEE TRANSACTIONS ON COMPUTATIONAL IMAGING, VOL. X, NO. X, AUGUST XXXX

is yet negligible. In this paper, we take a step towards
considering AI and present an adaptive and intelligent model
for digital image sampling and recovery. For a given number of measurements, we show the proposed sensing matrix
considerably increases the overall recovery performance, or
equivalently decreases the number of sampling points for
a specific recovery quality compared to Gaussian, random,
and dynamic sampling matrices employed by the state-ofthe-art CS methods [6, 10, 11, 13, 16, 19]. The number
of required samples in the spatial domain to reconstruct
image is adaptively determined based on space-frequencygradient information content of image. To this intent, a mixture
of uniform, random, and nonuniform sampling strategies is
suggested [2]. The devised intelligent sampling mechanism
donates the advantages of simple, fast, and efficient recovery
compared to complicated methods [7, 16, 19]. For reconstruction, we model the recovery as a Cellular Automaton
(CA) machine to iteratively restore the image with scalable
windows from generation to generation. To the best of the
authors’ knowledge, this is the first work for modeling of
image reconstruction issue via a CA. The convergence of the
proposed CA-based recovery algorithm is guaranteed after a
few generations, thus making it suitable for reconstructing
the present mega-pixel range images, whereas sophisticated
approaches such as ISP [16], TSCW-GDP-HMT [19], and
BCS [6] fail to process and store such high-dimensional
signals with general-purpose computers.
The adaptivity in signal measurements helps to automatically select enough sampling points based on local information
content of the scene under view. For instance, lowpass images
are sampled with low rates, whereas cluttered scenes are
acquired with more samples for an accurate reconstruction.
Generally, available methods in the literature do not have
such a degree of flexibility and for each signal are manually
adjusted at a fixed rate. This undesirable phenomenon may
lead to data-redundancy in smooth images or informationloss in textured ones. In practice, the measurement-adaptive
image sampling/recovery scheme includes various applications from intelligent compressive imaging-based acquisition
devices [20] to computer vision and graphics, and image
processing technology, e.g. inpainting, remote sensing missing
recovery, denoising of impulsive noise, and mesh-based image
representation [21]. For reproduction purposes, simulation
codes are available online1 .
Briefly, the main contributions of the paper are
‚

‚

‚

proposing a new measurement-adaptive image sampling
mechanism consisting of three uniform, random, and
nonuniform samplers which, in a synergistically manner,
incorporates local space-frequency-gradient information
content of the image,
introducing a novel nonuniform sampler based on directional gradients,
suggesting a novel cellular automaton model for image
recovery which its convergence is guaranteed after a few
generations, and

1 http://ee.sharif.edu/∼imat/

2

Uniform Sampler
Texture
Measure
QTs
Image
Patch
(R)

Uniform
Pattern
Generator

Pu

Random
Pattern
Generator

Pr

Nonuniform
Pattern
Generator

Pn

Random Sampler
DCT

Quantizer

∑

Mask
(M)

Nonuniform Sampler
Directional
Gradients

Sparsity
Locations
Localizer

Fig. 1: The proposed measurement-adaptive sparse sampler.
exposing various examinations to demonstrate the efficiency of the proposed sampling/recovery approach under
different settings and application areas.
The remainder of this paper is organized as follows. Section II describes the proposed sampler. In Section III, we
explain the CA-based recoverer. Section IV evaluates the
effectiveness of our framework and compares it to modern
approaches. The paper is finally concluded in Section V.
‚

II. I NTELLIGENT S AMPLING
Based on the content of local image patches, the proposed
sampling method performs flexibly at various rates to capture
required informative samples for recovery step. Figure 1
depicts the block diagram of the suggested measurementadaptive sampler. As shown, the system gets a gray-level
image patch and gives the related binary mask. Due to psychovisual considerations, we assume that patches are square with
the length b “ 8. The final sampling mask is generated
from the union of three uniform, random, and nonuniform
patterns obtained from spatial, frequency, and gradient spaces,
respectively. In the further sub-sections, we discuss about each
sampling strategy in detail.
A. Texture Measure
To adaptively sense the information content of each local
image patch, it may be required to measure a representative
criterion of signal’s complexity [22]. Here, we utilized Shannon’s joint entropy determined from Haralick’s co-occurrence
matrix [23, 24] as a quantitative metric to measure image block
texture. We calculate the texture percentage η P r0, 100s as
ˆ
˙
B B
100 ÿ ÿ
1
η“
p̂i,j log2
,
(1)
Hmax i“1 j“1
p̂i,j
p “ rp̂i,j s
where the gray-level co-occurrence matrix P
BˆB
is the joint Probability Mass Function (PMF) which is obtained by the quantized“ or scaled
version of the input image
‰
s
patch, namely Gs “ gi,j
.
The
joint probabilities can
bˆb
be defined in the horizontal, vertical, main diagonal, and
off-diagonal directions. We used` the horizontal neighbor
as
˘
s
s
defined in MATLAB by p̂i,j “ P gr,c
“ i, gr,c`1
“ j , where
1 ď i, j ď B, 1 ď r ď b, and 1 ď c ď b ´ 1. The
variables r and c represent row and column, respectively.
We considered the number of distinct gray levels in Gs as

IEEE TRANSACTIONS ON COMPUTATIONAL IMAGING, VOL. X, NO. X, AUGUST XXXX

Fig. 2: From left to right, patterns of the uniform sampler
defined by matrices PUVL , PULT , PUBT , PUHT , and PUVH .
s
B “ 8, hence gi,j
P r1, Bs , 1 ď i, j ď b. The normalizer
p which
Hmax in (1) denotes the maximum entropy of P,
obtains from the uniform distribution, i.e. p̂i,j “ B12 , @
řB řB
1 ď i, j ď B. Thus, Hmax “ i“1 j“1 B12 log2 pB 2 q, which
B “ 8 ñ Hmax “ 6.

B. Uniform Sampler
In the proposed uniform sampling strategy, we punch local
patches at certain regular locations based on their local texture
content in space domain. In fact, the punch operation makes
appropriate regular patterns for a stable recovery. As shown
in Fig. 1, the uniform sampler at first gets the image patch
R “ rri,j sbˆb . Afterwards, the textureness of the signal is
estimated. The number of punched points is proportional to
the calculated texture percentage. By this intuition, whenever
the textureness is high, we decorate the image with denser
regular points and vice versa. Therefore, the binary mask of
uniform pattern is determined by the following fuzzy rule
$
0 ď η ă 10
’
’ PUVL ,
’
’
& PULT , 10 ď η ă 25
PUBT , 25 ď η ă 45 ,
(2)
Pu “
’
’
P
,
45
ď
η
ă
70
’
UHT
’
%
PUVH , 70 ď η ď 100
in which the matrices PUVL , PULT , PUBT , PUHT , and PUVH
are assigned for very-low, low, bandpass, high, and veryhigh textures, respectively. Figure 2 shows the regular lattices
employed in the proposed uniform sampler. The number of live
cells (punched points) in the matrices PUVL , PULT , PUBT ,
PUHT , and PUVH are 20 , 21 , 22 , 23 , and 24 , respectively.
This reveals the relation between the texture and the number
of live cells is nonlinear with an increasing exponential form.
The configuration of punched points affects the recovery performance. Sub-optimal lattices should have minimum number
of sample points and maximum number of live cells in a
predefined neighboring of missing samples. They should also
consider patch boundaries and local correlation. We tried to
experimentally design patterns based on such constraints.
C. Random Sampler
The intermediate layer of Fig. 1 depicts the proposed
random sampler. Here, we first convert the data matrix R
from the spatial domain to the frequency space, Rf P Rbˆb ,
by using 2-D DCT transform as Rf “ TRTT , where the
matrix T represents the transformation kernel of DCT. In
DCT domain, samples are decorrelated and sparse. After an
adaptive quantization of DCT coefficients, we estimate the
rate of random sampling. Accordingly, a uniformly distributed
random mask is finally generated. In the following sections,
we discuss about creating this mask.

3

1) Random Sampling Rate Measure: To measure sparsity, one way is to threshold insignificant coefficients in a
transform domain [13]. However, such thresholding may be
inaccurate due to varying information of the signal under
process/acquisition. To overcome this problem, we suggest
adaptive quantization of DCT coefficients inspired from the
following observation.
Observation 1 (Adaptive quantization): In JPEG compression standard where DCT transform is utilized, the results
in [22] reveal that the number of zero entries of quantized
DCT coefficients is inversely proportional to both compression
quality level, l, and the image block texture, η. Now, let the
r f “ rr̃f s P Zbˆb be the quantized DCT matrix
matrix R
i,j
r f k “ b2 ´ ZpR
r f q,
for a b ˆ b image patch, then k R
`0
where k ¨ k`p denotes the `p -norm and the operator Zp¨q
counts the number of zero entries of a matrix/vector. From
rf k 91.
r f k 9 1 and k R
the above results, we have both k R
`0 η
`0
l
By assuming their equivalence, we have η9l. This shows the
texture percentage defined in (1) can be a good criterion for
adaptive quantization of information in the frequency domain
to measure the necessary random sampling rate.
Details: Based on Observation 1, we calculate the quantized
DCT matrix as
r f “ tRf c Qη s,
R
(3)
where the symbols c and t¨s denote the entry-by-entry diviand the nearest integer, respectively. The matrix Qη “
“sion
η ‰
η
qi,j
, @qi,j
P N, also represents the quantization table
bˆb
related to the texture η. For determining Qη , we exploited the
Quantization Tables (QTs) introduced by Independent JPEG
Group (IJG) [25]. Based on the value of η, the corresponding
quantization matrix is chosen and applied (See Appendix.).
Therefore, we adaptively derive sparsity number from the
number of nonzero entries of quantized DCT coefficients by
r f k . Now, we suggest measuring the rate of random
k“kR
`0
sampling as the following formula, which is consistent with
the theory of CS, i.e. c ¨ k ¨ log10 p nk q,
"
tc ¨ k ¨ log10 pd ¨ nk qs,
k‰0
Rrs “
,
(4)
0,
otherwise
in which n “ b2 , and c and d are tunable factors. The
parameter c changes the number of undersampled points,
whereas the coefficient d prevents from decaying sampling rate
on well-textured patches and saturates it at a relatively fixed
rate to control unnecessary information. By setting c “ 1.3 and
d “ 2.8, the function in (4) is a non-decreasing curve vs k. It
is important to note that, in practical situations, the number of
non-zero quantized DCT coefficients, k, may be a large value
even near the length of signal, n. This phenomenon occurs
especially in high quality level settings or textured regions.
Generally, high performance compressive sensing in such a
dense scenario is problematic because compressive sensing
theory for real-world applications imposes the sampling rate
to be at most 10% [26]. We tried to consider this limitation
for calculating the rate of random sampling in (4) by setting
c “ 1.3 and d “ 2.8, which covers both sparse and dense
cases.

IEEE TRANSACTIONS ON COMPUTATIONAL IMAGING, VOL. X, NO. X, AUGUST XXXX

4

NW d
2

ri , j

N

d1

1

W

d8 NE

x

g

2

d3

d7
ri , j

E

y

Fig. 3: A partitioning of 2-D frequency domain defined in
[22, 27] into 6 regions.
d4
SW

d6
SE

d5
S

2) Graduate Randomization Procedure: Instead of blind
random sampling adopted by researchers of compressive sensing field, we introduce a Graduate Randomization Procedure
(GRP). In this procedure, by receiving both the feedback of
the sampling rate in frequency domain and the generated
uniform mask, we add a certain randomness to the uniform
pattern. In other words, from low- to high- textured patches,
GRP gradually adds random nonuniformity to the uniform
base lattice. It is important to note that random samples are
selected and added to locations other than uniform marks. In
order to generate the random pattern, the GRP algorithm is
implemented as follows
‚ Find the number of zero elements in the uniform pattern,
i.e. ZpPu q, and their locations by storing in the set L fi
(ZpPu q
ppi , q i q i“1 , where pi and q i represent the row and
column of the ith location, respectively.
‚ Shuffle the location of zero entries and substitute the set
L with the result.
‚ Truncate (4) by
"
ZpPu q, Rrs ą ZpPu q
Rrs “
.
(5)
Rrs ,
otherwise

Fig. 4: The gradient directions corresponding to cardinal and
ordinal directions, and the directional decomposition of the
gradient vector g into two d1 and d8 adjacent directions.

D. Nonuniform Sampler

single DC mode D, the five subsets of AC modes exist, namely
S1 to S5 corresponding to low, horizontal, vertical, diagonal,
and high frequencies information, respectively. If nonzero
quantized DCT coefficients are sparsely located in horizontal,
vertical, or diagonal regions, the corresponding edge points are
sampled in space domain. For instance, in a block of wellstructured vertical edges, the nonzero coefficients are only
appeared in the partitions S1 and S2 . To detect the edge
maps, we first determine a D-directional decomposition of the
image patch gradient [28, 29]. For this purpose, we utilized
Sobel operator to calculate the gradient of image block,
Br
Br
g fi r Bxi,j , Byi,j sT , @i, j, and considered D “ 8 with four
cardinal directions N, E, S, W in addition to four ordinal directions NE, SE, SW, NW to decompose the block. In order to
preserve the image structure in the boundary of blocks, we also
acquired border pixels for computing the gradient. Let Fdγ “
d
rfi,jγ sbˆb , @γ P r1, Ds ” tN, NW, W, SW, S, SE, E, NEu ”
td1 , d2 , d3 , d4 , d5 , d6 , d7 , d8 u, be D-directional gradient features. To project the gradient vector g in px-yq plane onto
the two adjacent directions
$
& |p| BrBxi,j |´| BrByi,j |q|
, k g k `2 ‰ 0 ,
kgk`
α1 “
(6)
2
%
0,
otherwise
$ ?
& 2 minp| BrBxi,j |,| BrByi,j |q
, k g k`2 ‰ 0 ,
kgk`
(7)
α2 “
2
%
0,
otherwise

It is very important to preserve the structure of image
like edges in recovery side. Generally, edge preservation is a
difficult task in interpolation even for sophisticated algorithms.
Due to fuzziness nature of image edges, an edge-prior sampler
can reduce uncertainty in reconstruction phase. Although, the
GRP carry this task somewhat, in well-structured horizontal, vertical, and diagonal edges, or periodic configurations
which their recovery is perceptually important from subjective
evaluation viewpoint, more edge samples are required to
reduce uncertainty and unambiguously recover the missing
scattered samples. Hence, we suggest a sparsity location-aware
algorithm to capture horizontal, vertical, and diagonal edges
in sampling phase.
To implement the above idea, we address sparsity locations
in different frequency regions. Figure 3 depicts a partitioning
of 2-D frequency domain defined in [22, 27]. Except for the

the following rules are utilized.
Bri,j
Bri,j
Bri,j
d
‚ If Bx ě 0 ^ | By | ď | Bx |, then assign fi,j1 Ð α1 ;
Bri,j
Bri,j
Bri,j
d5
else if Bx ă 0 ^ | By | ď | Bx |, then fi,j
Ð α1 ; else
Bri,j
Bri,j
Bri,j
d7
if Bx ě 0 ^ | By | ą | Bx |, then fi,j Ð α1 ; else if
Bri,j
Bri,j
Bri,j
d3
Bx ă 0 ^ | By | ą | Bx |, then fi,j Ð α1 .
Bri,j
Bri,j
d
‚ If Bx ě 0 ^ By ě 0, then assign fi,j8 Ð α2 ; else if
Bri,j
Bri,j
Bri,j
d2
Bx ě 0 ^ By ă 0, then fi,j Ð α2 ; else if Bx ă
Br
Br
Br
d4
0^ Byi,j ă 0, then fi,j
Ð α2 ; else if Bxi,j ă 0^ Byi,j ě 0,
d6
then fi,j
Ð α2 .
Afterwards, we apply a threshold to the magnitude of
two directions with stronger gradients to suppress redundant
uninformative samples. Figure 4 illustrates the directional
decomposition of the gradient vector g into two adjacent
d1
d2
directions. For this example, we have fi,j
“ α1 , fi,j
“ 0,

‚

‚

‚

Generate Rrs integer values randomly without replacement in the range r1, ZpPu qs, namely v “
rv1 , v2 , ¨ ¨ ¨ , vRrs sT .
Define the matrix of random pattern, Pr “ rpri,j sbˆb , and
initialize it with Pr “ O, where O represents a b ˆ b
zero matrix.
For n “ 1 to Rrs , at first, assign pi, jq Ð Ltvn u, and
then pri,j Ð 1.

IEEE TRANSACTIONS ON COMPUTATIONAL IMAGING, VOL. X, NO. X, AUGUST XXXX

d8
d7
d6
d5
d4
d3
“ α2
“ 0, and fi,j
“ 0, fi,j
“ 0, fi,j
“ 0, fi,j
“ 0, fi,j
fi,j
based on the specified rules. The steps of nonuniform sampling
algorithm are implemented as follows.
‚ Initialize D-directional gradient matrices, i.e. Fdγ “ O,
@γ P r1, Ds.
‚ Extract D-directional gradient features as explained
above for the vector g.
‚ Obtain the normalized versions of Fdγ , @γ P r1, Ds,
between 0 and 1.
‚ Define the matrix of nonuniform pattern, Pn
“
rpni,j sbˆb , and initialize it as Pn “ O.
‚ For determining the structure of edges, initialize a binary
string with the length of the number of AC frequency
regions, Nfr “ 5, such as a “ “a1 a2 a3 a4 a5 ” “ “00000”.
f
‚ For s “ 1 to Nfr , obtain the set Bs fi tr̃i,j |pi, jq P Ss u.
If k Bs k`0 ‰ 0, set as Ð “1”.
‚ If the string a “ “11000” _ “10100” _ “10010” _
“11010” _ “10110”, which respectively represent vertical, horizontal, diagonal, vertical-diagonal,
horizontalřb řor
d
b
diagonal edges, at first, sort the set t i“1 j“1 fi,jγ uD
γ“1 .
Afterwards, for the two directions having stronger gradients, calculate its gradient magnitude, and then renormalize it between 0 and 1. Finally, for i “ 1 to b and j “ 1
to b, if the resulting magnitude of gradient exceeds the
predefined threshold 0 ď τ ď 1, set pni,j Ð 1.

E. Measurement-Adaptive Sampling Algorithm
As shown in Fig. 1, the ultimate binary mask is generated
by the union of uniform, random, and nonuniform patterns as
ď
M“
Pi .
(8)

5

Algorithm 1 The proposed adaptive sampling algorithm
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:

Input: The image patch R.
Measure the texture η as (1).
Generate the uniform sampling mask Pu by (2).
Determine the quantization table Qη in terms of η.
Transform the patch R to 2-D DCT domain Rf .
r f “ tRf c Qη s.
Quantize DCT coefficients via R
rf k .
Estimate the sparsity number k “ k R
`0
Calculate the rate of random sampling as (4).
Perform the GRP routine to obtain the random mask Pr .
Run the nonuniform sampler algorithm
Ť to generate Pn .
Combine the created masks as M “ iPtu,r,nu Pi .
Obtain the scattered sample points by S “ M d R.
Outputs: The mask M and the subsampled image S.

Fig. 5: Applying Algorithm 2 on a 32 ˆ 32 DC image. From
left to right, sampling matrices encompassed by the greencolored frame illustrate initial state, and the evolution at 1st ,
2nd , and 3rd generations, respectively.

uniformity, and synchronicity [30]. They are composed of
discrete cells equipped with special rules to be able to solve
sophisticated problems. Various applications have been found
for CA in image processing tasks such as noise removal, edge
detection, and forgery localization [31].

iPtu,r,nu

The regularity, randomness, and structure of patterns in the
above combined form create naturally a sampling mask
with chaotic behavior. For selecting the threshold τ in the
nonuniform sampler, a trade-off between subjective quality
enhancement
Ť and increased redundancy exists. If τ Ñ 1, then
M Ñ Pu Pr . Therefore, we empirically set τ “ 0.9 to
maintain sampling near an optimal sub-Shannon-Nyquist rate.
The downsampled image patch is constructed from
S “ M d R,

(9)

in which the symbol d denotes Hadamard product. Algorithm 1 summarizes the proposed adaptive intelligent image
sampling strategy. For the previously captured image I of the
dimension h ˆ w, the proposed block-wise algorithm can be
repeated on all nonoverlapping image patches to get the overall
sampling mask Ms “ rmsi,j shˆw and the subsampled image
Is “ risi,j shˆw .
III. C ELLULAR AUTOMATON - BASED I MAGE R ECOVERY
As mentioned in the literature, different approaches exist
for image recovery [6–9, 11–13, 16–19]. Here, we tailor a
novel technique to recover scattered samples obtained from
the intelligent sampling stage. We model the reconstruction
mechanism by a dynamic CA. Cellular automata are simple
intelligent agents with fundamental characteristics of locality,

A. Cellular Automaton Modeling for Recovery
In modeling, we consider a 2-state CA machine, in which
the states of dead and live cells represent the zero and
one values in the sampling mask Ms , respectively. Contrary to conventional fixed-neighbor CA models, the proposed
CA-based recoverer performs as an iterative method which
applies variable-scale windows to the sampled image Is .
The size of square window, Ω, increases at each generation
of CA. In summary, we define the elements of model as
MpMs , Is , Ω, σ, ζ, pI, Mg q, in which σ represents the standard
deviation of a Gaussian kernel, ζ ě 1 is a coefficient
for increasing σ during generations, and Mg and pI denote
the matrix at next generation of Ms and the recovered
image, respectively. The rule behind the suggested recovery
algorithm is simple yet efficient as follows. By considering
the correlation among adjacent samples, missing pixels are
reconstructed via a Gaussian-weighted averaging. We find live
cells of Is around the central dead cell msi,j and extract their
corresponding weights in pΩ ´ 1q-Moore neighborhood [30],
namely the dynamic-length vectors x “ rx1 , ¨ ¨ ¨ , xk sT and
ω “ rω1 , ¨ ¨ ¨ , ωk sT , respectively. If at least one live cell exists,
then, the dead cell msi,j is replaced with the weighted mean
of live cells in the subsampled image and that dead cell will
revive at next generation. At each generation of CA-based recoverer, the model elements are updated. Algorithm 2 presents

IEEE TRANSACTIONS ON COMPUTATIONAL IMAGING, VOL. X, NO. X, AUGUST XXXX

the pseudo-code of CA-based image recovery algorithm, in
which the symbols t¨u and r¨s represent the floor and ceiling
functions, respectively. We utilized the replication rule to cope
with boundary conditions of border cells.
After the last generation, we utilize a post processing stage
to alleviate possible blockiness artifacts due to the patch-based
nature of recovery. This phenomenon may be observed in plain
regions. Based on the coding style of measurement-adaptive
sampler, the number of live cells within a window can be
an appropriate criterion to discriminate plain regions from
textured ones in the recovery phase. Therefore, we apply a
conditional smoothing filter only for flat regions. To this intent,
at first, we initialize the post-processed result pIp “ rpipi,j shˆw
as pIp “ pI. Afterwards, for all i P r1, hs and j P r1, ws, we scan
the image pI “ rpii,j shˆw with a square window of the length
Ωf “ 3. If the number of live cells in the initial sampling
p0q
mask Ms inside the window are less than the threshold
tρ ¨ Ω2f s, the corresponding pixel values of pI in the window
σ
are smoothed via the Gaussian kernel Wσf “ rwi,jf sΩf ˆΩf ,
for which 0 ď ρ ď 1. Then, the result is placed into pipi,j . In
experiments, we set the coefficient ρ “ 0.3 and the standard
deviation of smoother σf “ 1.5.

Algorithm 2 The proposed CA-based recovery algorithm

B. Convergence Analysis

23:
24:
25:
26:
27:
28:

As illustrated in Fig. 6, the CA-based recoverer first starts
with exposing high-textured dense missing samples at a fine
scale. Then, by updating model parameters, this process continues and finally terminates by retrieving lowpass sparse samples at a coarse scale. The advised variable-scale windowing
in comparison to fixed window speeds up the convergence rate
of Algorithm 2 without sacrificing recovery quality. Also, this
mechanism controls better error propagation during recursion.
In Lemma 1, we prove the algorithm quickly converges after a
few generations for an image with arbitrary degree of texture.
Lemma 1 (Convergence guarantee): Let patches of the
image I “ rii,j shˆw be square of the length b “ 8. If Ng
denotes the number of CA generations, then, Algorithm 2 for
recovering the image I guarantees to converge at most at three
iterations, i.e. Ng ď 3.
Proof: To prove, we consider the worst case, i.e. the
image under reconstruction is a fully DC image as the possible
sparsest signal to determine the upper bound. In such a case,
for all r hb s ˆ r wb s patches, we have η “ 0 which results
in Pu “ PUVL , and the random and nonuniform patterns
equal Pr “ Pn “ O. Hence, M “ Pu for all patches. Let
the live and dead cells be represented by white- and blackcolored cells, respectively. For visualization, Fig. 5 illustrates
generations of the sampling mask matrix Ms for a given flat
DC image of the dimension 32 ˆ 32. It is important to note
that, if h and w are not multiples of b, we can use zero
pκq
padding. If the matrix Ms denotes the sampling mask at
th
generation κ , then it is shown in Fig. 5 that the number of
p3q
dead cells at the 3rd generation is Nd “ ZpMs q “ 0. This
implies that N g “ 3 for a completely DC image because the
distance between two live cells is 7 pixels in both horizontal
and vertical directions. For a given image with a certain degree

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:

29:
30:
31:
32:
33:
34:
35:
36:
37:

6

Inputs: The mask Ms and the subsampled image Is .
Initialize Ω “ 3, σ “ 1, ζ “ 1.05, pI “ Is , Mg “ Ms ,
p0q
and Ms “ Ms .
Calculate the number of dead cells by Nd “ ZpMg q.
while Nd ‰ 0 do
ωh fi t Ω2 u
σ
sΩˆΩ .
Obtain the Gaussian kernel Wσ “ rwi,j
for r Ð 1, h do
for c Ð 1, w do
if msr,c “ 0 then
ϑflag fi 0
t fi 1
for p Ð ´ωh , ωh do
for q Ð ´ωh , ωh do
if msr`p,c`q “ 1 then
xt fi isr`p,c`q
σ
ωt fi wp`r
Ω
Ω
2 s,q`r 2 s
ϑflag “ 1
tÐt`1
end if
end for
end for
if ϑflag “ 1 then
T
îr,c “ Σxk ωωk
mgr,c “ 1
end if
Release the vectors x and ω.
end if
end for
end for
ΩÐΩ`2
σ Ðζ ¨σ
Is “ pI
Ms “ Mg
Nd “ ZpMg q
end while
p0q
Apply the post processing on pI via Ms to obtain pIp .
Output: The recovered post-processed image pIp .

of texture, the density of live cells in Ms is always more than
or equal to the worst case, thus demonstrating Ng ď 3.
Example 1 (Sampling and recovery): Figure 6 shows the
original images of Baboon, Boat, and Lena for which the
sampling masks, subsampled images, the intermediate, and
final results of recovery step are visualized. These images
have different texture averages sorted in descending order, for
which the sampler automatically extracted 35.23, 25.16 and
18.12 percent pixels, and the recoverer estimated the images
with PSNRs of 29.38, 29.43 and 32.19 decibels, respectively.
C. Practical Considerations
In order to implement the proposed sampling and recovery
algorithms in practice, we can consider a mechanism similar
to the structure of image acquisition using Rice single-pixel
camera [20]. However, one of the main differences is to

IEEE TRANSACTIONS ON COMPUTATIONAL IMAGING, VOL. X, NO. X, AUGUST XXXX

7

CA-based recovery approach. The solution is to focus the
physical scene on a Digital Micro-mirror Device (DMD) via
a primary lens. The DMD itself is partitioned into nonoverlapping b ˆ b patches such as 8 ˆ 8 blocks. In a sequential
line-scan manner and at each time, only one of the analog
image patches is reflected to another lens by configuring mirror
positions. Other mirrors reflect the scene on a black absorbent
surface. This secondary lens focuses the light of that single
patch on an 8 ˆ 8 CCD or CMOS sensor. After amplifying
the signal and passing through an Analog to Digital Converter
(ADC), the digitized sub-image, R, is obtained. Afterwards,
the texture of the patch is measured. Then, based on the
proposed sampler, the sampling mask, M, and the sub-sampled
image, S, are determined. This process is repeated until the
generation of the binary mask as well as the sub-sampled
image pertaining to the last patch. The overall sampling mask,
Ms , and corresponding sampled image, Is , can be stored in a
memory card. In the decoding phase, the proposed recoverer
gets Ms and Is and finally restores the image scene, pIp .
IV. E XPERIMENTS
This section provides extensive experiments to support our
findings. Algorithms were implemented in MATLAB and run
on an Intel Core i7 2.2GHz laptop with 8GB RAM.

Fig. 6: From top to down, different images respectively show
the original images, sampling masks, subsampled images,
recovered results after 1st , 2nd and 3rd generations, and the
results of applying conditional smoother.

utilize a small array of sensors but not a single photo-diode.
This allows us to measure adaptively the local image content,
sample in 2-D space, and recover the scene by the suggested

A. Experimental Settings
In order to evaluate the performance of our sampling/recovery framework and comparing it with other related
approaches, we utilized various databases having different
statistical characteristics. Image sets include the standard
databases of NCID [32], CCID [33], UCID [34], Microsoft
Research Cambridge Object Recognition Image Database2 ,
and a collection of well-known test images such as Baboon,
Cameraman, Lena, etc. Beside, we gathered a set of IR test
images as well as a mega-pixel range images database called
RCID [35], to more precisely investigate the efficiency of our
algorithms. Some of these images are shown in Fig. 12. In
order to evaluate recovery performance, we employed standard objective criteria of PSNR, SSIM [36], and Normalized
Recovery Error (NRE) defined as NRE fi k pIp ´ I k`2 {k I k`2
as well as subjective evaluation. We compared our sampling
and recovery algorithms to related approaches including the
intelligent MAR sampler suggested in [10], the scattered data
recovery technique of spline interpolation and the state-of-theart methods of BCS [6], uHMT-NCS AMP [7], IMAT [13],
IMATI [15], ISP [16], and TSCW-GDP-HMT [19].
We adjusted the set of parameters so that, in the suggested
sampler, the PSNR is maximum for a given average rate.
For tuning parameters, we used images of UCID database as
validation set. To generate uniform matrices in the uniform
sampler, different regular and periodic lattices with various
punched points may be used. We examined a set of such
configurations and finally selected the patterns introduced in
Fig. 2, which optimize the aforementioned objective criterion.
We also set parameters of CA-based recoverer and its postprocessing stage to obtain a sub-optimal objective PSNR
2 https://www.microsoft.com/en-us/download/details.aspx?id=52644

IEEE TRANSACTIONS ON COMPUTATIONAL IMAGING, VOL. X, NO. X, AUGUST XXXX

8

31

40
The coefficient c variations
The threshold variations

35

30.8
30.6

Average PSNR (dB)

Average PSNR (dB)

30
25
20
15

30.4
30.2
30
29.8
29.6

10
29.4

5

29.2
29

0
0

5

10

15

20

25

30

35

40

0

45

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Average sampling rate (%)

Fig. 7: Average PSNR (dB) vs average sampling rate (%) on
images listed in Table I.

performance criterion in dB. The influence of post-processing
stage with appropriate parameters of σf and ρ in the CAbased recovery algorithm is to improve PSNR up to about
0.5dB. The standard deviation of smoother σf , meanwhile,
has negligible impact on performance. The number of punched
points in the uniform sampler, the coefficients c and d in
the random sampler, and the threshold τ in the nonuniform
sampler control the sampling rate in the suggested algorithm.
On one hand, for natural image patches, we have often η ď 45,
thus considering short intervals for such textures. This also
means the number of punched points is often less than or equal
to 4. On the other hand, the main usage of the coefficient d
is to make the formula acceptable for dense images, i.e. high
values of k. Therefore, the influential factors for increasing or
decreasing the adaptive sampling rate are the coefficient c and
the threshold τ .

Fig. 8: The influence of the coefficient ρ on recovery performance for images listed in Table I.

Fig. 9: PSNR for different methods in dB. From top to down
and left to right, the original image, recovered images obtained
by uHMT-NCS AMP (24.15), ISP (19.46), BCS (18.58), IMAT
(24.25), IMATI (24.95), TSCW-GDP-HMT (26.72), and the
proposed algorithm (27.33), respectively.
BCS
IMAT
IMATI
TSCW-GDP-HMT
Proposed method

106

B. Performance on Standard Databases
104

Run-time (s)

Here, we evaluate the performance of our sampling/recovery
method on 21 well-known test images of the size 512ˆ512 and
four standard databases. The adaptive sampling rate measured
by sampler and recovery quality on image processing test
images are tabulated in Table I. For an average 23.75%
dynamic sampling rate, the mean of PSNR was 30.48dB. For
variations 0 ď c ď 3 and 0 ď τ ď 1, the curves of the average
PSNR (dB) vs the average sampling rate (%) on 21 wellknown images are depicted in Fig. 7. The curves demonstrate
that the sensitivity of our algorithm is low for variations c and
τ . Figure 8 also shows average PSNR (dB) vs ρ. In this figure,
the coefficient ρ “ 0 means the recovery without applying
post-processing step and the optimal value is ρ˚ “ 0.3, thus
demonstrating about 0.5dB gain.
As an example, Fig. 9 illustrates recovery of Baboon image for compared state-of-the-art approaches. To be able to
compare different algorithms of various complexities on the
laptop, we reduced the image to the size 128 ˆ 128. For all
methods, the sampling rate was 36.62%, i.e. 6000 measurements. Based on PSNR, the proposed approach outperforms
other methods. For BCS, IMAT, IMATI, TSCW-GDP-HMT,
and the proposed algorithm, the average PSNR (dB) on all
images with the size of 128 ˆ 128 was 18.52, 23.34, 24.26,
28.54, and 27.99, respectively. Other methods such as ISP
failed to compute the image restoration process. The average

102

100

20

25

30

35

40

45

50

55

60

65

Sampling rate (%)

Fig. 10: The run-time (s) of different methods in terms of
sampling rate (%) on the image of Fig. 9.

sampling rate determined by our algorithm was 34.3%. For a
fair comparison, we set this rate for competing algorithms, too.
Although the average PSNR for TSCW-GDP-HMT algorithm
is 0.55dB more than the suggested method, its computational
complexity is very high. To clarify this issue, Fig. 10 compares
the run-time (s) of different methods in terms of sampling rate
(%) on the 128 ˆ 128 Baboon image. As seen, the proposed
method is slightly better than IMATI.
Table II reports the performance of our approach on NCID,
CCID, UCID, and RCID databases, which are sorted in
descending order of textural information; i.e. RCID and NCID
databases have minimum and maximum texture averages, respectively. The average sampling rate acquired by the proposed
dynamic sampler confirms this subject. Images of each image-

IEEE TRANSACTIONS ON COMPUTATIONAL IMAGING, VOL. X, NO. X, AUGUST XXXX

TABLE I: Results of the proposed sampling/recovery method
Performance metrics

Image number

Image name

Sampling rate (%)

PSNRpdBq

SSIM

NRE

1

Baboon

35.23

29.38

0.8814

0.0085

2

Barbara

28.86

26.8

0.8431

0.0143

3

Boat

25.16

29.43

0.8227

0.0093

4

Butterfly

24.01

30.27

0.8356

0.0082

5

Cameraman

15.39

32.77

0.9295

0.0074

6

Einstein

18.09

31.26

0.7867

0.0082

7

Elaine

19.56

31.26

0.7613

0.0106

8

Fruits

16.6

31.13

0.8605

0.0056

9

Goldhill

25.51

30.4

0.8199

0.0107

10

Jetplane

18.59

31.36

0.91

0.0078

11

Lake

25.67

29.21

0.8324

0.0155

12

Lena

18.12

32.19

0.8779

0.0068

13

Livingroom

29.04

29.27

0.8288

0.0116

14

Owl

34.48

28.53

0.819

0.0115

15

Peppers

18.97

31.19

0.8464

0.0129

16

Pirate

24.62

29.99

0.8488

0.0095

17

Shack

34.1

27.09

0.8045

0.0234

18

Walkbridge

38

27.38

0.8279

0.0119

19

Woman-blonde

21.73

29.73

0.8229

0.0115

20

Woman-darkhair

11.55

36.69

0.9265

0.0053

21

Zelda

15.41

34.72

0.8915

0.0076

23.75

30.48

0.8465

0.0104

Average

set have also intra-database variability of textural information.
However, we tuned parameters of our sampler so that dynamic
sampling rate for each image doesn’t approximately exceed
50%, based on the limitation imposed in compressive sensing
theory for real-world applications [26]. This leads to an
acceptable average PSNR with a high variance, so that the
range of PSNR (dB), [min(PSNR), max(PSNR)], is [17.39,
48.48], [21.91, 41.25], [18.41, 40.03], and [26.29, 44.35] for
NCID, CCID, UCID, and RCID databases, respectively. From
compression viewpoint, recovery performance in comparison
to sensed samples is promising for all databases.

9

this rate for competing algorithms. The performance of the
proposed method is the second best.
D. Influence of Uniform and Nonuniform Samples
In this experiment, we considered two scenarios to evaluate
the role of the proposed uniform and nonuniform samplers. In
the first scenario, we omitted our proposed uniform sampler
from the scheme of Fig. 1 and compared the resulting random/nonuniform sampler with the whole model. This can be
done by setting Mu “ O. In this case, for a fair comparison
by preserving samples at an equal rate, we can manually
increase the coefficient c in the random sampler or decrease τ
in the nonuniform sampler. Here, we only increased c from the
resulting random/nonuniform sampler. In this case, the average
PSNR on 21 well-known images listed in Table I was 29.78dB.
In comparison to the result reported in Table I, a decrease
of nearly 0.7dB in PSNR is seen which demonstrates the
importance of existing uniform sampler. It is noticeable that
performance reduction in low-pass and periodic images like
Cameraman and Barbara is much more than cluttered ones.
In the second scenario, the nonuniform sampler was neglected by setting τ “ 1. Similarly, we increased c to
equalize the rate of resulting uniform/random sampler with
the complete scheme. The average PSNR was 30.52dB, that
is 0.04dB more than the average of Table I. Although such an
objective evaluation shows a little performance improvement
in the absence of nonuniform sampler, different subjective
evaluations narrate a better reconstruction of image structure
such as edges in the presence of nonuniform sampler. For
instance, Fig. 11 visualizes this phenomenon for four images
in which reconstructed edges are sharper in the entire solution than the second scenario. Therefore, for approximately
equal sampling rate and objective PSNR metric, the combined
strategy has subjective superiority in performance.

C. Test on Infra-Red Imaging Systems
Infra-red imaging has a wide range of applications from
surveillance and Intelligent Transportation Systems (ITSs) to
medical imaging. These images are naturally sparse and of
interest in the compressive sensing community to fabricate
their low-cost CS-based imagers [37]. In this experiment, we
evaluate the performance of adaptive sampling/recovery framework on a test-set including 20 images grabbed by various
near IR cameras for ITSs applications such as License Plate
Recognition (LPR). Four representative IR images are shown
in the first row of Fig. 12. Specifically, for recovering images
1 to 4 in this figure, PSNRs(sampling rates) in terms of dB(%)
were 46.68(2.08), 35.89(11.58), 42.04(4.41) and 38.86(6.57),
respectively. Generally, all IR images were recovered by an
average PSNR of 39.28dB for the average adaptive sampling
rate 7.16%, thus demonstrating an excellent IR recovery
quality of our scheme for the high sample reduction.
For BCS, IMAT, IMATI, TSCW-GDP-HMT, and the proposed algorithm, the average PSNR (dB) on all images with
the size of 128 ˆ 128 was 28.44, 22.11, 26.34, 35.23, and
31.88, respectively. The average sampling rate determined by
our algorithm was 15.37%. For a fair comparison, we set

E. Recovery Performance under Fixed Sampling Structures
This section compares the performance of different recovery
approaches under fixed sampling structures. To this intent,
we considered two sampling scenarios including conventional
pure random and our measurement-adaptive samplers. For
recovering scattered samples obtained from the proposed intelligent sampling stage as well as random sampling structure,
various approaches such as scattered data interpolators can be
employed. Therefore, in addition to modern sparse recovery
techniques [6, 13, 15, 16, 19], we applied the prominent
scattered spline interpolator of spline to separately investigate
the efficiency of our CA-based recoverer under the considered
scenarios.
In this experiment, we address the problem of image reconstruction of different scales. To do this, we selected 4
mega-pixel images from RCID database with the original size
minph, wq ˆ maxph, wq “ 3456 ˆ 5184 or vice versa. These
images are shown in the second row of Fig 12. To evaluate
the efficiency of the mentioned approaches, by using the
bicubic interpolation, we generated 6 down-sampled versions
for each of which by a factor about Ó 2 to finally construct

IEEE TRANSACTIONS ON COMPUTATIONAL IMAGING, VOL. X, NO. X, AUGUST XXXX

10

TABLE II: Performance of the proposed image sampling/recovery framework on four standard databases
Database

Total images

Image dimensions

pµ ˘ σ 2 q of sampling rate (%)

pµ ˘ σ 2 q of performance metrics
PSNRpdBq

SSIM

NRE

0.8355 ˘ 0.0049

0.02 ˘ 9.56e-5

NCID [32]

5150

256 ˆ 256

37.29 ˘ 3.4

26.57 ˘ 26.04

CCID [33]

1096

576 ˆ 768

26.68 ˘ 1.35

29.95 ˘ 11.22

0.849 ˘ 0.0025

0.0091 ˘ 2.2e-5

UCID [34]

1338

384 ˆ 512

26.56 ˘ 0.91

28.58 ˘ 9.1

0.8654 ˘ 0.0024

0.0154 ˘ 0.63e-5

RCID [35]

208

3456 ˆ 5184

12.24 ˘ 0.54

36.49 ˘ 10.06

0.9282 ˘ 9.68e-4

0.0021 ˘ 1.3e-6

Fig. 12: From left to right, images in the first and the second
rows indexed by numbers 1 to 4, show representative infra-red
and RCID mega-pixel images [35], respectively.

Fig. 11: From left to right and top to down, the reconstruction
of Walkbridge, Livingroom, Jetplane, and Lake images without
and with nonuniform sampler, respectively.

a 7-level pyramid. Table III compares the proposed CAbased recovery algorithm to modern reconstruction methods
under random and adaptive sampling scenarios for original
and scaled images. From lowest to highest resolution, i.e. in
the direction of increasing correlation, the average dynamic
sampling rates (%) were 42.95, 38.63, 34.4, 30.21, 26.69,
22.37, and 15.59, respectively. For a fair comparison, we
set the random sampling rate equal to that of the average
dynamic rate determined by measurement-adaptive sampler.
The bolded values show the best performance. Except for very

low resolution levels, ISP, TSCW-GDP-HMT, BCS, and spline
methods failed to yield viable results due to computational
complexity of matrix operations or huge memory requirement.
Because of heavy blockiness effects, BCS algorithm fails
to reconstruct images in both the random measurement and
the adaptive scheme especially in lowpass images. Generally,
state-of-the-art results were obtained for configurations of
measurement-adaptive sampler+spline for low-resolution and
measurement-adaptive sampler+CA for high-resolution levels.
This demonstrates the power of measurement-adaptive sampling scheme.
It is noticeable that IMAT and IMATI recovery methods
reconstruct the whole image rather than block-wise processing.
The sub-optimal parameters for adaptive thresholding technique of IMAT were α “ 0.5, β “ 50, and Niteration “ 30.
In ISP algorithm, DCT was utilized as sparsifying transform
and hard thresholding function. We also set its sub-optimal
parameters as γ “ 0.4, σnoise “ 0, τf “ 5e-6, c “ 0.9,
I “ 3, and maxiter “ 300. TSCW-GDP-HMT method uses
Generalized Double Pareto (GDP) distribution for modeling
signal sparsity. An inherent limitation of this algorithm is that
the input image should be square with a length of multiples
of 8. To tackle this problem and compare different approaches
as fair as possible, we first considered the nearest length more
than or equal to dimensions reported in Table III, and then
rescaled the recovered image to the real size. BCS algorithm
utilizes Daubechies1 2-D wavelet at four decomposition levels.
We also used MATLAB built-in function to implement spline
interpolation.
F. Evaluation under Fixed Recovery Algorithms
In this experiment, we evaluate the performance of three
sampling methods including the pure random, the state-of-theart dynamic MAR sampling approach [10], and the suggested
sampler under fixed recovery algorithms. We utilized 21
well-known test images of the original size 512 ˆ 512 for
performance evaluation. As known in Table III, only IMAT,

IEEE TRANSACTIONS ON COMPUTATIONAL IMAGING, VOL. X, NO. X, AUGUST XXXX

11

TABLE III: Analyzing the performance of different recovery techniques in the case of two distinct sampling scenarios on
images shown in the second row of Fig. 12. For each scenario, the values in bold type show the best performance.
PSNRpdBq
Image number

Size at pyramid level

Random sampling scenario
ISP
TSCW-GDP-HMT
BCS

Spline

CA

IMAT

Proposed measurement-adaptive sampling scenario
IMATI
ISP
TSCW-GDP-HMT
BCS
Spline

20.11
F
F
F
F
F
F

27.7
29.56
F
F
F
F
F

7.93
8.77
F
F
F
F
F

25.71
26.44
F
F
F
F
F

25.64
26.02
25.07
25.36
27.81
31.49
35.33

21.2
21.03
20.27
20.15
19.87
18.59
15.58

20.27
21.73
21
22.36
23.26
24.19
26.48

18.56
F
F
F
F
F
F

25.54
28.55
F
F
F
F
F

8.43
8.94
F
F
F
F
F

27.49
28.37
F
F
F
F
F

26.56
27.15
26.88
27.93
30.27
33.3
36.7

19.7
23.72
24.36
26.22
27.74
29.61
31.51

18.26
F
F
F
F
F
F

25.96
27.3
F
F
F
F
F

8.75
7.63
F
F
F
F
F

25.62
26.06
F
F
F
F
F

24.63
25.27
25.7
26.7
28.1
30.27
32.77

23.86
23.9
23.22
22.44
21.49
19.19
14.73

19.4
23.87
24.43
25.05
25.39
25.25
25.44

20.5
F
F
F
F
F
F

24.03
26.23
F
F
F
F
F

7.01
8.79
F
F
F
F
F

27.18
27.84
F
F
F
F
F

25.78
26.63
27.3
28.17
29.78
32.19
35.16

21.99
23.21
23.99
25.07
26.47
27.36
24.37

20.77
22.8
26.07
28.43
31.77
35.47
35.9

16.42
F
F
F
F
F
F

25.82
28.7
F
F
F
F
F

11.19
10.35
F
F
F
F
F

24.68
25.55
F
F
F
F
F

23.51
24.33
25.82
27.48
30.43
34.41
36.76

22.84
22.33
19.87
17.29
16.04
14.91
13.94

21.13
23.64
19.51
19.16
20
21.73
26.33

17.23
F
F
F
F
F
F

24.25
27.47
F
F
F
F
F

10.24
9.24
F
F
F
F
F

26.68
27.55
F
F
F
F
F

24.85
26.51
27.87
29.43
32.3
35.71
37.06

23.11
22.97
22.72
22.35
21.97
21.74
20.14

19.12
23.64
24.49
24.54
24.03
24.33
25.15

20.01
F
F
F
F
F
F

27.84
28.35
F
F
F
F
F

7.88
8.26
F
F
F
F
F

26.91
26.29
F
F
F
F
F

24.95
25.4
25.36
25.25
25.07
25.88
27.3

22.11
21.14
20.67
20.27
19.64
18.57
16.32

19.64
22.61
22.85
23.04
23.14
23.23
23.39

20.84
F
F
F
F
F
F

25.81
27.55
F
F
F
F
F

9.71
8.34
F
F
F
F
F

27.58
27.68
F
F
F
F
F

26.08
26.55
26.92
27.34
27.59
28.99
31.2

IMAT

IMATI

1

56 ˆ 80
104 ˆ 160
216 ˆ 320
432 ˆ 648
864 ˆ 1296
1728 ˆ 2592
3456 ˆ 5184

22.79
22.88
22.61
22.82
23.84
24.74
22.45

20.27
24.1
24.02
25.46
27.82
31.57
35.31

2

56 ˆ 80
104 ˆ 160
216 ˆ 320
432 ˆ 648
864 ˆ 1296
1728 ˆ 2592
3456 ˆ 5184

21.92
22.82
23.22
23.64
24.22
24.58
21.92

3

56 ˆ 80
104 ˆ 160
216 ˆ 320
432 ˆ 648
864 ˆ 1296
1728 ˆ 2592
3456 ˆ 5184

4

56 ˆ 80
104 ˆ 160
216 ˆ 320
432 ˆ 648
864 ˆ 1296
1728 ˆ 2592
3456 ˆ 5184

35

30

CA

32
IMAT
IMATI
CA

IMAT
IMATI
CA

30
28

Average PSNR (dB)

Average PSNR (dB)

25

20

15

10

26
24
22
20
18
16

5

14
12

0
Random

MAR

Our

Sampling method

Fig. 13: The comparison of recovery performance for different
sampling approaches.
IMATI and the proposed CA-based recoveres can restore the
images of such dimensions. Hence, these recovery algorithms
were considered to be able to recover subsampled images.
MAR sensing matrix is basically generated from the lowresolution versions of images with the size 128 ˆ 128. For
a fair comparison with this method, we created subsampled
images from Hadamard product of MAR masks and original
images in recovery side. Figure 13 depicts average PSNR
criterion in dB for different sampling/recovery configurations.
The sampling rate for all samplers is the same and equal to
23.75%. The first and the second ranks are belonging to the
proposed measurement-adaptive sampler+CA recoverer and
MAR sampler+CA recoverer, respectively. It is noticeable that
the performance of IMAT and IMATI recovery algorithms
under the pure random sampling strategy are better than other
state-of-the-art samplers due to their nature.

0

0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

0.09

0.1

Noise variance

2
,
Fig. 14: Average PSNR (dB) vs the noise variance, σAWGN
for the case of noisy measurements.

model the noisy subsampled image, Is , as
Is “ Is ` pMs d Nq,

(10)

where matrices Is , Ms , and N denote the noise-free subsampled image, the binary sampling mask, and Additive
White Gaussian Noise (AWGN) with the distribution of
2
N p0, σAWGN
q, respectively. We then applied the normalized
version of Is between 0 and 1 to CA recoverer. Figure 14
2
plots average PSNR in dB vs the noise variance, σAWGN
, for
21 well-known test images of the size 512ˆ512. In this figure,
2
σAWGN
“ 0 means the noise-free case. The suggested method
shows better performance in terms of the noise variance than
IMAT and IMATI approaches. All curves exponentially decay
2
by increasing σAWGN
. As mentioned before, other algorithms
fail to recover these images due to data dimensionality.
H. Performance of Joint Sampling/Recovery Frameworks

G. Recovery Robustness in the Presence of Noise
In order to evaluate the robustness of the proposed CAbased recovery algorithm in the presence of noise, we can

In practice, design of sampling structure and recovery stage
are closely related to each other. Here, we jointly compare
the performance of the proposed sampling/recovery with other

IEEE TRANSACTIONS ON COMPUTATIONAL IMAGING, VOL. X, NO. X, AUGUST XXXX

Fig. 15: From top to down and left to right, test samples of
Microsoft Database indexed by numbers 1 to 27.
0

10

12

comparison to sophisticated techniques of ISP [16], TSCWGDP-HMT [19], and BCS [6], especially in reconstruction of
mega-pixel range imaging such as remote sensing. In CAbased recovery algorithm, updating rule for reviving dead
cells is done based on a simple weighted averaging. As a
future work, more precise predictors can be utilized for this
purpose. Also, the suggested sampling/recovery pipeline can
be generalized to other sparsifying transform domains like
wavelets. Extensive tests on standard image data-sets, infrared, and mega-pixel imaging show the capabilities of proposed
technique for practical implementations of compressively sampled imaging systems.

−1

10

NRE

A PPENDIX
−2

10

uHMT−NCS AMP (0.0755)
ISP (0.1788)
BCS (0.1728)
IMAT (0.0496)
IMATI (0.0492)
Proposed adaptive sampler + CA−based recoverer (0.0291)

−3

10

5

10

15
Image Number

20

25

Fig. 16: Normalized Recovery Error (NRE) of different approaches for the images labeled in Fig. 15.
state-of-the-art frameworks [6, 7, 13, 15, 16]. To do this, we
collected 27 randomly selected samples from Microsoft Object
Recognition Database as shown in Fig. 15, which were already
employed in [7]. At first, the randomly chosen images were
cropped for highlighting salient objects in the scene, and then
their gray-scale versions were resized to 128ˆ128 dimensions
using the bicubic interpolation to be able to compare different
algorithms of various complexity. Figure 16 plots NRE vs
image number as indexed by numbers 1 to 27 in Fig. 15 for
competing methods. The average NREs are given in legend
parentheses. As shown, the average NRE of the proposed
scheme is lower than other methods. BCS failed to recover
the smooth images 4 and 18, whereas our scheme retrieved
them with minimum recovery error. For a fair comparison, we
set the sampling rate of competing methods the same as the
average rate 32.44% obtained adaptively by our sampler.
V. C ONCLUSION
This paper suggested a measurement-adaptive sparse image
sampling and recovery framework. The proposed sampler
adaptively measures the required sampling rate and accordingly determines sample positions by a combined uniform,
random, and nonuniform sampling mechanism. Our idea originated from the fact that natural images may have lowpass,
bandpass, highpass, or sparse components depending on the
scene under view. Therefore, a measurement-adaptive mechanism was advised to trap informative samples. Unlike Gaussian
measurement sensing scenario, the proposed sparse coding
style does not need all samples in advance.
In the recovery phase, we modeled the problem by a cellular automaton machine. The suggested algorithm converges
always at a few generations. Low computational burden and
memory usage are two main advantages of this algorithm in

Based on JPEG standard [25], the quantization matrix
corresponding to the texture η can be formulated as
Qη “ max pt
1ďi,jďb

1
sη
Qr ` Ju, Jq,
100
2

(11)

in which the matrix J denotes a b ˆ b all-one matrix. Also,
the scaling factor, sη , and the reference quantization table, Qr ,
are determined by
"
5000
0` ď η ă 50
η ,
sη “
,
(12)
2p100 ´ ηq, 50 ď η ď 100
and

¨

˚
˚
˚
˚
˚
Qr “ ˚
˚
˚
˚
˚
˝

16
12
14
14
18
24
49
72

11
12
13
17
22
35
64
92

10
14
16
22
37
55
78
95

16
19
24
29
56
64
87
98

24
26
40
51
68
81
103
112

40 51
58 60
57 69
87 80
109 103
104 113
121 120
100 103

61
55
56
62
77
92
101
99

˛
‹
‹
‹
‹
‹
‹.
‹
‹
‹
‹
‚

ACKNOWLEDGEMENTS
This Postdoc research was jointly sponsored by Iran National Science Foundation (INSF) and ACRI of Sharif University of Technology under agreement numbers 95/SAD/47585
and 7000/6642, respectively. The authors would like to thank
Dr Z. Sadeghigol who provided TSCW-GDP-HMT codes for
comparison. We also thank Prof A. Amini, Mr A. Esmaeili,
and other researchers in Signal Processing and Multimedia
Lab for their priceless comments.
R EFERENCES
[1] Y. C. Eldar and G. Kutyniok, Compressed sensing:
Theory and applications. Cambridge University Press,
2012.
[2] F. Marvasti, Nonuniform sampling: Theory and practice.
Springer Science & Business Media, 2012.
[3] J. Romberg, “Imaging via compressive sampling,” IEEE
Signal Process. Mag., vol. 25, no. 2, pp. 14–20, 2008.
[4] D. L. Donoho, “Compressed sensing,” IEEE Trans. Inf.
Theory, vol. 52, no. 4, pp. 1289–1306, 2006.
[5] E. J. Candès, J. Romberg, and T. Tao, “Robust uncertainty principles: Exact signal reconstruction from

IEEE TRANSACTIONS ON COMPUTATIONAL IMAGING, VOL. X, NO. X, AUGUST XXXX

[6]

[7]

[8]

[9]

[10]

[11]

[12]

[13]

[14]

[15]
[16]

[17]

[18]

[19]

[20]

highly incomplete frequency information,” IEEE Trans.
Inf. Theory, vol. 52, no. 2, pp. 489–509, 2006.
S. Ji, Y. Xue, and L. Carin, “Bayesian compressive
sensing,” IEEE Trans. Signal Process., vol. 56, no. 6,
pp. 2346–2356, 2008.
B. Shahrasbi and N. Rahnavard, “Model-based nonuniform compressive sampling and recovery of natural images utilizing a wavelet-domain universal hidden markov
model,” IEEE Trans. Signal Process., vol. 65, no. 1, pp.
95–104, 2017.
Z. Wang and G. R. Arce, “Variable density compressed
image sampling,” IEEE Trans. Image Process., vol. 19,
no. 1, pp. 264–270, 2010.
M. L. Malloy and R. D. Nowak, “Near-optimal adaptive
compressed sensing,” IEEE Trans. Inf. Theory, vol. 60,
no. 7, pp. 4001–4012, 2014.
J. Yang, W. E. I. Sha, H. Chao, and Z. Jin, “High-quality
image restoration from partial mixed adaptive-random
measurements,” Multimed. Tools Appl., vol. 75, no. 11,
pp. 6189–6205, 2016.
J. A. Tropp and A. C. Gilbert, “Signal recovery from
random measurements via orthogonal matching pursuit,”
IEEE Trans. Inf. Theory, vol. 53, no. 12, pp. 4655–4666,
2007.
T. Blumensath and M. E. Davies, “Normalized iterative hard thresholding: Guaranteed stability and performance,” IEEE J. Sel. Topics Signal Process., vol. 4, no. 2,
pp. 298–309, 2010.
F. Marvasti, M. Azghani, P. Imani, P. Pakrouh, S. J. Heydari, A. Golmohammadi, A. Kazerouni, and M. Khalili,
“Sparse signal processing using iterative method with
adaptive thresholding (IMAT),” in IEEE ICT, 2012, pp.
1–6.
E. J. Candès and T. Tao, “Near-optimal signal recovery from random projections: Universal encoding strategies?” IEEE Trans. Inf. Theory, vol. 52, no. 12, pp. 5406–
5425, 2006.
A. I. Zayed and G. Schmeisser, New perspectives on
approximation and sampling theory. Springer, 2014.
M. Sadeghi and M. Babaie-Zadeh, “Iterative
sparsification-projection: Fast and robust sparse signal
approximation,” IEEE Trans. Signal Process., vol. 64,
no. 21, pp. 5536–5548, 2016.
H. Mohimani, M. Babaie-Zadeh, and C. Jutten, “A fast
approach for overcomplete sparse decomposition based
on smoothed `0 norm,” IEEE Trans. Signal Process.,
vol. 57, no. 1, pp. 289–301, 2009.
J. M. Duarte-Carvajalino and G. Sapiro, “Learning to
sense sparse signals: Simultaneous sensing matrix and
sparsifying dictionary optimization,” IEEE Trans. Image
Process., vol. 18, no. 7, pp. 1395–1408, 2009.
Z. Sadeghigol, M. H. Kahaei, and F. Haddadi, “Model
based variational Bayesian compressive sensing using
heavy tailed sparse prior,” Signal Process.: Image,
vol. 41, pp. 158–167, 2016.
M. F. Duarte, M. A. Davenport, D. Takbar, J. N. Laska,
T. Sun, K. F. Kelly, and R. G. Baraniuk, “Single-pixel
imaging via compressive sampling,” IEEE Signal Pro-

13

cess. Mag., vol. 25, no. 2, pp. 83–91, 2008.
[21] X. Li, “Anisotropic mesh adaptation for image representation,” EURASIP J. IMAGE VIDE., vol. 2016, no. 1,
p. 26, 2016.
[22] A. Taimori, F. Razzazi, A. Behrad, A. Ahmadi, and
M. Babaie-Zadeh, “Quantization-unaware double JPEG
compression detection,” J. Math. Imaging Vis., vol. 54,
no. 3, pp. 269–286, 2016.
[23] R. M. Haralick, K. Shanmugam, and I. H. Dinstein,
“Textural features for image classification,” IEEE Trans.
Syst., Man, Cybern., vol. SMC-3, no. 6, pp. 610–621,
1973.
[24] R. Jain, R. Kasturi, and B. G. Schunck, Machine vision.
McGraw-Hill, Inc., 1995.
[25] W. Luo, J. Huang, and G. Qiu, “JPEG error analysis and
its applications to digital image forensics,” IEEE Trans.
Inf. Forensics Security, vol. 5, no. 3, pp. 480–491, 2010.
[26] M. Lustig, D. L. Donoho, J. M. Santos, and J. M. Pauly,
“Compressed sensing MRI,” IEEE Signal Process. Mag.,
vol. 25, no. 2, pp. 72–82, 2008.
[27] R. Hasimoto-Beltrán, S. Baqai, and A. Khokhar, “Transform domain inter-block interleaving schemes for robust
image and video transmission in ATM networks,” J. Vis.
Commun. Image R., vol. 15, no. 4, pp. 522–547, 2004.
[28] Z.-L. Bai and Q. Huo, “A study on the use of 8directional features for online handwritten chinese character recognition,” in IEEE ICDR, 2005, pp. 262–266.
[29] L.-L. Huang, A. Shimizu, Y. Hagihara, and H. Kobatake,
“Gradient feature extraction for classification-based face
detection,” Pattern Recogn., vol. 36, no. 11, pp. 2501–
2511, 2003.
[30] J. L. Schiff, Cellular automata: A discrete view of the
world. John Wiley & Sons, 2011, vol. 45.
[31] P. Rosin, A. Adamatzky, and X. Sun, Cellular automata
in image processing and geometry.
Springer, 2014,
vol. 10.
[32] Q. Liu, A. H. Sung, and M. Qiao, “Neighboring joint
density-based JPEG steganalysis,” ACM Trans. Intell.
Syst. Technol., vol. 2, no. 2, 2011, article 16.
[33] A. Olmos and F. A. A. Kingdom, “A biologically inspired
algorithm for the recovery of shading and reflectance
images,” Perception, vol. 33, no. 12, pp. 1463–1473,
2004.
[34] G. Schaefer and M. Stich, “UCID: An uncompressed
color image database,” in SPIE SRMAM, vol. 5307, 2004,
pp. 472–480.
[35] A. Taimori, F. Razzazi, A. Behrad, A. Ahmadi, and
M. Babaie-Zadeh, “A novel forensic image analysis tool
for discovering double JPEG compression clues,” Multimed. Tools Appl., vol. 76, no. 6, pp. 7749–7783, 2017.
[36] Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli,
“Image quality assessment: From error visibility to structural similarity,” IEEE Trans. Image Process., vol. 13,
no. 4, pp. 600–612, 2004.
[37] H. Chen, N. Xi, B. Song, L. Chen, J. Zhao, K. W. C.
Lai, and R. Yang, “Infrared camera using a single nanophotodetector,” IEEE Sensors J., vol. 13, no. 3, pp. 949–
958, 2013.

