論文 / 著書情報
Article / Book Information
Title

A GPU-accelerated image reduction pipeline

AuthorMasafumi Niwano, Katsuhiro L. Murata, Ryo Adachi, Sili Wang, Yutaro
Tachibana, Yoichi Yatsu, Nobuyuki Kawai, Takashi Shimokawabe,
Ryosuke Itoh

Journal/Book namePublications of the Astronomical Society of Japan, Volume 73, Issue 1,
pp. 14-24

Issue date

2020, 10

Note

This is a pre-copyedited, author-produced version of an article acceptedfor publication in Publications of the Astronomical Society of Japanfollowing peer review. The version of record Publications of the
Astronomical Society of Japan, Volume 73, Issue 1, pp. 14-24 isavailable online at: https://doi.org/10.1093/pasj/psaa091.

Powered by T2R2 (Tokyo Institute Research Repository)Publ. Astron. Soc. Japan (2018) 00(0), 1–10
doi: 10.1093/pasj/xxx000

GPU-accelerated Image Reduction Pipeline
Masafumi N IWANO1 , Katsuhiro L. M URATA1 , Ryo A DACHI1 , Sili WANG1 ,Yutaro TACHIBANA1 , Yoichi YATSU1 , Nobuyuki K AWAI1 , Takashi
S HIMOKAWABE2 and Ryosuke I TOH1,3
1

Department of Physics, Tokyo Institute of Technology, 2-12-1 Ookayama, Meguro-ku, Tokyo
152-8551, Japan
2
Supercomputing Research Division, Information Technology Center, The University of Tokyo,
2-11-16 Yayoi, Bunkyo-ku, Tokyo 113-8658, Japan
3
Bisei Astronomical Observatory, 1723-70 Ookura, Bisei-cho, Ibara, Okayama 714-1411,
Japan∗ E-mail: niwano@hp.phys.titech.ac.jp

Received 2020 June 26; Accepted 2020 August 26

Abstract
We developed a high-speed image reduction pipeline using Graphics Processing Units (GPUs)
as hardware accelerators. Astronomers desire detecting EM counterpart of gravitational-wave
sources as soon as possible for sharing positional information to organize systematic followup observations. Therefore, high-speed image processing is important. We developed a newimage reduction pipeline for our robotic telescope system, which uses a GPU via a Python
package CuPy to achieve high-speed image processing. As a result, the processing speed
was increased by more than a factor of forty to that of the current pipeline, while maintaining
the same functions.
Key words: techniques: image processing, gravitational waves, gamma-ray burst: general

1 Introduction
Gamma-ray bursts (GRB) and gravitational waves (GW)2−3
are short-lived (<sec) and no one knows when or∼ 10
where they occur, and we call such phenomena as “transients”. In order to understand the physical processes oftransients, it is essential to perform multi-wavelength andcontinuous observations for before, during, and for a reasonable length of time after the occurrence. We aim forobserving their afterglow and perform the prompt followup observations. However, positional information about
the GW source obtained from laser interferometers such as
LIGO (LIGO Scientific Collaboration et.al. 2015), Virgo(Acernese et.al. 2014) and KAGRA (Akutsu et.al. 2019)has too large of an uncertainty (≈ 101−3 deg2 ) to performsystematic follow-up observations. Thus it is necessary to
immediately identify the EM counterpart and send morec 2018. Astronomical Society of Japan.
⃝

accurate information to the ground telescope network assoon as possible. To realize that, all processes from receiving an alert to finding the counterpart must be performed
in a short period (<∼ 0.1days), and therefore, high-speed
image processing becomes the key technology.
High-speed image processing is especially important in
the era of time-domain astronomy which requires highdata-rate observations consisting of a large number of pixels taken at high cadence. Wide-field surveys such as PanSTARRS (Chambers et.al. 2016), Subaru-HSC (Miyazaki
et.al. 2018), ZTF (Bellm et.al. 2019), and Tomo-e Gozen(Sako et.al. 2018) are typical examples. One of the solutions to deal with the high data-rate is CPU-based parallel
computing. For example, ZTF obtains very high-resolution
CCD images of about 600 Mpix every 40 sec, and its data
rate exceeds 200 Mbit/s. To process those extremely fastgrowing data, they perform parallel computing in a com-1

2

Publications of the Astronomical Society of Japan, (2018), Vol. 00, No. 0

Table 1. Performance comparison of CPU and GPU
Name
FP641
[TFLOPS]CPUIntel Xeon E5-2640v4
0.4
GPU
NVIDIA Tesla P100 (PCIe)
4.7!"#$%&'

Power2
[W]
90
250

*++

"=0>?@ (@09,2,909,
"=%&)

(%&)

,-.
/-0-

AB-:C:;D0,-E0B?@

1

Theoretical maximum number of 64-bit floating-point operations
performed per second. The CPU value is calculated with IPC=16(instructions per clock) and Clock=2.4GHz, and the GPU value isnominal.
2
The CPU value is the thermal design power, and the GPU value
is the maximum power consumption.

0123

FG9,:E-@ E;0C?;00124+-,HC:;D0,-E0B?@

0125)<-0CIB9</B@J
0126

A-/C2BK9<C1-:HB@J
puting system with 1192 physical CPU cores. As a result,they have achieved the real-time detection of transientsand moving objects (Masci et.al. 2019).
A Graphics Processing Unit (GPU) is a hardwareaccelerator for graphics processing. A GPU has thousandsof computing units called “shaders”, which correspond to
CPU cores, and the GPU achieves extremely high performance only for SIMD1 processing (e.g. matrix operation). Although that is achieved by specializing in graphics processing, this performance is sometimes utilized for
not only graphics processing, but also other purposes, andsuch GPU utilization is called GPGPU (General Purpose
computing on GPU). In the astronomy, GPGPU is applied for physical simulation, machine learning, analysis ofradio interferometer data, to name of few (Tobias et.al.
2009, Tuccillo et.al. 2018, Tazzari et.al. 2018). To build
the same-scale computing system as ZTF, 120 Xeon E52640v42 processors are required and these consume about
10 kW of power. However, with Tesla P100, the same theoretical performance can be achieved with only 10 modules,
and the power consumption will be less than 2.5 kW (table1). That advantage is limited to the case when the performance of the GPU can be fully exploited, but it is not
so diﬃcult because image processing is a speciality of theGPU.
The purpose of this research is to utilize GPU for astronomical image processing, and more specifically, to speed
up the image reduction pipeline of our robotic telescopes,
MITSuME-Akeno and MITSuME-Okayama (Kotani et.al.
2005, Yatsu et.al. 2007, Shimokawabe et.al. 2008).

2 Development
2.1 Image reduction at MITSuME
Both of MITSuME telescopes are equipped with three
CCD cameras to obtain simultaneously three-band images.
Each CCD generates a 16-bit integer image of 1074×1024&<BJ@19@00128!?L-//B@J

2

Single Instruction Multiple Data
CPU used at ZTF Science Data System,9:;<0Fig. 1. Current Pipeline overview: Rectangular objects in a frame labeled“HDD” (Hard Disk Drive) indicate FITS files, “raw data” is a set of unreduced data, “result” is a reduced data, and “tmp∗” are sets of a temporary
file. Arrows indicate data flows. Only the processing performed by IRAF isshown, and there are actually other processing.

pixels including an overscan area, which is stored as a FITS
file. Typical value of exposure time, overhead and data volume is 60 sec (for GW or GRB), 8 sec and 4 GB/night (for2 telescopes), respectively. To reduce the pixel randomization, we perform 9 points dithering. Therefore positions ofpoint sources on the detector diﬀers for each exposure evenobserving the same object. After derivation of pedestal
level and astrometry, FITS files are transferred to the dataserver in Tokyo Tech.
The current image reduction pipeline is implemented asa Python3 script which uses IRAF4 (Doug 1986, Doug
1993) via PyRAF5 , and its sequence is shown in figure 1.
Specifically, overscan cutout, bias subtraction, dark subtraction, flat fielding, bad pixel masking, alignment, andco-adding are performed by calling IRAF tasks. At thestage of bias subtraction, casting from 16-bit integer to
32-bit floating point is performed. In the overscan cutout,
overscans of 50 × 1024 pix and 2 pix from the edge areremoved, leaving the image of 1020 × 1020 pix. Becauseof the dithering, alignment is necessary before co-adding.
The pipeline uses bicubic-spline interpolation and sigmaclipped-mean algorithm for sub-pixel image shifting andco-adding, respectively. Only the processing performedby IRAF/PyRAF is shown here, and other processing
such as calculating the relative position of images using
WCSTools6 is not shown.
3

https://www.python.org/
http://ast.noao.edu/data/software
http://www.stsci.edu/institute/software hardware/pyraf
6
http://tdc-www.harvard.edu/wcstools/
4

1

0127

5

Publications of the Astronomical Society of Japan, (2018), Vol. 00, No. 0

3

2.2 New Reduction Pipeline

3.1 Consistency Evaluation

We created a new pipeline by using Python, because
Python is an easy-to-code language, and there are Python
packages available for the purpose of this research. The
new pipeline uses NumPy7 (Walt et.al. 2011), Astropy8(Astropy Collaboration et.al. 2018), PyWCS9 and CuPy10 .
CuPy is a Python package for numerical calculation on
GPU using CUDA11 . The most remarkable feature of
CuPy is the API compatibility with NumPy, thus a code
written with NumPy can be easily rewritten with CuPy.

We coded the bicubic-spline algorithm and sigma-clippedmean algorithm ourselves, because CuPy did not provide
corresponding functions. Therefore we had to confirm
whether the same result as the current pipeline can be
reproduced for these two processing methods. We defined∆Iˆ as a quantity for evaluating relative diﬀerence between
two images as follows,

The sequence of the new pipeline is described in figure 2. Specifically, the header and image data are readfrom FITS file as instances of astropy.io.fits.Header12and numpy.ndarray13 . At the same time, casting from 16bit integer to 32-bit floating point is also performed. The
header is passed to pywcs.WCS14 constructor for calculating
the relative position of the images, and the image data isconverted to cupy.ndarray15 while slicing overscan. Then,bias subtraction, dark subtraction, flat fielding, bad pixel
masking, alignment, and co-adding are performed on the
GPU with CuPy. The array of resulting image is converted to numpy.ndarray, and then written to a FITS filewith Astropy. The pipeline uses bicubic-spline interpolation and sigma-clipped-mean algorithm, the same as thecurrent pipeline. The image processing is performed on
the GPU-side because it is basically a SIMD operation,and the other processing is performed on the CPU-side because there are few benefits performing it on the GPU or
it cannot be performed on the GPU. The overscan cutout
is performed on the latter because the array-slicing operation is not SIMD, and transferring data not used for
calculation to the VRAM generates extra memory allocation time. astropy.wcs.WCS has the same functions aspywcs.WCS, but an initialization of the former takes ≈ 7ms
while the latter takes ≈ 200µs in our environment.

3 Verification and Performance Test
We evaluated consistency and performance of the new
pipeline. Information about the execution environment isshown in tables 2 and 3.

7

https://numpy.org/
https://www.astropy.org/index.html
9https://pypi.org/project/pywcs/
10
https://cupy.chainer.org/
11
https://developer.nvidia.com/cuda-zone
12
Python class provided by Astropy for storing contents of FITS header
13
Python class of multidimensional array provided by NumPy
14
Python class for converting coord between WCS and other system
15
Python class corresponding to numpy.ndarray in CuPy
8∆Iˆij ≡refIij − Iij
ref
Iij

(1)

ref
where Iij is the pixel value of the test image, Iijis thepixel value of the reference image, and i, j are indices of a
pixel in the image.

3.1.1 Shift Consistency
For the image processed up to bad-pixel masking in the
reduction, we created the following two sets of images with
ˆvarious shift amounts, and calculated ∆I.
A The image shifted by the new pipeline, using a shifted
image made by IRAF-imshift as a reference.
B The image that was shifted once and then shifted backto the original position by the new pipeline, using
the original image as a reference.
A is for evaluating the consistency with IRAF, and B isused to evaluating how much information of the original
image is conserved. We defined a coordinate system inwhich the lower left corner of the image is the origin, thehorizontal direction is the X-axis, and the vertical is the
Y-axis. Then, the shift amounts in the X, Y dimensions
are ∆X, ∆Y. The data used was obtained by applying
overscan cut out, bias & dark subtraction, flat fielding,and bad-pixel masking to four of the images acquired in
MITSuME-Akeno on May 14, 2018. These are shown by
the table 4 and figure 3.
Typical results are shown in figure 4 and figure 5 for A,
and figure 6 and figure 7 for B. Figure 4 and figure 6 areimages in which each pixel is appropriately colored withˆ and figure 5 and figure 7 are histograms of ∆I.
ˆ The
∆I,
−6
−5
diﬀerence from IRAF is about 10 − 10 . On the other
hand, the strange pattern seen in figure 4 suggests that
this diﬀerence is not random but systematic. Specifically,ˆ is small on the left sidewith the boundary X ≈ 240, |∆I|
and large on the right side. Around the locations where
the original image has point sources, ∆Iˆ is negative at nearside from the boundary and positive at far side. The diﬀerence from the original image is about 10−2 − 10−1 , and noparticular pattern as seen in figure 4 appears. The diﬀerence from IRAF is 3 to 4 orders of magnitude smaller thanthat from the original image. Therefore, the eﬀect of the

4

Publications of the Astronomical Society of Japan, (2018), Vol. 00, No. 0344!"#$%&'

1"#$2%&'"()*+, -,)./0/.)./&8)/+0(I9:"(=>?@A< *.57./A$,9:0(;,75//5(

/56
75)5

!9"(

KJ./8H5,
H9)A+9),9:0(;,75//5(<H90(;,75//5(
L-58A$A45/MA89N)/5H)-+,$A=B5)AO-.B7-,P"(C!@
DEFG%&E4.H!5BH /.B5)-J.A
0+8-)-+,8&8)/+0(/.89B)!#4&

L57A0-Q.BA:58M-,P&B-P,:.,)!+R577-,P,9:0(;,75//5(<=>?@A

H90(;,75//5(<,9:0(;,75//5(

Fig. 2. New pipeline overview: Rectangular objects on “HDD” indicate FITS files, “raw data” is a set of unreduced data, “result” is a reduced data. Arrows
indicate data flows. Roughly, “raw data” are read by Astropy, array calculation is performed by CuPy, and the result is written to “result” by Astropy.

Table 2. Machine Specification

M/B
CPU
RAM
GPU
SSD
HDD

Model
Supermicro X10SRA
Intel Xeon E5-1650 v4
A2ZEON RD4R16G44S2400 ×4
NVIDIA TITAN X (Pascal)Intel SSDSC2BB12
Seagate ST2000VN0001

Notes
C612
3.6-4.0GHz
DDR4 2400MHz Registered ECC, Quad Channel
3584 shaders, 1417-1531MHz, GDDR5X 12GB, PCIe 3.0 x16
SATA 6Gbps, OS and executable files are located
SATA 6Gbps, Script files and data are located

Table 3. Software versions
OS
Ubuntu 14.04 LTS x86 64
Python
2.7.6
CUDA
10.0
IRAF
2.16.1 NumPy
1.15.4
PyRAF
2.1.14 Astropy
2.0.9
SExtractor
2.8.6
PyWCS 1.12
WCSTools
3.8.7
CuPy
5.1.0

Image1

Image2

Image3

Image4∗

Python packages other than PyRAF were
installed with pip.

Table 4. Information of data used for verification of shiftconsistency
Label1
BG-LEVEL2 [ADU]
Image1
705
Image2
197
Image3
68
Image4699

BG-RMS3 [ADU]
23.74114.07213.51422.2671

An appropriate name to identify the data
Median of pixel values
3Standard deviation of pixel values∗
For all images, the shape is 1020 × 1020 and thenumerical format is a 32-bit float.
2

Fig. 3. Data used for verification of shift consistency

Publications of the Astronomical Society of Japan, (2018), Vol. 00, No. 0

5

Original : Image2, ∆X=0.3, ∆Y=0.2

Original : Image4, ∆X=0.6, ∆Y=0.7

1000

8

6

6
800

4
2

400

∆Î 10 −6

400

2
4

200

]
0

[

0

2

600

]

600

4

2
4

200

6

6
0

0

200

400

600

800

0

1000

Fig. 4. An image in which each pixel is appropriately colored with ∆Iˆ calculated for set A. Image2 is used, and (∆X, ∆Y) = (0.3, 0.2). With theˆ is small on the left side and large on the rightboundary X ≈ 240, |∆I|
side. Around the locations where the original image has point sources, ∆Iˆ
is negative at near side from the boundary and positive at far side.

0.5

[

800

∆Î 10 −2

1000

0

200

400

600

800

8

1000

Fig. 6. An image in which each pixel is appropriately colored with ∆Iˆ calculated for set B . Image4 is used, and (∆X, ∆Y) = (0.6, 0.7). No particularpattern appears other than the pattern around the vertical line in the original
image.

0.30

Original : Image2, ∆X=0.3, ∆Y=0.2

90% range: ±2.644e-06
95% range: ±3.507e-06

Original : Image4, ∆X=0.6, ∆Y=0.7

90% range: ±2.412e-02
95% range: ±2.894e-02

0.25
frequency (normalized)

frequency (normalized)

0.4
0.3
0.2

0.15
0.10
0.05

0.1
0.0

0.20

0.00
6

4

2

0
∆Î [10 −6 ]

2

4

6

Fig. 5. A histogram of ∆Iˆ calculated for set A. Image2 is used, and(∆X, ∆Y) = (0.3, 0.2). “90(95)% range” is the range in which 90(95)%of ∆Iˆ of all pixels fall. ∆Iˆ is about 10−6 − 10−5 .

shift operation itself on the result is more dominant thanthe eﬀect of the change from IRAF to the new pipeline.
We conducted additional experiments to determine theorigin of the 10−6 − 10−5 diﬀerence. First, we created animage in which two-dimensional cubic B-spline basis functions were arranged (figure 8 and figure 9). Since the cubic
B-spline basis function is exactly reproducible by cubicspline interpolation, a shifted image can be obtained without actually interpolating. Using the shifted image createdby this way as a reference image, we obtained ∆Iˆ for theimage shifted by IRAF-imshift and the image shifted by
new pipeline. The results are shown in figure 10 and 11.
For IRAF-imshift, ∆Iˆ ≈ 10−6 , and a line of discontinuity as shown in figure 4 appears. On the other hand, for

4

2

0
[ ]

∆Î 10 −2

2

4

Fig. 7. A histogram of ∆Iˆ calculated for set B. Image4 is used, and(∆X, ∆Y) = (0.6, 0.7). “90(95)% range” is the range in which 90(95)%of ∆Iˆ of all pixels fall. ∆Iˆ is about 10−2 − 10−1 .

the new pipeline, ∆Iˆ ≈ 10−7 , and there is no discontinuity.
Therefore, IRAF-imshift has a diﬀerence of about 10−6
from the accurate bicubic-spline interpolation. (It has notbeen tested in other environments, and it is possible thatthe bug occurs only in our environment.)
3.1.2 Co-adding Consistency
For the data set processed up to the alignment, we generated a co-added image by a function of the new pipeline,
and calculated ∆Iˆ using a image co-added by IRAFimcombine as the reference. The input images are acquired
at MITSuME-Akeno on May 14, 2018. Information about
the input-image sets are given in table 5. The co-added
images are shown in figure 12.

6Publications of the Astronomical Society of Japan, (2018), Vol. 00, No. 0

0.7

1.05
0.6
80
0.90
0.5

Y

0.60

0.3

∆Î [10 −7 ]

0.75

60
0.4

40
0.45

0.2

0.30
20

0.1

0.15

0

0.0
0.0

0.5

1.0

1.5

2.0

2.5

3.0

3.5

0.00
0

4.0

20

40

60

80

X

Fig. 8. Cubic B-spline basis function. This function satisfies the cubic-spline
constraint conditions, so it can be reproduced exactly with cubic-spline.

Fig. 11. An image in which each pixel is appropriately colored with ∆Iˆ calculated for new pipeline shifted image. (∆X,∆Y) = (0.3,0.2). ∆Iˆ is about
10−7 and there is no discontinuity.

Table 5. Information of data-sets used for verification ofco-adding consistency
Label1
Number of Images
Set1
35
Set2
81
Set3
97
Set4
44

Shape2 (X×Y) [pix]910 × 944
928 × 864
911 × 938
811 × 916

1

An appropriate name to identify the data-set
The reason why the shape of the image is diﬀerent fromthat before the alignment is that an area which does notoverlap with another image is cut out during the alignment.
∗
For all data-set, the numerical format is a 32-bit float.
2

Fig. 9. Two-dimensional cubic B-spline basis functions were arranged in thisimage. A shifted image of this can be obtained without actually interpolating.

Set1

Set2

Set3

Set4

1.25

1.00

80

0.75

0.25

40

∆Î [10 −6 ]

0.50

60

0.00

0.25

20

0.50

0.75

0
0

20

40

60

80

Fig. 10. An image in which each pixel is appropriately colored with ∆Iˆ calculated for IRAF shifted image. (∆X, ∆Y) = (0.3, 0.2). ∆Iˆ is about 10−6
and the line of discontinuity as seen in figure 4 appears.

Fig. 12. Co-added images of data sets used for verification of co-adding consistency

Publications of the Astronomical Society of Japan, (2018), Vol. 00, No. 0

7

Data-Set : Set3
1.0 1e790% range: ±2.464e-07
95% range: ±3.556e-07

Data set : Set3
8

800

0.8

6

frequency (normalized)

4
600

0
400

∆Î [10 −7 ]

2

2

0.6
0.4

4

0.2

200
6

8

0.0

0
0

200

400

600

800

ˆ DataFig. 13. An image in which each pixel is appropriately colored with ∆I.
set is Set3. No particular pattern appears.

Typical results are shown in figure 13 and figure 14.
Figure 13 is a image in which each pixel is appropriately
ˆ and figure 14 is a histogram of ∆I.
ˆ Incolored with ∆I,
figure 13, no particular pattern is seen. As can be seenfrom figure 14, ∆Iˆ is discrete, and is about 10−7 . 32-bit
floating point conforming to IEEE754 is assigned 1 bit fora sign, 8 bits for an exponent, and 23 bits for a significand.
In the case of a non-zero number, the maximum digit of
the significand is always 1. So, the significand is regarded
as having 24 bits by assuming that there is always 1 in the
24th digit (IEEE 2008). Therefore, the maximum numberof decimal digits that can be represented by 32-bit floating point is log10 224 ≈ 7. Since the last digit can easily
fluctuate depending on the order of calculations, the maximum number of significant digits is about 6, and this iscomparable to our result. Therefore, it can be said thatthe co-adding function of the new pipeline reproduces thesame function as IRAF-imcombine within the precision of
32-bit floating point.

3.2 Execution Speed
The procedure for measuring the execution time is described below. Using a set of images that have not been
reduced as input, we ran the pipeline 5 times while displaying the execution time with the Unix command time.
We considered an average of the five output values as theexecution time in this condition, and the standard deviation as the error. We define the number of inputs for thethree bands as n and the execution time as t.
3.2.1 Comparison with the current pipeline
We measured t while changing n by 30 (10 per band) forthe current pipeline, the new one, and the new one run-

6

4

2

0

∆Î

2

4

6
1e 7

ˆ Data set is Set3. “90(95)% range” is the range
Fig. 14. A histogram of ∆I.
in which 90(95)% of Diff of all pixels fall. The value of Diff is discrete, and isabout 10−7 .

ning entirely on the CPU without using the GPU. The
third is a code in which processing performed by CuPy isrewritten with NumPy to separate the eﬀect of performing
processing on the GPU from nonusage of IRAF/PyRAF.
Thereafter, we refer to the current pipeline as “Current”,
the new pipeline as “New(GPU)”, and the CPU version of
the new pipeline as “New(CPU)”. In order to evaluate thediﬀerence of the calculation time between the CPU andGPU as purely as possible, we also measured the breakdown of t for New(CPU) and New(GPU). This was performed by acquiring the UNIX time in each processing stepusing the Python module time16 .
We used a set of images acquired with MITSuMEAkeno at August 22, 2018. 188 images per band were
captured. In order to make a nice round number, the first
input was 180 images per band (i.e. n = 540).
Figure 15 is a plot of n and t for each pipeline. It showsthat t increases linearly with respect to n for n ≤ 540. Thisis because both the file I/O amount and the calculation
amount are approximately proportional to the data volume. Figure 16 is a plot of n and t-ratio for each pipeline.
The ratio diﬀers depending on n because t has a constantcomponent for n (e.g. starting up the interpreter). Whenn = 540, the ratio of t between Current and New(GPU)is 42. That is, New(GPU) runs 42 times faster. Table 6
shows the breakdown of t for New(CPU) and New(GPU).
Since FITS I/O and WCS operation do not use the GPU
for calculation even in New(GPU), the results are comparable between New(GPU) and New(CPU). On the other
hand, the other processing clearly shows diﬀerences among
the calculation speed of GPU and CPU, New(GPU) is
about 10-1000 times faster than New(CPU). In particu16

Python build-in module that provides time-related functions

8

Publications of the Astronomical Society of Japan, (2018), Vol. 00, No. 0%!

! ()*+,-./0123/4+5 678

$"

$!

#"

#!

"

!
!

%

&

'

#$

#"

!!! (9-4:+;20<2=/*+>75 6#!& =/*8

Fig. 15. Plot of n and t for each pipeline. Errors are not drawn because they
are comparable to the marker size. It is hard to see the plot of New(GPU)in the above figure, so the scale changed figure is shown below. For n ≤540, t of each pipeline increases almost linearly with n, and New(GPU) is
significantly faster than Current and New(CPU). New(CPU) is slower thanCurrent because our implementation of the bicubic-spline interpolation forimage alignment is less efficient than IRAF-imshift.

)!

!(%)*",-.
!(%)*+,-.

(!

'!

!*+,-./&!

%!

!"#$$%&'$!

!(%)*+,-.

#!

"!

!(%)*",-.
!"#$$%&'

!
!

"!!

#!!

$!!

%!!

&!!

'!!

" 012345+6/76.892-6.3,:5;<Fig. 16. Plot of n and the ratio of t. When n = 540, the ratio of t between
Current and New (GPU) is 42. That is, New(GPU) runs 42-times faster.

Table 6. Breakdown of t(n = 540)
New(GPU) [s]FITS I/O2
3.4
WCS operation
0.96
Bias,Dark,Flat
1.0 × 10−3
Bad-pixel mask
0.36
Alignment
0.48
Combining
0.14
1

New(CPU) [s]
2.7
0.94
0.74
10.55
4.0 × 102
12

Ratio1
0.81
0.97
720
29
830
91

New(CPU)/New(GPU)Including RAM-VRAM transfer time∗
Interpreter startup and module initialization time are not shown
2

Fig. 17. Plot of N × N and t. t has a downward convex change with respect
to N 2 . It seems that linear fitting can be performed up to N 2 ≤ 3 × 106 , but
the same fitting cannot be applied to all regions.

lar, since the calculational cost of bicubic-spline interpolation is very large, New(CPU) spends a very long timeon alignment. Therefore, if a less costly algorithm suchas nearest-neighbor or bilinear interpolation is used, t is
expected to be significantly reduced. In our environment,
NumPy uses only a single CPU core. However, depending
on how NumPy is built, multiple cores can be used, and
it is expected to run faster. We note that New(CPU) is
slower than Current by factor of 2. This is probably causedby the diﬀerent implementations of their bicubic-spline interpolations in the image alignment. In the bicubic-splineinterpolation, there is a process of solving simultaneous
linear equations for each column or row of the image. The
number of unknowns of each simultaneous equations isabout the same as N , the number of pixels in the eachcolumn or row (e.g. N = 1020 for MITSuME data), andthose simultaneous equations can be expressed as a tridiagonal matrix. The cost of the solution optimized for thetridiagonal matrix is O(N 2 ), and IRAF-imshift probably
uses such a solution. The new pipeline solves those equations using the inverse matrix which results in the cost ofO(N 3 ).
3.2.2 Application to various data
We also evaluated the execution time for the various image
sizes, assuming the application of this pipeline to observational data other than MITSuME. We generated a set ofdummy images in which random values were assigned to
N × N pixels, and measured the execution time t whilechanging N as long as the VRAM is suﬃcient. The number of images was fixed at 100 per band (i.e. n = 300). The
result is shown in figure 17.
Assuming that both the file-reading time and the cal-Publications of the Astronomical Society of Japan, (2018), Vol. 00, No. 0&!

! '()*+,-./012.3*4 567

%!

$!

#!

"!

!
!

&!!

"!!!

"&!!

#!!!

#&!!

$!!!

$&!!

%!!!

%&!!

" '8,39*:1/;1.0<,-1.3=>*64
Fig. 18. Plot of n and t. This is the same condition of New(GPU) in figure 15,
except that the input is dummy data. The dashed line is the model obtainedby linear-fitting the plot of New(GPU) in figure 15. At n > 1000, the plot
deviates significantly from the line and has a convex downward change.

culation time are basically proportional to the data volume, t is considered to change linearly with respect to N 2 .
However, in figure 17, t has a convex downward changewith respect to N 2 . The first conceivable cause is the calculation time of bicubic-spline interpolation. As mentioned
in section 3.2.1, our implementation of the bicubic-splineinterpolation has a calculational cost of O(N 3 ). Therefore,
this eﬀect does not depend on hardware. Figure 18 shows
the result of performing the same measurement while fixing N = 1020 and changing n as long as the VRAM is
suﬃcient. This is the same condition as the measurement
performed in section 3.2.1, except that the input data isdummy data. When n is large, t has a downward convex
change respect to n. The details have not been clarified
yet, but we believe that a non-linear eﬀect about the data
volume is aﬀecting the execution time. What we currently
suspect is the time required for the memory allocation andreleasing. However, we could not reproduce that eﬀect bysimply allocating or releasing the memory.

4 Summary
To speed up the image reduction processing of MITSuME
data, and to establish set of a low-cost, high-speed imageprocessing methods, we developed a GPU-accelerated
image-reduction pipeline. This pipeline uses CuPy, a
Python package for GPGPU, to perform image processing on the GPU. Because the library does not provide a
corresponding function, there is some processing which wascoded by ourselves. However, this pipeline reproduced thefunctions of the current pipeline with a relative diﬀerence−5of <∼ 10 . Furthermore, the new pipeline is 42 times faster9

than the current pipeline when processing 540 images.
Currently, we are performing real-time image reductionof MITSuME observational data with GPU-accelerated
pipeline, and using it to check results of prompt followup observation of GRBs. At present, the reduced images of this pipeline have been used in 6 reports to
GRB Coordinates Network Circular Service (Toma et.al.
2019; Adachi et.al. 2019; Niwano et.al. 2019; Adachi et.al.
2020; Oeda et.al. 2020; Hosokawa et.al. 2020).
We have developed a Python package based on our
new pipeline and released it on GitHub17 . Thus, GPUaccelerated image reduction can be performed with on datafrom other observatories.

Acknowledgments
This work was supported by JSPS KAKENHI Grant Numbers
JP16K13783, JP17H06362 and JP20K04011. This work wasalso supported by Optical and Near-Infrared Astronomy InterUniversity Cooperation Program of the MEXT of Japan, andby JSPS and NSF under the JSPS–NSF Partnerships forInternational Research and Education.

References
IEEE Standard for Floating-Point Arithmetic,” in IEEE Std
754-2008 , vol., no., pp.1-70, 29 Aug. 2008
The LIGO Scientific Collaboration, et.al. 2015, Class. Quantum
Grav.32 07400
F. Acernese, M. Agathos, K. Agatsuma, et.al. 2014,arXiv:1408.3978
T. Akutsu, M. Ando, K. Arai, et.al. 2019, Nature Astronomy,
3, pages35–40
K. C. Chambers, E. A. Magnier. N. Metcalfe, et.al. 2016,arXiv:1612.05560
E. A. Magnier, K. C. Chambers, H. A. Flewelling, et.al. 2016,arXiv:1612.05240
S. Miyazaki, Y. Komiyama, S. Kawanomoto, et.al. 2018, PASJ,
70, S1
H. Furusawa, M. Koike, T. Takata, et.al. 2018, PASJ, 70, S3
Eric C. Bellm, Shrinivas R. Kulkarni, Matthew J. Graham, et.al.
2019, PASP, 131, 018002
Frank J. Masci, Russ R. Laher, Ben Rusholme, et.al. 2019,
PASP, 131, 018003
S. Sako; R. Ohsawa, H. Takahashi, et.al. 2018, Proc. SPIE,
Volume 10702, id. 107020J 17 pp.
Tody Doug, 1986, in Crawford D. L., ed., Instrumentationin Astronomy VI,Vol.627 of Proc. SPIE, The IRAF Data
Reduction and Analysis System.p.733
Tody Doug, 1993, Astronomical Data Analysis Software andSystems II, ASP Conf. Ser.Series, Vol.52, 1993, R. J.
Hanisch, R. J. V. Brissenden, and Jeannette Barnes, eds.,
p.173.
T, Kotani, N. Kawai, K. Yanagisawa, et.al. 2005, Il Nuovo
17

https://github.com/MNiwano/Eclaire

10

Publications of the Astronomical Society of Japan, (2018), Vol. 00, No. 0

Cimento C, vol.28, Issue 4, p.755
Y. Yatsu, N. Kawai, T. Shimokawabe, et.al. 2007, Physica E,
vol.40, Issue 2, p.434
T. Shimokawabe, N. Kawai, T. Kotani, et.al. 2007, GAMMARAY BURSTS 2007:
Proceedings of the Santa Fe
Conference. AIP Conference Proceedings, vol.1000, pp.543
P. Tobias, V. Peter, P. Wolfgang, S. Johannes J, 2009, Journalof Computational Physics, Volume 228, Issue 12, p. 44684477.
D. Tuccillo, M. Huertas-Company, E. Decenciere, et.al. 2018,
Monthly Notices of the Royal Astronomical Society, Volume475, Issue 1, p.894-909
M. Tazzari, F. Beaujean and Leonardo Testi, 2018, Monthly
Notices of the Royal Astronomical Society, Volume 476,
Issue 4, p.4527-4542
Travis E. Oliphant, 2006, USA: Trelgol PublishingStéfan van der Walt, S. Chris Colbert and Gaël Varoquaux,
2011, Computing in Science & Engineering, 13, 22
Astropy Collaboration, Robitaille T. P, Tollerud E. J, et.al.
2013, Astronomy & Astrophysics, Volume 558, id.A33, 9pp.
Astropy Collaboration, Price-Whelan. A. M, Sipőcz. B. M,
et.al. 2018, The Astronomical Journal, Volume 156, Issue
3, article id. 123, 19 pp.
S. Toma, M. Niwano, R. Adachi, et.al. 2019, GRB Coordinates
Network, Circular Service, No.26019
R. Adachi, K. L. Murata, M. Oeda, et.al. 2019, GRB
Coordinates Network, Circular Service, No.26597
M. Niwano, K. L. Murata, M. Oeda, et.al. 2019, GRB
Coordinates Network, Circular Service, No.26776
R. Adachi, K. L. Murata, M. Niwano, et.al. 2020, GRB
Coordinates Network, Circular Service, No.26864
M. Oeda, R. Adachi, K. L. Murata, et.al. 2020, GRB
Coordinates Network, Circular Service, No.27236
R. Hosokawa,R. Adachi, K. L. Murata, et.al. 2020, GRB
Coordinates Network, Circular Service, No.27429