OUT-OF-CORE CONE BEAM RECONSTRUCTION USING MULTIPLE GPUS
Fumihiko Ino† , Yusuke Okitsu† , Taketo Kishi‡ , Syuhei Ohnishi‡ , Kenichi Hagihara††

Graduate School of Information Science and Technology, Osaka University
1-5 Yamadaoka, Suita, Osaka 565-0871, Japan‡
Analytical and Measuring Instruments Division, Simadzu Corporation
1 Nishinokyo-Kuwabara-cho, Nakagyo-ku, Kyoto 604-8511, Japan

ABSTRACTThis paper presents a graphics processing unit (GPU) based
method capable of accelerating cone-beam reconstruction oflarge volume data, which cannot be entirely stored in video
memory. Our method accelerates the Feldkamp, Davis andKress (FDK) algorithm in a multi-GPU environment. We
present how the entire volume can be efﬁciently decomposed
into small portions to reduce the usage of video memory oneach graphics card. Experimental results are also presented
to understand the reconstruction throughput on an nVIDIA
Tesla S1070 server. It takes approximately three minutes toreconstruct a 20483 -voxel volume from 720 20482 -pixel projections. The effective bandwidth of video memory reaches
137 GB/s per GPU, demonstrating a higher utilization of texture caches.

In this paper, we present how out-of-core reconstruction
can be efﬁciently realized for a large volume, which cannot be
entirely stored in video memory. Our algorithm decomposes
the volume into small portions, which are then processed using an in-core algorithm [2] in a multi-GPU environment. The
in-core algorithm is based on the Feldkamp, Davis and Kress(FDK) algorithm [5], which is widely used in many CBCT
systems. Our algorithm parallelizes the ﬁltering stage andthe backprojection stage, which are the time-consuming partof the FDK algorithm. The algorithm is implemented using
a general-purpose programming framework, called the compute uniﬁed device architecture (CUDA) [6].
The outline of this paper is as follows. Section 2 introduces previous work. Section 3 then describes the out-of-corealgorithm. Sections 4 and 5 present results and conclusions.

1. INTRODUCTION

2. RELATED WORK

Cone beam (CB) reconstruction is an imaging technique for
producing a three-dimensional (3-D) volume from a seriesof 2-D projections acquired by a CB computed tomography(CBCT) scan. This technique plays an important role in various ﬁelds, such as clinical diagnosis and nondestructive inspection of engineering products. For many years, CB reconstruction has been a compute-intensive task to us. However,
the graphics processing unit (GPU) has emerged as a powerfulaccelerator that achieves a 30-fold speedup over a CPU-based
method optimized using SSE instructions [1]. For example, areconstruction task of a 5123 -voxel volume can be completed
within several seconds using the commodity GPU.
However, most GPU-based methods [2–4] assume that thevolume is small enough to store the entire data in video memory. This assumption restricts the volume size by 5123 voxels
on current commodity GPUs with at most 1 GB of memory.
Though this limitation was not so serious problem, ﬂat panel
detectors in recent CBCT systems now have more than 10242
detector units. Since the data size reaches 4 GB for a 10243 voxel volume, an out-of-core algorithm is needed to reconstruct a large volume from high-resolution projections.

Many researchers [2–4] are trying to accelerate CB reconstruction on the GPU. Okitsu et al. [2] show that texture cache
optimization is the key technique to accelerate the backprojection stage on the CUDA-compatible GPU. Scherl et al. [3]propose another optimization strategy that minimizes the resource usage to exploit higher parallelism on the GPU. In contrast to these CUDA-based implementations, Yan et al. [4]propose a graphics-based method that implements the FDK
reconstruction using the OpenGL library. Their implementation takes only 5.2 seconds to reconstruct a 5123 -voxel volume from 360 5122 -pixel projections. The implementations
mentioned above are based on in-core algorithms that requires
the entire volume to be stored in video memory.
In contrast to these in-core algorithms, out-of-core algorithms can deal with larger data. Nesterets et al. [7] proposes
an efﬁcient algorithm that uses as little video memory for thereconstruction of each axial slice of the object as possible.
Their algorithm carries out reconstruction in a slice-by-slicemanner. It takes about 81 minutes on a GeForce GTX 260
card to reconstruct a 20483 -voxel volume from 1440 20482 pixel projections. A similar work is also done by Nöel etal. [8]. Their implementation runs on a GeForce 8800 GT
card and takes 61.3 seconds to reconstruct a 10243 -voxel volume from 106 10243 -pixel projections. Zhao et al. [9] im-This work was partly supported by JSPS Grant-in-Aid for Scientiﬁc Research (A)(2)(20240002) and the Global COE Program “in silico medicine”
at Osaka University.

978-1-4244-4126-6/10/$25.00 ©2010 IEEE

792

ISBI 2010

GPU #1

GPU #1…

MainmemoryY…

GPU #C

B/C
subvolumes

K filteredprojections

K/C projections…

plement an out-of-core algorithm using the OpenGL library.
It takes about 101.9 seconds on a Quadro FX 4600 card toreconstruct a 10243 -voxel volume from 720 10242 -pixel projections. In contrast to these out-of-core algorithms, we showthat our algorithm runs more efﬁciently with reaching the theoretical performance bound.

GPU #C

3. METHODSZ
XConcatenatedvolume…

Consider a reconstruction task that produces an N 3 -voxel volume V from a series of K M 2 -pixel projections using C
GPUs. The FDK algorithm [5] solves this problem by processing the ﬁltering stage and the backprojection stage. We
currently use the Shepp-Logan ﬁlter [10] at the ﬁltering stage.
3.1. Parallelization
The data dependence at the ﬁltering stage and at the backprojection stage can be summarized as follows.
Filtering stage: Different pixels on the same projection canbe simultaneously ﬁltered at this stage, because there isno data dependence between them. Furthermore, different projections can also be processed at the sametime, because they do not have data dependence.
Backprojection stage: The backprojection stage must beprocessed after the ﬁltering stage. However, different
voxels can be simultaneously processed at the backprojection stage. In addition, the backprojection operator
has an associativity property, so that projections can beprocessed in an arbitrary order. All projections must be
backprojected into every voxel in the volume.
According to the analysis mentioned above, we have decided to decompose the volume data V into multiple subvolumes V1 , V2 , · · · , VB , where B (> 1), because the volume
size usually limits the maximum reconstruction size. We determine the value of B such that (1) every subvolume is smallenough to store in video memory and (2) every GPU has atleast one assigned subvolume (B ≥ C). After this decomposition, the in-core algorithm [2] can be applied to each subvolume of N 3 /B voxels.
Figure 1 shows an overview of the proposed method,which exploits data parallelism in each stage. As shown inthis ﬁgure, each of C GPUs is responsible for K/C projections at the ﬁltering stage and for B/C subvolumes atthe backprojection stage. Both stages are processed by thein-core algorithm [2], which performs the ﬁltering of I projections at a time and performs the backprojection of I ﬁlteredprojections into a subvolume at a time, where 1 ≤ I ≤ K.
Therefore, the algorithm iteratively invokes the GPU program
to process all projections and subvolumes. The value of I isexperimentally determined as I = 5.
Notice here that all ﬁltered projections are sent back tomain memory before the backprojection stage. This intends
to gather projections and broadcast all of them to GPUs. The

793

K/C projections

Filtering stage

K filteredprojections

B/C
subvolumes

Backprojection stage

Fig. 1. Overview of proposed out-of-core algorithm.
Input: Projections P1 , P2 , . . . , PK , number B of subvolumes,
Output: Volume V
Algorithm Cone Beam Reconstruction()
1: Load P1 , P2 , . . . , PK from storage device;
2: for k = 1 to K do in parallel
3: Download Pk from main memory to video memory;
4: Qk ← Filtering(Pk );// Qk : k-th ﬁltered projection
5: Readback Qk to main memory;
6: endfor
7: for b = 1 to B do in parallel
8: Initialize subvolume Vb in video memory;
9: for k = 1 to K do
10:
Download Qk from main memory to video memory;
11:
Bind Qk as texture;
12:
Vb ← Backprojection(Qk , k, b);
13: endfor
14: Readback subvolume Vb to main memory;
15: endfor

Fig. 2. Pseudocode of proposed algorithm. The actual codeis optimized to process the ﬁltering and backprojection of I
projections at a time.
reason why we gather ﬁltered projections is that all K ﬁlteredprojections are needed at the backprojection stage on every
GPU, though each GPU is responsible for a part of the volume. Thus, our algorithm requires data transfer but prevents
us from ﬁltering all of K projections on every GPU. After thisparallel backprojection, B subvolumes are then concatenatedinto the ﬁnal volume.
Figure 2 shows a pseudocode of the proposed algorithm.
In this code, “download” represents the data transfer frommain memory to video memory. Similarly, “readback” represents the transfer in the opposite direction. Lines 2–6 and7–15 correspond to the ﬁltering stage and the backprojection
stage, respectively. In the backprojection stage, ﬁltered projections are stored in textures [6] to use texture units for interpolation of pixels on projections. Since texture units are separated from processing elements in the GPU, using textures isnecessary to maximize the entire performance by ofﬂoading
workloads from processing elements to texture units.
Note here that the code is optimized using an asynchronous execution technique [6]. This technique is applied

45
Throughput

4.5

4.5

4.0

4

Speedup

25

2.5

20

2.0

15

1.5

10

1.0

5
0

3.5

3.5

3

3.0

2.5

2.5

2

2.0

1.5

1.5

1

1.0

0.5

0.5

0.5

0.0

0

Throughput (pps)3.0
Speedup

3.5

30

C=2

0.0
C=1

C=4

(a)

to the ﬁltering stage to overlap data transfer (lines 3 and 5)
with GPU computation (line 4). In contrast, we cannot apply
this technique to the backprojection stage because textures
currently have to be sent in synchronous mode.
Finally, we explain how our algorithm decomposes thevolume. As shown in Fig. 1, our algorithm uses a 2-D
block decomposition scheme that separates the volume spacewith respect to the (x, y) coordinate. In other words, thisscheme does not divide the Z (axial) space, so that any
voxel V (x, y, z) located on the same (x, y) coordinate will beassigned to the same GPU. This allows us to maximize the effect of a data reuse technique [2], which reduces the amountof computation between different axial slices. Thus, data
decomposition depends on the data reuse technique because
data cannot be reused between different GPUs.

4.0

Speedup

35

C=1

4.5
Throughput

Speedup

40

Throughput (pps)

Table 1. Analytical results of time complexity and spacecomplexity. Parameter I represents the number of projections
processed at a time. We currently use I = 5.
Item
Original [2]Proposed
Download
O(KM 2 )O(BKM 2 /C)
FilteringO(KM 2 )O(KM 2 /C)
Readback (projection)—
O(KN 3 /C)
Backprojection
O(KN 3 )O(KN 3 /C)
3
Readback (volume)O(N )O(N 3 /C)
Concatenation—
O(B)Video memory usage O(N 3 + IM 2 ) O(N 3 /B + IM 2 )Main memory usage O(N 3 + KM 2 ) O(N 3 + KM 2 )

C=2

C=4

(b)

Fig. 3. Throughput and speedup ratio (a) for N = M = 1024
and (b) for N = M = 2048. The throughput is presented interms of the number of projections per second (pps).
Table 2. Breakdown of execution time T in seconds.
N = M = 1024
N = M = 2048
Breakdown
C=1C=2C=4C=1C=2C=4
Initialization
1.3
2.7
5.4
6.1
5.5
7.5
Download
4.4
2.4
1.3 255.8 190.7 103.6
Filtering
5.2
2.6
1.3 20.8 10.4
5.2
Readback (prj.) 1.1
0.6
0.3
3.8
3.5
1.6
Backprojection 39.7 19.9 10.0 235.3 117.5 59.3
Readback (vol.) 1.3
0.9
0.6 10.2
6.3
3.6
T (no stream)
53.0 29.1 18.9 532.0 333.9 180.8
T (stream)
51.5 28.4 18.0 523.8 328.3 178.4

3.2. Analytical Analysis
Table 1 shows an analytical analysis of the proposed algorithm and the original in-core algorithm [2] in terms of timecomplexity and space complexity. As shown in this table,
data decomposition allows us to reduce the space complexity
so that out-of-core reconstruction can be realized by selecting
the appropriate value for B.
With respect to the time complexity, our algorithm basically reduces the original complexity by a factor of 1/C because it fully exploits data parallelism using C GPUs. However, the time complexity of data download is the exception.
This can be explained by the pseudocode in Fig. 2, which has
a nested loop at line 9. The loop implies that the same projections can be iteratively sent to video memory for different
subvolumes. Though this redundancy increases the complexity by a factor of B, it saves the usage of video memory todeal with large volume data.
Note here that Table 1 shows results for per GPU. Therefore, the measured results can differ from them. For example,
we currently cannot send data to multiple GPUs at the sametime. This is due to the graphics driver, which sequentiallysends the data in a ﬁrst-in, ﬁrst-out (FIFO) manner. Furthermore, our algorithm is optimized using stream, which overlaps data transfer with program execution on the GPU.

794

4. RESULTS
We now show experimental results to evaluate the performance of the proposed method. A reconstruction of10243 /20483 -voxel volume from 720 10242 /20482 -pixel projections is carried out on an nVIDIA Tesla S1070 server(K = 720). The 10243 -voxel data and 20483 -voxel data
are decomposed into 4 portions (B = 4) and 64 portions(B = 64), respectively. The server is connected to a PC thathas a Xeon E5450 CPU and 16-GB main memory. It runs on
CentOS 5.3 with CUDA 2.2 [6] and driver version 185.18.08.
Figure 3 shows the reconstruction throughput in terms ofthe number of projections per second (pps). The throughput reaches 4 pps when using all four GPUs for N = 2048
(C = 4). On a single GPU, the throughput for 10243 -voxel
data and that for 20483 -voxel data are 14 pps and 1.4 pps, respectively. The speedup over a single GPU version is arounda factor of 3, which is slightly lower than the optimal speedup.
This inefﬁciency can be explained by Table 2, which shows
the breakdown of execution time T for each execution conﬁguration. Although the ﬁltering stage and the backprojection
stage are linearly accelerated using multiple GPUs, the GPU
initialization and data transfer emerge as performance bottle-

1
350
300

0.8

250
0.6
200
150

0.4

100
0.2
50
0

0
C=1

C=2

C=4

1.8
Effective arithmetic performance (Tflops)

400

Effective Memory bandwidth (GB/s)

Effective arithmetic performance (Tflops)

Arithmetic
Bandwidth

600
Arithmetic
Bandwidth

1.6

500
1.4
1.2

400

1
300
0.8
200

0.6
0.4

100
0.2

Table 3. Performance comparison with previous work.
Throughput (pps)
GPU speciﬁcation
Work
Normalized
Bandwidth Arithmetic(Measured)(GB/s)(Tﬂops)[8]3.9 ( 1.7)
44.8
0.504
[9]
10.8 ( 7.1)
67.2
0.336
This paper
14.0 (14.0)
102.0
0.933
[7]
0.3 ( 0.3)
111.9
0.715
This paper
1.4 ( 1.4)
102.0
0.933

Effective memory bandwidth (GB/s)

450

1.2

0

0
C=1

(a)

C=2

C=4

(b)

Fig. 4. Effective arithmetic performance and memory bandwidth of backprojection measured (a) for N = M = 1024
and (b) for N = M = 2048.
neck as we increase the number of GPUs. This might be dueto the limitation of current graphics drivers, which are notfully multithreaded yet. Thus, GPUs are serially initializedin our implementation. Similarly, graphics drivers currently
cannot simultaneously exchange data with multiple GPUs.
Figure 4 shows the effective arithmetic performance andmemory bandwidth at the backprojection stage. The theoretical bound of our GPU is 102 GB/s and 0.933 Tﬂops (teraﬂoating point number operations per second) for the memory
bandwidth and for the arithmetic performance, respectively.
As shown in this ﬁgure, we ﬁnd that the effective bandwidthfor N = M = 2048 reaches 137 GB/s per GPU, which exceeds the theoretical bound. This can be explained by texture caches, which saves the bandwidth of video memory. In
contrast, the arithmetic performance is lower than 0.4 Tﬂops,
which is equivalent to 41% efﬁciency. This lower efﬁciency
indicates that processing elements in the GPU have to wait fordata to be fetched from video memory. Therefore, we have tofurther increase the effective bandwidth to improve the effective arithmetic performance.
Finally, Table 3 shows a performance comparison withprevious algorithms [7–9]. Since different GPUs are used inexperiments, we have normalized the performance according
to the bandwidth of video memory. The proposed algorithm
achieves the highest efﬁciency among known algorithms.
5. CONCLUSION
We have presented an out-of-core reconstruction algorithm
for large volume data. The algorithm efﬁciently decomposes
the volume data into small portions. It also exploits data parallelism in the FDK algorithm to maximize the performancein a multi-GPU environment. In experiments, our algorithm
takes about 3 minutes to reconstruct a 20482 -voxel volume
from 720 20482 -pixel projections. We also show that the795

measured bandwidth exceeds the theoretical bandwidth forthe sake of texture cache utilization. Therefore, we think thatthe measured performance is close to the theoretical bound.
In the future, we plan to further minimize the amount of datatransfer between video memory and main memory.
6. REFERENCES[1] Marc Kachelrieß, Michael Knaup, and Olivier Bockenbach,“Hyperfast parallel-beam and cone-beam backprojection using
the cell general purpose hardware,” Medical Physics, vol. 34,
no. 4, pp. 1474–1486, Apr. 2007.
[2] Yusuke Okitsu, Fumihiko Ino, and Kenichi Hagihara, “Fastcone beam reconstruction using the CUDA-enabled GPU,”in Proc. 15th Int’l Conf. High Performance Computing(HiPC’08), Dec. 2008, pp. 108–119.
[3] Holger Scherl, Benjamin Keck, Markus Kowarschik, andJoachim Hornegger, “Fast GPU-based CT reconstruction using the common uniﬁed device architecture (CUDA),” inProc. Nuclear Science Symp. and Medical Imaging Conf.
(NSS/MIC’07), Oct. 2007, pp. 4464–4466.
[4] Guorui Yan, Jie Tian, Shouping Zhu, Yakang Dai, andChenghu Qin, “Fast cone-beam CT image reconstruction using
GPU hardware,” J. X-Ray Science and Technology, vol. 16, no.
4, pp. 225–234, Oct. 2008.
[5] L. A. Feldkamp, L. C. Davis, and J. W. Kress, “Practical conebeam algorithm,” J. Optical Society of America, vol. 1, no. 6,
pp. 612–619, June 1984.
[6] nVIDIA Corporation, “CUDA Programming Guide Version
2.2,” April 2009, http://developer.nvidia.com/cuda/.
[7] Ya. I. Nesterets and T. E. Gureyev, “High-performance tomographic reconstruction using graphics processing units,” inProc. 18th World IMACS/MODSIM Congress Modeling and
Simulation (MODSIM’08), July 2009, 7 pages.
[8] Peter B. Noël, Alan M. Walczak, Kenneth R. Hoffmann, Jinhui
Xu, Jason J. Corso, and Sebastian Schafer, “Clinical evaluationof GPU-based cone beam computed tomography,” in Proc.
High-Performance Medical Image Computing and Computer
Aided Intervention (HP-MICCAI’08), Sept. 2008.
[9] Xing Zhao, Jing jing Hu, and Peng Zhang, “GPU-based 3D
cone-beam CT reconstruction for large data volume,” Int’l J.
Biomedical Imaging, Article ID 149079, 2009, 8 pages.
[10] L. A. Shepp and B. F. Logan, “The fourier reconstruction of ahead section,” IEEE Trans. Nuclear Science, vol. 21, no. 3, pp.
21–43, June 1974.