15th Int Symp on Applications of Laser Techniques to Fluid Mechanics
Lisbon, Portugal, 05-08 July, 2010

Application of High Performance Computing on Volumetric Velocimetry
Processing
Thomas NonnDantec Dynamics A/S, DK-2740 Skovlunde, Denmark
Thomas.nonn@dantecdynamics.com

Abstract This paper examines the use of graphical processing units (GPU), a technology that has received
considerable fanfare in today’s scientific community, in processing and analysing data for volumetricvelocimetry applications. GPU’s are highly parallel devices that while lacking the data management
characteristics and flexibilities of CPU’s can execute parallel tasks at a far superior rate. CPU’s found inpersonal computers today typically contain 1 to 8 cores. GPU’s, on the other hand, can be found with up to512 cores, with the additional possibility of connecting several (up to four in a single pc) of these devices
together. Floating-point processing is often measured as FLOPS (floating-point operations per second). The
throughput of the fastest CPU is on the order of 30-50 GFLOPS, however, a single off-the-shelf GPU has
been measured to 1 Teraflop. The difference between floating-point and double-precision performance canbe considerable, therefore any algorithm must leverage performance with accuracy. GPU implementations ofvolumetric applications in general experience a 30-50-fold improvement over CPU implementations.

1. Introduction
Since the advent of Particle Image Velocimetry (PIV) in experimental fluids measurements there
has been a steady and sustained incline in the throughput capability and resolution of hardware
devices (i.e.; CMOS cameras) needed to acquire and transfer the copious amounts of image data.
With the introduction of tomographic volumetric measurement techniques the amount of datasuddenly increased by an order of magnitude. This demand is even further compounded by theinclusion of time-resolved volumetric techniques. While the development of hardware paces
reasonably well with the acquisition demand placed by current experiments, the ability ofcomputers and current algorithms to process and further reduce the data within a reasonable period
has fallen dramatically behind. Take for example a 3-camera tomographic application running at
one kilohertz for one second. If the cameras are running at full resolution, say 1K x 1K, the total
amount of raw data collected in one second will exceed (assuming 8-bit resolution) 3 gigabytes.
Processing this data on high performance computers will often take many hours or even days tocomplete.
Particle tracking techniques in themselves do not apply very well to large scale multiprocessor
architectures due to the inherent nature of tracking particle histories. Typically a large queue of
detected particles over a time T are recorded, and then in very much a serialised manner tracks aredetermined for those particles whose positions most likely fulfil a set of criteria: particle prehistory,particle distance, trajectory constraints, quality and so forth. A considerable amount of data sortingis necessary to reduce repetition. Some performance improvement is still possible throughmultithreading, though true parallel processing is not possible.
Least-Squares Matching (LSM) or Least-Squares Tracking (3D-LST), on the other hand, is ideal forparallel computation since the problem is dissected into many smaller, regularly sized cuboids thatcan be independently evaluated. Originally applied by Maas et al [4] for tracking structures in LIF(laser induced fluorescence) applications, LSM applies deformation and rotation to cuboids in an15th Int Symp on Applications of Laser Techniques to Fluid Mechanics
Lisbon, Portugal, 05-08 July, 2010

iterative fashion, using a least-squares adjustment procedure to determine the affine transformation
between a gray cuboid at time t1 and t2. The resulting affine transformation contains 12 parameters:
for translation, scale, rotation and shear.
x 2 � a 0 � a1 x1 � a 2 y1 � a 3 z1y 2 � b0 � b1 x1 � b2 y1 � b3 z1(1)z 2 � c0 � c1 x1 � c 2 y1 � c 3 z1The method requires that for two time-steps the cuboid g2 at time t2 can be represented as an affine
transformation of the cuboid g1 at time t1 and a noise expression:
g1(x,y, z)−e(x,y, z) = g2(x,y, z)(2)

The linearized observation equation can then be defined as:
g1(x,y, z)−e(x,y, z) = g02(x,y, z)+gxda0 +gxx02da1 +gxy02da2 +gxz02da3+
gydb0 +gyx02db1 +gyy02db2 +gyz02db3+gzdc0 +gzx02dc1 +gzy02dc2 +gzz02dc3(3)

2. Graphical Processing UnitsOriginally designed for games and visualisation applications, graphical processing units (GPU’s)
are recognised today as necessary tools for the processing of large amounts of data in parallel.
Computation is carried out in kernels, small fragment programs that operate independently of one
another on separate cores. A single GPU is made up of multiprocessors, where each multiprocessor
contains N cores. A single GPU can have in total up to 512 cores. Each core can access global
memory, local memory and shared memory, and in some cases system board memory.
GPU’s are used widely in scientific applications, with each successive generation of hardware
outdoing its predecessor by an order of magnitude. In figure 1 we can see a discernable deviation in
the performance improvements of GPU’s when compared to CPU’s.

Fig 1. Performance comparison between the CPU and GPU (Nvidia)

15th Int Symp on Applications of Laser Techniques to Fluid Mechanics
Lisbon, Portugal, 05-08 July, 2010

The fundamental difference between the CPU and GPU architectures lies not in how manytransistors are built-in, but how they are utilised. The majority of transistors on a GPU are focusedon computation rather than data caching and flow control, as shown in figure 2. While floatingpoint performance is high, it is not without restrictions. GPU’s were originally optimised to workwith single-precision floating-point numbers. With double precision the performance falls
dramatically, however, many operations in scientific applications can easily leverage precision,using double precision only when necessary. Other limitations include memory latency andcoalescence. While memory throughput is high on the GPU board, memory transfers from and tothe board are slow. To remedy this a large amount of RAM is installed on the board.

Fig 2. CPU and GPU architectures (NVidia)
Memory coalescence describes the order in which memory is accessed. Typically with CPU’s, notmuch thought is given as to what order elements in a memory block are read. The advanced flow
control and data caching ensure that performance is not appreciably affected. This is not the case forGPU’s, where the order of data access does matter and can severely hinder performance if care isnot taken. A simple example is the access of one element from a global memory block
simultaneously from multiple kernels or threads, creating a conflict. The flow of data is determined
ultimately by the programmer, who must design the code to minimise memory conflicts. The
current investigation focuses only on the LSM processing of data, however, it is clearly apparent
that the voxel reconstruction pre-processing of the data is also a candidate for parallel processing.

3. Experimental setup
To test the processing capabilities of the GPU’s, the results from another volumetric investigation[3] were employed as input data. A descending toroidal vortex was captured via three timeresolved cameras while illuminated by a volumetric light source. Two consecutive frames fromeach camera were used in the voxel reconstruction. The resulting reconstructed volume for thisapplication was 10 x 10 x 1 cm. The experimental setup consisted of a large quiescent tank,octagonal in form to ensure orthogonal orientations of the cameras, as seen in figure 3.
To quantify the actual data processing improvements several approaches were utilised:�� LSM processing carried out on CPU only machines, comparing the performance ofdifferent processor types (variations in the number of threads: 2, 4, 8)

15th Int Symp on Applications of Laser Techniques to Fluid Mechanics
Lisbon, Portugal, 05-08 July, 2010

�� LSM processing carried out on 2 different GPU cards: one with 128 cores (NVidia Quadro3700) and another with 256 cores (NVidia Tesla 1060C) .

Fig 3. Experimental layout.
Actual code implementation required usage of commercially available GPU libraries, namely
CUDA from NVidia. CUDA consists of an API for program flow, kernel execution and memory
management onboard the GPU. Processing takes place in kernel functions defined by theprogrammer and executed in parallel by the API. It may seem tempting to carry out entire
calculations on a kernel-by-kernel basis, but this is not practical. Conditional level program flowneeds to take place outside the kernel due to the operating system misinterpreting any delay from
hardware as a potential timeout. Furthermore, it is more effective to evaluate a three dimensional
interpolation of all the elements inside a cuboid using all available cores rather than using one core
to evaluate an entire cuboid.

Fig 4. LSM result of toroidal vortex
A number of kernel functions were developed for tackling tasks such as matrix operations, equation
solving, 3D interpolation, and others. It is not a straightforward task to convert multithreaded codeoriginally designed for CPU’s for GPU use. Computation tasks are often dissected into serial15th Int Symp on Applications of Laser Techniques to Fluid Mechanics
Lisbon, Portugal, 05-08 July, 2010

segments and parallel kernels. In reality, computation is often a mix of these two approaches and is
so termed heterogeneous computation. To minimise memory latencies, all processing, regardless ofwhether the processing was serial or parallel, took place onboard the GPU card. The CPU, now freeof the burden of processing, was simply used to orchestrate the process, allowing for other moremundane tasks to execute.

4. Analysis
The voxel volume is divided into equal sized volumetric chunks called cuboids. A cuboid at time t1
is termed the template cuboid, and at time t2 is termed the search cuboid. The search cuboid is
typically several times larger than the template cuboid. With each iteration the search cuboid would
be deformed (via an affine transformation) and compared against the template cuboid until asuitable match could be found.
To effectively evaluate the processing performance, cuboids would be selected from different
locations in the voxel volume. Regions of unidirectional or smooth flow converged rapidly after a
few iterations (5 <), whereas regions of high shear and rotation would take noticeable longer. The
number of iterations observed was dependent on:
�� Seeding density
�� Time between successive frames
�� The affine parameters enabled in the analysis
As in PIV, poor seeding densities would hinder convergence. Large frame intervals would meanlarger translations or rotations, thus extending the time of convergence. In some extreme cases thesolution would simply oscillate and be rejected.

5. Results
Preliminary tests were conducted using only multi-core CPU’s using the code originally developed
and then made multi-threaded. Varying the number of threads available a clear pattern could beobserved. A 5 x 5 x 5 cuboid block was tested, each cuboid being 25 x 25 x 25 voxels. A block wasplaced in the mean flow and in the region of maximum rotation in the vortex. The results were thanextrapolated for a larger region to give an estimation of the total processing time for a 100 x 100 x
100 vector field (one million vectors).
908070Time (s)60504030
201000246810

Num ber of CPU coresFig 5. Time to process 125 cuboids vs. number of CPU cores15th Int Symp on Applications of Laser Techniques to Fluid Mechanics
Lisbon, Portugal, 05-08 July, 2010

As indicated earlier, two GPU cards were investigated and compared.

6.00E-03

5.00E-03

Time in secs

4.00E-03

GPUCPU

3.00E-03

2.00E-03

1.00E-03

0.00E+00
050000 100000 150000 200000 250000 300000 350000 400000

Size of interpolation volume (count)Fig 6. Performance comparison between GPU and CPU for deformation processing of one cuboid
For the case of GPU’s, the performance improvements were dramatically better than CPU’s. The
bulk of the processing time for LSM takes place in the iterative deformation of the cuboids which
require efficient three-dimensional interpolation schemes. The processing of interpolations alone
were an impressive two orders of magnitude faster than on a typical workstation CPU. The fact thatthe time does not increase appreciably with increasing volume size points to the overhead of kernel
execution.
In total, a 128-core GPU could process a 25 x 25 x 25 cuboid 20-30 times faster than an 8-core
CPU. The 256-core GPU (Tesla) could interpolate a cuboid 40-50 times faster than an 8-core CPU.
In addition, other processes such as matrix multiplication, inversion, equation solving were all
markedly faster on GPU’s than on the CPU.

6. Conclusions
It is clear from the results that increasing the number of processing threads in the CPU improved
processing times for LSM. While the improvement is not entirely linear, it is nonethelessconsiderable. This discrepancy is most probably due to the fact that some aspects of the processing
are not parallel, and that there is considerable overhead in the validation segment of the processing.

15th Int Symp on Applications of Laser Techniques to Fluid Mechanics
Lisbon, Portugal, 05-08 July, 2010

As clearly seen in the results GPU implementation of LSM clearly reduces the overall time required
to process a large volumetric dataset. It is also clear that in flows that experience reasonably large
rotation and shear that a GPU solution would clearly reduce the workload experienced by the pc.
The lack of double-precision performance in the GPU did not affect the overall performance.
Further investigations are needed to optimise the existing code even further as the total processing
time is still too slow for the processing of many time-steps.

7. References[1] Kitzhofer, J., Westfeld, P., Pust, O., Maas, H.-G. & Brücker, Ch., 2010. Direct Estimation of 3D
Deformation and Rotation Rate Tensor from Gray Value Weighted Voxel Space via Least Squares
Matching. Accepted for publication in the Proceedings of the 15th International Symposium on
Applications of Laser Techniques to Fluid Mechanics, 5 – 8 July 2010, Lisbon.
[2] Kitzhofer, J., Brücker, Ch. (2009) Tomographic Particle Tracking Velocimetry, submitted toExperiments in Fluids[3] Kitzhofer J., Brücker Ch. & Pust O., 2009. Tomo PTV using 3D Scanning Illumination andTelecentric Imaging. Proceedings of the 8th International Symposium on Particle Image
Velocimetry, 25-28 August, Melbourne, Victoria, Australia.).
[4] Maas, H.-G., Westfeld, P., Putze, T., Bøtkjær, N., Kitzhofer, J., Brücker, C., 2009.
Photogrammetric techniques in multi-camera tomographic PIV. In: Soria, J., Atkinson, C. (eds.),
Proceedings of the 8th International Symposium on Particle Image Velocimetry, 25-28 August,Melbourne, Victoria, Australia, pp. 599-602.