This article has been accepted for publication in a future issue of this journal, but has not been

fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TII.2017.2724205, IEEE
Transactions on Industrial Informatics
IEEE TRANSACTIONS ON INDUSTRIAL INFORMATICS (ACCEPTED MANUSCRIPT: TII-17-0340) 1

Image Autoregressive Interpolation Model using

GPU-Parallel Optimization
Jiaji Wu, Member, IEEE, Long Deng, Student Member, IEEE, and Gwanggil Jeon, Member, IEEE

Abstract—With the growth in the consumer electronics indus- proposed. To preserve edge structures in interpolation, Li and
try, it is vital to develop an algorithm for ultra high definition Orchard [3] proposed to estimate the covariance of HR images
(HD) products that is more effective and has lower time com- from the covariance of LR images, and to then interpolate the
plexity. Image interpolation, which is based on an autoregressive
model, has achieved significant improvements compared with the missing pixels based on the estimated covariance. In recently
traditional algorithm with respect to image reconstruction, in- years, one of the most promising algorithms introduced was
cluding a better peak signal-to-noise ratio (PSNR) and improved autoregressive modeling image interpolation, as proposed by
subjective visual quality of the reconstructed image. However, the Zhang and Wu [6], in which a piecewise 2-D autoregressive
time-consuming computation involved has become a bottleneck in (PAR) model is used to recover the HR image block by block
those autoregressive algorithms. Because of the high time cost,
image autoregressive-based interpolation algorithms are rarely in a soft-decision way (SAI, soft-decision adaptive interpo-
used in industry for actual production. In this study, in order lation). This technique is also used in image compression
to meet the requirements of real-time reconstruction, we use and has yielded promising results [8]. Inspired by Zhang
diverse CUDA optimization strategies to make full use of the GPU and Wu, a number of optimization approaches have been
(NVIDIA Tesla K80), including a shared memory and register proposed, based on autoregressive modeling [9], [10]. In [9],
and multi-GPU optimization. To be more suitable for GPU-
parallel optimization, we modify the training window to obtain a the authors used the Gauss Seidel method to further optimize
more concise matrix operation. Experimental results show that, the SAI model. In [10], the SAI model was applied to arbitrary
while maintaining a high PSNR and subjective visual quality and factor image interpolation, and achieved promising results.
taking into account the I/O transfer time, our algorithm achieves In [11], the authors used an autoregressive model in sparse
a high speedup of 147.3 times for a Lena image and 174.8 times representation for image interpolation. But the core of all these
for a 720p video, compared to the original single-threaded C
CPU code with -O2 compiling optimization. autoregressive models was a time-cost approach, which makes
it difficult to meet the requirements of real-time reconstruction
Index Terms—Image interpolation, autoregressive model, par- in actual production. This is the main reason that GPU-parallel
allel optimization, CUDA, GPU
optimization is needed for those algorithms.
With the development of high-performance parallel com-
I. I NTRODUCTION puting technology, GPU-based acceleration has been used in
aspects of scientific computing [12-18] for their low cost and
Due to their high video quality, ultra high-definition (HD)
efficient power usage. Mussi and Daolio used the CUDA
products are becoming increasingly popular. However, ultra
computing platform to optimize particle swarms [12]. In [14],
HD resolution requires a more effective algorithm that has
the authors used a GPU-based model for the infrared atmo-
lower time complexity for real-time reconstruction. Interpola-
spheric sounding interferometer (IASI). In [17], CUDA was
tion is an important method in image processing and is a basic
used to optimize deep learning algorithms and also played
requirement for ultra HD products for using low resolution
an important role in computed tomography (CT) reconstruc-
(LR) images to generate high resolution (HR) images. This
tion [16]. With the increasing capacity of cores and memory
method is widely used in image magnification, deinterlacing,
in GPU architecture, it is now possible to solve more complex
and denoising. Conventional bilinear interpolation, cubic con-
computing problems. CUDA is a general purpose parallel
volution interpolation [1], and cubic spline interpolation [2]
computing platform and programming model that leverages
are widely used for real-time applications because of their
the parallel computing engine in NVIDIA GPUs to solve
relatively low complexity. However, these classical algorithms
many complex computational problems. As the pixels to be
can lead to interpolation artifacts such as blurring, ringing,
interpolated are independent from each other in diagonal or
jaggies, and zippering. To satisfy the demand to better dis-
vertical interpolations in autoregressive models, these models
play resolution, many excellent algorithms [3-7] have been
are suitable for optimization in parallel by CUDA.
Manuscript received July 30, 2016; revised February 23, 2017. This work is One of the key aspects in the autoregressive modeling
supported by National Natural Science Foundation of China (No. 61377011), image interpolation method is the use of a suitable training
Shaanxi Science and Technology Department (No. 2016KTCQ01-116), and window. To capture greater detail of the interpolated pixels
NRF Grant funded by the Korean Government (2015R1D1A1A01058171).
J. Wu and L. Deng are with the School of Electronic Engineering, Xidian while reducing the computation required, we can use the
University, Xi’an, China, 710071, China (email: wujj@[Link], training windows of Zhang and Wu [6], which are large
ldeng@[Link]). octagonal windows, with four pixels estimated simultaneously
G. Jeon is with the School of Electronic Engineering, Xidian University,
Xi’an, China, 710071, China and Department of Embedded Systems Engineer- in a larger slip window, and which is significant for serial
ing, Incheon National University, 22012, Korea (email: ggjeon@[Link]). computing. Considering the register, shared memory, and other

1551-3203 (c) 2017 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See [Link] for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TII.2017.2724205, IEEE
Transactions on Industrial Informatics
IEEE TRANSACTIONS ON INDUSTRIAL INFORMATICS (ACCEPTED MANUSCRIPT: TII-17-0340) 2

Fig. 1. Two interpolation steps based on the autoregressive model. The white
dots are the LR image pixels, and the dark and gray dots are the missing HR
pixels, interpolated in two steps.

Fig. 2. Spatial configuration of SAI. The white dots are existing LR pixels, the
resource limitations of GPUs, if the load of each thread is dark dots are missing HR pixels. − →
a is the 8-connected neighbors directional
too large, the concurrency will drop dramatically. To calculate −
→
coefficients of the pixels and b is the 4-connected neighbors directional
the interpolated pixels, a large number of matrix operations coefficients of the pixels. W is the training window.
must be performed, including multiple matrix transposition,
matrix multiplication, and matrix inversion. So, in this study,
we modified the training window to obtain a more concise autoregressive coefficient for the (m, n)th neighbor of pixel
matrix operation. Our experimental results show that while X(i, j). Let xi ∈ Il and yi ∈ Ih be the pixels of images
maintaining a high PSNR and subjective visual quality our Il and Ih , and yit (t = 1, 2, . . .) be the neighbors of pixel
algorithm achieves a high speedup in a GPU (NVIDIA Tesla location i in the HR image. The main idea of the image
K80). autoregressive interpolation model is that known pixels can
The remainder of this paper is organized as follows. In be described by unknown pixels, and unknown pixels can also
section II, we review the SAI model and describe the specific be described by known pixels according to geometric duality.
implementation of our algorithm, including some mathemat- The parameters α(m, n) remain constant or near constant in
ical optimizations for matrix operations. In section III, we a small locality. Structures such as edges and textures can
present GPU/CUDA parallel optimization, which is used in be learned by fitting samples of a local window. The image
many optimization strategies based on CUDA, and which autoregressive interpolation model has two parts: diagonal
we used in our experiment with the NVIDIA Tesla K80. In interpolation and vertical interpolation. Diagonal interpolation
section IV, we present our experimental results in terms of the is used to calculate the HR pixels whose neighboring pixels
PSNRs for the test images. In section V, we use our algorithm are known. In vertical interpolation, the pixels estimated in
for image interpolation with an arbitrary factor and for a 720p the diagonal interpolation are regarded as “known” in the HR
video. Finally, we draw our conclusions in section VI. image. First, we obtain the LR image from a HR image by
down-sampling. The white dots in Fig. 2 are the LR image
II. I MAGE AUTOREGRESSIVE I NTERPOLATION M ODEL pixels and the gray and dark dots are the missing HR pixels.
Then, the dark dots are interpolated by diagonal interpolation.
As the spatial details and edge sharpness of autoregressive Finally, the gray dots are interpolated by vertical interpolation.
interpolation model-based interpolated images are better than As vertical interpolation is performed in a similar way, we only
those produced by traditional methods [1], [2], [3], this model explain diagonal interpolation in the following.
produces results that are significantly better with respect to
Figure 2 shows the training window of our algorithm where
both PSNR and subjective visual quality. Fig. 1 shows the
the white dots are LR pixels denoted by x, and the dark dots,
interpolation model. Let Ih be the HR image to be estimated,
denoted by y, are the missing HR pixels to be interpolated by
Il be the observed LR image, which is a down-sampled version
vertical interpolation. The interpolation is formulated to follow
of the HR image by a factor of two. The arbitrary factor image
the square problem, where parameter λ is the Lagrangian
interpolation method is described in section V.
factor, as follows:
Similar to the method used in [6], in this study we propose
→
− (8)
X X
to model an image as follows: y = arg min{
c kyi − at xit k
X {y,a,b} i∈W t
X(i, j) = α(m, n)X(i + m, j + n) + υi,j , (1) X X (8)
X X (4)
(m,n)∈W + kxi − at yit k +λ kyi − bt yit k},
i∈W t i∈W t
where W is a local window shown in Fig. 2, containing 4 × 4 (2)
unknown pixels, υi,j is a random perturbation independent P P (4) P P (4)
of spatial location and the image signal. α(m, n) is the subject to i∈W kyi − t bt yit ≈ i∈W kxi − t bt xit k.

1551-3203 (c) 2017 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See [Link] for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TII.2017.2724205, IEEE
Transactions on Industrial Informatics
IEEE TRANSACTIONS ON INDUSTRIAL INFORMATICS (ACCEPTED MANUSCRIPT: TII-17-0340) 3

TABLE I
S AMPLE I MAGES SIZE ( FACTOR = 2)

Sample Image Lena Peppers Pentagon Airport

LR Size 256 × 256 256 × 256 512 × 512 512 × 512
HR Size 512 × 512 512 × 512 1024 × 1024 1024 × 1024

The matrix operation is very time-consuming, so we must

try to reduce the intermediate process. As matrix D is sparse,
we can rewrite D→ −
x explicitly, and define G = D→−x . We can
also decompose matrix C to an identity matrix I and a matrix
H, where I9 is a 9-by-9 identity matrix.
 
I
C= 9 (9)
H
Then, Eq. (8) is rearranged as in Eq. (10):
→
−
y = (I9 + H T H)−1 (G1 + H T G2 ),
c (10)
Fig. 3. Determination model of parameters −→a (8-connected neighbor) and
−
→ where
b (4-connected neighbor), estimated by known pixels.
 
a1 a2 0 a3 a4 0 0 0 0
→
− 0 a1 a2 0 a3 a4 0 0 0
In fact, these →
−
a and b parameters can also be estimated  
from the least squares method proposed by Li and Orchard [3], 0
H= 0 0 a1 a2 0 a3 a4 0 ,
0 0 0 0 a1 a2 0 a3 a4 
as shown in Fig. 3. This method is based on the geometric Hb
duality assumption that the image is piecewise stationary, so
a1 x1 + a2 x2 + a3 x5 + a4 x6
 
there is correspondence between the LR pixel covariance and
 a1 x2 + a2 x3 + a3 x6 + a4 x7 
the HR pixel covariance,
..
 
 
→
− .
a = arg min k→
−
x − A→
−  
c ak (3) 
..

→
−
a
 

 . 

 .. 
→
− →
− .
b = arg min k→
−
 
x −B bk (4)
c  
→
−
 .
..

b
   
G1  
G= = . , (11)
where →−x = [x1 , x2 , . . . , xn ]T is the data vector containing M 2 G2 

.
. 

pixels inside the local window. A is an M 2 × 4 matrix whose  .. 
(8)
 . 
ith row vector comprises the four 8-connected neighbors xit 
a x + a x + a x + a x 

2
of xi , and B is an M × 4 matrix whose i row vector th  1 11 2 12 3 15 4 16 
(4)
 x6 
comprises the four 4-connected neighbors xit of xi . We note  
 x7 
that t = {1, 2, 3, 4} and M = 4. Finally, we can obtain:  

 x10 

→
−
a = (AT A)−1 AT →
c −
x, (5)
 x 
11
0
→
−
b = (B T B)−1 B T →
−
x, (6)
c  
Hb = λ 0 −b1 0 −b3 1 −b4 0 −b2 0 . (12)
Then, Eq. (2) can be rewritten as Eq. (7) by substituting →
−
a The parameter λ in Eq. (12) is the Lagrangian factor in Eq. (2).
→
−
and b : The amount of computation required by Eq. (10) is less than
→
−
y = arg min kC →
c −
y − D→
−
x k, (7) that required by Eq. (8).
→
−
y The sample images in our experiments were Lena, Peppers,
Pentagon, and Airport, as shown in TableI. Table II compares
where → −
y = [y1 , y2 , . . . , y9 ] is the vector of nine unknown pix- the PSNR and CPU running times. We compared the running
els in the current window, as shown in Fig. 2. The pixel y5 is time of our algorithm with that of SAI. All simulations were
the estimated pixel in one block. Now, this has become a linear conducted on an Intel core (TM) i7-920, 2.66 GHz CPU.
least square optimization problem. The optimal estimation of The executable program of SAI was provided by Zhang and
→
−y is given by the following equation: Wu [6].
→
−
y = (C T C)−1 C T D→
c −
x, (8) Table II reveals that the CPU running times (single-thread)
of our method are slightly reduced compared to those of
The sizes of C and D are 14 × 9 and 14 × 14, respectively. SAI. This is because the matrixes become smaller, which

1551-3203 (c) 2017 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See [Link] for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TII.2017.2724205, IEEE
Transactions on Industrial Informatics
IEEE TRANSACTIONS ON INDUSTRIAL INFORMATICS (ACCEPTED MANUSCRIPT: TII-17-0340) 4

TABLE II
C OMPARISON OF PSNR AND RUNNING TIMES ON CPU

Sample PSNR (dB) Running Time

Images SAI Our Method SAI Our Method
Lena 34.67 34.59 2.834 s 0.682 s
Peppers 32.09 32.32 2.826 s 0.686 s
Pentagon 29.34 29.26 11.668 s 4.294 s
Airport 28.53 28.41 11.653 s 4.228 s

TABLE III
D EVICE INFORMATION

Tesla K80
Compute Capability 3.7
Framework Kepler
Number of Cores 2 × 2496
Memory 2 × 12 GB
Default Shared Memory per Block 49,152 B
Register File per Block 65,536

significantly reduces the amount of computation required. In

Fig. 4. Algorithm flowchart.
SAI, the sizes of C and D are 21 × 12 and 21 × 21, which
is much larger than the matrixes in Eq. (10) and Eq. (11).
In our method, we can only get one pixel in a slip window executed by the device [19], [20]. In CUDA programs, the
when preparing for GPU parallel computing. This means that GPU is responsible for the parallel aspect, called the kernel.
more data is overlapped between adjacent training windows. Today, parallel technology is used in all aspects of science and
However in SAI, four pixels are estimated simultaneously in technology, and the GPU-based image processing model [21-
a larger slip window, which reduces data overlap. Although 24] is growing in popularity. Parallel computing has shown
we have made mathematical improvements, the volume of enormous potential in image processing, which is significant
computation is still very large. So, in the next section, we with respect to some time-consuming methods. The CUDA
optimize the autoregressive interpolation model based on the code usually involves three steps: transfer the data into the
GPU. We use many CUDA optimization strategies to make memory of the GPU, run the CUDA kernel, and transfer the
full use of the GPU, each of which helps to effectively solve results from the GPU into the memory of the CPU. In some
these time-consumption problems. The PSNR and time will applications, the results may not need to be transferred back
be given in each optimization. to the CPU, but will be directly displayed on the GPU [25].
The diagonal interpolation process is as follows:
III. GPU- BASED O PTIMIZATION FOR I MAGE
1) Transfer the image data into the GPU and mirror repli-
AUTOREGRESSIVE I NTERPOLATION M ODEL
cation for edges.
A. Experimental Device Description 2) Generate a 4-connected neighbor pixels matrix and 8-
We performed GPU-based optimization for the image au- connected neighbor pixels in a slide training window
toregressive interpolation model on an NVIDIA Tesla K80, for training directional coefficients.
which comprise two GPUs, containing 2 × 2496 CUDA cores, 3) Construct matrixes H and G with the directional coef-
and a total of 26 SMX units. Compared with its previous GPU ficients of the previous step.
architecture, the number of registers and the shared memory 4) Calculate the pixel to be interpolated using Eq. (10).
in the Tesla K80 have been significantly improved to nearly 5) Deliver the results from the GPU to the memory of the
twice the original. GPUs are advantageous because of their CPU.
low cost and efficient power usage. In our experiments, the In the following optimization algorithms, the target image
serial C code ran on a CPU (Intel core (TM) i7-920, 2.66 data in our experiments were Lena, Peppers, Pentagon, and
GHz) with -O2 compiling optimization. Airport-the most classical image data in image processing. For
When running the CUDA program, the default total amount Lena, the HR image size is 512 × 512, so the size of the LR
of shared memory in every block in the Tesla K80 was 49,152 image down-sampled from the HR image is 256 × 256, as
B, and the total number of registers available per block was shown in Table I.
65,536. The specifications for our device are listed in Table III. There are two steps in our program and all processes are
carried out in parallel in the GPU. When diagonal interpolation
B. Algorithm Flow based on CUDA is complete, the image data is not delivered to the CPU but
CUDA C provides a simple path for users familiar with is prepared for vertical interpolation. In fact, before each
the C programming language to easily write programs to be interpolation, mirror replication for edges is needed. But this

1551-3203 (c) 2017 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See [Link] for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TII.2017.2724205, IEEE
Transactions on Industrial Informatics
IEEE TRANSACTIONS ON INDUSTRIAL INFORMATICS (ACCEPTED MANUSCRIPT: TII-17-0340) 5

limited space and bank conflict. The largest size of the shared
memory in every block is 49,152 bytes. Any read or write
requests made to memory with n addresses, and that fall
into n distinct memory banks, can therefore be accessed
simultaneously [19], [20], so it is important to avoid bank
conflicts as much as possible. The register memory is private
with respect to certain threads. That is, there is no bank conflict
in the register, so it is usually faster than shared memory. In
the Tesla K80, the L1 cache and shared memory share a 128-
KB configurable on-chip cache in each SMX. All SMXs share
an L2 cache with a maximum size of 1.5 MB. In our program,
the number of registers available per block is 65,536 and the
size of each register file is 32 bits, which is very valuable. In
the following optimization process, these operations on small
matrixes are completed in the register, as shown in Fig. 5. The
basic idea of the parallel scheme in our experiment is to use
Fig. 5. Correspondence between the data and CUDA threads. one thread to process one slide window. Since the L1 cache
can be used to deal with register spills, we use the command
“cudaFuncCachePerferL1” to set more of the L1 cache, which
part is relatively simple, almost no impact on the program, we means that the L1 cache will be set to 48 KB in each SMX.
will not focus on this part. As vertical interpolation involves The correspondence between the data and the thread is shown
the same process as diagonal interpolation, here we introduce in Fig. 5. Coordinating the shared memory with the register
only the diagonal interpolation case, as shown in Fig. 4. is a classic optimization technique in the CUDA program.
As the simple method, such as bicubic interpolation, suffice Through this effort, we obtained a significant speedup of
to interpolate smooth area. If the local variances estimated 79.3 for Lena.
from LR pixels is above a threshold, it will be interpolated In this version, there are 8 × 8 threads in one block, so
based on autoregressive. If the local variance is below a we must deal with 8 × 8 unknown pixels in a block. For
threshold, it will be interpolated based on bicubic, which each unknown pixel, we need to build a 6 × 6 slide window,
can reduce the amount of computing, without sacrificing the so the image data in one block is then 13 × 13 with the
effectiveness. And we empirically find that the determinants boundary overlap, as shown in Fig. 5. Since we do not know
of the training Matrix A in Eq. (3) and the training matrix B the C code details of SAI and have an executable program
in Eq. (4) play a key role. The above is the implementation only, here we only compare our single-threaded CPU program
details of this paper. But the core problem of this optimization with our GPU program in this GPU-parallel optimization
has not changed, we have to solve the Eq. (10). section. The parallel matrix optimization has many excellent
algorithms [26], [27], [28]. Considering various matrix sizes
in our algorithm, including 4 × 4, 5 × 9, and 1 × 14, methods
C. Optimization based on CUDA based on a fine-grained model are too complicated for our
1) Optimization with Shared Memory and Register: GPU method, as the matrix size is not fixed, causing the fine-grained
memory plays an important role in the CUDA optimization thread division becomes very difficult, but will be considered
processes. The latency associated with reading and writing in our follow-up work. Table IV shows time comparisons to
data from global memory is a bottleneck in some massively single-threaded C CPU code after using the shared memory
CUDA-dependent programs. Compared with the global mem- and register with the 8 × 8 block size.
ory, the shared memory and register have a higher memory Considering various matrix sizes in our algorithm, including
bandwidth and lower access latency. For faster access, the 4×4, 5×9, and 1×14, methods based on a fine-grained model
original LR image data down-sampled from the HR image are too complicated for our method, but will be considered in
is first translated to shared memory. The correspondence our follow-up work. Table III shows time comparisons after
between the data and CUDA threads are shown in Fig. 5. We using the shared memory and register with the 8 × 8 block
copy the data in each block to shared memory, as shown in size.
Fig. 5. The image data in each block will firstly be copied to 2) Finding the Most Suitable Block Size: From subsec-
__shared__ Image[13][13]. According to these data, tion III.C.1, we know that the resources of shared memory and
64 pixels will be estimated with the block size of 8 × 8 in register in the Tesla K80 are scarce and expensive. As such,
this version. In order not to form a boundary, there is a data the number of threads per block must be reasonable. If the
overlap between blocks. Coordinating the shared memory with block size is too small, resources will not be fully utilized.
the register is a classic optimization technique in the CUDA However, if the block size is too large, this will reduce the
program. degree of parallelism. Usually, the CUDA program achieves
All threads in the same block can communicate with each good performance when the number of threads per block is a
other, so all threads can share the LR data in one block. multiple of 16. But we must also consider that the optimal
However, we must take into account the shared memory’s block size is not always the same, and is closely related

1551-3203 (c) 2017 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See [Link] for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TII.2017.2724205, IEEE
Transactions on Industrial Informatics
IEEE TRANSACTIONS ON INDUSTRIAL INFORMATICS (ACCEPTED MANUSCRIPT: TII-17-0340) 6

TABLE IV TABLE V
S PEEDUP FROM USING SHARED MEMORY AND REGISTER WITH THE 8×8 S PEEDUP USING ADT
BLOCK SIZE
Sample PSNR (dB) Running Time Speed
Sample PSNR (dB) Running Time Speed Images CPU GPU CPU GPU Up
Images CPU GPU CPU GPU Up
Lena 34.59 34.59 0.682 s 7.72 ms 88.3
Lena 34.59 34.59 0.682 s 8.60 ms 79.3 Peppers 32.32 32.32 0.702 s 8.11 ms 86.5
Peppers 32.32 32.32 0.702 s 9.12 ms 76.9 Pentagon 29.26 29.26 4.294 s 46.55 ms 92.2
Pentagon 29.26 29.26 4.294 s 51.93 ms 82.6 Airport 28.41 28.41 4.228 s 46.21 ms 91.4
Airport 28.41 28.41 4.228 s 51.87 ms 81.5

Fig. 9. Comparison of single GPU and multi-GPU.

Fig. 6. Speedup comparison between block sizes for sample image data.

to all the GPU resources. To obtain the best performance, engine). With this architecture, ADT can be managed through
we conducted numerous experiments with respect to speedup streams, where a stream is a sequence of commands that are
using different block sizes. As shown in Fig. 6, a block size executed in order. In this version, the LR image data is divided
of 128 achieved the best performance in our program, with a into four pieces and four streams are created accordingly, with
speedup of 86.4 for Lena, 83.9 for Peppers, 87.6 for Pentagon, one stream responsible for one piece of data. Data transfer and
and 87.2 for Airport. data computing can be performed simultaneously with the two
3) Further Optimization Using Asynchronous Data Trans- different engines. If the volume of data is large but the required
fer: In our program, LR image data must be copied from the computing is not complicated, the effect of optimization by
host device to the GPU device, and then the HR image data ADT will be significant, as shown in Fig. 7. In this version,
is stored back in the host memory, which involves a relatively four streams are created in the GPU.
large overhead when the volume of image data is large. However, ADT is not suitable for use with all programs.
Using asynchronous data transfer (ADT), we can effectively If the transmission time between the host and the device is
hide the host-device transmission time in our program. There insufficient to hide the time management for the streams, ADT
are two memory copy engines in the Tesla K80: a host-to- is useless. In our method, the data transfer time is smaller than
device engine (H2D engine) and a device-to-host engine (D2H the kernel execution time and when streams are being used the
new edges must be mirrored. In this case, the speed is only
slightly improved to an 88.3 speedup for Lena, as shown in
Table V.
4) Optimization using a Multi-GPU: There are two GPUs
in the Tesla K80. A multi-GPU is suitable for single-node
desktop systems with multiple computing devices. So we can
fully use the hardware resources afforded by multi-GPU tech-
nology. As the ADT strategies introduced in subsection III.C.3
did not achieve the desired performance, we cannot use this
optimization approach to simplify the programming. However,
since one GPU relies on one stream, we create only two
streams in this version. The core code for a multi-GPU is
Fig. 7. Asynchronous data transmission using the Tesla K80. shown in Fig. 8. Each GPU handles half of the image data,

1551-3203 (c) 2017 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See [Link] for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TII.2017.2724205, IEEE
Transactions on Industrial Informatics
IEEE TRANSACTIONS ON INDUSTRIAL INFORMATICS (ACCEPTED MANUSCRIPT: TII-17-0340) 7

Fig. 8. Core code for a multi-GPU.

TABLE VI
S PEEDUP USING A MULTI -GPU

Sample PSNR (dB) Running Time Speed

Images CPU GPU CPU GPU Up
Lena 34.59 34.58 0.682 s 4.63 ms 147.3
Peppers 32.32 32.32 0.702 s 4.90 ms 143.2
Pentagon 29.26 29.26 4.294 s 26.89 ms 159.6
Airport 28.41 28.41 4.228 s 26.82 ms 157.6

and different devices execute each instruction in parallel. After

the final optimization, the speedup performance improved (a) (b)
significantly: up to 147.3 times for the Lena image, as shown
in Table VI. After the final optimization, PSNRs almost have
no change, only with a decrease of 0.01dB for Lena. Fig. 9
shows a comparison of a single GPU (without ADT) and a
multi-GPU. Compared to a single GPU, a multi-GPU has the
obvious advantage.

IV. E XPERIMENTAL R ESULTS

In section III, we explained how we achieved high speedup
performance compared to a serial C code running on a
CPU. To evaluate the proposed interpolation algorithm, we (c) (d)
conducted numerous experiments, which are described in this
section. Figure 10 shows the test images we used in this Fig. 10. Four sample images in the test set. (a) Lena, (b) Peppers, (c)
Pentagon, and (d) Airport.
simulation, and Fig. 11 shows the HR images of these test
images after using our interpolation method. Figs. 10-12 show
that both methods (SAI and our proposed method) achieved
good visual effect, and there is almost no difference in visual TABLE VII
quality between them. Table VII shows the PSNR results of PSNR C OMPARISON BETWEEN TWO METHODS : SAI AND OUR METHOD
two interpolation methods (SAI and the proposed method with
the final optimization) for the test images. As the target CPU Image SAI Our method Gain
have 4 cores, we also give the comparison between multi- Lena 34.67 34.58 -0.09dB
thread on CPU and K80 with our method, shown in Table VIII. Peppers 32.09 32.32 +0.23dB
From Table VIII, we can know that when four-thread is Pentagon 29.34 29.25 -0.09dB
used in program, it does not make full use of 4-core CPU. Airport 28.53 28.4 -0.13dB
So we create eight threads to take full advantage of the CPU

1551-3203 (c) 2017 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See [Link] for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TII.2017.2724205, IEEE
Transactions on Industrial Informatics
IEEE TRANSACTIONS ON INDUSTRIAL INFORMATICS (ACCEPTED MANUSCRIPT: TII-17-0340) 8

TABLE VIII
C OMPARISON BETWEEN MULTI - THREAD ON CPU AND K80 WITH OUR METHOD

Running Time Speed Up

Sample
Images CPU (one CPU (four CPU (eight GPU vs GPU vs
GPU
-thread) -thread) -thread) One-thread Eig-thread
Lena 0.682 s 0.254 s 0.211 s 4.63 ms 147.3 45.5
Peppers 0.702 s 0.265 s 0.218 s 4.90 ms 143.2 44.4
Pentagon 4.294 s 1.603 s 1.176 s 26.89 ms 159.6 43.7
Airport 4.228 s 1.543 s 1.154 s 26.82 ms 157.6 43

(a) (b) (a) (b)

(c) (d) (c) (d)

Fig. 11. Reconstructed images using our proposed method. (a) Lena, PSNR Fig. 12. Reconstructed images using SAI. (a) Lena, PSNR = 34.67 dB, (b)
= 34.58 dB, (b) Peppers, PSNR = 32.32 dB, (c) Pentagon, PSNR = 29.25 dB, Peppers, PSNR = 32.09 dB, (c) Pentagon, PSNR = 29.34 dB, and (d) Airport,
and (d) Airport, PSNR = 28.40 dB. PSNR = 28.53 dB.

core. Table VI and Table VIII show that the speedup of our
method is maintained at about 140 times compared to the
single-threaded and about 40 times compared to eight-threaded
C code running time on the CPU. According to the results,
performance improvement based on eight-thread compared to
one-thread on CPU is very obvious, but its performance is
still far less than multi-core GPU. From Table VII, we can
see that the PSNR results of our method are slightly worse
than those of SAI, as it misses a few details in the smaller
training window. But for the Peppers image, we achieved a
higher PSNR result, as the SAI may result in over fitting in this
image. A more reasonable approach is to select the training
window adaptively. In different small localities, the size of the Fig. 13. Image interpolation (scaling factor of 2.5) based on our proposed
training window can be changed adaptively according to the method.
LR pixel covariance. But if the size of the training window
is unfixed, it will lead to an imbalance in the thread load
in the GPU-parallel optimization. Further research is needed
to develop a GPU-based fine-grained implementation for the
autoregressive interpolation algorithm that balances the thread load.

1551-3203 (c) 2017 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See [Link] for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TII.2017.2724205, IEEE
Transactions on Industrial Informatics
IEEE TRANSACTIONS ON INDUSTRIAL INFORMATICS (ACCEPTED MANUSCRIPT: TII-17-0340) 9

(a) (b) (c)

Fig. 14. (a) Original image 1280 × 720 in size; (b) Interpolated image
2560 × 1440 in size by GPU; (c) Interpolated image 2560 × 1440 in size
by CPU.

V. GPU-PARALLEL AUTOREGRESSIVE I NTERPOLATION

A LGORITHM FOR I MAGE WITH A RBITRARY FACTOR AND
FOR A 720 P V IDEO

In the above analysis of our proposed method, the inter-

Fig. 15. Speedup (on the left) and FPS (on the right) for a 720p video with
polation factor was two. For a natural image, arbitrary factor GPU-parallel autoregressive interpolation algorithm on a Tesla K80.
interpolation is applied more widely. For an arbitrary scaling
factor (not an integer), the SAI-based algorithm can also work
well [10], by first using the SAI algorithm to expand the input improves. When the GPU resources are exhausted, the speedup
image by an integer time, and then applying a conventional will remain stable, as shown in Fig. 13. For this 720p video, we
image interpolation algorithm. In [10], SAI combined with achieved a speedup of 174.8 times and an FPS of 25.27, which
bilinear interpolation achieved promising results. In this sec- meet the requirements of real-time reconstruction. For typical
tion, we apply the GPU-parallel autoregressive interpolation deinterlacing videos, the amount of computing is smaller so
algorithm in this arbitrary factor image interpolation. Fig. 13 the frame rate can be further improved.
shows two examples with a scaling factor of 2.5 for Lena (for
a size from 256 × 256 to 640 × 640) and Airport (for a size
from 512 × 512 to 1280 × 1280), and Table IX shows the
speedup for four sample images using our proposed method. VI. C ONCLUSION
From Table IX, since the volume of computation for linear
interpolation is far less than for autoregressive interpolation, Image interpolation based on autoregressive models has
there is almost no change in the speedup values from those in achieved significant improvement in visual quality. But it is
Table VI, with the speedup for Lean, Peppers, Pentagon and difficult to apply in actual ultra HD products as it is complex
Airport remaining above 40 times for eight-threaded C CPU and time-consuming. In our proposed approach, we modified
code. the training window of the SAI, which is suitable for GPU-
We also make an experiment about color image with Lena parallel optimization. In addition, we added some mathemat-
(scaling factor of 2), shown in Table 10. As the color image ical optimizations to the matrix operations. To optimize our
has three channels of R, G and B, we decompose those three algorithm, we used various CUDA optimization strategies and
channels of the LR image firstly. Then, estimate unknown used them in a Tesla K80. For example, we used the shared
pixels in HR image based on our method. Finally, those results memory and register to reduce the latency of the data access.
can be combined into color images. We used ADT to hide the host-device transmission time. To
Image processing is the foundation of videos, and many make full use of the Tesla K80, we also used multi-GPU
excellent image processing algorithms are also used in videos. technology. Compared with SAI, the PSNR of our algorithm
To improve the quality of video reconstruction, taking the is similar. Experimental results show that, compared with a
processing speed into account, we applied the GPU-parallel single-threaded CPU counterpart in a 2.97-GHz Intel Core i7
autoregressive interpolation algorithm to video processing. To host computer, our algorithm achieves a high speedup of up
test the video processing performance, we adopted a 720p to 147.3 times for the Lena image. At the same time, our
video for our experiment (Fig. 14). Figure 15 shows the algorithm maintains a high PSNR and high subjective visual
speedup results compared to single-threaded CPU and frames quality.
per second (FPS). With respect to image interpolation with a scaling factor
Each frame in the video is a separate image. In Table V, of 2.5, we achieved a speedup of 144.6 for the Lena image.
we showed that the processing time for Lena with a size of For a 720p video, we achieved a speedup of 174.8 times and a
512 × 512 (scaling factor of 2) is less than 5 milliseconds. FPS above 24, which meet the requirements of real-time image
For this 720p video we also achieved a significant effect, reconstruction. These are promising results that support the use
with the FPS remaining above 24. If the FPS remains above of autoregressive interpolation model-based image processing
24, the video is played without interruption. The advantage algorithms. To further expand the degree of concurrency, fine-
of a GPU is the ability to process large amounts of data in grained implementation will be addressed in our follow-up
parallel. With increasing amounts of data, the speedup slowly research.

1551-3203 (c) 2017 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See [Link] for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TII.2017.2724205, IEEE
Transactions on Industrial Informatics
IEEE TRANSACTIONS ON INDUSTRIAL INFORMATICS (ACCEPTED MANUSCRIPT: TII-17-0340) 10

TABLE IX
S PEEDUP FOR IMAGE INTERPOLATION ( SCALING FACTOR OF 2.5) BASED ON OUR PROPOSED METHOD

Sample Running Time Speed Up

Image CPU (one CPU (eight GPU vs GPU vs
GPU
-thread) -thread) One-thread Eig-thread
Lena 0.687 s 0.218 s 4.75 ms 144.6 45.8
Peppers 0.708 s 0.226 s 5.04 ms 140.4 44.8
Pentagon 4.312 s 1.214 s 27.16 ms 158,7 44.6
Airport 4.248 s 1.187 s 27.08 ms 156.8 43.8

TABLE X
S PEEDUP FOR RGB IMAGE OF L ENA

Running Time PSNR (dB) Speed Up

Image CPU (One CPU (Eig GPU vs GPU vs
GPU CPU GPU
-thread) -thread) One-thread Eig-thread
Lena 2.129 s 0.640 s 14.62 ms 32.96 32.95 145.6 43.7

R EFERENCES fast feature embedding,” in Proc. ACM International Conference on

Multimedia, 2014, pp. 675–678.
[1] R. G. Keys, “Cubic convolution interpolation for digital image process- [18] J. Wu, L. Deng, G. Jeon, and J. Jeong, “GPU-parallel interpolation using
ing,” IEEE Trans. Acoustics, Speech and Signal Processing, vol. 29, no. the edge-direction based normal vector method for terrain triangular
6, pp. 1153–1160, 1981. mesh,” Journal of Real-Time Image Processing, pp. 1–10, 2016.
[2] H. S. Hou and H. Andrews, “Cubic splines for image interpolation and [19] NVIDIA, CUDA Programming Guide, USA, 2015.
digital filtering,” IEEE Trans. Acoustics, Speech and Signal Processing, [20] NVIDIA, CUDA C Best Practices Guide, USA, 2015.
vol. 26, no. 6, pp. 508–517, 1978. [21] Z. Yang, Y. Zhu, and Y. Pu, “Parallel image processing based on CUDA,”
[3] X. Li and M. T. Orchard, “New edge-directed interpolation,” IEEE in Proc. International Conference on Computer Science and Software
Trans. Image Processing, vol. 10, no. 10, pp. 1521–1527, 2001. Engineering, 2008, vol. 3, pp. 198–201.
[4] L. Zhang and X. Wu, “An edge-guided image interpolation algorithm [22] I. K. Park, N. Singhal, M. H. Lee, S. Cho, and C. W. Kim, “Design
via directional filtering and data fusion,” IEEE Trans. Image Processing, and performance evaluation of image processing algorithms on GPUs,”
vol. 15, no. 8, pp. 2226–2238, 2006. IEEE Trans. Parallel and Distributed Systems, vol. 22, no. 1, pp. 91–
[5] X. Wu and X. Zhang, “Image interpolation using texture orientation map 104, 2011.
and kernel fisher discriminant,” in Proc. ICIP, 2005, vol. 1, pp. I–49–52. [23] M. Kraus, M. Eissele, and M. Strengert, “GPU-based edge-directed
image interpolation,” in Proc. SCIA 2007, pp. 532–541.
[6] X. Zhang and X. Wu, “Image interpolation by adaptive 2-D autore-
[24] Y. Kui-Ying, J. Lin, Y. Jun-Peng, and X. Lu-Ping, “Processing piece-
gressive modeling and soft-decision estimation,” IEEE Trans. Image
wise autoregressive model image interpolation algorithm on GPU with
Processing, vol. 17, no. 6, pp. 887–896, 2008.
CUDA,” in Proc. International Conference on Wireless Communications
[7] C. Lee, S.-W. Kim, and C. Yoo, “VADI: GPU virtualization for an
and Signal Processing, 2011, pp. 1–4.
automotive platform,” IEEE Trans. Industrial Informatics, vol. 12, no.
[25] J. Wu, L. Deng, and A. Paul, “3D terrain real-time rendering method
1, pp. 277–290, Feb. 2016.
based on CUDA-OpenGL interoperability,” IETE Technical Review, vol.
[8] X. Wu, X. Zhang, and X. Wang, “Low bit-rate image compression via 32, no. 6, pp. 1–8, 2015.
adaptive down-sampling and constrained least squares upconversion,” [26] G. Sharma, A. Agarwala, and B. Bhattacharya, “A fast parallel Gauss
IEEE Trans. Image Processing, vol. 18, no. 3, pp. 552–561, 2009. Jordan algorithm for matrix inversion using CUDA,” Computers &
[9] K. Tang, O. C. Au, L. Fang, Z. Yu, and Y. Guo, “Image interpolation Structures, vol. 128, pp. 31–37, 2013.
using autoregressive model and gauss-seidel optimization,” in Proc. [27] L. B. With Jensen, A. Kjær-Nielsen, K. Pauwels, J. B. Jessen, M.
ICIG, 2011, pp. 66–69. Van Hulle, and N. Krüger, “A two-level real-time vision machine
[10] K. Tang, O. C. Au, Y. Guo, J. Pang, J. Li, and L. Fang, “Arbitrary factor combining coarse-and fine-grained parallelism,” Journal of Real-Time
image interpolation by convolution kernel constrained 2-D autoregres- Image Processing, vol. 5, no. 4, pp. 291–304, 2010.
sive modeling,” in Proc. ICIP, 2013, pp. 996–1000. [28] L. Chen, O. Villa, and G. R. Gao, “Exploring fine-grained task-
[11] W. Dong, L. Zhang, R. Lukac, and G. Shi, “Sparse representation based based execution on multi-GPU systems,” in Proc. IEEE International
image interpolation with nonlocal autoregressive modeling,” IEEE Trans. Conference on Cluster Computing, 2011, pp. 386–394.
Image Processing, vol. 22, no. 4, pp. 1382–1394, 2013.
[12] L. Mussi, F. Daolio, and S. Cagnoni, “Evaluation of parallel particle
swarm optimization algorithms within the CUDAT M architecture,”
Information Sciences, vol. 181, no. 20, pp. 4642–4657, 2011.
[13] J. Wu, Z. Song, and G. Jeon, “GPU-parallel implementation of the edge-
directed adaptive intra-field deinterlacing method,” Journal of Display
Technology, vol. 10 no. 9, pp. 746–753, 2014.
[14] B. Huang, J. Mielikainen, H. Oh, and H.-L. A. Huang, “Development of
a GPU-based high-performance radiative transfer model for the Infrared
Atmospheric Sounding Interferometer (IASI),” Journal of computational
Physics, vol. 230, no. 6, pp. 2207–2221, 2011.
[15] N. Zhang, Y.-s. Chen, and J.-l. Wang, “Image parallel processing based
on GPU,” in Proc. ICACC, 2010, vol. 3, pp. 367–370.
[16] X. Jia, Y. Lou, R. Li, W. Y. Song, and S. B. Jiang, “GPU-based fast
cone beam CT reconstruction from undersampled and noisy projection
data via total variation,” Medical Physics, vol. 37, no. 4, pp. 1757–1760,
2010.
[17] Y. Jia, E. Shelhamer, J. Donahue, S. Karayev, J. Long, R. Girshick,
S. Guadarrama, and T. Darrell, “Caffe: Convolutional architecture for

1551-3203 (c) 2017 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See [Link] for more information.