Acceleration of Variance of Color
Differences-Based Demosaicing Using CUDA
Muhammad Ismail Faruqi∗ , Fumihiko Ino† , and Kenichi Hagihara‡
Graduate School of Information Science and Technology, Osaka University
1-5 Yamadaoka, Suita, Osaka 565-0871, Japan
Email: {i faruqi∗ , ino† , hagihara‡ }@ist.osaka-u.ac.jp

Abstract—Image demosaicing algorithms are used to reconstructa full color image from the incomplete color samples output(RAW data) of an image sensor overlaid with a Color Filter
Array (CFA). Better demosaicing algorithms are superior interms of acuity, dynamic range, signal to noise ratio, and artifact
suppression, which make them suitable for high quality delivery
such as theatrical broadcast.
In this paper, we present our efforts in examining the feasibilityof exploiting the Graphics Processing Unit (GPU) as an emerging
accelerator to create an on-the-fly implementation of Variance ofColor Differences (VCD) demosaicing, a state-of-the-art heuristic
demosaicing algorithm developed to eliminate false-color artifacts
in texture region of images.
Our efforts in this paper are 1) implementing the algorithm as
several kernels to separate the bottleneck portion of the algorithm
from the rest and to minimize idle threads and 2) reducing
I/O between shared and global memory when performing green
channel interpolation by separating the input RAW data into fourchannels. We then compare the implementation featuring both
acceleration methods with a single kernel implementation. Basedon experimental results, our proposed acceleration methods
achieved per-frame processing time of 343 ms on an nVidia GTX
480, which translates into 2.95 fps. Additionally, our proposed
methods were also able to accelerate the kernel time and theeffective memory bandwidth by a factor of 2.1x compared with
its single kernel counterpart.
Keywords—Parallel processing; image demosaicing; CUDA; GPU

I. I NTRODUCTION
Image demosaicing [1] is an integral part of color imaging
pipeline. It is the first step of the pipeline, in which theluminance data known as RAW data in each photosite isexpanded into RGB values. One method to evaluate the qualityof demosaicing algorithms is to measure how effective they
approximate the remaining color values while not introducing
artifact known as moiré. The less moiré introduced in the demosaiced images, the better quality the demosaicing algorithm
is.
However, sophisticated demosaicing algorithms tend to becomputationally expensive and impractical to be implemented
on almost all digital cameras. Therefore, demosaicing algorithms used in many digital cameras are generally favoring
execution speed over quality, resulting pictures with moiré.
Hence, to enable users accessing higher quality images,

higher-end cameras offer a feature to directly record RAW data
into storage without going through the entire imaging pipeline.
Such data is then processed by some external processors
capable executing demosaicing in shorter time.
In another spectrum, users of video cameras equipped with
Bayer filter [2] have suffered by the long time and huge spacerequired to acquire high quality files from their camera. To
process the files, users have to wait for a long time until
the demosaicing algorithm finishes processing the frames, andthen saving the demosaiced file into storage. Here lie thechallenges this paper is going to solve: to enable users acquire
high quality demosaicing result from high resolution video
files quickly without having to store large demosaicing result
into storage.
In this paper, we propose an acceleration of Variance of Color
Differences (VCD) [3]-based demosaicing, a high quality
demosaicing algorithm specifically developed to combat moiréin texture region of images, using Compute Unified Device
Architecture (CUDA) [4]. The objective of this implementation is to demosaic video RAW files on-the-fly as fast as
possible, so that the video editing workflow will be accelerated
and the storage requirement to work on demosaicing result
can be eliminated. To achieve this, we first introduce thewavefront processing as the base method of the algorithm
parallelization. We then propose two implementation methods,
which are 1) implementing the algorithm as multiple kernels
to separate the bottleneck portion of the algorithm and tominimize idling threads, and 2) reducing input and output
transfer between global and shared memory [4] in the green
channel interpolation phase by separating the input RAW data
into separate channels.
II. R ELATED W ORK
Several works have been dedicated to implement demosaicing
using GPU. For example, McGuire [5] accelerated MalvarHe-Cutler [6] image demosaicing algorithm using OpenGL in
real-time speed. Fung et al. [7] show two examples of CUDAbased demosaicing based on bilinear and Lanczos [8] method.
However, these algorithms have inferior moiré suppression
compared to Chung’s algorithm. The first commercial application known to deliver real-time preview and grading for 4K
RAW was IRIDAS [10] Speedgrade XR which was launchedin 2009. However, based on its high speed, we suspect it
to perform a demosaicing algorithm with low complexity,
although the exact algorithm is closed.
Meanwhile, recent high-quality demosaicing algorithms include Mairal et al. [12] algorithm. This algorithm along with
VCD has a high CPSNR figure. However, Menon et al. [13]benchmarks recent demosaicing algorithms including Mairal’sand Chung’s algorithm, and we found Chung’s algorithm has
the best CIELab [14] figure which closely resembles humaneye evaluation result.
III. VARIANCE OF C OLOR D IFFERENCES (VCD)
D EMOSAICING (C HUNG ’ S A LGORITHM )This section explains VCD-based demosaicing. VCD is ademosaicing algorithm developed by Chung and Chan [3] tosuppress artifact known as moiré when demosaicing Bayer
CFA. It extends Adaptive Color Plane Interpolation (ACPI)
algorithm [9] by heuristically determining the interpolation
direction of the green channel. In this paper, we describe
Chung’s algorithm as an algorithm consisting of four steps
explained from section III-A to section III-D.
For the rest of the explanation of Chung’s algorithm, we will
denote the value of a sampling position from the input RAW
data at (x, y) coordinate as Px,y . Each Px,y is overlaid by acolor filter element whose channel type is either red, green, or
blue. The output of the algorithm is an image with complete
information in each sampling position, where the red, green,
and blue channel values at (x, y) will be addressed as rx,y ,gx,y , and bx,y respectively.
The problem of a demosaicing algorithm for Bayer CFA is tointerpolate two missing channel values in each Px,y . It needs
to interpolate gx,y and bx,y at red CFA sampling positions,rx,y and bx,y at green CFA sampling positions, and rx,y andgx,y at blue CFA sampling positions. By completing all of
those channels, we construct a full color image from a RAW
image with Bayer CFA.

D
ĝx,y
=

Px−1,y + Px+1,y + Px,y−1 + Px,y+1
4
4Px,y − Px−2,y − Px+2,y − Px,y−2 − Px,y+2
+
.(3)
8

In the actual algorithm, those values are computed on-thefly while interpolating green channel. However, since thosevalues are used frequently in green channel interpolation
phase, moving those values calculation into precomputation
will reduce its execution steps from O(w + h) steps into O(1)steps.
B. Interpolation of Green Channel Values at Red and Blue
CFA Sampling Positions
The next step, which is the contribution of Chung and Chanis to heuristically interpolate green channel at red and blue
CFA sampling positions. First, they classify whether a CFA
sampling position belongs to a texture or to an edge. Thisis achieved by computing an edge classifier ex,y , which is
obtained from summing the difference of the neighboring
pixels in horizontal (LH ) and vertical (LV ) direction. The
values of ex,y , LH , and LV are obtained from (4), (5), and(6) respectively.
(
ex,y = max

LH =



∑

As the first step, the algorithm precalculates green channel values for each pixel on blue and red CFA sampling positions inHVDhorizontal (ĝx,y), vertical (ĝx,y), and diagonal (ĝx,y) direction.
They are obtained as follows:
2Px,y − Px−2,y − Px+2,yPx−1,y + Px+1,y
H
+
,
ĝx,y
=
2

V
ĝx,y
=

4

(1)

Px,y−1 + Px,y+1
2Px,y − Px,y−2 − Px,y+2
+
, (2)
2
4

,

(4)

|Px+dx,y+dy − Px,y+dy |,−2≤dx≤2,x̸=0

(5)

LV =



∑



∑



−2≤dx≤2

A. Precalculation of Three Green Channel Values at Red andBlue CFA Sampling Positions

)



∑



−2≤dy≤2

LH LV
,
LV LH

|Px+dx,y+dy − Px+dx,y |.

−2≤dy≤2,y̸=0(6)
The CFA sampling position is classified as an edge if ex,y ≥T , where T is a predefined threshold by user. Its green channel
value gx,y can then be interpolated according to (7).
{gx,y =ĝ,Hĝ,Vif LH < LV ,otherwise.

(7)

On the other hand, the CFA sampling position is classified as atexture if ex,y < T . Chung and Chan stated that for a texture,
three variance of color differences values can additionallysupplied to increase the accuracy of the interpolation direction
for this CFA sampling position. Those values are variance ofcolor differences in horizontal, vertical, and diagonal directions which are obtained by computing (8)–(10).

∑1∑((dx+i,y ) −(dx+j,y )),
9
9

(8)

∑
1∑
((dx,y+i ) −(dx,y+j )),
9
9

(9)

1
2
H σx,y =

i∈Ψ

y

y

j∈Ψ

1
2
V σx,y =

i∈Ψ

x

x

j∈Ψ

∑
1∑
((fx+i,y ) −(fx+j,y ))
2 9
9
i∈Ψ
j∈Ψ
1∑
1∑
+
((fx,y+i ) −(fx,y+j ))),
9
9

1 1
2
D σx,y = (

i∈Ψ

(a)

(10)

j∈Ψ

2
2
2
where H σx,y
, V σx,y, and D σx,ydenote variance of color
differences in horizontal, vertical, and diagonal directions
respectively. The range of those summations is defined by
Ψ = {0, ±1, ±2, ±3, ±4}．Furthermore, dx+i,y , dx,y+i ,
fx+i,y , and fx,y+i are defined with (11)–(14).

 Px+i,y − gx+i,y ,H,
Px+i,y − ĝx+i,y
dx+i,y = dx+i−1,y +dx+i+1,y2

 Px,y+i − gx,y+i ,VPx,y+i − ĝx,y+i,
dx,y+i = dx,y+i−1 +dx,y+i+12

 Px+i,y − gx+i,y ,D,
Px+i,y − ĝx+i,y
fx+i,y = fx+i−1,y +fx+i+1,y2

 Px,y+i − gx,y+i ,DPx,y+i − ĝx,y+i,
fx,y+i = fx,y+i−1 +fx,y+i+1
2,

if i = −4, −2,
if i = 0, 2, 4,
if i = ±1, ±3,

(11)

(12)

,

if i = −4, −2,
if i = 0, 2, 4,
if i = ±1, ±3,

(13)

,

if i = −4, −2,
if i = 0, 2, 4,
if i = ±1, ±3,,

if i = −4, −2,
if i = 0, 2, 4.
if i = ±1, ±3.

(14)Finally, the green channel value gx,y can be interpolated by
(15).
 H
 ĝ,
ĝ V
gx,y =
 ,D
ĝ,

2
2
2
2
if H σx,y= min(H σx,y,V σx,y,D σx,y),2222
if V σx,y= min(H σx,y,V σx,y,D σx,y),222
2
if D σx,y= min(H σx,y,V σx,y,D σx,y).

(b)

Figure 1. (a) The location of green CFA sampling position surrounded by red
CFA horizontally and blue CFA sampling positions vertically. (b) The locationof green CFA sampling position surrounded by blue CFA horizontally and red
CFA sampling positions vertically.

hx,y = gx,y +Px−1,y − gx−1,y + Px+1,y − gx+1,y
,
2

(16)

vx,y = gx,y +Px,y−1 − gx,y−1 + Px,y+1 − gx,y+1
.
2

(17)

Those two values are then assigned as interpolated blue orgreen channel value, according to the location of the sampling
position. If it is located as Fig. 1(a) shows, the missing red
and blue values will be assigned as rx,y = hx,y , bx,y = vx,y ,otherwise if it is located as Fig. 1(b) they will be assigned asbx,y = hx,y , rx,y = vx,y .
D. Interpolation of Red Channel Values on Blue CFA Sampling Positions and Blue Channel Values on Red CFA Sampling Positions
Finally, we interpolate the red channel value at blue CFA
sampling positions, and the blue channel value at red CFA
sampling positions. The equation for red channel value interpolation is described in (18). The value for the blue channel
at red CFA sampling positions is also computed in the same
fashion.

rx,y = gx,y +

1 ∑
4

∑

(Px+dx,y+dy −gx+dx,y+dy ) (18)dy=±1 dx=±1(15)

C. Interpolation of Red and Blue Channel Values on Green
CFA Sampling Positions
Using the result of the interpolated green channel, we cancompute the values of red and blue channels on green CFA
sampling positions. To do this, we interpolate hx,y value
horizontally and vx,y value vertically at (x, y) according to(16) and (17).

IV. PARALLELIZATION S TRATEGIES OF C HUNG ’ S
D EMOSAICING A LGORITHM
In this chapter, we will describe the parallelization strategy
we have gone through in accelerating Chung’s demosaicing
algorithm. We choose CUDA as the platform to accelerate
Chung’s algorithm since it has the data parallelism, where the
GPU has potential to perform it faster than the CPU.
As an overview, we start by explaining the wavefront method
which we believe is the natural parallelization method for

m=3
X-4

X-3

x-2

x-1

x

x+1

x+2

X+3

X+4

Block
(1,0)

y-4
y-3
y-2

Block
(0,1)

y-1

n= 3

y

Block
(m-1,n-1)

y+1
y+2
y+3

Figure 3. The execution area P is divided into m × n blocks. Here m =3, n = 3, α = 4, and β = 3. A block wave Dk consists of all blocks withthe same k = m + n. For example, block (0, 1) and (1, 0) belong to D1 .

Y+4

(a)

(b)
Figure 2. (a) An illustration of data dependencies when interpolating thegreen channel at Px,y . The light green box indicates interpolated green
channel values, and the arrows indicates the direction of the dependency.
(b) An example of wavefront processing. Here, one thread is responsible tointerpolate the green channel at one sampling position. Threads in the sameanti diagonal, i.e. thread wave, are executed in parallel. Meanwhile, eachthread waves are executed sequentially.

data dependency introduced by (11)–(14). We then propose
two strategies to accelerate Chung’s algorithm, which are 1)
separating data-dependent portion from the rest of computation
by implementing Chung’s method as multiple kernels and 2)
reducing idling threads during I/O transfer by separating input
data into several channels.
A. Parallel Processing by Wavefront Method
In the green channel interpolation phase of Chung’s algorithm,(11) and (13) require previous green channel interpolation
results of Px−2,y and Px−4,y , meanwhile (12) and (14) requireprevious interpolation result of Px,y−2 and Px,y−4 . Thisimplies a form of data dependency which is illustrated inFig. 2(a). As a result, Px,y is dependent to the green channelinterpolation results at Px−2,y , Px−4,y , Px,y−2 , and Px,y−4 .
Therefore, we hypothesize that the natural parallelization
method for this phase while guaranteeing the required data
dependency is by a wavefront processing method [15] which
is illustrated by Fig. 2(b).
To conform with CUDA’s hierarchical execution model [4],the implementation of the wavefront method for Chung’salgorithm in this paper consists of two levels: the grid-level
and the block-level. Fig. 3 illustrates the grid-level wavefront

processing. At the grid-level, the RAW data P with width
w and height h is divided into m × n blocks. Each blockhhas width α and height β, hence m = ⌈ w
α ⌉ and n = ⌈ β ⌉.
All blocks where k = m + n are grouped into a block wave
Dk . The application processes the block waves sequentiallystarting from top left to bottom right, where block wave
Dk−1 is processed before block wave Dk . Each block wave
is mapped into a CUDA kernel launch with one dimensional
grid size. In CUDA, synchronization between block waves isdone automatically between kernel launches.
Inside a block, if a pixel Px,y has ex,y < T , it is mapped intoa thread. All threads with thread index (tx , ty ) are grouped asa thread wave Tv , where v = tx + ty . The kernel consists of aloop in which each thread decides to process a pixel or not in
the wave depending on its thread index. The thread waves
are processed sequentially, where Tr−1 is executed before
Tr . Between each thread wave execution, synchronization is
performed once. Fig. 2(b) illustrates the block-level wavefront
processing. In CUDA, synchronization between thread waves
is accomplished manually by __synchthreads() function
provided by the CUDA Software Development Kit (SDK).
The computation of thread waves inside a block is performed
in α + β − 1 steps, while the computation of the entire blockwaves is performed in m + n − 1 steps. Hence, this CUDAbased wavefront method takes a total w + h + m(β − 1) +n(α − 1) − (α − β − 1) steps to interpolate the green channelat blue and red CFA sampling positions. In other words, it hasO(w + h) steps which will potentially become the bottleneckof the application.
B. Separation of Data-dependent Portion from the Rest ofComputation
Fig. 4 illustrates that although in the green channel interpolation step Px,y depends on the previous interpolation resultof its neighbors, the rest of steps are independent and can be
accelerated by an embarrassingly parallel method. Hence, if allof four Chung’s algorithm’s steps are implemented as a single
monolithic kernel, the performance of the phases without any
data dependency will be hindered by the bottleneck mentioned

Single kernel implementaon
Monolithic kernel

Fully parallelMulple kernel implementaon
Kernel 1

Kernel 2

Fully parallel

Data dependentmax(Agrn,Avar)

P

Q

R

max(Agrn,Avar)

Data dependentw

Fully parallelhKernel 3

Fully parallel

Avar

Avar Agrn

Agrn

Figure 4. When the algorithm is implemented as a single monolithic kernel,
the performance of data-independent portions of calculation will be reducedby the bottleneck portion of the algorithm.

C. Classification Based on Differences of Kernel Execution
Area

TABLE IF INAL KERNEL LIST.
Step
N/A
—
1
2
—
2
3&4
—

N/A

Remarks
Transfers data from host to global
memory.
Converts input from unsigned short
to float.
Precomputes temporary green
channels.
Performs edge detection.
Performs the channel separation of
input RAW data.
Calculate variance values and interpolate G channel value.
Interpolate the rest of the channels.
Combine per-channel float-typed
output data into unsigned short
RGBA data.
Transfers data from global to hostmemory.

Figure 5. Computation areas. R denotes the execution area of input loading
phase. Q denotes the computation area of precomputation of green channel
values phase, where P denotes the computation area of other phases.

Area
R

Kernel
N/A

R

1

Q

2

P
R

3
4

P

5

P
R

6
7

R

N/A

before. On the other hand, if the algorithm is implemented asseveral kernels, the bottleneck portion can be isolated into its
own kernel without hindering other steps’ performance.
Based on this fact, we propose to classify the computation
steps in Chung’s algorithm based on their data dependency.
The classification policies are:1) The computation phase with data dependency is separated from the rest and is implemented as one kernel.
2) The computation phases without any data dependency
are united together into one kernel.
Based on those policies, because there is only one step withdata dependency, we logically need to implement Chung’salgorithm as three CUDA kernels: 1) a kernel that performs
all computations before green channel interpolation phase, 2)
a kernel to perform the green channel interpolation step, and
3) a kernel that handles the rest of computations. However,
due to optimization schemes introduced in Sec. IV-C–IV-D,
the implementation is finally implemented as seven kernels asshown in Table I.

Besides the data dependency-based separation policies explained in Section IV-B, another factor that must be consideredwhen implementing a CUDA kernel is that if each kernel handles different computation areas, there will be idling threads
which will lead into kernel ineffectiveness. A computation
area is defined as an area in which a kernel will store its
computation result into.
Fig. 5 illustrates different computation areas of each step. The
H
V
computation area for step 1 which computes ĝx,y, ĝx,y, andDĝx,y is denoted by Q. Based on (11) and (13), dx,y and fx,yHDcalculations in step 2 require ĝx,yand ĝx,yin positive direction
with maximum range Avar = x + 4. Additionally, (12) andVD(14) implies that step 2 also requires ĝx,yand ĝx,yin positive
direction with maximum range Avar = y + 4. Hence, Q has
width w + Avar and height h + Avar .
Meanwhile, based on (1)–(3), to compute Q, we require valuesof neighboring pixel within a maximum radius of Agrn inhorizontal and vertical directions as input. In [3], Chungdefines Agrn = 2. Hence, the required data that must be loaded
into global memory [4] and converted to floating point typeto accommodate all calculation phases exists in area R, which
has width w + Avar + Agrn + max(Avar , Agrn ) and heighth + Avar + Agrn + max(Avar , Agrn ).
As mentioned above, if different areas are processed by thesame kernel, there will be idling threads when computing thephase with least area. Fig. 6 shows an example where threads
are idling in block (0, 0) and in block (m − 1, n − 1) when
computing P . This problem happens because the kernel withm × n thread blocks tries to process P , Q, and R in different
steps.
To solve this problem, we propose to do further classificationof the computation steps based on the following policies:
1) Separate the kernel which contains different computation
areas into multiple kernels according to its computationarea.
2) Combine any subsequent computation steps which have

Idle threads when compu!ng P

TABLE IIS PECIFICATION OF THE EXPERIMENT MACHINE .

0,0

P Q R
Idle threads when compu!ng P
m-1,n-1

Figure 6. Idle threads that will happen if P, Q, and R are processed by a same
kernel. Each cell represents a thread. Gray-colored cells are busy threads, andwhite-colored cells are idle threads.

the same computation area into one kernel.
Based on the policies above, since kernel 1, 2, and 3 on Table
I processes R, Q, and P , we decide to separate those process
into three kernels. Similarly, since kernel 6 performs two steps
described in Section III-C and III-D in P , we can combine
them into one kernel.
D. Reducing Input Loading Time of the Green Channel Interpolation Phase by Separating Input Data
We noticed that in (11)–(14), the only neighbors accessed in
input data and output data by Px,y only Px+α,y+α , where
α = ±2, ±4. Furthermore, all neighboring CFAs of Px+α,y+αhave the same CFA as the CFA being computed, and other
CFAs are not used.
Based on this fact, in order to save computation time fromloading unnecessary CFAs, we propose to separate all values in
input data at the same CFA into their own matrix in a separate
kernel. Specifically, the input RAW data P is separated into
four matrices R, B, GR, and GB in place. Matrix R will holdvalues from red CFA, matrix GR will hold green CFAs which
are surrounded horizontally by two red CFAs, matrix B willhold blue CFAs, and finally matrix GB will hold green CFAs
which are surrounded horizontally by two blue CFAs. For P
with width w and height h, the width and height of matrices
R, B, GR and GB will be w/2 and h/2, respectively.
Since kernel 5 in Table I that handles green channel interpolation phase is executed multiple times by the wavefront
processing, we hypothesize that this optimization will help
reduce its total execution time. In this paper, we implement
the input separation phase as a kernel which is launched after
the edge detection kernel is executed. This kernel is denoted as
Kernel 4 in Table I. As the final result of those optimizations,the algorithm is implemented as 7 kernels.
V. E XPERIMENTAL R ESULTS AND A NALYSIS
In this section, we present the experimental results of our
proposed CUDA-based implementation of Chung’s algorithm.
In order to evaluate our proposed acceleration methods, we
set several objectives for our experiments below.
1) Measuring the kernels’ performance by comparing execution times among our CUDA-based implementations
and a CPU-based of our optimization.

Item
Processor
MemoryGPU
GPU Memory
Operating System
CUDA Version
IDE

Value
Intel Xeon X5450, 3 GHz, 4 cores
DDR3 8 GB
nVidia GeForce GTX 480
GDDR5 1.5 GB
Windows XP 64-bit XP3
4.0
Visual Studio 2005

2) Measuring the efficiency of optimization methods interms of execution time, floating point operation persecond (FLOPS), and memory bandwidth.
The specifications of the machine used for executing theexperiments are shown in Table II.
As we do not have access to any video camera capable ofoutputting a RAW image with 4K resolution, we chose tosimulate it by using a RAW image from a still digital camera
with nearly similar resolution. The RAW image frame used as
input for our experiment has size of 4608 × 3072 pixels, anduses a Bayer type CFA.
To determine the optimum block width α, width β, andregister count r combination for each kernel, we performed
the following experiment. First, for each kernel we choose
several (α, β) whose static shared memory size [4] fits in a
device. Next, using each (α, β) we ran the kernel with anarbitrary number of registers using CUDA Visual Profiler [11]to obtain its static memory size. We then input the obtained
static memory size, α, and β into the CUDA Occupancy
Calculator [16] to get the biggest register count ropt that results
the highest occupancy number. This ropt will give the fastest
result for its corresponding α and β. Finally, we compare theresults of each (α, β, ropt ) combination and chose the fastestone.
A. Performance Evaluation
We used two CUDA implementations, which are a monolithic
kernel version and a multiple kernels version which implemented our acceleration methods. The first version performs
wavefront processing on all interpolation steps of Chung’salgorithm, while the multiple kernels version only performs
wavefront processing on the green channel interpolation step,and processes other steps using embarassingly parallel methods. We will refer the monolithic kernel version as the single
kernel version hereafter.
We start the evaluation by presenting Table III, which shows
the total execution time of our GPU implementations while
compared to the CPU implementation which uses a single
core. In this table, the total execution time consists of CPU
time and GPU time. Furthermore, the GPU time consists ofhost to device (D → H) memory transfer time, kernel time(s),and the device to host (H → D) memory transfer time.
Based on Table III, the total execution times of our single

TABLE III
D ETAIL OF TOTAL EXECUTION TIME FOR EACH IMPLEMENTATION .
Execution times (ms)
CUDA (multiple kernel) CUDA (single kernel) CPU
CPU time
84.5
136.9
5742.0
D→H
13.5
13.5
—
Kernel time
178.9
370.5
—
H→D
65.9
65.9
—
Total time
342.8
586.8
5742.0TABLE IV
M EMORY BANDWIDTH COMPARISON BETWEEN CUDA- BASED AND
CPU- BASED VCD DEMOSAICING IMPLEMENTATION .

Detailkernel and multiple kernel implementations were 587 ms
and 343 ms respectively, which are 9.9x and 16.8x fasterthan its CPU implementation, respectively. As for the CPU
implementation, with execution time of 5,742 ms per frame itis definitely not fast enough to demosaic a RAW video stream
consisting hundreds of thousands of 4K frames. Even if allof four X5450 cores is used to perform Chung’s algorithm
in parallel, the execution time would be still slower than our
CUDA-based implementations.
B. Efficiency Analysis of CUDA-based Implementations
Two metrics known as memory bandwidth and floating point
operations per second (FLOPS) are widely used to measure
the effectiveness of CUDA kernels. Given a measured kernel
bandwidth bM , peak GPU bandwidth bT , measured kernel
FLOPS performance fM , and peak GPU FLOPS fT , one cancalculate kernel’s bandwidth effectiveness EB = bM /bT andkernel’s FLOPS performance effectiveness EF = fM /fT .
When EB < EF the kernel is said to be memory-bound,
otherwise the kernel is said to be arithmetic-bound. Also,
when EB ≥ 1, the kernel is said to use the shared memory
effectively. The nVidia GTX 480 used for this experiment has
a bT of 177.4 GB/s and a fT of 1344.96 GFLOPS.
To evaluate the efficiency of the proposed acceleration methods, we chose to analyze and compare the memory bandwidthof both CUDA-based implementations and CPU-based implementations.

Kernel
1
2
3
4
5
6
7

Total time (ms)Multi Single CPU
0.9
6.2
5.5
1.6 370.5 5742.0
152.6
9.3
2.8

Bandwidth (GB/s)
EB (%)
Multi Single CPU Multi Single
144.1
81.2
122.5
69.0
363.4
200.4
101.4 7.3 42.6 57.1 4.1
15.2
8.5
90.3
50.9
131.8
74.2

efficiently, which are indicated by their high EB percentages.
For example, the EB of kernel 3 achieved 200.4%. We
observed that such high EB corresponds with the high memory
access count inside the kernel. According to (5) and (6), thetotal access count for each Px,y is 80. Since those accesses
are performed by reading data from the fast shared memory,a high EB is obtained. This indicates that kernel 3 uses theshared memory effectively.
On the other hand, kernel 2 and 6 only achieved about halfof GTX 480 peak memory bandwidth. According to (1)–(3),HV
kernel 2 refers 19 cells for each Px,y to compute ĝx,y, ĝx,y, andDĝx,y , and kernel 6 which minimizes global memory transfer
for (16)–(18) also refers 19 cells for each Px,y to interpolate
the rest of channels. Theoretically, both kernels should achieve
memory bandwidth about a quarter of kernel 3, since their total
memory access count is a factor of 0.24 times of kernel 3. Bothkernels’ EB confirmed this hypothesis, where the numbers areabout a quarter of kernel 3’s EB .
Finally, while the CPU implementation does not seem to be
fast enough, it is actually very efficient by looking at its
bandwidth number. A memory bandwidth of 42.62 GB/s which
far exceeds peak DDR3 bandwidth of 10.6 GB/s by about
400% indicates that 1) the CPU-based implementation uses
CPU cache pretty well and 2) performing Chung’s algorithm
on a 4K frame is very memory intensive.

Table IV presents the further breakdown of kernel times for the
CUDA-based multi kernels implementation. In this table, kernel 5 contains the isolated code of green channel interpolation
step. Based on results shown in Table IV, by isolating green
channel interpolation in one kernel and eliminating unused
input data in that kernel, we are able to reduce the total kernel
time of single kernel implementation from 370.5 ms to 178.9
ms, which represents a 2.1x speedup. At the same time, we
are also able to increase its memory bandwidth from 7.3 GB/s
to 15.2 GB/s which represents a 2.1x increase.

We have presented a CUDA-based acceleration of Chung’s
demosaicing algorithm. The acceleration strategy we proposedconsists of two methods, which are 1) implementing thealgorithm as multiple kernels to separate the bottleneck portionof the algorithm from the rest and to minimize idle threads,
and 2) reducing I/O between shared and global memory when
performing the green channel interpolation step by separating
the input RAW data.

Meanwhile, since the wavefront processing on kernel 5 performs at O(w + h) steps, it only yields EB of 8.5% evenafter we applied our optimizations. This kernel actually takes
about 85.2% of the total kernel execution time in our multiple
kernels version. Therefore, it is still able to benefit from any
further memory-based optimizations.

Our experimental results show that methods 1) and 2) succeeded in accelerating the kernel time by a factor of 2.1x
faster than its single kernel counterpart. Additionally, we
also discover that the green channel interpolation step whichbecomes the bottleneck of the implementation spends 85.2%of the total kernel time.

We are also able to let other computation phases to perform

As future work, we plan to extend the implementation to

VI. C ONCLUSION AND F URTHER W ORK

support multiple GPUs usage. We also plan to further increase
the performance of green channel interpolation kernel by
reducing idle threads during wavefront computation.
ACKNOWLEDGMENT
This work was partly supported by JSPS Grant-in-Aid forYoung Researchers (B)(23700057) and Scientific Research(B)(23300007). The authors would like to thank the anonymous reviewers for their helpful comments to improve thequality of the paper.
R EFERENCES[1] B.K. Gunturk, J. Glotzbach, Y. Altunbask, R.W. Schafer, and R.M.
Mersereau, “Demosaicing: color filter array interpolation,” IEEE Signal
Processing Magazine, vol. 22, no. 1, pp. 44–54, (2005).
[2] B.E. Bayer, “Color imaging array,” U.S. Patent 3 971 065, July 1976.
[3] K.H. Chung, and Y.H. Chan, “Color Demosaicing Using Variance ofColor Differences,” IEEE Transactions on Image Processing, vol. 15,
Issue 10, pp. 2944–2955 (2006).
[4] NVIDIA Corporation, “CUDA PROGRAMMING GUIDE 4.0” (2011).
[5] M. McGuire, “Efficient, High-Quality Bayer Demosaic Filtering on
GPUs,” Journal of Graphics, GPU, and Game Tools, vol. 13, issue 4,
pp. 1–16 (2008).
[6] H.S. Malvar, L.W. He, and R. Cutler, “High quality linear interpolation
for demosaicing of Bayer-patterned color images,” Microsoft Research(2004).
[7] J. Stam and J. Fung, “Image De-Mosaicing,” GPU COMPUTING GEMSEMERALD EDITION, Morgan Kaufmann, pp. 583–398 (2011).
[8] C.E. Duchon, “Lanczos filtering in one and two dimensions,” Journal ofApplied Meteorology and Climatology, vol. 18, pp. 1016–1022 (1979).
[9] J.F. Hamilton and J.E. Adams, “Adaptive color plane interpolation insingle sensor color electronic camera”, U.S. Patent 5 629 734 (1997).
[10] IRIDAS, http://www.iridas.com (2009).
Corporation:“CUDA
VisualProfiler”,[11] NVIDIA
http://developer.nvidia.com/nvidia-visual-profiler (2011).
[12] J. Mairal, M. Elad, and G. Sapiro, “Sparse Representation for Color
Image Restoration,” IEEE Transactions on Image Processing, vol. 17,
issue 1, pp. 53–69 (2011).
[13] D. Menon and G. Calvagno G, “Color Image Demosaicking: An
overview,” Signal Processing Image Communication (2011).
[14] K. McLaren, “The development of the CIE 1976 (L*a*b*) uniform
color-space and colour-difference formula,” Journal of the Society ofDyers and Colourists, vol. 92, pp. 338–341 (1976).
[15] M.
Snir,“Resourceson
ParallelPatterns”,
http://www.cs.uiuc.edu/homes/snir/PPP/ (2011).
Corporation:“CUDA
Occupancy
Calculator”,[16] NVIDIA
http://developer.download.nvidia.com/compute/DevZone/docs/html/
C/tools/CUDA Occupancy Calculator.xls (2011).