DOI:10.1587/transinf.2016EDP7178
Publicized:2016/12/12
This advance publication article will be replaced by the finalized version
after proofreading.

IEICE TRANS. ??, VOL.Exx–??, NO.xx XXXX 200x

1

PAPER

Cache-aware, In-place Rotation Method for Texture-based
Volume Rendering
Yuji MISAKI† , Nonmember, Fumihiko INO† , and Kenichi HAGIHARA† , Members

SUMMARY We propose a cache-aware method to accelerate texture-based volume rendering on a graphics processing unit
(GPU) that is compatible with the compute unified device architecture. The proposed method extends a previous method suchthat it can maximize the average rendering performance whilerotating the viewing direction around a volume. To realize this,the proposed method performs in-place rotation of volume data,
which rearranges the order of voxels to allow consecutive threads(warps) to refer to voxels with the minimum access strides. Experiments indicate that the proposed method replaces the worst
texture cache (TC) hit rate of 42% with the best TC hit rateof 93% for a 10243 -voxel volume. Thus, the average frame rateincreases by a factor of 1.6 in the proposed method comparedwith that in the previous method. Although the overhead of inplace rotation slightly decreases the frame rate from 2.0 frames
per second (fps) to 1.9 fps, this slowdown occurs only with a few
viewing directions.
key words: Cache optimization, volume rendering, in-place algorithm, GPU, CUDA.

1.

Introduction

Volume rendering [1] is a visualization technique thatcan generate a 2-D projection of 3-D volume data withan arbitrary viewing direction. For example, volume
rendering visually supports the understanding of timevarying fluid simulation [2–4] and computer-aided diagnosis in clinical fields [5, 6]. A typical visualization
technique is ray casting [7], which propagates rays fromthe viewpoint to every pixel on the screen, as shown in
Fig. 1. A pixel value can be computed by accumulating
the values of penetrated voxels on the ray at regular intervals. Ray casting can be easily parallelized because
each pixel computation is independent in terms of data
dependence. However, ray casting is a memory-boundapplication because of the low locality of references;
penetrated voxels can only be reused between neighboring rays.
Therefore, many volume renderers [8] are accelerated by a graphics processing unit (GPU) [9], which
provides not only one magnitude higher memory bandwidth over the CPU but also thousands of processing
cores within a chip. For example, NVIDIA’s Maxwellarchitecture [9] provides 336 GB/s of memory band-†

The authors are with the Graduate School of Information Science and Technology, Osaka University, Suita-shi,
565–0871, Japan.

zzVolumexVolumex

ScreenScreenViewpoint

y

y

(a)

= 270

=0

= 90

(b)

Fig. 1 Geometric relations between the volume axes and thescreen: (a) rays are propagated from the viewpoint, which canbe rotated around the volume; (b) θy represents the viewing angle around the y-axis. Angles θy = 90 and 270 are unfavorable
whereas angle θy = 0 is favorable.

width with 3072 processing cores. These rich resources
use millions of GPU threads to independently compute pixel values. Furthermore, the GPU provides
hardware-accelerated trilinear interpolation. This acceleration is realized by texture units [10] that are separated from the processing cores. Thus, interpolation forevery sampling point on rays can be oﬄoaded from theprocessing cores to the texture units. To take advantage of this technique, existing volume renderers [11–14]
typically store volume data in texture memory, whichis then accessed by the processing cores to rapidly fetch
interpolated texels through a texture cache (TC). Consequently, maximizing the TC hit rate is the key toachieving high performance for memory-bound volume
rendering on a GPU.
Previous studies have accelerated volume rendering by maximizing the TC hit rate. Sugimoto et al.
[13] presented a cache-aware volume rendering methodthat dynamically selects the best shape (i.e., width andheight) of thread blocks [10], according to the viewing
direction. Their dynamic method attempts to minimize the memory access stride for warps [10], i.e., aseries of consecutive threads that execute the same instruction every clock cycle. To achieve this, they analyzed the memory access stride, which varies according
to the direction of the neighboring voxels, i.e., the x-,y-, and z-directions, as shown in Fig. 2. For example,
their method selects a horizontally long thread block(i.e., horizontally long warps) when the x-axis of thevolume, which has the smallest stride between neighboring voxels, is rendered horizontally on the screen.
With an NVIDIA GeForce 580 GTX GPU, this minimization increases the worst frame rates by a factor ofc 200x The Institute of Electronics, Information and Communication Engineers
Copyright ⃝IEICE TRANS. ??, VOL.Exx–??, NO.xx XXXX 200x2z
32 pixels rendered bya horizontally long warp

Volumex
Voxels accessed
simultaneously bya warpScreenViewpointy
32 pixels rendered bya vertically long warp

Fig. 2 Principles of thread block shaping [13]. Thirty-two
threads in a warp are assumed to simultaneously access voxels
on a facing place (the xy-plane, in this example). The shape of
thread blocks is determined such that the memory access stride
between simultaneously accessed voxels is minimized. In thisexample, a horizontally long shape rather than a vertically longshape is selected because neighboring voxels along the x-axis havea smaller stride than those along the y-axis.

2.2 when projecting a 10243 -voxel volume onto a 10242 pixel screen. The main drawback of this method is thatthe minimization is insuﬃcient for unfavorable viewing
directions where the x-axis is parallel to the penetrating rays as shown in Fig. 1(b). In general, threads
in the same warp simultaneously access voxels on theyz-plane in this case, and this access pattern causes
relatively large strides because the simultaneously accessed voxels are not oriented along the x-axis. In fact,
the best frame rate was 2.7 times higher than the worst
frame rate.
In this study, we present a cache-aware method toaccelerate texture-based volume rendering on a GPU.
The proposed method extends a previous method [13]such that it can increase the average frame rate obtained while rotating the volume using an in-place rotation algorithm that dynamically rearranges the order of voxels to allow unfavorable viewing directions
to be interpreted as favorable. Here, an in-place algorithm is defined as an algorithm that requires O(1)memory space (except the original volume data). The
in-place procedure is activated when the viewing direction enters from a favorable (unfavorable) region to an
unfavorable (favorable) region. This data rearrangement takes place at run-time; therefore, its overhead
can drop the worst rendering performance. However,after changing the order of voxels, unfavorable directions can be processed as favorable directions rapidly;
thus, the average performance can be significantly increased. Our renderer was implemented on a GPU thatis compatible with the compute unified device architecture (CUDA) [10].
Our in-place implementation is based on a tilebased matrix transpose method [15] that utilizes (onchip) shared memory [10] for optimization of oﬀ-chip
memory access. Because the volume rotation cannot
be realized with matrix transpose, we extend this previous method for volume rotation. We show that (1) atile-based approach is useful for in-place rotation and(2) 3-D tiles rather than 2-D tiles realize eﬃcient rotation with reduced (oﬀ-chip) memory access.

The remainder of this paper is organized as follows. Related studies are discussed in Section 2. Section 3 summarizes the previous method [13], which isthe basis of the proposed cache-aware method. Section
4 describes the proposed method, and Section 5 gives
experimental results. Section 6 concludes this paper
and provides suggestions for future studies.
2.

Related Work

Weiskopf et al. [16] presented a data structure thatsplits a volume into smaller subvolumes called bricks.
Each brick comprises 43 voxels, and the bricks are oriented in diﬀerent directions. Consequently, for any
viewing direction, half the bricks are oriented along theviewing direction and the remaining bricks are oriented
in the perpendicular direction. Thus, the performancevariance for diﬀerent viewpoints is averaged. In contrast to this averaging strategy, the proposed method
accepts some overheads at specific viewing directions
to obtain the maximum performance for other viewing
directions.
A similar averaging strategy was presented by
Wang et al. [14, 17], who presented a sampling strategy called warp marching to realize constant rendering
performance while rotating the volume by exploiting
the data parallelism not only between rays but also oneach ray. In other words, a warp is responsible for a
ray so that threads in the warp can simultaneously access voxels along the ray. This is useful for unfavorable
viewing directions where the x-axis is parallel to therays. With an NVIDIA GeForce GTX TITAN GPU,
their method increased the worst frame rates from 15
frames per second (fps) to 30 fps when projecting a10243 -voxel volume onto a 5122 -pixel screen. However,the best frame rate dropped from 60 fps to 45 fps because the best performance was obtained while facing
the xy-plane. In other words, the warps responsible forthe same ray can suﬀer from large strides when accessing voxels neighboring along the z-axis.
Jönsson et al. [18] presented a data structure that
exploits inter-ray coherence by sharing cached data
called marching cache between neighboring rays. Theirmethod stores a marching cache in shared memory [10],which can be rapidly accessed as software managed
cache. However, one drawback of shared memory isa lack of hardware-accelerated interpolation. Compared with a texture-based implementation, the rendering performance decreased by roughly 71% using
trilinear interpolation. Consequently, marching caches
are useful for high-quality volume rendering that uses
a complex filter, such as the Catmull-Rom spline [19]or a B-spline [20]. A similar approach was presented
by Mensmann et al. [21], who stored the subvolume inshared memory.
In general, cache-aware methods are useful for
accelerating memory-bound applications. Zheng et

MISAKI et al.: CACHE-AWARE, IN-PLACE ROTATION METHOD FOR TEXTURE-BASED VOLUME RENDERING

3
z

z
z

6

x

y

1

2
y

y

(a)
Fig. 3 Memory access strides between neighboring voxels in a
3-D texture. The access stride ratio of the x-, y-, and z-axis is1:2:6 [13] if the texture is stored in the Morton order [25].

al. [22] presented a cache-aware memory-schedulingscheme for cone beam backprojection [23, 24], which isan inverse volume rendering problem. In other words,
a series of 2-D projections are backprojected into 3D space to reconstruct volume data. To increase thecache hit rate during backprojection, their method rotates the volume such that warps can access a small region on the projections. This idea is similar to our data
rearrangement strategy, but we employ an in-place rotation algorithm for an inverse problem wherein warps
primarily access a small region in the volume (3-D texture) rather than the projections (2-D textures).
3.

x

x

Principles of Warp-Level Cache Optimization

The basic idea of the previous method [13] is the warplevel cache optimization for volume rendering. This
warp-level strategy is eﬀective on the GPU because
threads in the same warp issue the same memory transaction simultaneously. To achieve this warp-level optimization, Sugimoto et al. [13] analyzed the memory access strides along the volume axes. They focused on the
Morton order [25] and showed that the memory stride
ratio along the x-, y-, and z-axis of a 3-D texture isaround 1:2:6, as shown in Fig. 3. This non-uniform ratio indicates that the rendering performance can vary
according to the viewing direction. In fact, rays canpropagate along the x-, y-, or z-axis. We assume thata thread is responsible for computing a pixel value on
the screen. Here, let θy be the viewing angle aroundthe y-axis, as shown in Fig. 1.
Voxels are sampled at regular intervals from theviewpoint; thus, a warp simultaneously accesses voxels
on the surface of a sphere. For simplicity, the previous
method assumes that this spherical surface can be approximated with a plane. Under this approximation, awarp simultaneously accesses voxels on a plane that areparallel to the screen (Fig. 2). For such a facing plane,
there are three cases, i.e., the xy-, yz-, and xz-plane.
Moreover, there are two variations for each case, i.e.,the vertical (horizontal) axis of the plane has smaller
access strides than the horizontal (vertical) axis. Thus,
any viewing direction can be classified into these 3 × 2
groups.

(b)

Fig. 4 Favorable and unfavorable viewing directions: (a) favorable viewing direction θy = 0; (b) unfavorable viewing direction
θy = 90 (the facing plane for the former is the xy-plane and thatfor the latter is the yz-plane, which does not include the x-axis).

Given a viewing direction, the previous method
first identifies the facing plane and selects an axis fromthe two axes that constitute the facing plane such thatthe selected axis has smaller access strides. For example, the facing plane in Fig. 4(a), where θy = 0,
is the xy-plane and the plane axis with smaller accessstrides is the x-axis, which can be seen horizontally on
the screen. In this case, the previous method selects
horizontally long thread blocks (i.e., horizontally longwarps) to allow warps to access voxels with the minimum access strides. In contrast, for the viewing direction presented in Fig. 4(b), where θy = 90, the method
selects vertically long thread blocks because the vertical axis has smaller access strides than the horizontal
axis (i.e., the y-axis rather than the z-axis).
Although the previous method selects the best
thread block shape for the viewing direction, this selection is insuﬃcient when the facing plane comprises
the y- and z-axis. Such an unfavorable viewing direction is shown in Fig. 4(b). In this example, the x-axis
is parallel to the viewing direction; thus, warps fail toaccess voxels with the minimum access strides. To address this issue, a ray must be assigned to a warp (i.e.,
multiple threads), as proposed by Wang et al. [14]. Analternative solution is to rearrange the order of voxels,
which we present in the next section.
4.

Proposed Method

The basic idea of the proposed method is to rearrange
the order of voxels to allow unfavorable viewing directions to be interpreted as favorable. For example, theviewing direction in Fig. 4(b) can be considered as thatin Fig. 4(a) if the volume is rotated 90◦ around the yaxis prior to rendering. After this data rearrangement,
the facing plane changes from the yz-plane to the xyplane, which can be rapidly accessed by warps with theminimum access strides.
Note that this rotation is intended to change theorder of voxels in the texture rather than update theview matrix. Consequently, a naive method might distinguish the output volume from the input volume during rotation. However, this method consumes twice thevideo memory, which reduces the maximum volume size

IEICE TRANS. ??, VOL.Exx–??, NO.xx XXXX 200x

4(b)Fig. 6 Tiled approach for rotating an xz-slice: (a) before rotation and (b) after rotation. Four data-dependent tiles T0,0 , T0,3 ,T3,3 , and T3,0 , each containing t × t voxels, are rotated with acyclic permutation processed by a thread block. Reads and writes
to a tile are parallelized by t × t threads within the thread block.

Invoke the in-placerotation kernel

End

Fig. 5 Flowchart of the proposed volume rendering method.
The “Start” and “End” states are associated with a CPU rendering process.

that can be rapidly rendered without swapping datafrom the video to the main memory. Therefore, theproposed method realizes in-place rotation of the volume, which integrates the input volume with the output
volume to limit memory consumption.
Figure 5 shows an overview of the proposed cacheaware volume rendering method. The proposed method
first identifies whether the viewing direction is on aboundary of an unfavorable region. If a viewing angle
θy is within the range of [45, 135) or [225, 315), the viewing direction is unfavorable; otherwise, it is favorable.
After this identification process, the proposed method
invokes an in-place rotation kernel if the viewing direction enters from a favorable (unfavorable) region toan unfavorable (favorable) region. Finally, the method
invokes a rendering kernel to produce the final image
for the viewing direction. This rendering kernel is invoked with an appropriate shape of the thread blocks toperform warp-level cache optimization [13]. Note thatthe in-place rotation kernel must be separated from therendering kernel because barrier synchronization is required between the two kernels.
4.1

x

T3,3

Yes

T3,0

Invoke therendering kernel

x
(a)

T0,3

T0,3

Update theviewing direction

No

x
T0,0

No

Boundary of
unfavorable region?

T3,3

z
Yes

T0,0

T3,0

Data padding
if needed
Finishrendering?

z

z

Start

In-place rotation

Similar to the matrix transpose method [15], the proposed method is based on a tile-based approach thatfacilitates in-place rotation of the volume. We assume
that xz-slices of the volume are isometric, i.e., the volume size Nx along the x-axis is the same as Nz along
the z-axis. For non-isometric data, dummy transparent
voxels are padded to make the data isometric prior torendering. This padding operation can be processed as
an initialization phase, which occurs only once after thedataset to be rendered has been loaded. The proposed
method avoids rendering such dummy voxels using an
empty space skipping technique [26].

Here, let N (= Nx = Nz ) be the volume size. Let
t×t be the thread block size of the in-place rotation kernel. Figure 6 shows a data arrangement example thatrotates an xz-slice 90◦ around the y-axis. As shown,the proposed method partitions the xz-slice into squaretiles of t×t voxels. Note that, for simplicity, we describe
our method using 2-D tiles of t × t voxels; however, ourimplementation uses 3-D tiles of t × t × k voxels, wherek (> 0) represents the tile size along the y-axis. Here,let Ti,j be a tile of the i-th row and the j-th column,
where 0 ≤ i, j < N/t. We assume that N is a multipleof t for simplicity.
A series of data-dependent tiles is then assigned toa thread block. For example, the series of tiles T0,0 ,T0,3 , T3,3 , and T3,0 in Fig. 6 must be processed by thesame thread block because the i-th tile moves to the((i + 1) mod 4)-th tile, where 0 ≤ i < 4. In otherwords, the rotation of these four tiles can be considered a cyclic permutation, which must be processed insequence. However, voxels in a tile can be moved inparallel so that tile movement can be parallelized by athread block. More formally, a series of data-dependent
tiles can be given by Ti,j , TN/t−j−1,i , TN/t−i−1,N/t−j−1
and Tj,N/t−i−1 , where 0 ≤ i, j < N/t.
Algorithm 1 shows the pseudocode of the in-placerotation kernel, which is processed by GPU threads ina single-instruction multiple-thread manner [10]. We
assume that the local index of a thread is given by(o1 , o2 ), where 0 ≤ o1 , o2 < t. Given the original volume V and the global coordinates (x, y, z) of the local
origin of the first tile Ti,j , the GPU threads output therotated volume using shared memory. To do this, eachthread iterates the following steps four times:
1. Tile load (lines 1 and 7). Each thread accesses
video memory to load a voxel of the current tile ina register. Threads in the same warp simultaneously access contiguous memory addresses during
this step.
2. Tile rotation (line 9). Loaded voxels are written
into shared memory such that the shared memory
holds the rotated tile.

MISAKI et al.: CACHE-AWARE, IN-PLACE ROTATION METHOD FOR TEXTURE-BASED VOLUME RENDERING

5
Table 1

Algorithm 1 In-place rotation kernel
Input: Original volume V
Input: The global coordinates (x, y, z) of the local origin of thefirst tile Ti,jInput: The local index (o1 , o2 ) of a thread, where 0 ≤ o1 , o2 < t
Output: Rotated volume V
1: r1 ← V [x + o1 ][y][z + o2 ]; ▷ load a voxel of the first tile in aregister
2: for i = 0 to 3 do
▷ for each data-dependent tile3:
tmp ← z;▷ compute the coordinates of the next tile4:z ← x;
5:x ← N/t − tmp − 1;
6:
if i ̸= 3 then7:
r2 ← V [x + o1 ][y][z + o2 ]; ▷ load a voxel of the next
tile in a register8:
end if
9:
shared[o2 ][t − o1 − 1] ← r1 ;▷ store the rotated tile inthe shared memory
10:
syncthreads();
11:
V [x + o1 ][y][z + o2 ] ← shared[o1 ][o2 ]; ▷ store the rotatedtile in the video memory12:
r1 ← r2 ;▷ prepare for the next iteration
13:
if i ̸= 3 then14:
syncthreads();
15:
end if
16: end for3. Tile store (line 11). Each thread stores the rotated tile from the shared memory to the video
memory. Similar to the first step, threads in thesame warp simultaneously access contiguous memory addresses during this step.
Note that t-strided memory accesses are required
to rotate a tile by threads. Such noncontiguous memory
accesses cannot be handled eﬃciently with (oﬀ-chip)video memory because memory transactions are performed per warp [10]. Therefore, our kernel uses sharedmemory (i.e., on-chip memory) to avoid strided accesses
to the video memory. However, synchronization is required prior to loading shared data, which can be stored
asynchronously by other threads.
As mentioned above, our method can be regarded
as an extension of the matrix transpose method [15].
However, this method cannot be directly used for volume rotation, which requires cyclic permutations instead of simple swap operations. Each permutationconsists of a sequence of shift operations applied to fourtiles as illustrated in Fig. 6. Consequently, our methodassigns these four tiles to a thread block whereas thematrix transpose method assigns two tiles to a thread
block.
4.2

Issues with writable textures

As mentioned in Section 1, a volume is typicallystored in a 3-D texture to take advantage of hardwareaccelerated trilinear interpolation.
However, theCUDA currently prohibits GPU threads from writing
to textures [10]. Consequently, the texture data cannot
be directly rotated by GPU threads.

Experimental environment.

Item
CPU
Main memory capacity
GPU
Video memory capacity
OS
CUDA version
Driver version

Specification
Intel Core i7-3770K 3.5 GHz
16 GB
NVIDIA GeForce GTX TITAN X
12 GB
Windows 7
7.0
353.30To address this issue, we adopt a two-way approachto access the volume data. This approach uses texture
memory and surface memory [10], which is writablefrom GPU threads but does not provide hardwareaccelerated interpolation. In other words, the CUDA
array [10] allocated for the volume is bound to a texture object and a surface object. The in-place rotation kernel, which never requires interpolated values,
then writes the volume data through the surface object. In addition, the rendering kernel accesses the volume data through the texture object to take advantageof hardware-based interpolation.
Notice here that only a single CUDA array is required to store the entire volume with padded area.
The presented two-way approach never requires a copyof the volume data. Consequently, the overall memory
usage can be approximated by the number of voxels in
the original volume and padded region.
5.

Experimental Results

We evaluated the proposed cache-aware method bycomparing it to a previous method [13] and a naive
method distributed as CUDA sample code [10]. Comparisons were performed in terms of the TC hit rate andthe frame rate. Similar to the previous method [13], acolor map table, which associates color and opacity values with each voxel, was stored in shared memory to
avoid perturbation of TC behavior. All methods activated early ray termination [26] to skip accumulationthat did not significantly contribute to the pixel.
According to [27], at least four over-samplings are
required to reconstruct the ray integral with suﬃcientaccuracy. Consequently, the sampling distance along
the rays was set to max(Nx , Ny , Nz )/4, so that the renderer used at most 4×max(Nx , Ny , Nz ) samples per ray,starting from the first penetrated voxel in the volume.
We determined experimentally that t = 32 achieves thehighest rendering performance on the deployed GPU.
Details of the experimental environment are listed inTable 1.
We used four volume datasets: head magnetic resonance imaging (MRI), Hazelnut, Flower, and Porsche,
as shown in Fig. 7. The head MRI and Porsche datasets
are non-isometric; the Hazelnut and Flower datasets areisometric. The isometric datasets diﬀer with respect tovolume size. All voxels had an 8-bit scalar value, andthe volume was stored in a 3-D texture in unsigned

IEICE TRANS. ??, VOL.Exx–??, NO.xx XXXX 200x

6

(a)

(b)

(c)

(d)

1.0

1.0

0.8

0.8
Cache hit rate

Cache hit rate

Fig. 7 Experimental datasets rendered on a 1024 × 1024-pixel screen: (a) head MRI (256 × 256 × 124 voxel); (b) Hazelnut(512 × 512 × 512 voxel); (c) Flower (1024 × 1024 × 1024 voxel); and (d) Porsche (559 × 1023 × 347 voxel). Head MRI
and Porsche datasets were obtained from http://www.gris.uni-tuebingen.de/edu/areas/scivis/volren/datasets/new.html andhttp://www9.informatik.uni-erlangen.de/External/vollib/, respectively. Hazelnut and Flower datasets were obtained fromhttp://www.ifi.uzh.ch/vmml/research/datasets.html.

0.6
0.4
0.2

0

45

90 135 180 225 270
Viewing angle (degree)

0.4
0.2

Our method
Previous method
CUDA sample

0.0

0.6

Our method
Previous method
CUDA sample

0.0
315

360

0

45

1.0

1.0

0.8

0.8

0.6
0.4
0.2

0

45

90 135 180 225 270
Viewing angle (degree)

360

315

360

0.6
0.4
0.2

Our method
Previous method
CUDA sample

0.0

315

(b)

Cache hit rate

Cache hit rate(a)

90 135 180 225 270
Viewing angle (degree)

Our method
Previous method
CUDA sample

0.0
315

360

(c)

0

45

90 135 180 225 270
Viewing angle (degree)(d)Fig. 8 TC hit rates of the proposed method, previous method [13], and CUDA sample
code [10] for (a) head MRI, (b) Hazelnut, (c) Flower, and (d) Porsche datasets.

char format. The screen size was set to 10242 pixels
during the experiments. A pair of movie files capturingreal-time visualization of the Flower dataset is available from http://www-hagi.ist.osaka-u.ac.jp/
research/movie/misaki-{proposed,previous}.mp4.
5.1

Rendering performance

Figure 8 shows the TC hit rates of the rendering kernelmeasured while rotating the viewing direction aroundthe y-axis, which significantly varies the rendering performance. As for the rotation around the x-axis or
z-axis, where the volume is always rendered from favorable viewing directions, the rendering performance
was kept high [13]. Consequently, these results are notpresented here. The TC hit rates were obtained using
the CUDA toolkit nvprof profiling tool [28].
As shown in Fig. 8(a), there was no significant gap
between the TC hit rate of the proposed method andthat of the previous method when the smallest head

MISAKI et al.: CACHE-AWARE, IN-PLACE ROTATION METHOD FOR TEXTURE-BASED VOLUME RENDERING

7

45

Our method
Previous method
CUDA sample

Our method
Previous method
CUDA sample

40
35

150

Frame rate (fps)Frame rate (fps)

200

100

50

30
25
20
15
10
5

0

0
0

45

90 135 180 225 270
Viewing angle (degree)

315

360

0

45

(a)
30

120

360

315

360

Our method
Previous method
CUDA sample

100
Frame rate (fps)Frame rate (fps)

315

(b)

Our method
Previous method
CUDA sample

25

90 135 180 225 270
Viewing angle (degree)

20
15
10
5

80
60
40
20

0

0
0

45

90 135 180 225 270
Viewing angle (degree)

315

360

(c)

0

45

90 135 180 225 270
Viewing angle (degree)(d)Fig. 9 Frame rates of the proposed method, previous method [13], and CUDA sample
code [10] for (a) head MRI, (b) Hazelnut, (c) Flower, and (d) Porsche datasets.

MRI dataset was rendered. Obviously, such a small
dataset can be rendered eﬃciently with many cache hits
because the largest stride between neighboring voxels
is given by Nx × Ny = 65, 536. Consequently, 64 KB ofmemory is suﬃcient to store an xy-slice of the volume.
On the other hand, the deployed GPU has 1152 KB of
L1 texture cache; the GeForce GTX TITAN X has 24
streaming multiprocessors, each having 2 × 24 KB of
L1 texture cache [29, 30].
However, as shown in Figs. 8(b) and 8(c), the gap
between the proposed method and the previous method
increased for larger datasets that comprised at least
5123 voxels. Both large datasets show similar behavior
in terms of the TC hit rate. The proposed method
significantly increased the TC hit rates for unfavorable viewing directions within the ranges of [45, 135)and [225, 315). For example, with the Flower dataset(Fig. 8(c)), the worst TC hit rates at θy = 90 and
270 increased from 42% to 93%. Moreover, this increased TC hit rate was almost equal to the best hit rateachieved by the previous method, i.e., 93% at θy = 0
and 180. In other words, the proposed method successfully translated unfavorable viewing directions into
favorable viewing directions.
A similar periodical behavior was observed with
the non-isometric Porsche dataset. Note that the proposed method decreased the TC hit rate at four viewing

regions: [46, 59), [126, 136), [226, 236), and [302, 316)
(Fig. 8(d)). This side eﬀect was due to our rotation approach, which exchanges Nx with Nz . In other words,
this non-isometric dataset satisfies Nx > Nz (559 >
347) before rotation; however, rotation changes this to
Nx < Nz . Thus, the number of non-dummy voxels
neighboring along the x-axis was reduced after rotation. Therefore, the TC hit rate decreased within theabovementioned four viewing regions. For the same
reason, the peaks at θy = 90 and 270 are more pronounced with the proposed method than those of theprevious method and the CUDA sample code at θy = 0
and 180.
Next, we analyzed the rendering performance,
which includes the in-place rotation kernel run-timeoverhead. Figure 9 shows the frame rates observed
while rotating the viewing direction around the y-axis.
With the smallest dataset, i.e., head MRI, there wasno significant gap between the frame rates of the proposed method and those of the previous method (Fig.
9(a)). However, similar to the TC hit rate behavior, theproposed method successfully increased frame rates forlarge datasets, as shown in Figs. 9(b) and 9(c). Thus,
maximizing the TC hit rate is critical to maximizing
rendering performance for large datasets.
Compared with the previous method, the worst
frame rate for the Flower dataset was reduced from

IEICE TRANS. ??, VOL.Exx–??, NO.xx XXXX 200x

160

35

140

30

120

25

Frame rate (fps)Frame rate (fps)

8

100
80
60
40

15
10
5

20
0

20

Our method

Previous method

0

CUDA sample

Our method

(a)

12

Frame rate (fps)Frame rate (fps)

14
10
8
6
4
2
Our method

Previous method

CUDA sample

(b)

16

0

Previous method

CUDA sample

(c)

90
80
70
60
50
40
30
20
10
0

Our method

Previous method

CUDA sample(d)Fig. 10 Box plots of frame rates of the proposed method, previous method [13], and CUDA sample code [10] for (a)head MRI, (b) Hazelnut, (c) Flower, and (d) Porsche datasets. Maximum, 75% percentile, 50% percentile (median), 25%
percentile, and minimum rates are presented.

approximately 2.0 fps to 1.9 fps with the proposed
method. These dropped frame rates were observed
when θy = 45, 135, 225, or 315, where the in-place rotation kernel was invoked to rearrange the order of voxels. Similarly, the overhead of the rotation kernel was
ignorable for the Hazelnut dataset as shown at θy = 45,
135, 225, and 315 in Fig. 10(b). In contrast, the overhead slightly decreased the frame rate for non-isometric
datasets: the head MRI and Porsche datasets as canbe seen in Figs. 10(a) and 10(d), respectively. Thiswas due to padded voxels. Although these transparent
voxels were skipped during rendering, the frame rate
slightly dropped at these angles. However, the invocation overhead allowed unfavorable regions [45, 135) and[225, 315) to be processed as rapidly as the remaining
favorable regions. Consequently, compared to the previous method, the execution times spent from θy = 0 to
θy = 359 were reduced by 13%, 38%, and 18% for theHazelnut, Flower, and Porsche datasets, respectively;
as for the smallest MRI Head dataset, the execution
time increased by 5%. Thus, the rendering performance
was primarily dominated by texture access.
Box plots of the measured frame rates are shown in
Fig. 10. As volume size increased, the performance gainof the proposed method was evident. The most significant improvement can be seen in the median of themeasured performance, which increased from 2.5 fps to
9.8 fps for the Flower dataset (Fig. 10(c)). The average frame rate also increased from 5.7 fps to 9.3 fps. A
similar significant increase was observed with the nonisometric Porsche dataset. Consequently, the proposed
method can increase the average rendering performanceat the cost of four viewing directions, where in-place ro-tation occurs (θy = 45, 135, 225, and 315).
Note that with the previous method, there was a
relatively large gap between the median and average
frame rates; the gap was 3.2 (= 5.7 − 2.5) fps withthe previous method, whereas it was 0.5 (= 9.8 − 9.3)
fps with the proposed method. This large gap implies
that the previous method suﬀers from many low frame
rates below the median. Such low frame rates wereprimarily observed when the volume was rendered withunfavorable viewing directions. Similar behavior wasobserved with the CUDA sample code.
As presented above, our method increased theframe rates for many viewing directions at the expenseof few directions (slowdowns). Because volume rendering is designed for visualization with an arbitrary
viewing direction, we think that these slowdowns are
acceptable if the rendering performance increases formany viewing directions. In fact, when rotating the
Flower dataset, we faced little (4%) slowdowns only for
1.1% of viewing directions, whereas the frame rates increased by ×3.5 for 50% of viewing directions. No significant performance improvement (>2%) was observed
for the remaining (49%) directions. In this case, the influence of the worst frame rate was limited, and thus,
demonstrating the practicality of our method.
5.2 Eﬃciency of in-place rotation
We also analyzed the eﬃciency of the in-place rotation
kernel because rendering performance can be limited by
its run-time overhead. As shown in Fig. 9, the worst
frame rates at θy = 45, 135, 225, and 315 decreasedslightly with the proposed method.

MISAKI et al.: CACHE-AWARE, IN-PLACE ROTATION METHOD FOR TEXTURE-BASED VOLUME RENDERING

Execution time (ms)

9
50
45
40
35
30
25
20
15
10
5
0

47

6.
23.2
14.9

14.1

14.2
N/A

1

2

4
8
16
k: tile size along the y-axis

32

Fig. 11 Execution time of the in-place rotation kernel with
diﬀerent tile sizes t × t × k, where t = 32 for Flower dataset. Dueto the limited size of thread blocks (i.e., the maximum numberof threads in a thread block), execution with k = 32 resulted ina failure.

Figure 11 shows the execution times of the in-placerotation kernel with diﬀerent tile sizes t × t × k, wheret = 32 and 1 ≤ k ≤ 32. The execution time was minimized when k = 16. With k = 16, the eﬀective memory
throughput reached 231 GB/s, which was 93% of theeﬀective oﬀ-chip memory bandwidth (248 GB/s) measured by the CUDA software development kit bandwidthTest program. Because the performance of inplace rotation is limited by oﬀ-chip memory access, we
believe that the achieved performance is the best resultwith the deployed GPU.
The execution times at k = 1 and k = 2 were,respectively, 4.1 and 2.2 times longer than the shortest execution time (k = 16). According to the nvprof
profiling tool, execution with k = 1 (k = 2) incurred
four (two) times more oﬀ-chip memory accesses thanwith k = 4. In contrast, oﬀ-chip memory access wasnearly constant when k ≥ 4. This behavior indicatesthe internal structure of surface memory, which uses
byte addressing [10]. Owing to this byte addressing,surface memory might be optimized for 3-D surface objects of depth of at least four. Thus, 3-D tiles (k > 1)
rather than 2-D tiles (k = 1) realize eﬃcient rotation
with reduced oﬀ-chip memory access.
Finally, the execution time of in-place rotation wasmuch shorter than that of data transfer between the
CPU and GPU. We found that it took at least 270 ms
to swap out the Flower dataset from the video memory to the main memory, whereas the in-place rotation
took 9.3 ms, as shown in Fig. 11. Therefore, without the proposed in-place algorithm, the worst rendering performance dropped from 1.9 fps to 1.3 fps dueto the data transfer required to swap out the volumedata from the video memory to the main memory. In
contrast, our in-place algorithm avoided dropping theworst rendering performance because it took 9.3 ms tocomplete the rotation kernel; the rotation kernel canbe processed at 108 (=1000/9.3) fps even though it isinvoked after every frame.

Conclusion

In this study, we have presented a cache-aware method
for accelerating texture-based volume rendering on aGPU. The proposed method maximizes the TC hit rateby changing the order of voxels at specific viewing directions. This data rearrangement is performed by anin-place algorithm, which requires O(1) memory space(except the volume data). Thus, a large volume can be
rapidly rendered without swapping data from the video
memory to main memory. We integrated the proposed
method into a previous cache-aware method [13] thatis capable of selecting the best thread organization according to the viewing direction.
We compared the proposed method with the previous method and the CUDA sample code. The experimental results show that the proposed method successfully increased the worst TC hit rate (42%) to 93%.
Accordingly, the average frame rate increased from 5.7
fps to 9.3 fps, achieving a 1.6× speedup over the previous method. We also found that the performance ofthe proposed in-place rotation kernel was limited byoﬀ-chip memory bandwidth.
In the future, we plan to develop a sophisticatedmechanism to dynamically activate the in-place rotation kernel according to the movement history of theviewing directions.
Acknowledgments
This study was supported in part by the Japan Society for the Promotion of Science KAKENHI Grant
Numbers 15K12008, 15H01687, and 16H02801, and theJapan Science and Technology Agency CREST program, “An Evolutionary Approach to Construction ofa Software Development Environment for MassivelyParallel Computing Systems.” We are also grateful tothe anonymous reviewers for their valuable comments.
References[1] R.A. Drebin, L. Carpenter, and P. Hanrahan, “Volume
rendering,” Computer Graphics (Proc. SIGGRAPH’88),
vol.22, no.3, pp.65–74, Aug. 1988.
[2] H.W. Shen, L.J. Chiang, and K.L. Ma, “A fast volume rendering algorithm for time-varying fields using a time-spacepartitioning (TSP) tree,” Proc. 10th IEEE Visualization
Conf. (VIS’99), pp.371–377, Oct. 1999.
[3] E.B. Lum, K.L. Ma, and J. Clyne, “Texture hardware assisted rendering of time-varying volume data,” Proc. 12th
IEEE Visualization Conf. (VIS’01), pp.263–270, Oct. 2001.
[4] D. Nagayasu, F. Ino, and K. Hagihara, “A decompression
pipeline for accelerating out-of-core volume rendering oftime-varying data,” Computers and Graphics, vol.32, no.3,
pp.350–362, June 2008.
[5] P.S. Calhoun, B.S. Kuszyk, D.G. Heath, J.C. Carley, and
M.E.K. Fishman, “Three-dimensional volume rendering ofspiral ct data: Theory and method,” Radiographics, vol.19,
no.3, pp.745–764, May 1999.

IEICE TRANS. ??, VOL.Exx–??, NO.xx XXXX 200x

10[6] A. Takeuchi, F. Ino, and K. Hagihara, “An improved
binary-swap compositing for sort-last parallel rendering on
distributed memory multiprocessors,” Parallel Computing,vol.29, no.11/12, pp.1745–1762, Nov. 2003.
[7] M. Levoy, “Display of surfaces from volume data,” IEEE
Computer Graphics and Applications, vol.8, no.3, pp.29–
37, May 1988.
[8] J. Beyer, M. Hadwiger, and H. Pfister, “State-of-the-art
in GPU-based large-scale volume visualization,” Computer
Graphics Forum, vol.xx, no.xx, pp.xx–xx, Sept. 2015.
[9] NVIDIA Corporation, “NVIDIA GeForce GTX 980,” Nov.
2014.
[10] NVIDIA Corporation, “CUDA C Programming Guide Version 7.0,” March 2015.
[11] I. Boada, I. Navazo, and R. Scopigno, “Multiresolution volume visualization with a texture-based octree,” The VisualComputer, vol.17, no.3, pp.185–197, May 2001.
[12] J. Krüger and R. Westermann, “Acceleration techniques for
GPU-based volume rendering,” Proc. 14th IEEE Visualization Conf. (VIS’03), pp.287–292, Oct. 2003.
[13] Y. Sugimoto, F. Ino, and K. Hagihara, “Improving cache
locality for GPU-based volume rendering,” Parallel Computing, vol.40, no.5/6, pp.59–69, May 2014.
[14] J. Wang, F. Yang, and Y. Cao, “Cache-aware samplingstrategies for texture-based ray casting on GPU,” Proc.
4th IEEE Symp. Large Data Analysis and Visualization(LDAV’14), pp.19–26, Nov. 2014.
[15] M. Harris, “Optimizing CUDA,” SC’07 Tutorial, Nov.
2007.
http://gpgpu.org/static/sc2007/SC07_CUDA_5_
Optimization_Harris.pdf.
[16] D. Weiskopf, M. Weiler, and T. Ertl, “Maintaining constantframe rates in 3D texture-based volume rendering,” Proc.
21st Computer Graphics Int’l (CGI’04), pp.604–607, June
2004.
[17] J. Wang, F. Yang, and Y. Cao, “Cache-aware samplingstrategies for texture-based ray casting on GPU,” Proc.
8th IEEE Pacific Visualization Symp. (PacificVis’15), April
2015.
[18] D. Jönsson, P. Ganestam, A. Ynnerman, M. Doggett, andT. Ropinski, “Explicit cache management for volume raycasting on parallel architectures,” Proc. 12th Eurographics
Symp. Parallel Graphics and Visualization (EGPGV’12),
pp.31–40, May 2012.
[19] E. Catmull and R. Rom, A Class of Local Interpolating
Splines, pp.317–326, Academic Press, New York, NY, 1974.
[20] S. Lee, G. Wolberg, and S.Y. Shin, “Scattered data interpolation with multilevel B-splines,” IEEE Trans. Visualization and Computer Graphics, vol.3, no.3, pp.228–244, July
1997.
[21] J. Mensmann, T. Ropinski, and K. Hinrichs, “An advanced
volume raycasting technique using GPU stream processing,” Proc. 5th Int’l Conf. Computer Graphics Theory andApplications (GRAPP’10), pp.190–198, May 2010.
[22] Z. Zheng and K. Mueller, “Cache-aware GPU memory
scheduling scheme for CT back-projection,” Proc. Nuclear
Science Symp. and Medical Imaging Conf. (NSS/MIC’10),
pp.2248–2251, Nov. 2010.
[23] Y. Okitsu, F. Ino, and K. Hagihara, “High-performancecone beam reconstruction using CUDA compatible GPUs,”
Parallel Computing, vol.36, no.2/3, pp.129–141, Feb. 2010.
[24] Y. Lu, F. Ino, and K. Hagihara, “Cache-aware GPU optimization for out-of-core cone beam CT reconstruction ofhigh-resolution volumes,” IEICE Trans. Information and
Systems, vol.E99-D, no.12, pp.xx–xx, Dec. 2016.
[25] G.M. Morton, “A computer oriented geodetic data base anda new technique in file sequencing,” tech. rep., IBM Ltd,
Ottawa, Ontario, Aug. 1966.

[26] M. Levoy, “Eﬃcient ray tracing of volume data,” ACM
Trans. Graphics, vol.9, no.3, pp.245–261, July 1990.
[27] S. Roettger, S. Guthe, D. Weiskopf, T. Ertl, and
W. Strasser, “Smart hardware-accelerated volume rendering,” Proc. 5th Eurographics/IEEE TCVG Symp. Visualization (VisSym’03), pp.231–238, May 2003.
[28] NVIDIA Corporation, “Profiler User’s Guide Version 7.0,”
March 2015.
http://docs.nvidia.com/cuda/pdf/CUDA_
Profiler_Users_Guide.pdf.
[29] X. Mei, K. Zhao, C. Liu, and X. Chu, “Benchmarking thememory hierarchy of modern GPUs,” Proc. 11th IFIP Int’l
Conf. Network and Parallel Computing (NPC’14), pp.144–
156, Sept. 2014.
[30] X. Mei, K. Zhao, C. Liu, and X. Chu, “Dissecting GPU
memory hierarchy through microbenchmarking,” IEEE
Trans. Parallel and Distributed Systems, http://arxiv.
org/pdf/1509.02308.pdf.

received the B.E. deYuji MISAKI
gree in information and computer sciencesfrom Osaka University, Osaka, Japan, in
2015. He is currently working toward the
M.E. degree at Osaka University. His current research interests include computer
graphics and high performance computing.

received the B.E.,
Fumihiko INO
M.E., and Ph.D. degrees in informationand computer sciences from Osaka University, Osaka, Japan, in 1998, 2000, and2004, respectively. He is currently an Associate Professor in the Graduate Schoolof Information Science and Technology at
Osaka University. His research interests
include parallel and distributed systems,
software development tools, and performance evaluation.

received the
Kenichi HAGIHARA
B.E., M.E., and Ph.D. degrees in information and computer sciences from Osaka University, Osaka, Japan, in 1974,
1976, and 1979, respectively. From 1994
to 2002, he was a Professor in the Department of Informatics and Mathematical Science, Graduate School of Engineering Science, Osaka University. Since 2002,
he has been a Professor in the Departmentof Computer Science, Graduate School ofInformation Science and Technology, Osaka University. From1992 to 1993, he was a Visiting Researcher at the University of
Maryland. His research interests include the fundamentals andpractical application of parallel processing.