GPU DEVICES AND
PARALLELISM OF IMAGE PROCESSING
TIAN XUWEN
2013/09/05
此内容仅代表个人观点
AGENDA
Steps for parallel programming
Introduction of GPUs
Image Processing
1. Affine transformation
2. Convolution
1. STEPS FOR PARALLEL PROGRAMMING
Finding concurrency
Algorithm Structure
Supporting Structures
Implementation Mechanisms
Finding Concurrency
Dependency AnalysisDecomposition
1. STEPS FOR PARALLEL PROGRAMMING
Algorithm Structure
Supporting Structures
Implementation Mechanisms
Task Decomposition
Data Decomposition
Group Tasks
Order Tasks
Data Sharing
Design Evaluation
2. INTRODUCTION OF AMD GPUS
Processing Elements (PE)
Data Parallel (SIMD)
• Each PE computed independently
• No dependencies between “threads”
• High Arithmetic Intensity (compute to
bandwidth ratio)
• Hide memory latency
2. INTRODUCTION OF AMD GPUS
Simplified Diagram of PE
2. INTRODUCTION OF AMD GPUS
VLIW 5
TABLE OF CONTENTS
�Steps for parallel programming
� Introduction of GPUs
� Image Processing
‒Affine transformation
‒Convolution
�Questions and Answers
A REPRESENT FORM OF NATURAL IMAGES
sample grid space quantization intensity quantization digital image
2-D AFFINE TRANSFORMATION
3-D AFFINE TRANSFORMATION
3-D AFFINE TRANSFORMATION
�′
�′
�′
�′
=
��� �� ��
���
� � � �
� �
�
�
���
�
��
�
�
�
��
���
�′
�′
�′
�′
A =
5 0 0 0
0 1 0 0
0
0
0
0
1 0
0 1
B =
3/2 −1/2 0 0
1/2 3/2 0 0
0
0
0
0
1 0
0 1
C =
1 0 2 0
0 1 3 0
0
0
0
0
1 0
0 1
D =
1 0 0 2
0 1 0 1
0
0
0
0
1 2
0 1
H = ���� =
5 3/2 −1/2 2 2
5/2 3/2 3 1
0
0
0
0
1 2
0 1
GENERAL PROJECTIVE TRANSFORMATION
Spatial
Transformation
�a’b’c’d’: Block to be transformed
Previous frame
a’ b’
c’ d’
(u,v)
� abcd: Predicted block
Current frame
a
c d
b
(x,y)
GENERAL PROJECTIVE TRANSFORMATION
……
float4 m4x4[4], vi, vo;
= m4x4[0].x * + m4x4[0].y * + m4x4[0].z * + m4x4[0].w * ;
= m4x4[1].x * + m4x4[1].y * + m4x4[1].z * + m4x4[1].w * ;
= m4x4[2].x * + m4x4[2].y * + m4x4[2].z * + m4x4[2].w * ;
= m4x4[3].x * + m4x4[3].y * + m4x4[3].z * + m4x4[3].w * ;
……
The OpenCL code of Projective Transformation
Note：：：：
These lines can be further abbreviated by using fused multiple add (FMA) on AMD GPU/APU.
EXAMPLE 1
Image Transformation
EXAMPLE 2
Geometric Correct Remote Sensing
BILINEAR INTERPOLATION DIAGRAM
BICUBIC INTERPOLATION DIAGRAM
BICUBIC INTERPOLATION
If the values of a function f(x) and its derivative are known at x=0 and x=1, then the function can be
interpolated in the interval [0,1] by a third degree polynomial. This is called cubic interpolation. The third
degree polynomial and its derivative are:
� � =a�
+ b� + c� + d
�′ � =3a� + 2b� + c
The values of the polynomial and its derivative at x=0 and x=1:
� 0 = d
� 1 = a + b + c + d
�′ 0 = c
�′ 1 = 3a + 2b + c
� = 2� 0 - 2 � 1 + �′ 0 + � 1
� = −3� 0 + 3� 1 - 2�′ 0 - �′ 1
� = �′ 0
� = � 0
Often the derivative of the function is unknown. Smoother curves can be found by using the slope of a line
between the previous and the next point as the derivative at a point. This polynomial is called a Catmull-Rom
spline. Suppose the values p0, p1, p2 and p3 at respectively x=-1, x=0, x=1, and x=2, and assign the values of
f(0), f(1), f'(0), and f'(1), use the formulas below to interpolate between p1 and p2.
BICUBIC INTERPOLATION
� 0 = ��
� 1 = �
�′ 0 = � − �� 2⁄
�′ 1 = �
− �� 2⁄
� = −
�
�� +
�� −
� +
�
�
� = �� −
5
2
�� + 2� −
1
2
�
� = −
�
�� +
�
�
� = ��
� ��, ��, � , �
, � = −
�
�� +
�� −
� +
�
�
�
+ �� −
"
�� + 2� −
�
�
�
+ −
�
�� +
�
� �
+ ��
BICUBIC INTERPOLATION POLYNOMIAL
BICUBIC INTERPOLATION CODE
double a00, a01, a02, a03, a10, a11, a12, a13, a20, a21, a22, a23, a30, a31, a32, a33;
double value; // the value of the pixel to be interpolated
double x2 = x * x;
double x3 = x2 * x;
double y2 = y * y;
double y3 = y2 * y;
a00 = p[1][1];
a01 = *p[1][0] + *p[1][2];
a02 = p[1][0] - *p[1][1] + 2*p[1][2] - *p[1][3];
a03 = *p[1][0] + *p[1][1] - *p[1][2] + *p[1][3];
a10 = *p[0][1] + *p[2][1];
a11 = *p[0][0] - *p[0][2] - *p[2][0] + *p[2][2];
a12 = *p[0][0] + *p[0][1] - p[0][2] + *p[0][3] + *p[2][0] - *p[2][1] + p[2][2] - *p[2][3];
a13 = *p[0][0] - *p[0][1] + *p[0][2] - *p[0][3] - *p[2][0] + *p[2][1] - *p[2][2] +
*p[2][3];
a20 = p[0][1] - *p[1][1] + 2*p[2][1] - .5*p[3][1];
a21 = *p[0][0] + .5*p[0][2] + *p[1][0] - *p[1][2] - p[2][0] + p[2][2] + .25*p[3][0] - .25*p[3][2];
a22 = p[0][0] - *p[0][1] + 2*p[0][2] - .5*p[0][3] - *p[1][0] + *p[1][1] - 5*p[1][2] + *p[1][3] +
2*p[2][0] - 5*p[2][1] + 4*p[2][2] - p[2][3] - .5*p[3][0] + *p[3][1] - p[3][2] + .25*p[3][3];
a23 = *p[0][0] + *p[0][1] - *p[0][2] + .5*p[0][3] + *p[1][0] - *p[1][1] + *p[1][2] -
*p[1][3] - p[2][0] + 3*p[2][1] - 3*p[2][2] + p[2][3] + .25*p[3][0] - .75*p[3][1] + .75*p[3][2] - .25*p[3][3];
a30 = *p[0][1] + *p[1][1] - *p[2][1] + .5*p[3][1];
a31 = .25*p[0][0] - .25*p[0][2] - .75*p[1][0] + .75*p[1][2] + .75*p[2][0] - .75*p[2][2] - .25*p[3][0] + .25*p[3][2];
a32 = *p[0][0] + *p[0][1] - p[0][2] + .25*p[0][3] + *p[1][0] - *p[1][1] + 3*p[1][2] - .75*p[1][3] -
*p[2][0] + *p[2][1] - 3*p[2][2] + .75*p[2][3] + .5*p[3][0] - *p[3][1] + p[3][2] - .25*p[3][3];
a33 = .25*p[0][0] - .75*p[0][1] + .75*p[0][2] - .25*p[0][3] - .75*p[1][0] + *p[1][1] - *p[1][2] +
.75*p[1][3] + .75*p[2][0] - *p[2][1] + *p[2][2] - .75*p[2][3] - .25*p[3][0] + .75*p[3][1] - .75*p[3][2] +
.25*p[3][3];
value = (a00 + a01 * y + a02 * y2 + a03 * y3) + (a10 + a11 * y + a12 * y2 + a13 * y3) * x + (a20 + a21 * y + a22 *
y2 + a23 * y3) * x2 + (a30 + a31 * y + a32 * y2 + a33 * y3) * x3;
EXAMPLE 2 RECONSTRUCT MR IMAGES
Cartesian Scan Data Spiral Scan Data
Gridding
Least-Squares (LS)
Spiral scan data + LS
Superior images at expense of significantly more computation
kx
ky
FFT
IMAGE CONVOLUTION
IMAGE CONVOLUTION
IMPLEMENTATION ISSUES
SEQUENTIAL CODE OF IMAGE CONVOLUTION
kCenterX = kCols / 2;
kCenterY = kRows / 2; // find center position of kernel (half of kernel size)
for(i=0; i < rows; ++i) // rows
{
for(j=0; j < cols; ++j) // columns
{
for(m=0; m < kRows; ++m) // kernel rows
{
mm = kRows - 1 - m; // row index of flipped kernel
for(n=0; n < kCols; ++n) // kernel columns
{
nn = kCols - 1 - n; // column index of flipped kernel
// index of input signal, used for checking boundary
ii = i + (m - kCenterY);
jj = j + (n - kCenterX);
// ignore input samples which are out of bound
if( ii >= 0 && ii < rows && jj >= 0 && jj < cols )
out[i][j] += in[ii][jj] * kernel[mm][nn];
}
}
}
}
CONVOLUTION KERNEL
__kernel void Convolve_Float4If(
const __global float * pInput,
__constant float * pFilter,
__global float * pOutput,
const int nInWidth,
const int nFilterWidth
)
{
const int nWidth = get_global_size(0);
const int xOut = get_global_id(0);
const int yOut = get_global_id(1);
const int xInTopLeft = xOut;
const int yInTopLeft = yOut;
float4 sum4 = 0;
Note: from AMD Developer Central
opencl/
for (int r = 0; r < nFilterWidth; r++){
const int idxFtmp = r * nFilterWidth;
const int yIn = yInTopLeft + r;
const int idxIntmp = yIn * nInWidth + xInTopLeft;
int c = 0;
int c4 = 0;
while (c <= nFilterWidth-4){
float4 filter4 = vload4(c4, pFilter+idxFtmp);
float4 in4 = vload4(c4, pInput +idxIntmp);
sum4 += in4 * filter4;
c += 4;
c4++;
}
int cMod = nFilterWidth – c;
if (cMod == 1){
int idxF = idxFtmp + c;
int idxIn = idxIntmp + c;
+= pFilter[idxF]*pInput[idxIn];
}
else if (cMod == 2)
{ //Use float4 here to further optimize the kernel
int idxF = idxFtmp + c;
int idxIn = idxIntmp + c;
+= pFilter[idxF]*pInput[idxIn];
+= pFilter[idxF+1]*pInput[idxIn+1];
}
else if (cMod == 3)
{ //Use float4 here to further optimize the kernel
int idxF = idxFtmp + c;
int idxIn = idxIntmp + c;
+= pFilter[idxF]*pInput[idxIn];
+= pFilter[idxF+1]*pInput[idxIn+1];
+= pFilter[idxF+2]*pInput[idxIn+2];
}
} //for (int r = 0…
const int idxOut = yOut * nWidth + xOut;
pOutput[idxOut] = + + + ;
}
Q & A
