Exp Fluids (2006) 41:933–947

DOI 10.1007/s00348-006-0212-z

R E S E A R C H A RT I C L E

Tomographic particle image velocimetry

G. E. Elsinga Æ F. Scarano Æ B. Wieneke Æ
B. W. van Oudheusden

Received: 13 April 2006 / Revised: 4 September 2006 / Accepted: 19 September 2006 / Published online: 11 October 2006

Springer-Verlag 2006

Abstract This paper describes the principles of a strates the capability and potential of the proposed
novel 3D PIV system based on the illumination, system with four cameras. The capability of the tech-
recording and reconstruction of tracer particles within nique in real experimental conditions is assessed with
a 3D measurement volume. The technique makes use the measurement of the turbulent flow in the near
of several simultaneous views of the illuminated par- wake of a circular cylinder at Reynolds 2,700.
ticles and their 3D reconstruction as a light intensity
distribution by means of optical tomography. The
technique is therefore referred to as tomographic
particle image velocimetry (tomographic-PIV). The 1 Introduction
reconstruction is performed with the MART algo-
rithm, yielding a 3D array of light intensity discretized The instantaneous measurement of the 3D velocity
over voxels. The reconstructed tomogram pair is then field is of great interest to fluid mechanic research as it
analyzed by means of 3D cross-correlation with an enables to reveal the complete topology of unsteady
iterative multigrid volume deformation technique, coherent flow structures. Moreover, 3D measurements
returning the three-component velocity vector distri- are relevant for those situations where the flow does
bution over the measurement volume. The principles not exhibit specific symmetry planes or axes and sev-
and details of the tomographic algorithm are discussed eral planar measurements are required for a sufficient
and a parametric study is carried out by means of a characterization. In particular this applies to flow tur-
computer-simulated tomographic-PIV procedure. The bulence, which is intrinsically 3D and its full descrip-
study focuses on the accuracy of the light intensity field tion therefore requires the application of measurement
reconstruction process. The simulation also identifies techniques that are able to capture instantaneously its
the most important parameters governing the experi- 3D structure, the complete stress tensor and the vor-
mental method and the tomographic algorithm ticity vector. The advent of PIV and its developments
parameters, showing their effect on the reconstruction (stereo-PIV, Arroyo and Greated 1991; dual-plane
accuracy. A computer simulated experiment of a 3D stereo-PIV, Kähler and Kompenhans 2000) showed the
particle motion field describing a vortex ring demon- capability of the PIV technique to quantitatively visu-
alize complex flows. Several different methods were
also proposed to achieve a 3D version of the technique
(scanning light sheet, Brücker 1995; holography,
G. E. Elsinga (&)
F. Scarano
B. W. van Oudheusden
Department of Aerospace Engineering, Hinsch 2002; 3D PTV, Maas et al. 1993).
Delft University of Technology, Kluyverweg 1, A novel system for 3D velocity measurements based
Delft 2629 HS, The Netherlands on tomographic reconstruction of the 3D particle dis-
e-mail: g.e.elsinga@tudelft.nl tribution is investigated in the present study. Record-
B. Wieneke ings of particle images from an illuminated volume
LaVision GmbH, Göttingen, Germany taken from several viewing directions simultaneously

123
934 Exp Fluids (2006) 41:933–947

are used to reconstruct the 3D light intensity distribu- sensor (Digital-Holographic-PIV, Coëtmellec et al.
tion. The method is therefore referred to as tomo- 2001). In that case the light intensity distribution in the
graphic particle image velocimetry (tomographic-PIV). measurement volume is evaluated numerically, usually
Provided that two subsequent exposures of the particle by solving the Fresnel diffraction formula on the
images are obtained, the measurement technique re- hologram (near-field diffraction, Pan and Meng 2002).
turns the instantaneous velocity field within the mea- CCD sensors, however, have a very limited spatial
surement volume by means of 3D particle pattern resolution compared to the photographic plate ret-
cross-correlation. urning about 2 to 3 orders less particle images and
The motivation for such a 3D PIV system is pre- velocity vectors. Moreover, the large pixel pitch
sented in the next section together with an overview of requires that the recording is obtained at a relatively
current alternative techniques. Then the working small angle (a few degrees between reference beam
principle and the tomographic reconstruction method and scattered light) in order to resolve the interference
are discussed in detail, followed by a numerical simu- pattern, hence strongly limiting the numerical aperture
lation study assessing the effect of the relevant exper- and depth resolution (Hinsch 2002).
imental and reconstruction parameters and the The scanning-PIV technique is directly derived from
full-scale capabilities of tomographic-PIV. The study standard 2C or stereo PIV with the light sheet scanning
concludes with the application of the proposed tech- through the measurement volume (Brücker 1995). The
nique to a cylinder wake flow in the turbulent regime. volume is sliced by the laser sheet at sequential depth
The measurement of 3D coherent structures within the positions where the particle image pattern is recorded.
wake is discussed and assessed. The second recording at that depth position can be
taken either directly after the first or after the complete
scan of the volume. The procedure returns planar
2 Current 3D PIV techniques velocity fields obtained slightly shifted in space and
time, which can be combined to return a 3D velocity
Among the different 3D velocimetry techniques pres- field. The strong point of this technique is the high in-
ently available holographic-PIV has received most plane spatial resolution and the fact that the analysis of
attention (Hinsch 2002; Chan et al. 2004). It uses the the recordings is straightforward. However, scanning-
interference pattern of a reference light beam with light PIV requires high-repetition systems (kHz) to ensure
scattered by a particle, which is recorded on a holo- that the complete volume recording is almost simul-
gram, to determine the particle location in depth. The taneous. The underlying hypothesis of scanning-PIV is
in-plane position in principle is given by the position of that the volume scanning time needs to be small if
the diffraction pattern in the image. Illumination of the compared with the characteristic time scale of the
hologram with the reference light beam reproduces the investigated flow structure. Cameras with high a
original light intensity field in the measurement volume recording rate are thus required, which is not a prob-
at the time of recording, the intensity being highest at lem in low speed flow as shown by Hori and Sakaki-
the original particle location. The reconstructed inten- bara (2004) measuring a turbulent jet in water.
sity field is scanned by a sensor, e.g. a CCD, to obtain a However, the technique is unsuited for air flows and in
digital intensity map, which can be used for cross-cor- particular in high speed flows. Moreover, high repeti-
relation yielding the velocity field. So far holographic- tion rate lasers are expensive and provide relatively
PIV has shown a great potential in terms of a high data low pulse energy. It should also be remarked that the
yield. However, its drawbacks are that the recording experimental setup is significantly more complicated
medium is a holographic film requiring wet processing, by the addition of a scanning mechanism.
which makes the process time consuming and somehow As an alternative to scanning, dual plane stereo PIV
inaccurate due to misalignment and distortion when re- (Kähler and Kompenhans 2000) records the particle
positioning the hologram for the object reconstruction. images in different planes simultaneously using light
The technique was successfully applied to measure a polarization or different colors to distinguish the scat-
vortex ring in air, the wake of a mixing tab in water (Pu tered light from the two planes. In principle, mea-
and Meng 2000), a cylinder wake flow in air and a free surements can be performed over more than two
air nozzle flow (Herrmann et al. 2000) returning large planes with each plane requiring a double-pulse laser,
numbers of vectors (up to 92,000 using individual par- however separation by polarization is the most com-
ticle pairing, Pu and Meng 2000). monly adopted solution and is possible only over two
Instead of recording on a photographic plate, the planes. Furthermore, using different colors complicates
hologram can also be captured directly by a CCD the optical arrangement.

123
Exp Fluids (2006) 41:933–947 935

Relatively recent options are 3D PTV (Maas et al. 3 Working principle of tomographic-PIV
1993) and defocusing PIV (Pereira et al. 2000). These
techniques rely on the identification of individual par- The working principle of tomographic-PIV is sche-
ticles in the PIV recordings. The exact position of the matically represented in Fig. 1. Tracer particles im-
particle within the volume is given by the intersection mersed in the flow are illuminated by a pulsed light
of the lines of sight corresponding to a particle image in source within a 3D region of space. The scattered light
the recordings from several viewing directions (typi- pattern is recorded simultaneously from several view-
cally three or four). The implementation of the particle ing directions using CCD cameras similar to stereo-
detection and location varies with the methods. In PIV, applying the Scheimpflug condition between the
comparison with the previous two methods, the 3D- image plane, lens plane and the mid-object-plane. The
PTV approach offers the advantage of being fully particles within the entire volume need to be imaged in
digital and fully 3D without the requirement for mov- focus, which is obtained by setting a proper f/#. The 3D
ing parts. The velocity distribution in the volume is particle distribution (the object) is reconstructed as a
obtained from either particle tracking or by 3D-cross- 3D light intensity distribution from its projections on
correlation of the particles pattern (Schimpf et al. the CCD arrays. The reconstruction is an inverse
2003). However, the procedure for individual particle problem and its solution is not straightforward since it is
identification and pairing can be complex and as it is in general underdetermined: a single set of projections
common for planar PTV several algorithms have been can result from many different 3D objects. Determining
proposed, which significantly differ due to the prob- the most likely 3D distribution is the topic of tomog-
lem-dependent implementation. The main limitation raphy (Herman and Lent 1976), which is addressed in
reported in literature is the relatively low seeding
density to which these techniques apply in order to
keep a low probability of false particle detection and of
overlapping particle images. Moreover, the precision
of the volume calibration or in the description of the
imaging optics is finite. This means the lines of sight for
a particle almost never truly intersect and an inter-
section criterion is needed. Consequently the maxi-
mum seeding density in 3D PTV is kept relatively low.
Maas et al. (1993) report a seeding density of typically
0.005 particles per pixel for a three camera system.
The development of the proposed tomographic-PIV
technique is motivated by the need to achieve a 3D
measurement system that combines the simple optical
arrangement of the photogrammetric approach with a
robust particle volume reconstruction procedure,
which does not rely on particle identification. As a
consequence the seeding density can be increased, with
respect to 3D PTV, to around 0.05 particles per pixel
(as will be shown later), which is not far from the
particle image density used in planar experiments. As a
consequence the flow field within the 3D domain can
be represented with a number of velocity vectors
comparable or slightly higher than obtained in planar
PIV. Furthermore, the proposed technique features the
instantaneous flow field measurement, as opposed to
scanning PIV, opening the possibility to perform 3D
measurements in several conditions irrespective of the
flow velocity. Finally the introduction of high-repeti-
tion rate PIV hardware is expected to further extend
the measurement technique capability to a 3D time-
resolved flow diagnostic tool as already demonstrated
in a recent study (Schröder et al. 2006). Fig. 1 Principle of tomographic-PIV

123
936 Exp Fluids (2006) 41:933–947

the following section. The particle displacement (hence system. In the present approach the measurement
velocity) within a chosen interrogation volume is then volume containing the particle distribution (the object)
obtained by the 3D cross-correlation of the recon- is discretized as a 3D array of cubic voxel elements in
structed particle distribution at the two exposures. The (X, Y, Z) (in tomography referred to as the basis
cross-correlation algorithm is based on the cross cor- functions) with intensity E(X, Y, Z). A cubic voxel
relation analysis with the iterative multigrid window element has a uniform non-zero value inside and zero
(volume) deformation technique (WIDIM, Scarano outside and its size is chosen comparable to that of a
and Riethmuller 2000) extended to 3D. pixel, because particle images need to be properly
The relation between image (projection) coordi- discretized in the object as it is done in the images.
nates and the physical space (the reconstruction vol- Moreover, the interrogation by cross-correlation can
ume) is established by a calibration procedure common be easily extended from a pixel to a voxel based object.
to stereo PIV. Each camera records images of a cali- Then the projection of the light intensity distribution
bration target at several positions in depth throughout E(X, Y, Z) onto an image pixel (xi, yi) returns the pixel
the volume. The calibration procedure returns the intensity I(xi, yi) (known from the recorded images),
viewing directions and field of view. The tomographic which is written as a linear equation:
reconstruction relies on accurate triangulation of the X
views from the different cameras. The requirement for wi;j EðXj ; Yj ; Zj Þ ¼ Iðxi ; yi Þ; ð1Þ
a correct reconstruction of a particle tracer from its j2Ni

images sets the accuracy for the calibration to a frac-

tion of the particle image size (see Sect. 5). Therefore, where Ni indicates the voxels intercepted or in the
a technique for the a-posteriori correction for the sys- neighborhood of the line of sight corresponding to the
tem misalignment can significantly improve the accu- ith pixel (xi,yi) (shaded voxels in Fig. 2). The weighting
racy of reconstruction (Wieneke 2005). The mapping coefficient wi,j describes the contribution of the jth
from physical space to the image coordinate system can voxel with intensity E(Xj, Yj, Zj) to the pixel intensity
be performed by means of either camera pinhole I(xi, yi) and is calculated as the intersecting volume
model (Tsai 1986) or by a third-order polynomial in x between the voxel and the line of sight (having the
and y (Soloff et al. 1997). cross sectional area of the pixel) normalized with the
voxel volume. The coefficients depend on the relative
size of a voxel to a pixel and the distance between the
4 Tomographic reconstruction algorithm voxel center and the line of sight (distance d in Fig. 2).
Note that 0 £ wi,j £ 1 for all entries wi,j in the 2D array
The novel aspect introduced with tomographic-PIV is W and that W is very sparse, since a line of sight
the reconstruction of the 3D particle distribution by intersects with only a small part of the total volume.
optical tomography. Therefore, a separate section is The weighting coefficients can also be used to account
devoted to the tomographic reconstruction problem for different camera sensitivities, forward or backward
and algorithms for solving it. scatter differences or other optical dissimilarities be-
By considering the properties of the measurement tween the cameras. Alternatively the recorded images
system, it is possible to select a-priori the reconstruc- can be pre-process in an appropriate way to compen-
tion method expected to perform adequately for the sate for these effects.
given problem. First the particle distribution is dis- In the above model, applying geometrical optics, the
cretely sampled on pixels from a small number of recorded pixel intensity is the object intensity E(X, Y,
viewing directions (typically 3 to 6 CCD cameras) and Z) integrated along the corresponding line of sight. In
secondly it involves high spatial frequencies. In these that case the reconstructed particle is represented by a
conditions algebraic reconstruction methods are more 3D Gaussian-type blob, which projection in all direc-
appropriate with respect to analytical reconstruction tions is the diffraction spot.
methods, such as Fourier and back-projection methods A range of algebraic tomographic reconstruction
(Verhoeven 1993). For this reason, only the former algorithms is available to solve these equations. How-
class of methods is considered for further evaluation. ever, due to the nature of the system, the problem is
underdetermined (Sect. 3) and the calculation may
4.1 Algebraic methods converge to different solutions, which implies that
these algorithms solve the set of equations of Eq. 1
Algebraic methods (Herman and Lent 1976) iteratively under different optimization criteria. A detailed dis-
solve a set of linear equations modeling the imaging cussion and analysis of these optimization criteria is

123
Exp Fluids (2006) 41:933–947 937

Fig. 2 Representation of the

imaging model used for
tomographic reconstruction.
In this top-view the image
plane is shown as a line of
pixel elements and the
measurement volume is a 2D
array of voxels. The gray level
indicates the value of the
weighting coefficient (wi,j) in
each of the voxels with Z
respect to the pixel I(x1,y1) d

Y
X w
θ
line-of-sight

camera 1 I(x1,y1) I(x2,y2) camera 2

beyond the scope of the present study and can be end loop 1
found in Herman and Lent (1976). Instead the per-
where l is a scalar relaxation parameter, which for
formance of two different tomographic algorithms is
ART is between 0 and 2 and for MART must be £1. In
evaluated by numerical simulations of a tomographic-
ART the magnitude of the correction depends on the
PIV experiment focusing the evaluation upon the P
residual Iðxi ; yi Þ  j2Ni wi;j EðXj ; Yj ; Zj Þ multiplied by
reconstruction accuracy and convergence properties.
a scaling factor and the weighting coefficient, so that
The comparison is performed between the additive and
only the elements in E(X, Y, Z) affecting the ith pixel
multiplicative techniques referred to as ART (alge-
are updated. Alternatively in MART the magnitude of
braic reconstruction technique) and MART (multipli-
the update is determined by the ratio of the measured
cative algebraic reconstruction technique), respectively
pixel intensity I with the projection of the current object
(Herman and Lent 1976). Starting from a suitable ini- P
tial guess (E(X, Y, Z)0 is uniform, see next section) the j2Ni wi;j EðXj ; Yj ; Zj Þ The exponent again ensures that
only the elements in E(X, Y, Z) affecting the ith pixel
object E(X, Y, Z) is updated in each full iteration as:
are updated. Furthermore, the multiplicative MART
1. for each pixel in each camera i: scheme requires that E and I are definite positive.
2. for each voxel j:
4.2 2-D numerical assessment
ART: EðXj ; Yj ; Zj Þkþ1
P A domain with reduced dimensionality was adopted to
Iðxi ; yi Þ  wi;j EðXj ; Yj ; Zj Þk
j2Ni evaluate the performances of the two methods: a
¼ EðXj ; Yj ; Zj Þk þ l P 2 wi;j
wi;j 50 · 10 mm2 2D slice of the 3D volume. The 2D par-
j2Ni ticle field contains 50 particles, which images are re-
ð2Þ corded by three linear array cameras with 1,000 pixel.
For the ART reconstruction the relaxation parameter
is set at 0.2 and the initial condition is a uniform zero
MART: EðXj ; Yj ; Zj Þkþ1 intensity distribution, while for MART relaxation and
, !lwi;j
X k
initial condition are both 1.
k
¼ EðXj ; Yj ; Zj Þ Iðxi ; yi Þ wi;j EðXj ; Yj ; Zj Þ The reconstructed particle fields after 20 iterations
j2Ni are shown in Fig. 3. The maximum intensity is 75
ð3Þ counts and values below 3 counts are blanked for
readability. The ART algorithm leaves traces of the
particles along the lines of sight, while the MART
end loop 2 reconstruction shows more distinct particles. The

123
938 Exp Fluids (2006) 41:933–947

additive ART scheme works similar to an OR-opera- numerical simulations. The number of iterations, the
tor: in order to have a non-zero intensity in the object it number of cameras, the viewing directions, the particle
is sufficient to have a particle at the corresponding image density, the calibration accuracy and the image
location in one of the PIV recordings. Still the intensity noise are considered.
is highest at the actual particle locations. The multi- In the simulation the order of the problem is re-
plicative MART scheme behaves as AND-operator: duced from a 3D volume with 2D images to a 2D slice
non-zero intensity only at locations where a particle with 1D images, which simplifies the computation and
appears in all recordings. The suitability of MART to the interpretation of the results without loosing gen-
reconstruct objects with sharp gradients or spikes has erality on the results. In fact the 2D volume can be
been confirmed by Verhoeven (1993). In conclusion, seen as a single slice selected from a 3D volume and
the artifacts in the ART reconstruction are undesir- similarly the 1D image as a single row of pixels taken
able, especially considering that the signal has to be from a 2D image.
cross-correlated for interrogation to provide the de- Tracer particles are distributed in a 50 · 10 mm2
sired velocity (displacement) information. The result slice, which is imaged along a line of 1,000 pixels from
from the MART algorithm better suits tomographic- different viewing directions h (Fig. 2) by cameras
PIV. Moreover, Fig. 4 shows that the individual parti- placed at infinity, such that magnification and viewing
cles reconstructed with MART are reconstructed at the direction are constant over the field of view and the
correct position with an intensity distribution slightly magnification is identical for all views and close to 1.
elongated in depth (dz ~ 1.5dx for the present case). Furthermore, the entire volume is assumed to be
Besides the sequential update algorithms of Eqs. 2 within the focal depth. Given the optical arrangement
and 3 schemes that update the reconstructed object the particle location in the images is calculated and the
simultaneously at every pixel exist, using all equations in application of diffraction (particle diameter is 3 pixels,
a single step, such as the conjugate gradient method which is justified by the large f/# required) results in the
(additive scheme) as well as several implementations of synthetic recordings. The particle image intensity is
the simultaneous MART algorithm (Mishra et al. 1999). assumed uniform. The 2D particle field is recon-
Those methods return similar results and are potentially structed from these recordings at 1,000 · 200 voxel
computationally more efficient. The investigation of resolution using the MART algorithm described in the
those methods goes beyond the scope of the present previous section. Unless stated otherwise the following
study and may be the subject of further investigations. (baseline) experimental parameters are assumed: three
cameras at h = {–20, 0, 20} degrees, 50 particles
(0.05 particles per pixel), 5 reconstruction iterations
5 Parametric study and no calibration errors or image noise.
To quantify the accuracy of the reconstruction pro-
This section discusses the effect of experimental and cess, the reconstructed object E1(X, Z) is compared
algorithm parameters on the reconstruction based on with the exact distribution of light intensity E0(X, Z)

Fig. 3 Particle field ART

reconstructed using ART 10
60
(top) and the same field
reconstructed using MART
z (mm)

40
(bottom). The actual particle 5
positions are indicated by 20
circles. The gray level
represents the intensity level 0 0
-25 -20 -15 -10 -5 0 5 10 15 20 25
x (mm)

MART
10
60
z (mm)

40
5
20

0 0
-25 -20 -15 -10 -5 0 5 10 15 20 25
x (mm)

123
Exp Fluids (2006) 41:933–947 939

MART experimental results show that with additional itera-

4 tions the returned vector field changes only within the
60
noise level.
3.5 50 Figure 6 shows the dependence of the reconstruc-
40 tion quality Q on the most relevant experimental
z (mm)

3 parameters, namely the number of cameras, the view-

30
ing angle, the particle image density (in particles per
20 pixel), and the calibration error. The diagrams show
2.5
10
clear trends, which provides the experimentalist with a
first indication of the optimum experimental arrange-
2 0
-9 -8.5 -8 -7.5 -7 ment and the limitations of the system. A correlation
x (mm) coefficient of 75% with respect to the exact spatial
distribution (Q = 0.75 corresponding to an expected
Fig. 4 Detail of the MART reconstruction (Fig. 3) showing
cross-correlation coefficient of Q2 = 0.56) is used as a
individual particles
cut-off value, above which the reconstruction should be
where the particles are represented by a Gaussian considered sufficiently accurate.
intensity distribution of three voxel diameter. The The effect of the number of cameras is clear: adding
reconstruction quality Q is defined as the normalized camera gives additional information on the object,
correlation coefficient of the exact and reconstructed which increases reconstruction accuracy. A 2-camera
intensity distribution, according to: system (h = –20 and 20) is largely insufficient,
whereas Q rapidly increases going to three and four
P
E1 ðX; ZÞ
E0 ðX; ZÞ cameras (h = –20, 0, 20 and 40) and approaches
X;Z
Q ¼ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
ffi: ð4Þ unity with five cameras (h = –40, –20, 0, 20 and
P 2 P 2
E1 ðX; ZÞ
E0 ðX; ZÞ 40).
X;Z X;Z The viewing angles are changed maintaining the
symmetric camera arrangement. The angle indicated in
A direct estimate of the correlation coefficient to be Fig. 6 is the angle between the outer cameras and the
expected in the cross-correlation analysis of the z-axis. The graph shows an optimum near 30. For
reconstructed objects is therefore given by Q2 in the smaller angles the depth resolution decreases resulting
assumption of perfect cross-correlation of the corre- in elongation of the reconstructed-particle in depth.
sponding Gaussian particle fields and uncorrelated For larger angles the intercepted length of the line-of-
reconstruction noise. sight increases, which causes a larger number of par-
First the convergence of the residual and the ticle to be formed with respect to those actually present
reconstruction is considered. The residual (Fig. 5, left) in the illuminated volume. Such extra particles will be
decreases monotonically in the tomographic iterative referred to as ‘‘ghost particles’’ (Maas et al. 1993). This
process. Although this behavior does not guarantee is a problem of ambiguity, which increases with the
that the reconstructed object converges to the exact number of particles, the particle diameter and the
solution (Watt and Conery 1993), for the present length of the line of sight in the volume. The latter
reconstruction algorithm and noise free recordings, the increases with the viewing angle in the present con-
solution converges correctly as evident from Fig. 5, figuration, hence the increase in ghost particles. The
right. A small tendency to divergence after four to five configurations returning an optimum have a viewing
iterations was observed only for a noise level in the angle in the range of 15–45.
recordings in excess of 50% with respect to the particle An increased particle density produces a larger
peak intensity. Such divergence can be countered by percentage of ghost particles consequently decreasing
applying a smaller relaxation (l = 0.2) as suggested by the reconstruction quality. On the other hand, the
Mishra et al. (1999), although that requires more iter- larger number of particles allows a higher spatial
ations to reach a similar reconstruction quality Q. In sampling rate of the flow, returning a potentially higher
the present parametric study and in the experiments spatial resolution. Therefore, a high particle density is
presented in Sect. 7 the reconstruction process is desirable. Based on the simulation results (Fig. 6) the
stopped after five iterations since both the residual and maximum imaged particle density yielding an accept-
reconstruction quality Q do not change significantly able reconstruction quality is 0.075 and 0.15 particles
performing further iterations. Moreover, any diver- per pixel for the three and four camera system,
gence phenomenon is avoided. Furthermore, the respectively. In a 3D system this would approximately

123
940 Exp Fluids (2006) 41:933–947

Fig. 5 Convergence of the 100 1

normalized residual (left) and
the reconstruction quality 0.9
(right)

residual [-]

Q [-]
10-1 0.8

0.7

0.6
10-2
0 5 10 15 20 0 5 10 15 20
iteration iteration

Fig. 6 Reconstruction quality 1 1

Q as a function of number of 0.9 0.9
cameras (top-left), viewing 0.8 0.8
direction (top-right), particle
0.7 0.7
density (bottom-left) and
Q [-]

Q [-]
calibration accuracy (bottom- 0.6 0.6
right) 0.5 0.5
0.4 0.4
0.3 0.3
2 3 4 5 0 10 20 30 40 50 60
number of cameras θ [deg]

1 1
0.9 0.9
0.8 0.8
0.7 0.7
Q [-]

Q [-]
0.6 0.6
0.5 0.5
2 cameras
0.4 3 cameras
4 cameras
0.4
0.3 5 cameras 0.3
0 0.05 0.1 0.15 0.2 0 0.2 0.4 0.6 0.8 1
N_p [ppp] calibration error [pixels]

translate to 0.025 and 0.05 particles per pixel assuming Stereo-PIV (Scarano et al. 2005; Wieneke 2005), which
each particle image spans three pixel lines. Moreover, can reduce calibration errors to less than 0.1 pixel.
Fig. 6-lower-left shows that additional cameras (in the Finally two types of image noise are considered:
configuration mentioned above) allow a higher seeding random noise added to the recordings and background
density for a given reconstruction quality. particles, which are located outside the reconstruction
In actual experiments calibration errors may occur, volume. Figure 7-left shows the effect of random noise
which result in a dislocation of the lines of sight in the fluctuations, which range is expressed as a percentage
reconstruction. As shown by Watt and Conery (1993) of the particle peak intensity. The image noise deteri-
this reduces the accuracy of the reconstruction. To orates the reconstructed particle shape and increases
quantify the error and to find the necessary calibration the number of ghost particles, as noise is mistaken for
accuracy, a calibration error is introduced after particles. As a result the reconstruction quality Q
recording the images displacing the center-camera to- strongly deteriorates with increasing random image
wards the right and the right-camera to the left by an noise. However, image pre-processing can be applied
equal amount (being the calibration error). Lines of to reduce the number of ghost particles. Subtracting a
sight from the different cameras that originally inter- sliding average from the recoding using a window of
sected now form a triangle in the reconstruction vol- 61 pixels significantly improves the reconstruction re-
ume. As seen from Fig. 6 a calibration error of sults (Fig. 7, left). Random image noise up to 25% of
0.4 pixel is the maximum acceptable for an accurate the particle peak intensity has only a small effect after
reconstruction. Stereo PIV calibration methods meet- pre-processing. For higher noise level (50%) the
ing this requirement have been developed for planar change of the particle image shape is important and

123
Exp Fluids (2006) 41:933–947 941

Fig. 7 Effect on the 1 1

reconstruction of white image 0.9 0.9
noise (left) and unmatched 0.8 0.8
particles (right)
0.7 0.7

Q [-]

Q [-]
0.6 0.6
0.5 0.5
outside volume
0.4 original images 0.4 inside volume
subtract sliding average
0.3 0.3
0 10 20 30 40 50 0 10 20 30 40 50
noise level [%] added particles [%]

cannot be recovered by the pre-processing. The effect 0.75. The reconstruction accuracy can be further
of added background particles (located outside the assessed by counting the number of intensity peaks or
volume) is compared with the same amount of particles particles. Each reconstructed object contains 24,400
added inside the volume in Fig. 7-right. This situation actual particles (peaks within 1.5 voxel radius in the
occurs when the reconstructed domain does not in- correct position) and 91,600 ghost particles considering
clude all the illuminated particles either because the only peaks values above 10 in arbitrary units. Even
laser light sheet has a Gaussian profile or because of though the number of ghost particles exceeds signifi-
stray-light and uncontrolled reflections illuminating cantly the number of actual tracer particles, their peak
particles outside of the measurement volume. The intensity Ip is lower as shown by the probability density
particle located outside the volume have a stronger function of peak intensity (Fig. 8). The expected peak
effect compared to the ones inside, therefore the intensity for an actual particle is 70 against 23 for a
reconstruction should include the entire illuminated ghost. Therefore, the contribution of the ghost particles
volume to yield the best results. to the cross-correlation map, hence velocity measure-
ment, is limited.
The displacement field is obtained from the recon-
6 Synthetic 3D experiment structed object using a 3D extension of the WIDIM
algorithm (Scarano and Riethmuler 2000) with
The numerical assessment of tomographic-PIV in a 3D 413 voxels interrogation volumes at 75% overlap.
configuration is performed by simulation of the particle Therefore, each interrogation volume contains on
motion field around a ring vortex. The vortex core is average 25 particle tracers. The measured vector field
located in the center plane z = 0 mm and forms a circle contains 66 · 66 · 10 vectors shown in Fig. 9-left
of 10 mm diameter. The analytical expression of the where the overall motion pattern is well captured. The
displacement field (in voxel units) d is given by: surface corresponding to a vorticity magnitude of 0.13
  voxels/voxel returns the expected torus. Figure 9-right
u presents the vectors in the cross section at x = 0.25 mm
  8R ðRÞ
d¼ 
v¼ l e
l ; ð5Þ with the corresponding vorticity in x-direction (con-
w
tours), where the flow symmetry and the two vortex

where R is the distance to the voxel-center ring and

l = 2 mm is a length scale that defines the width of the 0.4
vortex. The maximum displacement is 2.9 voxels. Inside actual particles
the 35 · 35 · 7 mm3 measurement volume 24,500 par- 0.3
ghost particles

ticles are distributed, which are imaged by four cameras

PDF [-]

placed at infinity. Each camera has a 30 viewing angle 0.2

(to the left, to the right, upward and downward) and
records the images onto 700 · 700 pixels. Consequently 0.1
the particle image density is 0.050 particles per pixel.
Furthermore, the magnification is 1 and the particle 0
0 50 100 150
image diameter is 3 pixels.
I [-]
The measurement volume is reconstructed at 700 · p

700 · 140 voxel resolution performing ten iterations Fig. 8 PDF of particle peak intensity for actual and ghost
and l = 1. The returned reconstruction quality Q is particles in the reconstructed objects

123
942 Exp Fluids (2006) 41:933–947

sections are clearly visible. The normalized cross-cor- diameter D is 8 mm and the free stream velocity is
relation peak value is 0.56 and the mean signal-to-noise 5 m s–1.
ratio exceeds 5 indicating a high confidence level for The flow is seeded with 1 lm droplets to a particle
the measurement in such configuration. image density of approximately 0.05 ppp. The illumi-
The measurement accuracy is presented in the form nation source is a Quanta Ray double cavity Nd:YAG
of 2D scatter plots of the displacement error (Fig. 10), laser from Spectra-Physics with a pulse energy of
where 90% of the vectors have an absolute error 400 mJ. A knife-edge slit is added in the path of the
smaller than 0.10 voxels in u and v and less than laser light sheet to cut the low intensity side lobes from
0.16 voxels in w (the slightly larger uncertainty in the light profile and create a sharply defined illumi-
depth direction is due to the 30 viewing direction). nated volume of 8 mm thickness. The cylinder axis is
The asymmetric scatter plot of the error on the w oriented vertically, parallel to the laser propagation
component (Fig. 10, right) appearing as a bias error direction in order to match the largest dimensions of
results from the limited spatial resolution (modulation the measurement volume with the streamwise and
error) in combination with the asymmetric w distribu- spanwise coordinates. This choice still allows to mea-
tion in the flow field. Such bias also appears when sure the Kármán vortices along the wake thickness
cross-correlating the exact distribution of light inten- since the measurement volume is one diameter deep.
sity E0(X, Y, Z) and is therefore not associated to the The time separation between exposures is 35 ls
reconstruction process. yielding a free stream particle displacement of
From the synthetic 3D experiment it can be con- 0.18 mm corresponding to 3.2 voxels.
cluded that tomographic-PIV is capable of the instan- A four-camera system (Fig. 11) is used to record
taneous measurement of flow structures at a good 12-bit images of the tracer particles at 1,280 ·
resolution. 1,024 pixels resolution from different directions. The
cameras are equipped with Nikon objectives set at f/8.
Scheimpflug adapters are used to align the mid-plane
7 Application to a cylinder wake flow of the illuminated area with the focal plane. The
effective field of view common to all cameras is
The experimental validation of the technique is per- 50 · 50 mm2. Table 1 lists the properties specific to
formed in a low-speed open-jet wind tunnel with a each camera, such as the magnification M, the focal
0.40 · 0.40 m2 square cross section. The wake behind a length f and the viewing angle in x and y-direction hx
circular cylinder is measured at ReD = 2,700 where the and hy, respectively.

Fig. 9 Velocity field and iso-

vorticity surface (0.13 voxels/
voxel) from a simulation of a
15
vortex ring (left) and velocity
vectors in the cross section at
X = 0.25 mm (right) the 15 10
contours show the vorticity in
x-direction (in voxels/voxel)
10 5
y [mm]

5 Vorx [-]
0
0.25
y [mm]

0 0.20
-5 0.15
0.10
-5
0.05
-10 0.00
-10 -0.05
-0.10
-15
-0.15
-15 -15 -0.20
-10
-5 -0.25
0 -2 0 2
5
x [mm 10 -2 ]
] 15 2 0 mm z [mm]
z[

123
Exp Fluids (2006) 41:933–947 943

Fig. 10 Scatter plots of 0.3 0.3

the particle displacement
error. Left: x and y 0.2 0.2
displacement. Right: x and

∆w [voxels]
∆v [voxels]
z displacement 0.1 0.1

0 0

-0.1 -0.1

-0.2 -0.2

-0.2 -0.1 0 0.1 0.2 -0.2 -0.1 0 0.1 0.2

∆u [voxels] ∆u [voxels]

Fig. 11 The optical

arrangement for the
tomographic-PIV
experiments in the low speed
wind tunnel

The imaging system is calibrated by scanning a plate In real experimental conditions an assessment of the
with 15 · 12 marks (crosses) through the volume in reconstruction quality is not straightforward, since the
depth direction in steps of 4 mm over a total range of actual particle distribution is unknown. It is possible,
8 mm. In each of the three calibration planes the however, to estimate the number and intensity distri-
relation between the physical coordinates (X, Y, Z) bution of the ghost particles by comparing the light
and image coordinates is described by a third order intensity reconstructed outside and inside the illumi-
polynomial fit. The calibration accuracy is approxi- nated region (representing reconstruction noise and
mately 0.2 pixels. Linear interpolation is used to find signal plus noise, respectively). Note that the light
the corresponding image coordinates at intermediate sheet position can be clearly identified within the
z-locations. reconstructed volume (Fig. 12) due to the application
The intensity distribution is reconstructed in a of a slit in the light path. The results confirm that the
40 · 40 · 10 mm3 volume discretized with 730 · 730 · expected peak intensity for the ghost particles is lower
184 voxels using the MART algorithm with five itera- than for the actual tracer particles (as in Fig. 8).
tions and with the relaxation parameter l = 1. How- Moreover, the ratio of actual particles over ghost
ever, the light sheet thickness covers approximately particles is 2 considering only the range of peak
150 voxels in depth. The reconstruction process is intensities corresponding to the actual particles. For
improved by means of image pre-processing with further detail on the assessment of reconstruction
background intensity removal, particle intensity equal- volumes in real experiments the reader is referred to
ization and a Gaussian smooth (3 · 3 kernel size). Elsinga et al. (2006).

Table 1 Optical arrangement

Camera # Camera type Pixel pitch (lm) M (–) f (mm) hy () hx ()

1 Imager Intense 6.45 0.13 60 +20 +22

2 Imager Intense 6.45 0.12 50 –10 +22
3 SensiCam 6.7 0.12 50 –10 –22
4 SensiCam QE 6.45 0.14 60 +20 –22

123
944 Exp Fluids (2006) 41:933–947

Fig. 12 Top view of the

reconstructed volume
showing the light intensity
integrated in y-direction. The
green lines indicate the
position of the light sheet

The cross-correlation analysis returned 64 · 64 · 30 To improve the visualization of the structural orga-
velocity vectors using an interrogation volume size of nization of the flow the span wise and the combination
41 · 41 · 21 voxels with 75% overlap. Data validation of stream wise and z component of vorticity are color-
based on a signal-to-noise ratio threshold of 1.2 and on coded in Fig. 14. The four uncorrelated snapshots show
the normalized median test with maximum threshold different phases of the vortex shedding cycle. The top-
of 2 (Westerweel and Scarano 2005) returns 4% spu- left snapshot (corresponding to Fig. 13) contains parts
rious vectors. The average signal-to-noise ratio and of four Kármán vortices; two from the lower cylinder
normalized correlation coefficient are 3.8 and 0.6, surface (green) at x/D = 2.8 and 6, and two from the
respectively. upper surface (cyan) at x/D = 2.5 and 4.5. The nor-
An example of an instantaneous velocity distribu- malized vorticity level in these vortices xyD/u¥ = 2.2
tion is presented in Fig. 13. In the plot the y-axis cor- agrees fairly well with the recent planar PIV mea-
responds to the cylinder axis and the x-coordinate is surements in the Reynolds number range 2,000 to
the distance from the cylinder axis in flow direction. 10,000 from Huang et al. (2006), who report an average
Low velocity and flow reversal is observed for x/D < 3 normalized peak vorticity of 2.1 at x/D = 3.
followed by a recovery of the flow velocity to approx- The secondary vortex structures are also clearly
imately 80% of its free stream value at x/D = 6. A visible in Fig. 14 (blue and red depending on the ori-
large swirling motion due to a Kármán vortex shed entation of vorticity in stream wise direction) and ap-
from the lower surface of the cylinder is observed pear to be organized in counter rotating pairs.
around x/D = 2.8 and z/D = –0.2 (labeled as A in Figure 15 shows in detail how the secondary vortices
Fig. 13), characterized by a vorticity peak. Besides the curl in between the primary Kármán vortices. Their
Kármán vortices the vorticity iso-surface in Fig. 13 also effect is to first distort the primary rollers, as observed
reveals a number of secondary streamwise vortex in the upper two snapshots of Fig. 14 at x/D = 4.5, and
structures (indicated by B), which interact with the finally they cause the breakup of the primary vortices.
primary rollers. At the present Reynolds number the The two snapshots of Fig. 14-bottom show a more
shear layers separating from the cylinder are transi- regular span wise organization of the secondary vorti-
tional and three-dimensionality on the scale of the ces yielding a quasi periodic behavior. From a visual
Kármán vortices is expected (Williamson 1996). inspection the normalized spatial wave length ky/D is

Fig. 13 Instantaneous
Z
velocity distribution showing B
one every third vector and Y
vorticity magnitude iso-
surface (|x| = 2.3 · 103 s–1) X

0.5 2 u [m/s]
5
z/D

0 1 4
3
-0.5 0 2
D
2 y/ 1
3 -1 0
4 -1
x/D 5 -2 -2
A
6

123
Exp Fluids (2006) 41:933–947 945

Fig. 14 Instantaneous
vorticity iso-surfaces
(x = 1.4 · 103 s–1). Color
coding: green xy < 0; cyan
xyffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
p > 0; blue
2 þ x2
signðx Þ\0; red
pxffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
x z x
x2x þ x2z
signðxx Þ[0: The
cylinder is shown in yellow

noise or reconstruction artifacts, when present, are in

fact expected to propagate particle (displacement)
information along viewing lines ultimately resulting in
a quasi-uniform particle displacement along that line.
Clearly such effects do not play a dominant role in the
current experiment.
Further confidence in the measurement technique is
given by recent quantitative comparisons between
tomographic- and stereo-PIV, which showed a good
agreement of the returned flow statistics (mean and
RMS) along profiles in the cylinder wake (Elsinga et al.
2006) and velocity distribution in a vortex ring (Wie-
neke and Taylor 2006).

Fig. 15 Detail of the interaction between primary and secondary

vortices taken from upper left snapshot of Fig. 14 8 Conclusions

Tomographic-PIV was introduced as a novel technique

estimated at 1.2, which is again in good agreement with for 3D velocity measurements. The 3D particle distri-
Huang et al. (2006) and slightly in excess with the re- bution within a volume is reconstructed by optical
ported values at lower Reynolds number (Williamson tomography from 2D particle image recordings taken
1996). A more detailed discussion of the cylinder wake simultaneously from several viewing directions.
flow is beyond the scope of the present study, but can Velocity information results from 3D particle pattern
be found in Scarano et al. (2006). correlation of two reconstructions obtained from sub-
The flow velocity exhibits a large variation along sequent exposures of the particle images. The tech-
the viewing directions of the cameras, as shown nique makes use of digital imaging devices as opposed
above, which allows to conclude that independent to holographic-PIV, and it allows a typical spatial
velocity information is measured along the z-coordi- resolution that is in between 3D-PTV (low seeding
nate corresponding with the volume depth. Possible density) and scanning-PIV (high seeding density).

123
946 Exp Fluids (2006) 41:933–947

Moreover, the technique measures the particle motion flow measurements using bacteriorhodopsin (bR). Meas Sci
instantaneously in the 3D observation domain, which Technol 15:647–655
Coëtmellec S, Buraga-Lefebvre C, Lebrun D, Özkul C (2001)
makes it suitable for the analysis of flows in several Application of in-line digital holography to multiple plane
regimes and irrespective of the flow speed. velocimetry. Meas Sci Technol 12:1392–1397
A numerical simulation of a tomographic-PIV exp- Elsinga GE, Van Oudheusden BW, Scarano F (2006) Experi-
eriment showed that the multiplicative algebraic mental assessment of tomographic-PIV accuracy. In: 13th
international symposium on applications of laser techniques
reconstruction techniques (MART) are the most suited to fluid mechanics, Lisbon, Portugal, paper 20.5
to particle reconstruction returning distinct particles Herman GT, Lent A (1976) Iterative reconstruction algorithms.
with limited artifacts. Furthermore, it was shown that Comput Biol Med 6:273–294
the reconstruction algorithm converges to the desired Herrmann S, Hinrichs H, Hinsch KD, Surmann C (2000)
Coherence concepts in holographic particle image veloci-
solution for noise free particle image recordings. Based metry. Exp Fluids 29(7):S108–S116
on a parametric study the number of cameras, viewing Hinsch KD (2002) Holographic particle image velocimetry. Meas
direction, particle density, calibration accuracy and Sci Technol 13:R61–R72
image noise were identified as critical factors deter- Hori T, Sakakibara J (2004) High-speed scanning stereoscopic
PIV for 3D vorticity measurement in liquids. Meas Sci
mining the quality of the 3D particle reconstruction. Technol 15:1067–1078
The calibration must be accurate within 0.4 pixel. Im- Huang JF, Zhou Y, Zhou T (2006) Three-dimensional wake
age pre-processing (e.g., removing background inten- structure measurement using a modified PIV technique. Exp
sity) reduces the effect of random image noise. Adding Fluids 40:884–896
Kähler CJ, Kompenhans J (2000) Fundamentals of multiple
a camera to the system provides extra information, plane stereo particle image velocimetry. Exp Fluids 29:S70–
which can be used to increase measurement resolution S77
(increasing the particle density) or accuracy. Maas HG, Gruen A, Papantoniou D (1993) Particle tracking
A 3D simulation of a ring vortex flow showed that velocimetry in three-dimensional flows. Exp Fluids 15:133–
146
the 3D measurement configuration with four cameras Mishra D, Muralidhar K, Munshi P (1999) A robust MART
can yield 66 · 66 · 10 vectors in a 35 · 35 · 7 mm3 algorithm for tomographic applications. Num Heat Transfer
volume, with a typical measurement error of 0.1 pixels. Part B 35:485–506
The application of the technique to a turbulent Pan G, Meng H (2002) Digital holography of particle fields:
reconstruction by use of complex amplitude. Appl Opt
cylinder wake flow at a Reynolds number ReD = 2700, 42:827–833
provided an experimental assessment of the capability Pereira F, Gharib M, Dabiri D, Modarress D (2000) Defocusing
of tomographic PIV. The measurement volume was digital particle image velocimetry: a 3-component 3-dimen-
aligned with the cylinder axis returning the stream- sional DPIV measurement technique. Exp Fluids 29(7):S78–
S84
wise evolution of the wake as well as its spanwise Pu Y, Meng H (2000) An advanced off-axis holographic particle
organization within a depth of one cylinder diameter image velocimetry (HPIV) system. Exp Fluids 29:184–197
(8 mm). The instantaneous velocity returned the Scarano F, Riethmuller ML (2000) Advances in iterative
Kármán street with counter-rotating vortices alterna- multigrid PIV image processing. Exp Fluids 29:S51–S60
Scarano F, David L, Bsibsi M, Calluaud D (2005) S-PIV
tively shed and cross-linked by secondary vortex comparative assessment: image dewarping + misalignment
structures aligned in the streamwise-binormal direc- correction and pinhole + geometric back projection. Exp
tion. The results give a further confirmation that the Fluids 39:257–266
tomographic-PIV technique is suited to the study of Scarano F, Elsinga GE, Bocci E, Van Oudheusden BW (2006)
Investigation of 3-D coherent structures in the turbulent
complex 3D flows. cylinder wake using Tomo-PIV. In: 13th international
symposium on applications of laser techniques to fluid
Acknowledgments The authors wish to thank Prof. Dave Watt mechanics, Lisbon, Portugal, paper 20.4
(University of New Hampshire) and Dr. Dirk Michaelis (LaVi- Schimpf A, Kallweit S, Richon JB (2003) Photogrammetric
sion GmbH) for the valuable suggestions on optical tomography particle image velocimetry. In: 5th international symposium
and for the software development, respectively. on particle image velocimetry, Busan, Korea, paper 3115
Schröder A, Geisler R, Elsinga GE, Scarano F, Dierksheide U
(2006) Investigation of a turbulent spot using time-resolved
tomographic PIV. In: 13th international symposium on
References applications of laser techniques to fluid mechanics, Lisbon,
Portugal, paper 1.4
Arroyo MP, Greated CA (1991) Stereoscopic particle image Soloff SM, Adrian RJ, Liu ZC (1997) Distortion compensation
velocimetry. Meas Sci Technol 2:1181–1186 for generalized stereoscopic particle image velocimetry.
Brücker Ch (1995) Digital-particle-image-velocimetry (DPIV) in Meas Sci Technol 8:1441–1454
a scanning light-sheet: 3D starting flow around a short Tsai RY (1986) An efficient and accurate camera calibration
cylinder. Exp Fluids 19:255–263 technique for 3D machine vision. In: Proceedings of IEEE
Chan VSS, Koek WD, Barnhart DH, Bhattacharya N, Braat conference on computer vision and pattern recognition,
JJM, Westerweel J (2004) Application of holography to fluid Miami Beach, FL, USA, pp 364–374

123
Exp Fluids (2006) 41:933–947 947

Verhoeven D (1993) Limited-data computed tomography algo- Wieneke B (2005) Stereo-PIV using self-calibration on particle
rithms for the physical sciences. Appl Opt 32:3736–3754 images. Exp Fluids 39:267–280
Watt DW, Conery DG (1993) Computational aspects of optical Wieneke B, Taylor S (2006) Fat-sheet PIV with computation of
tomography for fluid mechanics and combustion. In: 3rd full 3D-strain tensor using tomographic reconstruction. In:
ASME symposium on experimental and numerical flow 13th international symposium on applications of laser
visualization, New Orleans, LA, USA, pp 241–251 techniques to fluid mechanics, Lisbon, Portugal, paper 13.1
Westerweel J, Scarano F (2005) A universal detection criterion Williamson CHK (1996) Vortex dynamics in the cylinder wake.
for the median test. Exp Fluids 39:1096–1100 Annu Rev Fluid Mech 28:477–539

123