RGBA Packing for Fast Cone Beam Reconstruction on
the GPU
Fumihiko Inoa , Seiji Yoshidaa , and Kenichi Hagiharaa
a Graduate School of Information Science and Technology, Osaka University,

1-3 Machikaneyama, Toyonaka, Osaka 560-8531, Japan
ABSTRACTThis paper presents a fast cone beam reconstruction method accelerated on the graphics processing unit (GPU).
We implement the Feldkamp, Davis, and Kress (FDK) algorithm on the OpenGL graphics pipeline, which allows
us to exploit the full resources and capabilities available on the GPU. The proposed method diﬀers from previous
GPU-based methods in having an RGBA packing scheme capable of directly dealing with projections without
rebinning. It also reduces the amount of computation by using a data reuse scheme, which is useful to savethe memory bandwidth for this memory-intensive problem. Both schemes contribute to reduce the number ofrendering passes, namely the number of kernel invocations on the GPU, realizing fast cone beam reconstruction.
We show some experimental results obtained on a desktop PC with an nVIDIA GeForce 8800 GTX card. As aresult, the proposed method takes 8.1 seconds to reconstruct a 5123 -voxel volume from 360 5122 -pixel projection
images. This execution time is equivalent to a 15.6-fold speedup over a CPU implementation, showing 10%
higher performance as compared with a previous OpenGL-based method that requires the single-slice rebinningof projections for acceleration. With respect to non-rebinned data, our method provides approximately threetimes higher performance than the previous method.
Keywords: RECON, CTCB, ALG, GPU, high performance computing

1. INTRODUCTION
1Tomographic reconstruction is the process of reconstructing three-dimensional (3-D) objects from their integral
projections. In general, cone beam X-ray computed tomography (CT) scanners are required to reconstruct volumedata within a short time, because they usually assist surgical procedures. However, cone beam reconstruction is
a compute- and memory-intensive problem, which takes O(IN 3 ) time, where N represents the volume size andI represents the number of projections obtained during a scan. Although optimized CPU implementations thatuse SSE instructions2 can complete a reconstruction task in approximately two minutes,3 it is not suﬃcient toproduce a volume at an interactive rate. Therefore, we need some acceleration techniques to develop high-speed
imaging systems capable of producing a volume immediately after a scan.
One acceleration technique is to use accelerators such as the graphics processing unit (GPU),4–7 Cell Broadband Engine,3 and ﬁeld programmable gate array (FPGA).8 To the best of our knowledge, a GPU-based method4
demonstrates the highest performance of 8.9 seconds on an nVIDIA GeForce 8800 GTX card. The GPU9 here
is originally specialized to accelerate graphics tasks, but this commodity chip is rapidly increasing performancewith more ﬂexible programmability, allowing us to accelerate non-graphics tasks for interactive use. Although
the fastest performance is demonstrated by the GPU, the previous method4 requires the single-slice rebinning of
projections. If this requirement is not satisﬁed, it takes 2–3 times longer time for reconstruction.4 Furthermore,
the single-slice rebinning causes artifacts if applied to real projections obtained by a cone beam CT scanner.
Such a rebinning procedure can be applied to phantom data, where the cone angle can be easily set as almost
zero degree during projection, but not to real data, where the cone beam penetrates multiple slices due to a widecone angle.

Further author information: (Send correspondence to Fumihiko Ino)
Fumihiko Ino: E-mail: ino@ist.osaka-u.ac.jp

In this paper, we present a GPU-accelerated cone beam reconstruction method, aiming at developing highspeed cone beam CT scanners for intraoperative imaging. Our OpenGL-based method accelerates a widelyused ﬁltered backprojection algorithm, namely the Feldkamp, Davis, and Kress (FDK) method,1 using bothprogrammable processing units and hardwired components on the GPU. Our method has two main contributions
as compared with previous GPU-based methods. The ﬁrst one is an RGBA packing scheme that realizes fast
backprojection by packing four scalar data into an RGBA data. As compared with a previous method,4 ourscheme does not require rebinning before reconstruction. The proposed method makes the packing scheme more
ﬂexible with a 2-fold speedup over an unpacked implementation. The other one is a data reuse scheme thatomits computation by sharing computational results between diﬀerent voxels in the volume. This data reuse
scheme reduces the computational amount to achieve further acceleration. The two schemes mentioned above
reduce the number of rendering passes to 1/16, minimizing the CPU overhead such as texture switching neededat every rendering pass.
The paper is structured as follows. Section 2 presents related work and Section 3 summarizes the graphics
pipeline in the GPU. Section 4 describes our method with the acceleration schemes. Section 5 shows experimental
results obtained using a commodity graphics card. Finally, Section 6 concludes the paper with future direction.

2. RELATED WORKTo the best of our knowledge, Cabral et al.10 ﬁrstly present a GPU-based method that implements the FDK
algorithm using texture mapping hardware. Mueller et al.4, 5 also implement the FDK algorithm on the GPU.
They propose two methods,5 called MP-GPU and AG-GPU, which focus on programmable processing unitsof the GPU, such as vertex processors (VPs) and fragment processors (FPs). MP-GPU uses only FPs while
AG-GPU uses both FPs and VPs to achieve load balancing11 between processing units. AG-GPU improves
the performance of MP-GPU from 14.5 projections per second (pps) to 40.4 pps. They also show that earlyfragment kill (EFK) is useful to reduce the reconstruction time if region of interest (ROI) can be speciﬁed before
reconstruction.
With respect to RGBA packing, two schemes have been proposed in the past. The ﬁrst scheme12 packs four
orthogonal projections into a single RGBA projection. This scheme can be easily applied to phantom data but
it is not easy for real scanners to obtain such precise data due to the mechanical limitation. The other scheme4
packs adjacent four volume slices into a single RGBA slice. However, it cannot be directly applied to real data,
where the cone beam penetrates multiple slices due to a wide cone angle.
Schiwietz et al.7 propose a reconstruction pipeline running on the GPU. Riabkov et al.6 show that there
are two diﬀerent policies with respect to data assignment of the video memory. However, these methods are notfaster than AG-GPU.
Li et al.13 use compute uniﬁed device architecture (CUDA)14 to implement the FDK algorithm on the
nVIDIA GPU. Their method takes 3.6 seconds to perform backprojection (without ﬁltering of projections).
One drawback of CUDA-based methods is that they cannot be used with optimization techniques for graphics
applications, such as EFK.

3. GRAPHICS PIPELINE
Figure 1 shows an overview of the graphics pipeline in the GPU.9 This pipeline has two programmable processing
units, namely VPs and FPs, which are interconnected by the rasterizer. Each unit can access data stored on thevideo memory. Given vertex data that represents the shape and location of objects, the pipeline outputs pixels
of rendered images. Oﬀ-screen rendering techniques15 can be used to perform general purpose computation on
the GPU.
In the pipeline mechanism mentioned above, VPs are responsible for processing vertex data. Given vertex
data from VPs, the rasterizer generates fragments from the region speciﬁed by the vertex data. Fragment data
here corresponds to pixel data and has color and coordinate information linearly interpolated from the vertex
data. That is, the rasterizer converts per-vertex data into per-fragment data. The remaining programmableunits, namely FPs, take the responsibility for processing fragment data in parallel. The vertex program and thefragment program must be written as per-vertex operations and per-fragment operations, respectively.

MainmemoryCPUGPU
Vertex processors(VPs)
Rasterizer

Fragment processors
(FPs)

Video
memoryFigure 1. Graphics pipeline in the GPU.

As we mentioned earlier, fragment data has interpolated values of vertex data, so that linear operations forfragments can be replaced with operations for vertices. That is, some linear operations in the fragment program
can be moved to the vertex program. This code motion11 achieves load balancing between processing units. It
also reduces the time complexity if the number of vertices is smaller than that of fragments.
Note here that EFK is an optimization technique that limits the drawing region to reduce the number offragments passed from the rasterizer to FPs. Since FPs invoke the fragment program for every fragment, theeﬀects of EFK increases with the number of eliminated fragments. This reduction contributes to less amount ofcomputation as well as less amount of data transfer between the video memory and FPs.
The GPU is allowed to access data only on the video memory. Therefore, the data must be downloaded fromthe main memory to the video memory in advance of rendering. In contrast, the data must be readback to themain memory if the CPU requires results computed by the GPU.
A typical implementation strategy for GPGPU applications is to store the input/output data as twodimensional (2-D) textures on the video memory such that texels on output textures can be computed in parallel.
Since the GPU is optimized for 2-D textures, it is better to convert 3-D data into 2-D textures. RGBA textures
can have four-component data as a single texel data. Each component here can have a 32-bit ﬂoating point
number.

4. METHODS
1

The FDK algorithm consists of two processing stages: the ﬁltering stage and the backprojection stage. Our
method accelerates both stages on the GPU. However, we describe here only the backprojection stage, which isthe performance bottleneck of this algorithm.
Figure 2 illustrates the coordinate system used in our method. There is the volume F and a projection Qi
taken from the X-ray source at distance d0 with angle θi , where 1 ≤ i ≤ I. Notation di represents the distance
between the X-ray source and the volume origin. Let F (x, y, z) be the value of voxel (x, y, z) in the volume F .
Given I ﬁltered projections Q1 , Q2 , . . . , QI , the backprojection stage generates the volume F according to thefollowing equation:I
1 ∑F (x, y, z) =W2 (x, y, i)Qi (u(x, y, i), v(x, y, z, i)),
(1)2πI i=1

vzuy
X-ray source
θi
x
di
d’

Figure 2. Coordinate system used in our method.

v

u
z

z
Projections

Projection(a)Volume slice(b)Figure 3. (a) RGBA packing scheme used in our method. Our scheme packs four ﬁltered projections into a single RGBA
texture. (b) Previous scheme4 packs four succeeding pixels along the v-axis into a single RGBA pixel. Therefore, it
requires the single-slice rebinning of projections and assumes that projections have the same resolution as the volume.
Both schemes also pack four succeeding volume slices perpendicular to the z-axis.

whereu(x, y, i)v(x, y, z, i)
W2 (x, y, i)d0 (−x sin θi + y cos θi ),di − x cos θi − y sin θid0 z=,
di − x cos θi − y sin θi
)2
(
di
.
=
di − x cos θi − y sin θi
=

(2)
(3)
(4)

Since the GPU architecture is optimized to 2-D data structure, a naive method decomposes the volume F
into slices F0 , F1 , . . . , FN −1 , each stored in a 2-D texture. Similarly, each of projections is stored in a 2-D texture.
In this case, we have to render IN times to backproject I projections into N slices. The proposed method isdesigned to reduce this iteration number to IN/16 by packing input/output textures.

4.1 RGBA Packing Scheme
There is no data dependence between backprojection of an image and that of another image. Therefore, we
can simultaneously process the backprojection of multiple images. To do this, we pack four ﬁltered images
Qi , Qi+1 , Qi+2 , and Qi+3 into a single RGBA texture. Figure 3(a) shows how our RGBA packing scheme
works at the backprojection stage. In our scheme, an RGBA texel (u, v) in the packed texture contains scalartexels at the same coordinates: Qi (u, v), Qi+1 (u, v), Qi+2 (u, v), and Qi+3 (u, v). This allows us to utilize vector
operations to accelerate computation, as shown in Fig. 4. For example, we can vectorize the computationof u, v, and W2 for packed texels Qi (u, v), Qi+1 (u, v), Qi+2 (u, v), and Qi+3 (u, v) at lines 18 and 20 in Fig.
4. That is, a vector of u(x, y, i), . . . , u(x, y, i + 3) and that of W2 (x, y, i), . . . , W2 (x, y, i + 3) are computed at

Input: Projections P1 , P2 , . . . , PI and parameters d0 , d1 , d2 , . . . , dI , θ1 , θ2 , . . . , θI
Output: Volume slices F0 , F1 , . . . , FN −1
1: ProposedReconstruction()
2: begin
3: Download projections P1 , . . . , PI ;4: (Q1 , . . . , QI ) ← Filtering(P1 , . . . , PI );5: for n = 0 to N/4 − 1 do begin6:
z = 4n;7:
for m = 0 to I/4 − 1 do begin8:
i = 4m + 1;
9:
(Fz , Fz+1 , Fz+2 , Fz+3 )← Backprojection((Qi , Qi+1 , Qi+2 , Qi+3 ), (di , di+1 , di+2 , di+3 ),(θi , θi+1 , θi+2 , θi+3 ), d0 , z); // Rendering
10:
end;
11:
Readback volume slices Fz , . . . , Fz+3 ;12: end;
13: end
14: function Backprojection((Qi , Qi+1 , Qi+2 , Qi+3 ), d, θ, d0 , z)
15: begin
16: // x and y is given by the GPU
17: Compute d − x cos θ − y sin θ;
18: Compute u and W 2 ; // Eqs. (2) and (4)19: for n = 0 to 3 do begin
20:
Compute v with z = z + n; // Eq. (3)21:Convert u and v to texture coordinates coord1 , . . . , coord4 ;22:
for j = 1 to 4 do begin
23:
valj ← Interpolated texel at coordj in Qi+j ; // bilinear interpolation
24:
end;
25:
val ← W 2 · val;
26:
outn+1 ← val1 + val2 + val3 + val4 ;27: end;
28: current ← (Fz (x, y), Fz+1 (x, y), Fz+2 (x, y), Fz+3 (x, y));29: return out + current ;30: endFigure 4. Pseoudocode of the proposed method. An RGBA value that contains four scalar values a1 , . . . , a4 is denoted
as a in bold font. The number of rendering passes is reduced from IN to IN/16 in our method. A hardware accelerated
interpolation mechanism is employed to obtain a texel at line 23.

a time. Similarly, we can apply vectorization to computation of fi (x, y, z), . . . , fi+3 (x, y, z) at line 25, wherefi (x, y, z) = W2 (x, y, i)Qi (u(x, y, i), v(x, y, z, i)).
The RGBA packing scheme allows us to backproject four projections per rendering pass. Therefore, we canreduce the number of rendering passes by 75%. These fewer passes contribute to reduce the CPU overhead, suchas texture switching incurred before each rendering.

4.2 Data Reuse Scheme
Equations (2)–(4) indicate that some terms are independent of z: u(x, y, i), W2 (x, y, i), and the denominator ofv(x, y, z, i). Therefore, such z-independent terms must be computed for every line (x, y) rather than every voxel(x, y, z). On the other hand, our method can process every four slices of the volume at a time if the slices arepacked into an RGBA texture. Therefore, as we did for the input texture, we also pack the output texture inorder to reuse the z-independent terms within four volume slices.
The data reuse scheme further reduces the number of rendering passes, because we have 75% fewer output
textures. Thus, we have to invoke the rendering kernel (fragment program) only IN/16 times, which is shownby the loops at lines 5–7 in Fig. 4. In contrast, we have more iterations at lines 19 and 22 in the invoked kernel.
This means that we have more computation per rendering pass. Such a coarse-grained rendering task leads to aTable 1. Implementations used for experiments.

Implementation
A
B
C
D
E

Data reuse
No
Yes
No
Yes
Yes

RGBA packing
No
No
Yes
Yes
Yes

Load balancing5, 11
No
No
No
No
Yes

25
135

A

Execution Time (s)

20

15
C
10

D

5

0
CPU3

MPGPU5AGGPU5

Rapid
CT4

Naive w/ RGBA w/ data
method packing reuse

Figure 5. Execution time for reconstruction of 5123 -voxel volume from 360 5122 -pixel projections. The right three bars
correspond to our methods and the remaining bars represent previous methods.3–5 The results for previous methods arequoted from their papers, because it is measured on the same GPU.

higher eﬃciency because each rendering task is parallelized by the GPU. In particular, it is useful to improve theeﬃciency if the kernel takes a relatively shorter time than its preprocessing overhead such as texture switching.

5. RESULTS
To evaluate the proposed method in terms of performance and image quality, we have implemented the method
using the C++ language, the OpenGL library, and the Cg toolkit. We also have implemented diﬀerent versionsof the proposed method to analyze the eﬀect of the data reuse scheme, the RGBA packing scheme, and the loadbalancing scheme.5, 11 Table 1 summarizes the diﬀerences between implementations.
We perform experiments using a desktop computer equipped with a Core 2 Extreme CPU running at 2.93
GHz clock speed. This computer has an nVIDIA GeForce 8800 GTX card, running on driver version 162.01. We
use Windows XP for the operation system. Reconstruction is done for the Shepp-Logan phantom and clinicalbrain data. The output volume is stored in frame buﬀer object (FBO)15 to perform oﬀ-screen rendering.
Figure 5 shows the timing results. With respect to our method, the execution time contains times for RGBA
packing, download, ﬁltering, backprojection, readback, and RGBA unpacking, as shown in Table 2. As a result,the proposed method takes 8.1 seconds to reconstruct a 5123 -voxel volume from 360 5122 -pixel projections. This
execution time is equivalent to a 15.6-fold speedup over an SSE-optimized CPU implementation.3 It also achieves
10% higher performance as compared with the fastest previous method.4 Thus, the proposed method is usefulto produce volume data immediately after a scan of real objects.
We also analyze the eﬀectiveness of our method in terms of the memory bandwidth and the ﬂoating point
operations per second (FLOPS). Table 3 shows that the eﬀective bandwidth reaches 70.8 GB/s out of 86.4 GB/s.
With respect to the ﬂoating point performance, it provides 231.6 GFLOPS out of 518.4 GFLOPS. Thus, theperformance is limited by the memory bandwidth. Therefore, we have to reduce the amount of data being fetchedfrom video memory in order to achieve further acceleration. Such a reduction can be realized by EFK,5 whichis available in the graphics pipeline but not available in CUDA. Actually, we observe that EFK further reduces
the reconstruction time to 7 seconds if 75% of volume region is skipped during backprojection.

Table 2. Breakdown of execution time (s). Total time here includes the synchronization overhead between the CPU andthe GPU. Such an overhead incurs only when obtaining breakdown information.

Breakdown
RGBA packing
Download
Filtering
Backprojection
Readback
RGBA unpacking
OthersTotal time
Total time w/o synchronization

A
0.3
0.4
5.6
16.3
0.5
1.0
3.2
27.2
20.7

Implementation
B
C
D
0.3
0.6
0.6
0.4
0.5
0.5
5.6
2.5
2.4
7.4
6.8
3.7
0.4
0.5
0.4
0.9
0.9
0.6
4.0
1.3
1.4
18.9 13.0
9.6
17.0 10.8
8.3

E
0.6
0.5
2.4
4.0
0.4
0.6
1.9
10.4
8.1

Table 3. Eﬀective performance and bandwidth of each implementation. Note that implementation E reduces the timecomplexity by load balancing. This slightly reduces the reconstruction time but decreases the eﬃciency in terms of ﬂoating
point performance.

Metrics
Processing unit
Texture unit
Total
Memory bandwidth (GB/s)
Computation(GFLOPS)

A
117.0
33.9
150.9
28.7

Eﬀective performance
B
C
D
112.0
62.6
99.3
69.4
70.6 132.3
181.5 133.1 231.6
37.1
59.8
70.8

E
61.3
128.3
189.6
68.6

Theoretical
performance
345.6
172.8
518.4
86.4

With respect to the image quality, GPU-based methods cannot reconstruct volumes identical to those generated by CPU-based methods. This is due to the fact that the GPU has a diﬀerent architecture that does notprovide the same instruction set as the CPU. For example, a multiply-add operation is executed by a single
instruction on the GPU while it is executed by two instructions on the CPU. In addition, the GPU has larger
errors than the CPU. Despite of these diﬀerences, we ﬁnd that this problem is not critical in practical situations.
Actually, there is no visible diﬀerence in visualization results, as shown in Figs. 6(a) and (b). In addition, peaksignal-to-noise ratio (PSNR) values of the raw volume are more than 100 dB, which indicates that it is hard tosee image degradation.

6. CONCLUSION
We have presented a GPU-accelerated reconstruction method for cone beam CT scanners. Our method makes
it possible to pack projections into RGBA data without rebinning. It also reuses data to reduce the amountof computation and to save the memory bandwidth for this memory-intensive problem. Both RGBA packing

(a)

(b)

Figure 6. Maximum intensity projection (MIP) of the brain data generated (a) on the GPU and (b) on the CPU.

and data reuse schemes contribute to minimize the CPU overhead, such as texture switching needed at every
rendering pass.
As a result, the proposed method takes 8.1 seconds to reconstruct 5123 -voxel volume, which is 10% fasterthan a previous method. In addition, our method does not assume rebinned data. With respect to the eﬃciency,
the eﬀective bandwidth reaches a peak of 70.8 GB/s, which is 82% of the theoretical bandwidth. Therefore,
we think that further acceleration will be achieved if reducing the data amount being fetched from the video
memory.
Future work includes the support of large volume such as 10243 voxel data. Some data decomposition
techniques are needed to deal with such large size on current graphics cards.

ACKNOWLEDGMENTS
This work was partly supported by JSPS Grant-in-Aid for Scientiﬁc Research (A)(2)(20240002), Young Researchers (B)(19700061), and the Global COE Program “in silico medicine” at Osaka University. We are also
grateful to Dr. Kazuyoshi Nishino, Dr. Yoshihiro Mukuta, and Dr. Yukio Mishina from Shimadzu for theirvaluable discussions as well as clinical data.

REFERENCES[1] Feldkamp, L. A., Davis, L. C., and Kress, J. W., “Practical cone-beam algorithm,” J. Optical Society ofAmerica 1, 612–619 (June 1984).
[2] Klimovitski, A., “Using SSE and SSE2: Misconceptions and reality,” in [Intel Developer Update Magazine ],(Mar. 2001).
[3] Kachelrieß, M., Knaup, M., and Bockenbach, O., “Hyperfast parallel-beam and cone-beam backprojection
using the cell general purpose hardware,” Medical Physics 34, 1474–1486 (Apr. 2007).
[4] Mueller, K., Xu, F., and Nephytou, N., “Why do commodity graphics hardware boards (GPUs) work sowell for acceleration of computed tomography?,” in [Proc. SPIE Medical Imaging 2007 ], (Feb. 2007).
[5] Xu, F. and Mueller, K., “Real-time 3D computed tomographic reconstruction using commodity graphics
hardware,” Physics in Medicine and Biology 52, 3405–3419 (June 2007).
[6] Riabkov, D., Xue, X., Tubbs, D., and Cheryauka, A., “Accelerated cone-beam backprojection using GPUCPU hardware,” in [Proc. 9th Int’l Meeting Fully Three-Dimensional Image Reconstruction in Radiology
and Nuclear Medicine (Fully 3D ’07)], 68–71 (July 2007).
[7] Schiwietz, T., Bose, S., Maltz, J., and Westermann, R., “A fast and high-quality cone beam reconstruction
pipeline using the GPU,” in [Proc. SPIE Medical Imaging 2007 ], 1279–1290 (Feb. 2007).
[8] Gac, N., Mancini, S., and Desvignes, M., “Hardware/software 2D-3D backprojection on a SoPC platform,”in [Proc. 21st ACM Symp. Applied Computing (SAC’06)], 222–228 (Apr. 2006).
[9] Owens, J. D., Luebke, D., Govindaraju, N., Harris, M., Krüger, J., Lefohn, A. E., and Purcell, T. J., “Asurvey of general-purpose computation on graphics hardware,” Computer Graphics Forum 26, 80–113 (Mar.
2007).
[10] Cabral, B., Cam, N., and Foran, J., “Accelerated volume rendering and tomographic reconstruction using
texture mapping hardware,” in [Proc. 1st Symp. Volume Visualization (VVS’94) ], 91–98 (Oct. 1994).
[11] Ikeda, T., Ino, F., and Hagihara, K., “A code motion technique for accelerating general-purpose computationon the GPU,” in [Proc. 20th IEEE Int’l Parallel and Distributed Processing Symp. (IPDPS’06) ], (Apr. 2006).
10 pages (CD-ROM).
[12] Xu, F. and Mueller, K., “Accelerating popular tomographic reconstruction algorithms on commodity PC
graphics hardware,” IEEE Trans. Nuclear Science 52, 654–663 (June 2005).
[13] Li, M., Yang, H., Koizumi, K., and Kudo, H., “Fast cone-beam CT reconstruction using CUDA architecture,” Medical Imaging Technology 25, 243–250 (Sept. 2007). (In Japanese).
[14] nVIDIA Corporation, “CUDA Programming Guide Version 2.0,” (July 2008).
[15] OpenGL Extension Registry, “Gl ext framebuﬀer object,” (2006).
http://oss.sgi.com/projects/
ogl-sample/registry/EXT/framebuffer object.txt.