日本ソフトウェア科学会第 32 回大会 (2015 年度) 講演論文集

MRI Image Processing with OpenCL
Juan Luis Soler Ferrer

Kiminori Matsuzaki

Functional magnetic resonance imaging (fMRI) is now widely used in brain science, and it takes brain
images in four dimensions, usual three dimension images over time. The goal of this project is to obtainhigher spatio-resolution images from the original ones, making use of super-resolution techniques in two
dimensional, three dimensional and four dimensional scopes. In this study, we employ the open standard
for parallel programming of heterogeneous systems OpenCL. With this approach of accelerated computing
the code can be executed on diverse CPUs and GPUs, this is usually a great scenario for image processing, benefiting of the multiple cores and the balancing of the workload to get a higher performance of the
hardware. We report the experiment results comparing the sequential implementation with the parallel
implementation on a PC with two GPU devices, as well as lessons obtained through this project.
近年，脳科学において機能的磁気共鳴画像化法（fMRI）が広く利用されている．fMRI で撮像される脳画像は 4 次
元のデータである，すなわち，通常の 3 次元画像の時系列である．本プロジェクトの目標は，超解像手法を 2 次元，
3 次元および 4 次元データに適用することで，通常の fMRI 画像から高空間解像度の画像を得ることである．本研
究では，ヘテロジニアス環境における並列プログラミングのための標準のひとつである OpenCL を適用する．このアクセラレータを用いた計算手法により，プログラムコードをさまざまな CPU や GPU 上で動作させることができる．これは，画像処理に適しており，ハードウェアの複数のコアを利用することで高いパフォーマンスを得ることが期待できる．2 つの GPU デバイスを持つ計算機上で逐次プログラムと並列プログラムを用いた実験の結果ならびに，本プロジェクトにおいて得た知見について報告する．

1 Introduction
Functional magnetic resonance imaging (fMRI)
is now widely used in brain science. In fMRI, we
often use a imaging method that takes an image ofthe (whole) brain in short time, 2–3 seconds. However, compared to the high temporal-resolution,
the spatio-resolution of fMRI images are low: typical fMRI images have voxels of 3mm each. We
are studying to cope with this rather low spatioresolution of fMRI images [9, 10], especially basedon the so-called super-resolution technique [12].
In this paper, we improved the performance of∗ OpenCL を用いた MRI 画像処理

This is an unrefereed paper. Copyrights belong to
the authors.
Juan Luis Soler Ferrer, 高 知 工 科 大 学, Universitat
Politècnica de València, Kochi University of Technology, Universitat Politècnica de València.
Kiminori Matsuzaki, 高知工科大学, Kochi Universityof Technology.

the super-resolution program by using the generalpurpose GPU computing (GPGPU) technique. Especially, we developed a program for heterogeneous
GPU environment with OpenCL [6]. We will report
how we can develop a program for super-resolving
fMRI images with OpenCL and the experiment results.
The rest of the paper is organized as follows. In
Section 2, we introduce some important technologies used in this study. In Section 3, we show how
we developed a parallel program with OpenCL. In
Section 4, we report the experiment results conducted on a laptop PC with two GPU devices. In
Section 5, we discuss how we learned in this projectand our future directions.

2 Preliminaries
2. 1 Functional Magnetic Resonance Imaging
Recently functional magnetic resonance imaging(fMRI) is widely used in brain science. Using

fMRI for brain science has two important advantages: fMRI is non-invasive and we can take images of interior portion. In clinical use of MRI,
we often take structural images with high spatioresolution of 0.5–1 mm, which take several minutes
for the whole brain image. In fMRI, we take functional images with the time-resolution of 2–3 seconds but the image has low spatio-resolution of 3mm. Since the diﬀerence of signals is so small (often
less than several percents) that we need statisticalpost-processing to get an activation map and judge
the brain activities.
Images taken by fMRI have rather low spaceresolution. Increasing the space-resolution is helpful for the growth of brain science. One approach is to apply post-processing to obtain images
with higher space-resolution, which is the superresolution technique [12]. There have been a lot of
studies for the super-resolution technique applied
to the structural images (MRI) [4, 11]. However,there have been only a few studies for the functional images (fMRI) [1, 5, 8].
2. 2 Super-resolution Images
It is known as super-resolution (SR) [12] the setof techniques and algorithms designed to increase
the spatial resolution of an image, typically froma sequence of lower resolution images. They diﬀer
from traditional techniques of image scaling thatthe latter ones only use an image for the resolution
increase, focusing its objective to maintain sharpedges, without the appearance of new details. In
contrast, in super-resolution the case is to merge information from several images taken from the same
scene, to represent details that initially are not signiﬁcant in the original images.
In this project we use the reconstruction approach that it takes a set of low-resolution images with sub-pixel shift and estimates the highresolution image so that the error between the input
images and simulated images is minimized utilizing
the MAP algorithm.
Maximum a posteriori (MAP) Estimation [2] isa widely-used algorithm for reconstruction-based
super-resolution. In this MAP estimation, we reconstruct the high-resolution image by maximizing
the posterior probabilities based on the Bayesianestimation.
Here, we denote the high-resolution image x, thenumber of observed images n, the observed images
yi , and the transformation matrices Ai = Di Bi Mi .
A typical model in MAP estimation assumes thatthe noises follow the normal distribution and thatthe image is somehow smooth. Under this assumption, it minimizes the following objective function.

n||yi − Ai x|| + α||Hx||
E(x) =
i=1
In this equation, α is a parameter that control
the smoothness of the image, H is a high-pass ﬁlter. The ﬁrst term decreases if the estimated image
is close to observed images. The second term decreases if the estimated image is smooth.
2. 3 OpenCL
Open Computing Language (OpenCL) [6] isa framework for writing programs that executeacross heterogeneous platforms consisting of CPUs,
GPUs, DSPs, FPGAs and other processors.
OpenCL applications consist of two parts. The
OpenCL host program is a pure software routine
written in standard C/C++ that runs on any sortof microprocessor. At a certain point during the execution of this host software routine, there is likely
to be a function that is computationally expensive
and can beneﬁt from the highly parallel acceleration on a more parallel device: a CPU, GPU,FPGA, etc. This function to be accelerated is referred to as an OpenCL kernel. These kernels are
written in standard C; however, they are annotatedwith constructs to specify parallelism and memory
hierarchy. A Kernel is executed in the computedevices that typically consists of several computeunits, which in turn comprise multiple processing
elements (PEs). A single kernel execution can run
on all or many of the PEs in parallel.
SR needs to apply an algorithm per pixel in multiple iterations, this scenario is perfect for parallelizing using the OpenCL technology.
2. 4 OpenGL and SFML
Open Graphics Library (OpenGL) [7] is a crosslanguage, multi-platform application programming
interface (API) for rendering 2D and 3D vector
graphics. The API is typically used to interact
with a graphics processing unit (GPU), to achieve
hardware-accelerated rendering. For this project it
is used because two reasons. The ﬁrst one because
can work directly with OpenCL, for example they

can share texture or image buﬀers without using
any intermediate command. And the second one,for better user interface, being capable of checking
the ﬁnal results with an image, or an animation ora 3d model voxel rendering improves the debugging
process.
Simple and Fast Multimedia Library (SFML) [3]
is a cross-platform software development library
designed to provide a simple interface to various
multimedia components in computers. Also handles window creation and input as well as the creation and management of OpenGL contexts. In
this project it is used as auxiliary library for dealing with particularities of the diﬀerent OS like thetime measurement or the OpenGL handling thatcan diﬀer totally depending in the system that is
running.

3 Implementation
We took the fMRI images under the following
conditions. The subject was a man in his thirties.
The task consisted of 15 sets of 30 seconds restingstate and 30 seconds right-hand tapping motion, in
turn. The whole brain was scanned every three seconds, and each scan took 36 slice images (64 × 64
resolution, 12-bit depth).
In total we have 310 3D images of 64 × 64 × 36,
making it 11160 slices that will be processed using
the techniques previously explained for the purpose of improving the resolution from 64 × 64 to
128 × 128.
We have the algorithm implemented two times,
one in C++ in sequential code and the other in C99
for OpenCL. Both implementations are practically
identical with a few exceptions in their respective
idiosyncrasies. For the sequential code implementation it is straight forward, and does not need any
kind of preparation.
For the OpenCL implementation needs some
preparation prior to execution. Before creating
the context for OpenCL it is necessary to know inwhat platform (usually vendor dependent) to create the context. First, it is mandatory to checkthe OpenCL compatible platforms in the host computer. When it is decided the platform to use, thenwe can create the working context. A platform consists of one or multiple devices, where our code will
be executed, such as CPU, GPU, FPGA, DSP, etc.

We compile the program that is going to be executed in the device. It is important to stress that
this part is done on runtime.
Fig. 1 is the function that we call in the host program for the SR algorithm. The input is a vector of
4 low-resolution images and the output is the image
processed.
Line 1: The class i2D is a container of 2d image
in the format used by fMRI.
Lines 2–4: Create 3 buﬀers, 2 read-only for input
data and 1 write-only for output data.
Lines 6–10: Create the kernel that manage thecalls to the function ”superResolution” (Fig. 2) inthe devices, also put the buﬀers previously createdas arguments for the kernel.
Lines 12–13: Create the command queue that it
will our medium to sending commands to the device. All the commands processed by the device
they may past through the command queue.
Lines 15–18: Create two new images that will be
used for intermediate processing by the SR algorithm, and other auxiliary variables for the kernel
execution.
Lines 20–23: We write the vector with the 4 lowresolution images to one of the input buﬀers, andoptimize using only one concatenating the 4 images.
Line 25: Start the iteration loop.
Lines 26–27: Write the image in the buﬀer that
we have worked on.
Lines 28–29: Execute the kernel program.
Lines 30–31: Copy the output buﬀer directly tothe input buﬀer for the iteration process.
Lines 34–35: Read the last output buﬀer as aresult image.
Lines 37–41: Release OpenCL resources from thedevice.
Line 43: Return the processed image and ﬁnish
the function.

4 Results
The results obtained were gathered from a MSI
GE60 2PE Apache Pro Laptop Computer, thiscomputer has an Intel Core i7 (8 cores) and twographic cards, a Low consumption Intel graphic
card and a high performance NVidia GTX860M
graphic card. That implies that we have 3 OpenCL
compatible devices in one computer for testing.
Table 1 shows the experiment results where wemeasured the execution time and calculated the rel-

Table 1

Experiment results

time (s)DeviceIntel CPU (Sequential)
Intel GPU (OpenCL)
NVidia GPU (OpenCL)

Time (s)
535.5
223.9
108.5

Speedup
—
2.39
4.93

600
500
400
300
200
100
0

References

Intel CPU Intel GPU NVidia GPU(Sequential) (OpenCL) (OpenCL)

Speedupof the drivers. Also some vendors like NVidia put
OpenCL implementation under the CUDA technology, a similar but proprietary technology.
This project can be expanded in directions: forexample we may get better performance using an
Altera or Xilinx FPGA; we should analyze the results expanding with more criteria like power consumption and power eﬃciency; and conclude which
technology is more appropriate for end-user usage.

Fig. 3

Process timeIntel CPU
(Sequential)

Intel GPU
(OpenCL)

6
5
4
3
2
1
0

Fig. 4

NVidia GPU(OpenCL)

Speedup

ative speedup with respect to the sequential code
executed on CPU. Figures 3 and 4 plot the execution time and speedup. From this results, we cansee that the OpenCL version runs 2–5 times fasterwith a GPU.

5 Discussion
The results of the project are very satisfactory,
the utilization of the OpenCL technology in thistype of case is a great scenario of parallel computing. Having to process each pixel/voxel individually without dependencies with each other makes
it a very optimizable algorithm in this kind of context.
A negative part is that, although OpenCL isa good technology, in the end it depends on thevendors or manufacturers of the devices and theirimplementations. The closed-source nature of thedrivers makes it hard to work with it in Linux, whilein the Windows counterpart there is a newer version[ 1 ] Ben-Eliezer, N., Goerke, U., Ugurbil, K.,
and Frydman, L.: Functional MRI using superresolved spatiotemporal encoding, Magnetic Resonance Imaging, Vol. 30, No. 10(2012), pp. 1401–
1408.
[ 2 ] Bishop, C.: Pattern Recognition and Machine
Learning, Springer, 2006.
[ 3 ] Gomila, L.: Simple and Fast Multimedia Library(SFML), http://www.sfml-dev.org/.
[ 4 ] Greenspan, H.:
Super-resolution in medical imaging, The Computer Journal, Vol. 52,
No. 1(2009), pp. 43–63.
[ 5 ] Joshi, S., Marquina, A., Osher, S., Dinov, I.,
Darrell, J., Horn, V., and Toga, A.: Image resolution enhancement and its applications to medical image processing, Technical Report TechReport
CAM08-62, Department of Mathematics, UCLA,
2008.
[ 6 ] Khronos group: OpenCL, https://www.khronos.
org/opencl/.
[ 7 ] Khronos group: OpenGL, https://www.khronos.
org/opengl/.
[ 8 ] Kornprobst, P., Peeters, R., Nikolova, M., Deriche, R., Ng, M., and Hecke, P. V.: A superresolution framework for fMRI sequences and its impact
on resulting activation maps, Medical Image Computing and Computer-Assisted Intervention (MICCAI03), 2003, pp. 117–125.
[ 9 ] Matsuzaki, K. and Miyazaki, R.: Evaluation ofSuper-Resolution for fMRI Images, Proceedings of
5th International Symposium of Frontier Technology (ISFT 2015), 2015.
[10] Matsuzaki, K. and Miyazaki, R.:
SuperResolution for fMRI Images and its Evaluation,
Kochi University of Technology Research Bulletin,Vol. 12, No. 1(2015), pp. 131–138.
[11] Reeth, E., Tham, I., Tan, C., and Poh, C.:
Super-resolution in magnetic resonance imaging: A
review, Concepts in Magnetic Resonance, Vol. 40A,
No. 6(2012), pp. 306–325.
[12] Simpkins, J. and Stevenson, R.: An introduction
to super-resolution imaging, Mathematical Optics:
Classical, Quantum, and Computational Methods,
CRC Press, 2012.

1
2
3
4

i2D superResolutionOpenCL(vector<i2D> imas) {cl_mem _in_imas = CL_CHECK_ERR(clCreateBuffer(OCL.context, CL_MEM_READ_ONLY, sizeof(usint)*64*64*4, NULL, &_err));
cl_mem _in_res = CL_CHECK_ERR(clCreateBuffer(OCL.context, CL_MEM_READ_ONLY, sizeof(usint)*128*128, NULL, &_err));
cl_mem _out_next = CL_CHECK_ERR(clCreateBuffer(OCL.context, CL_MEM_WRITE_ONLY, sizeof(usint)*128*128, NULL, &_err));5
6

cl_kernel kernel;
kernel = CL_CHECK_ERR(clCreateKernel(OCL.program, "superResolution", &_err));
CL_CHECK(clSetKernelArg(kernel, 0, sizeof(_in_imas), &_in_imas));
CL_CHECK(clSetKernelArg(kernel, 1, sizeof(_in_res), &_in_res));
CL_CHECK(clSetKernelArg(kernel, 2, sizeof(_out_next), &_out_next));

7
8
9
10
11
12

cl_command_queue queue;
queue = CL_CHECK_ERR(clCreateCommandQueue(OCL.context, OCL.devices[0], 0, &_err));

13
14
15

i2D res;res.setSize(128,128,2047);
i2D next;next.setSize(128,128);cl_event kernel_completion;size_t global_work_size[1] = { 128*128 };

16
17
18
19
20

CL_CHECK(clEnqueueWriteBuffer(queue, _in_imas, CL_TRUE, 0*sizeof(usint)*64*64, sizeof(usint)*64*64, &imas[0].datos[0], 0, NULL, NULL));
CL_CHECK(clEnqueueWriteBuffer(queue, _in_imas, CL_TRUE, 1*sizeof(usint)*64*64, sizeof(usint)*64*64, &imas[1].datos[0], 0, NULL, NULL));
CL_CHECK(clEnqueueWriteBuffer(queue, _in_imas, CL_TRUE, 2*sizeof(usint)*64*64, sizeof(usint)*64*64, &imas[2].datos[0], 0, NULL, NULL));
CL_CHECK(clEnqueueWriteBuffer(queue, _in_imas, CL_TRUE, 3*sizeof(usint)*64*64, sizeof(usint)*64*64, &imas[3].datos[0], 0, NULL, NULL));

21
22
23
24
25

for (int i=0;i<50;i++) {CL_CHECK(clEnqueueWriteBuffer(queue, _in_res, CL_TRUE, 0, sizeof(usint)*128*128, &res.datos[0], 0, NULL, &kernel_completion));
CL_CHECK(clWaitForEvents(1, &kernel_completion));CL_CHECK(clReleaseEvent(kernel_completion));
CL_CHECK(clEnqueueNDRangeKernel(queue, kernel, 1, NULL, global_work_size, NULL, 0, NULL, &kernel_completion));
CL_CHECK(clWaitForEvents(1, &kernel_completion));CL_CHECK(clReleaseEvent(kernel_completion));
CL_CHECK(clEnqueueCopyBuffer(queue,_out_next,_in_res,0,0,sizeof(usint)*128*128,NULL,NULL,&
kernel_completion));
CL_CHECK(clWaitForEvents(1, &kernel_completion));CL_CHECK(clReleaseEvent(kernel_completion));
}

26
27
28
29
30
31
32
33
34

CL_CHECK(clEnqueueReadBuffer(queue, _out_next, CL_TRUE, 0, sizeof(usint)*128*128, &next.datos[0], 0, NULL, &kernel_completion));
CL_CHECK(clWaitForEvents(1, &kernel_completion));CL_CHECK(clReleaseEvent(kernel_completion));

35
36
37

clReleaseCommandQueue(queue);clReleaseKernel(kernel);clReleaseMemObject(_in_imas);clReleaseMemObject(_in_res);clReleaseMemObject(_out_next);

38
39
40
41
42
43
44

return next;}

Fig. 1The C++ function that manages the OpenCL algorithm implementation

1
2
3
4
5
6

__kernel void superResolution(__global unsigned short int *_in_imas, __global unsigned short int*_in_res, __global unsigned short int *_out_next){
int i = get_global_id(0);int x=i%128;
int y=i/128;
int alpha = 4;78

int dd[4][2] = {{0,0}, {0,1}, {1,0}, {1,1}};9
10

int dEdp=0;

1112

for (int l=0;l<4;l++) {int dx = dd[l][0];int dy = dd[l][1];
int newx = (x + dx) / 2; if (newx < 0) continue; if (newx >= 64) continue;
int newy = (y + dy) / 2; if (newy < 0) continue; if (newy >= 64) continue;
int pl=0;

13
14
15
16
17
18
19

for ( int lx = newx * 2 - dx; lx < (newx + 1) * 2 - dx; lx++) {for (int ly = newy * 2 - dy; ly < (newy + 1) * 2 - dy; ly++) {if (lx < 0 || lx >= 128 || ly < 0 || ly >= 128) {//pl += 0;} else {pl += _in_res[lx+(ly*128)];}}}
pl /= 4;
int ql = _in_imas[newx+(newy*64)+(l*64*64)];dEdp += pl - ql;

20
21
22
23
24
25
26
27
28
29
30
31
32

if (x > 0 && y > 0 && x < 128 - 1 && y < 128 - 1) {int pn = _in_res[(x-1)+(y*128)] + _in_res[(x+1)+(y*128)] + _in_res[(x)+(y*128)-128] +_in_res[(x)+(y*128)+128];dEdp += alpha * (4 * _in_res[x+y*128] - pn) / 4;}
int pix;
if (dEdp > 0) {pix = _in_res[x+y*128] - (dEdp + 2) / 4;
} else {pix = _in_res[x+y*128] - (dEdp - 2) / 4;}
_out_next[x+y*128]=clamp(pix,0,2047);

33
34
35
36
37
38
39
40
41
42
43

}
if (x == 0 || y == 0 || x == 127 || y == 127) {_out_next[x+y*128]=0;
}

44
45
46
47

}

Fig. 2The OpenCL algorithm implementation