Measurement of Three-Dimensional Coherent Fluid Structure in
High Reynolds Number Turbulent Boundary Layers
A dissertation submitted for the degree of
Doctor of Philosophy
at the University of Cambridge by
T. H. Clark
Trinity Hall, Cambridge
February 2012

Originality and length
This thesis describes work carried out between October 2007 and January
2012 in the Division A Fluids Group, at Cambridge University Engineering
Department. The work contained herein was completed under the combined
supervision of Doctor T. B. Nickels and Professor H. Babinsky.
Except where stated otherwise, this thesis is the result of my own work and
contains nothing which is the outcome of work done in collaboration except
where specifically indicated in the text. This thesis has not been submitted in
whole or in part for any degree or diploma at this or any other university.
The length of this thesis does not exceed 65,000 words, including appendices,
bibliography, footnotes, tables and equations. The number of figures is less
than 150.
T. H. Clark
Cambridge, February 2012

To Tim, who helped make Fluid Dynamics wonderful

Acknowledgements
Firstly, I would like to thank my supervisors, Prof. Holger Babinsky and Dr.
Tim Nickels. Tim encouraged me to do a PhD and instigated the research
within this thesis; he patiently and enthusiastically supervised it until Octo-
ber 2010. Following the sadness of Tim’s death, Holger unhesitatingly took
on the supervisory role; providing valuable insight and motivation for me to
finish the work. Tim and Holger: I will always be grateful to you both.
Two colleagues in particular deserve a great deal of thanks. Dr. Nicholas
Worth patiently taught me about Tomographic PIV, answered endless ques-
tions and shared his TomoPRO codes, which provided a foundation for the
Tomographic PIV Toolbox described in Chapter 3. Working with Jose´ Cardesa-
Duen˜as has yielded many of the most enjoyable moments of the past few
years; Jose´ has become a friend, a resource and a ’partner-in-crime’ in un-
ending - and always edifying - discussions on the finer points of experimen-
tal method, turbulence, coherent structure and flow topology.
Elsewhere in the Engineering Department, Peter Benie offered unfailing com-
puting support on issues far beyond my own skill. My thanks extend to Mick
Underwood, Rob Leroy and the rest of the Hydraulics Laboratory technicians
for their help and support, spanning three years of experimental work.
My friends and family have facilitated this work; offering support and en-
couragement throughout its challenges. Particular thanks are due to my par-
ents, step-sister Claire, Arielle, Susie, Greg, Catherine, Tom, Natalia, Tim,
Caroline, Chris, Andi, Wolfgang, Hatem and Ellie.
Sincere thanks also go to the Green-Tide Turbines team; for support far be-
yond the call of duty, for patience, encouragement and for the financial spon-
sorship required to complete this work.
Finally, I wish to thank the Master, Dean, Fellows and Staff of Trinity Hall for
providing a wonderful community within which to study - and for providing
partial funding via the Trinity Hall Brockhouse Scholarship and travel grants.

Abstract
The turbulent boundary layer is an aspect of fluid flow which dominates the
performance of many engineering systems - yet the analytic solution of such
flows is intractable for most applications. Our understanding of boundary
layers is therefore limited by our ability to simulate and measure them.
Tomographic Particle Image Velocimetry (TPIV) is a recently developed tech-
nique for direct measurement of fluid velocity within a 3D region. This al-
lows new insight into the topological structure of turbulent boundary layers.
Increasing Reynolds Number increases the range of scales at which turbu-
lence exists; a measurement technique must have a larger ’dynamic range’
to fully resolve the flow. Tomographic PIV is currently limited in spatial dy-
namic range (which is also linked to the spatial and temporal resolution) due
to a high degree of noise. Results also contain significant bias error.
This work proposes a modification of the technique to use more than two ex-
posures in the PIV process, which (for four exposures) is shown to improve
random error by a factor of 2 to 7 depending on experimental setup parame-
ters. The dynamic range increases correspondingly and can be doubled again
in highly turbulent flows. Bias error is reduced by up to 40%.
An alternative reconstruction approach is also presented, based on applica-
tion of a reduction strategy (elimination of coefficients based on a first guess)
to the tomographic weightings matrixWij. This facilitates a potentially sig-
nificant increase in computational efficiency.
Despite the achieved reduction in error, measurements contain non-zero di-
vergence due to noise and sampling errors. The same problem affects visual-
isation of topology and coherent fluid structures. Using Projection Onto Con-
vex Sets, a framework for post-processing operators is implemented which
includes a divergence minimisation procedure and a scale-limited denoising
strategy which is resilient to ’false’ vectors contained in the data.
Finally, developed techniques are showcased by visualisation of topological
information in the inner region of a high Reynolds Number boundary layer
(δ+ = 1890, Reθ = 3650). Comments are made on the visible flow structures
and tentative conclusions are drawn.
vii

Contents
Abstract vii
1 Introduction 1
2 Literature Review 7
2.1 Tomographic Particle Image Velocimetry . . . . . . . . . . . . 7
2.2 Tomographic PIV in Turbulent Boundary Layers . . . . . . . . 13
2.3 Post-processing for Tomographic PIV results . . . . . . . . . . 19
3 Experimental Methodology 31
3.1 Design of Experiments . . . . . . . . . . . . . . . . . . . . . . . 31
3.2 Experimental Facility and Boundary Layer Characteristics . . 35
3.3 Experimental Setup . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.4 Processing Software . . . . . . . . . . . . . . . . . . . . . . . . . 53
4 Enhancement of Tomographic PIV 57
4.1 Correlation Tracking Enhancement . . . . . . . . . . . . . . . . 57
4.2 Acceleration of Tomographic PIV by Weighting Reduction . . 68
4.3 Image Pre-processing for Tomographic PIV . . . . . . . . . . . 85
5 Effect of the Correlation Tracking Enhancement (CTE) 99
5.1 Experimental Results . . . . . . . . . . . . . . . . . . . . . . . . 99
5.2 Numerical Accuracy Study . . . . . . . . . . . . . . . . . . . . 113
5.3 Review of Results . . . . . . . . . . . . . . . . . . . . . . . . . . 122
6 Constrained Restoration of Tomographic PIV Data 125
6.1 Projection Onto Convex Sets . . . . . . . . . . . . . . . . . . . . 126
6.2 A Divergence Reduction Operator in 3 Dimensions . . . . . . 131
6.3 Wavelet-Based Band Limiting: The P4 Operator . . . . . . . . 135
6.4 The RELAX Algorithm and Variants . . . . . . . . . . . . . . . 141
6.5 Application of POCS to Tomographic PIV Data . . . . . . . . . 146
6.6 Review of Operators . . . . . . . . . . . . . . . . . . . . . . . . 165
7 Coherent Structures in High Reynolds Number Boundary Layers 167
7.1 Review of Boundary Layer and Experimental Parameters . . . 167
ix
x Contents
7.2 Comparison of Standard and Enhanced Techniques . . . . . . 172
7.3 Observation of Structural Elements in the Flow . . . . . . . . . 177
7.4 Preliminary Conclusions on Flow Structure . . . . . . . . . . . 189
8 Conclusions 191
8.1 On Software . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 191
8.2 On Enhancements to Tomographic PIV . . . . . . . . . . . . . 192
8.3 On the Accuracy of Tomographic PIV . . . . . . . . . . . . . . 193
8.4 On Constrained Restoration of Tomographic PIV Data . . . . . 195
8.5 On High Reynolds Number Turbulent Boundary Layers . . . 196
8.6 Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 197
References 201
List of Figures
1.1 Horseshoe vortex [Theodorsen, 1952] . . . . . . . . . . . . . . . . 1
1.2 Conditional averaging and visualisation of topology [Schro¨der et al.,
2011; Elsinga et al., 2010] . . . . . . . . . . . . . . . . . . . . . . . . 3
2.1 Steps of Tomographic PIV [Elsinga et al., 2006] . . . . . . . . . . . 7
2.2 Graphical explanation of bias error . . . . . . . . . . . . . . . . . . 9
2.3 Performance of Multiplicative First Guess [Worth & Nickels, 2008] 12
2.4 Coherent structures visualised using DNS data [Adrian, 2007; Adrian
& Liu, 2002] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.5 ’Lambda’ vortex structure [Perry & Chong, 1982] . . . . . . . . . . 16
2.6 Structure in the Boundary Layer [Elsinga et al., 2010] . . . . . . . 18
3.1 Elongation of reconstructed particles in the Z direction . . . . . . 33
3.2 C.U.E.D. Boundary Layer Water Tunnel . . . . . . . . . . . . . . . 35
3.3 General arrangement of the Tomographic PIV apparatus . . . . . 36
3.4 Estimated tomographic reconstruction volume . . . . . . . . . . . 37
3.5 Reynolds Stresses with wall distance [Schlichting & Gersten, 2003] 39
3.6 Power density spectrum for a turbulent boundary layer [David-
son et al., 2006] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.7 Alignment and general arrangement of laser beam optics. . . . . . 45
3.8 Laser beam expander optics . . . . . . . . . . . . . . . . . . . . . . 46
3.9 Mounting detail for SA1.1 cameras. . . . . . . . . . . . . . . . . . . 47
3.10 Setup of bird bath and goggles to view through the free surface . 48
3.11 Wide-grooved aluminium calibration plate. . . . . . . . . . . . . . 49
3.12 TomoPIV Toolbox script: Image pre-processing for self calibration 50
3.13 Histograms of calibration error . . . . . . . . . . . . . . . . . . . . 51
3.14 Disparity maps showing final calibration error . . . . . . . . . . . 52
3.15 Workflow chart for the Tomographic PIV Toolbox . . . . . . . . . 55
4.1 Windowing pattern for Correlation Tracking Enhancement . . . . 62
4.2 PDF of signal to noise ratio using VODIM and CTE . . . . . . . . 65
4.3 Isosurfaces showing correlation peaks for VODIM and CTE . . . 67
4.4 Surface plots showing correlation peaks for VODIM and CTE . . 68
4.5 Input images for Tomographic PIV . . . . . . . . . . . . . . . . . . 71
xi
xii List of Figures
4.6 Idealised voxel and pixel arrays for a tomographic reconstruction 72
4.7 Elimination of redundant weighting coefficients . . . . . . . . . . 73
4.8 Final elimination step for reduced weightings matrix . . . . . . . 74
4.9 Reduction of a weighting matrix for an idealised setup . . . . . . 75
4.10 Sparsity pattern of Wrs for a 3003 domain size experiment . . . . . 76
4.11 Least Squares Conjugate Gradient matrix solver . . . . . . . . . . 79
4.12 SMART algorithm implemented for a reduced weighting matrix . 80
4.13 Sparsity patterns of Wrs for a four-camera setup . . . . . . . . . . 82
4.14 Quality factor and Normalised Median Residual plotted against
Time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
4.15 Central 32 voxels3 of a domain reconstructed using MART and Wrs 85
4.16 Reconstruction of a single particle at a range of characteristic sizes. 87
4.17 Isosurfaces of performance in parameter space . . . . . . . . . . . 91
4.18 Function prototype for image optimisation in the TomoPIV Toolbox 94
4.19 Convergence of a Nelder-Mead simplex search to an optimal im-
age processing parameter set . . . . . . . . . . . . . . . . . . . . . 95
4.20 Vector field slices before and after image optimisation . . . . . . . 96
5.1 Convergence of mean flow velocity components and turbulent in-
tensities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
5.2 Convergence of velocity field derivatives and topological invari-
ants . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
5.3 Clauser plot for determination of skin friction velocity . . . . . . . 103
5.4 Boundary layer mean profile measured with CTE and VODIM . . 105
5.5 Crossover of measured and analytic velocity profiles . . . . . . . 107
5.6 Bias error as a function of wall-normal distance . . . . . . . . . . . 108
5.7 Bias error as a function of particle displacement . . . . . . . . . . 108
5.8 Turbulent intensities . . . . . . . . . . . . . . . . . . . . . . . . . . 111
5.9 Variation of eX with Z and particle density . . . . . . . . . . . . . 117
5.10 Variation of eY with Z and particle density . . . . . . . . . . . . . 118
5.11 Variation of eZ error with Z and particle density . . . . . . . . . . 119
6.1 Projection onto convex sets in Hilbert Space . . . . . . . . . . . . . 127
6.2 Effect of set concavity on the convergence of a POCS algorithm . 128
6.3 Requirement for orthogonality of projections in Hilbert space . . 133
6.4 Application of homogeneous Dirichlet boundary condition for so-
lution of the Poisson equation . . . . . . . . . . . . . . . . . . . . . 135
List of Figures xiii
6.5 Application of the P1 operator to a velocity field . . . . . . . . . . 136
6.6 Filtering a 2D scalar field using a wavelet transform . . . . . . . . 137
6.7 Filtering a 2D scalar field using a spectral-space top-hat filter . . . 138
6.8 Denoising a flow field variable using wavelet filtering . . . . . . . 140
6.9 Application of the P4 operator to a velocity field . . . . . . . . . . 141
6.10 Function prototype for a projection operator . . . . . . . . . . . . 143
6.11 Flow of data within the *RELAX family of POCS algorithms . . . 144
6.12 Decrease in convergence rate with tangency angle between sets . 144
6.13 Probability Density Functions for normalised divergence values . 148
6.14 Velocity, vorticity magnitude, and wall-normal vorticity compo-
nent data for a turbulent boundary layer . . . . . . . . . . . . . . . 149
6.15 Isosurface plots of the 2nd invariant Λ2 from raw results and fol-
lowing application of RELAX34 . . . . . . . . . . . . . . . . . . . . 151
6.16 Isosurface plot of the 2nd invariant Λ2 following application of a
gaussian smoothing filter . . . . . . . . . . . . . . . . . . . . . . . 152
6.17 Isosurface plots of the 2nd invariant Λ2 shown in 3D, from raw
results and following application of RELAX34 . . . . . . . . . . . 153
6.18 Instantaneous stream lines for visualisation of flow topology . . . 154
6.19 Change in velocity components following use of UNIRELAX3 . . 155
6.20 Improved topological visualisation in the Q-R plane . . . . . . . . 156
6.21 Magnitude of difference between original and filled-in vectors . . 158
6.22 Restoration of field with 20% missing data . . . . . . . . . . . . . 160
6.23 Restoration of field with 50% missing data . . . . . . . . . . . . . 161
6.24 Median error in restoration of a field with varying amount of miss-
ing data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162
6.25 Restoration of data using MULTIRELAX with a signal to noise ra-
tio based relaxation factor . . . . . . . . . . . . . . . . . . . . . . . 162
6.26 Restoration of data using MULTIRELAX with a divergence-based
relaxation factor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163
6.27 Unsteady convergence of MULTIRELAX. . . . . . . . . . . . . . . 164
7.1 Comparison of TPIV measurement volumes . . . . . . . . . . . . . 169
7.2 CTE technique applied for window size 243 . . . . . . . . . . . . . 171
7.3 Comparison of field divergence . . . . . . . . . . . . . . . . . . . . 173
7.4 Quiver plots for standard and enhanced velocity fields . . . . . . 174
7.5 Isosurfaces of Q for standard and enhanced fields . . . . . . . . . 175
xiv List of Figures
7.6 Skin friction lines for standard and enhanced fields . . . . . . . . 176
7.7 Trains of hairpins and low-speed zones . . . . . . . . . . . . . . . 178
7.8 Small-scale hairpin train distorted by larger structure . . . . . . . 180
7.9 Small-scale hairpin train distorted by larger structure . . . . . . . 181
7.10 Skin friction lines (whole field) . . . . . . . . . . . . . . . . . . . . 183
7.11 Skin friction lines (train) . . . . . . . . . . . . . . . . . . . . . . . . 184
7.12 Skin friction lines (isolated hairpins) . . . . . . . . . . . . . . . . . 185
7.13 Vortex rings in the flow . . . . . . . . . . . . . . . . . . . . . . . . . 186
7.14 An isolated vortex ring . . . . . . . . . . . . . . . . . . . . . . . . . 187
7.15 A train of vortex rings . . . . . . . . . . . . . . . . . . . . . . . . . 188
List of Tables
2.1 Boundary layer parameters from Tomographic PIV experiments . 17
2.2 Experimental parameters from Elsinga et al. [2006] . . . . . . . . . 27
2.3 Experimental parameters from Elsinga et al. [2007] . . . . . . . . . 27
2.4 Experimental parameters from Schro¨der et al. [2008] . . . . . . . . 27
2.5 Experimental parameters from Schro¨der et al. [2011] . . . . . . . . 28
2.6 Experimental parameters from Elsinga et al. [2010] . . . . . . . . . 28
2.7 Experimental parameters from Worth et al. [2010] . . . . . . . . . 28
2.8 Experimental parameters from Ku¨hn et al. [2011] . . . . . . . . . . 29
3.1 Boundary layer parameters measured in the C.U.E.D. tunnel [Tee
& Nickels, 2008] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
4.1 Parameters for the two cases used to assess performance of the
Wrs technique. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
4.2 Value ranges for image processing parameters. . . . . . . . . . . . 90
5.1 Boundary layer parameters measured in the CUED tunnel in-sync
with TPIV acquisition . . . . . . . . . . . . . . . . . . . . . . . . . . 104
5.2 Percentage improvement of CTE bias compared to VODIM bias
for varying displacements . . . . . . . . . . . . . . . . . . . . . . . 109
5.3 Variation of true / ghost particle ratio with number of particles
per pixel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116
7.1 Comparison of experimental parameters . . . . . . . . . . . . . . . 168
xv

Nomenclature
Symbols
A+ Boundary layer inner region van Driest profile constant.
C Boundary layer log-law relationship constant defined by Coles
[1956].
Dp Particle diameter. For a gaussian approximation, Dp = 4σ.
E Reconstructed 3D intensity field (voxels domain).
F(Z) Window-filtered analytical boundary layer profile.
I Image 2D intensity field (pixels domain).
Ncams Number of cameras used in a tomographic reconstruction.
Q Topological invariant as defined by Perry & Chong [1982].
Q− crit. Criterion for swirling flow defined by Cucitore et al. [1999].
Q− f ac. Quality Factor of tomographic reconstruction process.
R Topological invariant as defined by Perry & Chong [1982].
Reθ Momentum-thickness Reynolds Number.
U∞ Free-stream speed.
W, Wij Weightings matrix for full tomographic reconstruction prob-
lem.
Wrs Weightings matrix for reduced tomographic reconstruction
problem; also Weightings Reduction Scheme describing the
procedure for reduction of Wij to Wrs.
Λ ’Lambda’ vortex structure described by Perry & Chong [1982].
Λ2 Second invariant of the deformation tensor (usually denoted
λ2 but upper case version used here to prevent conflict).
xvii
xviii Nomenclature
Ω Field domain to which POCS operators are applied.
δ¯i Measured mean displacement in the ith direction.
u¯i Mean velocity component in the ith direction.
δ+ Skin-friction Reynolds Number.
δt Time lapsed between successive image frames.
δ99 Boundary layer thickness (at 99% of U∞).
δir Binary mask used to relate pixel indices in matrix Wij to the
reduced pixel indices in Wrs.
δjs Binary mask used to relate voxel indices in matrix Wij to the
reduced voxel indices in Wrs.
ei Error in the ith velocity component.
u′iu
′
j
uτ Turbulent Reynolds Stress described for all i, j ∈ [XYZ].
γ External (1-padded) boundary grid surrounding Ω.
λ, λi, λiXYZ Denotes the relaxation used in the *RELAX algorithm for the
ith operator at position [X, Y, Z].
A, B, C, D Letter denoting reconstructed 3d particle intensity field at a
point in time e.g. field B occurs at time δt after A.
uA⊗B 3D3C velocity field u formed by cross-correlating reconstructed
intensity field A with field B.
C0 Intersection of applied constraint sets Cm.
Cm POCS constraint represented by a convex set in Hilbert Space
H where m ≥ 1.
H Hilbert Space.
N Dimensionality of the POCS domain Ω.
Pm Operator function for application of a POCS constraint to
input data where m ≥ 1.
R A subspace ofH.
Nomenclature xix
µ Relaxation parameter for MART and SMART algorithms.
σ Standard deviation of a gaussian approximation to a parti-
cle. In Chapter 4 only, σ denotes standard deviation of a
Gaussian blurring kernel.
θ Boundary Layer Momentum Thickness.
c Column index (in image plane).
k+ Summated Turbulent Intensity.
r Row index (in image plane).
uτ Skin friction velocity.
ui Velocity component in the ith direction.
u′i Fluctuating velocity component in the i
th direction.
X, Y, Z Streamwise, cross-stream and wall normal directions (respec-
tively).
Superscripts
+ Indicating lengthscale normalised by wall-units or velocity
scale normalised by skin friction velocity uτ.
c Corrector (e.g. corrector field used to apply constraint to in-
put).
p Dimensionality of Sobolev space W . p = 2 results in a
Hilbert SpaceH.
s Degree of differentiability within Sobolev spaceW or Hilbert
SpaceH.
Subscripts
i Relates to a cartesian component X, Y, or Z (i = 1 : 3) of po-
sition, displacement or velocity. Exceptions: in section 4.2 i
denotes the pixel index used in the full tomographic recon-
struction. In context of the *RELAX algorithm, λi denotes
relaxation factor used for the ith operator in a POCS recon-
xx Nomenclature
struction.
j Relates to a cartesian component X, Y, or Z (i = 1 : 3) of po-
sition, displacement or velocity except in section 4.2 where j
denotes the voxel index used in the full tomographic recon-
struction.
p Relating to particle position.
Acronyms
CG Conjugate Gradient (matrix solution technique).
NNLS Non-Negative Least Squares (matrix solution technique).
PCG Preconditioned Conjugate Gradient (matrix solution tech-
nique).
ppp Number of Particles Per Pixel (particle density within im-
ages).
*RELAX RELAXation based algorithm used to apply POCS constraints.
3D3C 3-Dimensional, 3-Component [velocity measurements].
AEM Attached Eddy Model.
CFD Computational Fluid Dynamics.
CLM Closed Linear Manifold.
CTE Correlation Tracking Enhancement.
CUED Cambridge University Engineering Department.
DNS Direct Navier Stokes.
GPU Graphical Processing Unit.
HPIV Holographic Particle Image Velocimetry.
HWA Hot Wire Anemometry.
LoS Line of Sight.
MART Multiplicative Algebraic Reconstruction Technique.
Nomenclature xxi
MFG Multiplicative First Guess.
MLOS Multiplicative Line Of Sight [first guess technique].
MPEG Moving Picture Experts Group.
MTE Motion Tracking Enhanced.
PDF Probability Density Function.
PIV Particle Image Velocimetry.
POCS Projection Onto Convex Sets.
PTV Particle Tracking Velocimetry.
PVR Pixel to Voxel Ratio.
SMART Simultaneous Multiplicative Algebraic Reconstruction Tech-
nique.
SNR Signal to Noise Ratio (in the PIV correlation volumes).
TIFF Tagged Image File Format.
TPIV Tomographic Particle Image Velocimetry.
VLSM Very Large Scale Motion.
VODIM VOlume-Deforming Iterative Method.
WRS Weightings Reduction Scheme describing the procedure for
reduction of Wij to Wrs.

1 Introduction
The turbulent boundary layer is an aspect of fluid flow which dominates
the performance of many engineering systems - from aircraft to pipelines,
system efficiencies are governed by the drag exerted on the fluid at the wall.
However, the complex nature of turbulence renders the analytic solution of
such flows intractable for many (most) applications.
FIGURE 1.1 The first conceptual sketch of a ’horseshoe’ (or ’hairpin’) vortex - from
Theodorsen [1952].
Throughout a century of research by the scientific community, models to de-
cipher and describe the turbulent boundary layer have improved from basic
observations of boundary layer behaviour, through descriptions of the types
of turbulence found and conceptualisation of its structural form [Theodorsen,
1952] to the point of implementing advanced predictive models [Marusic
et al., 2010b].
Within the field of study, strong focus is now placed on the concept of coher-
ent structures:
’A three-dimensional region of the flow over which one fundamental flow variable exhibits
significant correlation with itself or with another variable over a range of space and/or time
that is significantly larger than the smallest scales of the flow
- Robinson [1991] -
1
2 Introduction
However, a limiting aspect in our understanding of these flows is our abil-
ity to simulate and measure their properties. Computational simulations
are strongly limited by Reynolds Number and the boundary conditions ap-
plied, whilst experimental measurements (such as dye visualisation, Hot-
Wire Anemometry and Particle Image Velocimetry) are all bounded by their
accuracy, their spatio-temporal resolution and limited dimensionality (i.e.
ability to make volumetric measurements rather than at points or in slice
planes through a flow).
Recently, measurement techniques have been introduced [Elsinga et al., 2006]
which address the limit on dimensionality. Amongst contenders such as
Holographic PIV and Particle Tracking Velocimetry, Tomographic Particle
Image Velocimetry (TPIV or TomoPIV) allows ’3D3C’ measurement of flow
fields - all three components (’3C’) measured within a three-dimensional vol-
ume (’3D’).
Among the ’3D3C’ techniques, Particle Tracking Velocimetry (PTV) is limited
by the density (quantity) of particle tracers (’seeding’) in the flow - strongly
limiting the spatial resolution achievable. Holographic PIV provides excep-
tional spatial resolution but the volume which can be measured is consid-
erably smaller than that available with TPIV or PTV. TPIV offers a trade-off
between volume size and spatial resolution which allows many flows of in-
terest to be measured - such as the turbulent boundary layer.
A variety of studies have been undertaken which use TPIV in a turbulent
boundary layer. In particular, the work of Schro¨der et al. [2011] and Elsinga
et al. [2010] utilise the ability of TPIV to provide topological information:
these authors are able to use analyses not possible with 1D or 2D data (such
as HWA or PIV measurements) in evaluating structural information within
the flow. The literature review (section 2.2) gives a full discussion of differ-
ent experiments of this nature, showing that investigations of the near-wall
region in boundary layers have not been conducted at Reynolds Numbers
above δ+=800. As Reynolds Number increases, the separation (in size) be-
tween the smallest and largest turbulent eddies in a flow increases. Thus, for
a given largest eddy (set by the boundary layer size), increasing the Reynolds
Number requires a finer spatial resolution to capture the smallest features
(vortices, or coherent structures) within the flow. Experiments conducted at
higher Re than δ+=800 have sacrificed spatial resolution (i.e. intentionally
3FIGURE 1.2 The use of Tomographic PIV has facilitated analyses and visualisation
previously only possible with computational, rather than
experimental, results. (left) Conditional averaging based on
topological invariants from Schro¨der et al. [2011] and (right)
visualisation of instantaneous structures within a boundary layer
from Elsinga et al. [2010].
not measured fine-scale information) in order to observe just the large-scale
structures in the flow.
In the same way that increasing Reynolds Number reduces the size of small-
est scale features in the flow, the time-scales of the flow are also reduced
(smaller structures ’turn-over’ in a reduced time). The temporal resolution of
a measurement technique must also improve to capture the smallest scales.
The spatio-temporal resolution of which a measurement technique is capable
is characterised by its ’dynamic range’ - the ratio between largest and small-
est size or time scale which can be measured.
Aspects of spatial and temporal resolution are linked. For example, in a PIV
interrogation, if particles convect by a greater displacement than the win-
dow size then particles cannot be matched in the cross-correlation. If par-
ticles displace too little, then the measurement is within the random error
of the procedure. In the latter case, users must adjust the δt used between
correlation frames to ensure that the degree of displacement is acceptable
for a given window size and level of noise - thereby affecting the temporal
resolution of the experiment. To improve spatial resolution, window size is
usually reduced - decreasing the number of tracer particles within a window
and therefore increasing random error in the experiment. With an increase
4 Introduction
in random noise, the value of δt must increase to ensure valid measurement:
thus a trade-off exists between spatial and temporal resolution. Improving
the accuracy of the cross correlation (i.e. reducing noise) therefore allows a
larger dynamic range to be measured. The dynamic range of TPIV is there-
fore bounded by the experimental error, which is considerably higher than
for conventional 2D PIV - see Worth & Nickels [2007]; Worth et al. [2010] and
Atkinson et al. [2011].
The purpose of this work is to facilitate further analyses of turbulent bound-
ary layers - at higher Reynolds Numbers, with finer spatial resolution (i.e.
improved accuracy with respect to existing literature). Particular attention
is paid to resolving the buffer and lower logarithmic regions of turbulent
boundary layer flow.
Since Tomographic PIV is a relatively new technique, the ’state of the art’ is
still developing. The literature review (chapter 2) considers not only aspects
of theory relating to boundary layers but also the limitations of the tomo-
graphic PIV technique itself - the spatial and temporal resolution, error asso-
ciated with the process, computational cost and post-processing techniques.
Existing experiments using TPIV are summarised.
In order to investigate flows at higher Reynolds Numbers, a higher dynamic
range in the experimental technique is required. Chapter 3 looks in close de-
tail at aspects relating to the dynamic range of experiments in a turbulent
flow, with reference to existing experimental datasets and accuracy studies
discussed in the literature review. Chapter 3 moves on to formulate exper-
imental parameters for a Tomographic PIV setup allowing better resolution
of the near-wall flow and (potentially) an improved dynamic range.
On the basis of considerations in chapter 3, chapter 4 exposes three novel
techniques for improving (primarily) the accuracy and (secondly) the speed
and robustness of Tomographic PIV, with the ultimate aim of improving dy-
namic range of the technique (which can be traded-off for improved spatial
and/or temporal resolution).
Chapter 5 revisits studies relating to the accuracy of Tomographic PIV. Work
is repeated and extended to include and validate the techniques of chapter 4
within both an experimental and numerical framework.
5Chapter 6 extends the scope of the work by introducing new methods for
post-processing the ’3D3C’ Tomographic PIV data. In particular, the method
of Projection Onto Convex Sets is used for restoration of Tomographic PIV
results to meet a-priori constraints. Novel approaches for smoothing, band-
limiting and gap-filling of datasets are demonstrated.
Chapter 7 highlights results from the experimental work, combining the in-
dividual methods of previous chapters in assessing the net performance of
the improvements.
Chapter 8 reviews the completed work and draws conclusions.

2 Literature Review
2.1 TOMOGRAPHIC PARTICLE IMAGE VELOCIMETRY
Tomographic Particle Image Velocimetry (TPIV, or TomoPIV) was first de-
scribed by Elsinga et al. [2006], as a novel technique for ’3D3C’ fluid flow
measurement (measurement within a 3d region or volume of all three com-
ponents of velocity). Their work describes the underlying principle, in which
a region of flow is seeded and illuminated (typically using a laser, as in stan-
dard digital PIV). Multiple cameras are used to image the volume, then a 3D
particle field is reconstructed using tomography. Two particle fields sepa-
rated by time δt are then cross correlated to obtain velocity components on a
grid within the volume. The process is shown diagramatically in figure 2.1.
FIGURE 2.1 Basic steps of Tomographic PIV from Elsinga et al. [2006].
7
8 Literature Review
The tomographic reconstruction is performed by projecting views from dif-
ferent cameras into a discretised 3D volume. Elsinga et al. explore issues
relating to accuracy of the reconstruction using numerical studies, param-
eterising camera layouts and dictating aspects of experimental setup (such
as camera separation angle). In addition, the authors present experimental
data, having used TPIV to visualise turbulent flow behind a cylinder at low
Reynolds Number.
The Tomographic PIV technique benefits from a great deal of the literature
originally intended for 2D PIV. In particular, the multigrid iterative interroga-
tion approach [Scarano & Riethmuller, 2000; Scarano, 2002], the normalised
median outlier test [Westerweel & Scarano, 2005], lessons regarding the accu-
racy of different interrogation techniques [Piirto et al., 2005] and the selection
of optimal sub-pixel interpolation techniques [Lourenco & Krothapalli, 1995;
Roesgen, 2003; Nobach, 2004] all generalise directly to Tomographic PIV -
the result being that many performance aspects of TPIV are documented and
understood better than might be expected for such a new technique.
Tomographic PIV is particularly sensitive to calibration error, since inaccu-
rate calibration causes the reconstruction quality to decrease substantially.
Following Wieneke [2005] which describes calibration procedures for Stereo-
PIV, Wieneke [2008] proposes a method for self-calibration of cameras imag-
ing a volume, which reduces calibration error to an acceptable level through-
out the volume (found by Elsinga et al. to be approximately 0.1 pixels).
At the time of writing, self calibration, multigrid iterative interrogation, win-
dow deformation (using cardinal and/or cubic type interpolation), subpixel-
accurate correlation peak location (using Whittaker/cardinal reconstruction)
and the normalised median test are considered well-proven and have been
incorporated into validated commercial TPIV software.
2.1.1 Bias Error
An aspect of TPIV performance which is increasingly well understood (but
which has not yet been fully addressed) is that of bias error. The phenomenon
is described by Elsinga et al. [2011] - it occurs when so-called ’ghost particles’
(artefacts, or ’shadows’, in the reconstruction process) convect with flow out-
side of each interrogation window (as shown in figure 2.2), causing measure-
2.1 Tomographic Particle Image Velocimetry 9
ments to be biased toward the mean in the presence of velocity gradients
(especially where the direction along which the gradient is taken is normal
to the nominal ’line of sight’ of the camera group). In cases with high shear
in the volume, the entire solution can be biased toward the mean flow.
Ghost particle (blue) moves 
with outer flow, true particles 
(red) move with local flow
A 'Ghost' particle is a 
shadow of multiple true 
particles, where camera 
Lines-of-Sight intersect 
Shear profile in fliud 
(eg when measuring a 
boundary layer)
Time t Time t + dt
Window A x Window B
Motion of ghost causes an 
unwanted additional peak in the 
cross-correlation. Where the 
unwanted peak merges with 
the true peak, 'bias' error ocurs.
FIGURE 2.2 Diagram giving a graphical explanation of bias error.
Biasing error is discussed further by Atkinson et al. [2011], who utilise nu-
merical studies to investigate the reconstruction accuracy. While the authors
also consider noise (see below), the flow field investigated is the mean profile
of a turbulent boundary layer - i.e. a flow in which bias is expected to feature
strongly as a source of error. Atkinson et al. show bias error up to 1.6 voxels
in the worst case.
Novara et al. [2010] propose a scheme - Motion Tracking Enhanced Multi-
plicative Algebraic Reconstruction Technique (MTE-MART) - in which the
tomographic reconstruction is iteratively updated by convecting the recon-
structed intensity field between the two timesteps. Application of the tech-
nique is shown to substantially improve the quality of the reconstruction.
The degree of bias error may be affected to some extent, although no studies
have been reported in the literature to date to confirm or .
10 Literature Review
2.1.2 Random Error
In addition to the bias error (above), a principal source of error in Tomo-
graphic PIV is that of random measurement noise. Several studies [Worth &
Nickels, 2007; Wieneke & Taylor, 2006; Atkinson et al., 2011] have reported
high levels of noise in TPIV data, typically 0.2 voxels (in the measurement
’plane’) and 0.3 voxels (out of ’plane’) although the error can exceed these
values depending on experimental parameters: the noise is shown by Atkin-
son et al. to increase strongly with seeding density.
It should be noted that random errors reported by Atkinson et al. (especially
for cases with window sizes 323) are expected to be lower than in ’real-life’,
since random error typically increases with a reduced number of particles per
window. To ensure comparability between their results, the authors had to
enlarge window sizes in the cross-stream and streamwise directions - thereby
maintaining a consistent but artificially high particle density (i.e. number
of particles per window). Using flat ’tile’ windows in this way is usually
inappropriate for general applications.
The MTE-MART approach of Novara et al. [2010] shows significant promise
in reduction of random error, but as discussed above (in section 2.1.1) the
advantage remains unconfirmed. Even if an expected improvement in ac-
curacy1 is confirmed experimentally, the authors report that convergence of
the approach can take up to 10 iterations of the TPIV process - increasing the
computational demand of the process by a considerable factor.2
Using existing approaches, the random error expected from TPIV measure-
ments is of the order 0.2 to 0.3 voxels (displacement error). In contrast,
the random error typically associated with single-camera 2D planar PIV is
around 0.1 pixels [Raffel et al., 2007; Nobach & Bodenschatz, 2009], substan-
tially lower than from TPIV (or stereo-PIV) measurements. The difference be-
tween these two values comes from several sources: the need to align multi-
ple cameras accurately, the additional step of reconstructing 3D particle fields
and (typically) a reduced number of tracer particles per windows.
1Improvements in measurement accuracy are expected based on the based on the promis-
ing improvements in reconstruction quality demonstrated by the Novara et al. [2010].
2On top of the large number of iterations, the sub-pixel interpolation scheme used to
convect particle fields must be of a high order, and is usually computationally expensive.
Application of the VODIM algorithm in an iterative manner is also a considerable burden.
2.1 Tomographic Particle Image Velocimetry 11
Depending on the flow field (e.g whether or not significant shear is present
in the flow, affecting the bias error), the random error associated with TPIV
can be more or less significant than the bias. However, it can be stated that
random error is typically of a similar order to bias error thus these sources of
error are both significant.
2.1.3 Computational Cost
The tomographic reconstruction process effectively solves the equation Ax =
b, where x represents the unknown intensity distribution in a discretised vol-
ume (illuminated particles resulting in a bright region in the volume). The
matrix A is extremely sparse, but nevertheless is too large to fit in memory
(even for state-of-the-art workstations at the time of writing). In order to
solve the problem, coefficients are calculated on-the-fly during the solution
(although some methods invoke some precomputation in order to accelerate
this process). As a result, direct solution methods are not appropriate and an
iterative approach required.
Elsinga et al. [2006], Petra et al. [2009] and Mishra et al. [1999] discuss a va-
riety of techniques for solution of the problem. Having a limited number
of views (from arbitrary locations) restricts the approach to a real-space to-
mography method (as opposed to fourier based techniques frequently used
in medical imaging). The Multiplicative Algebraic Reconstruction Technique
(MART) was found to be superior by both of the above works - MART is
originated by Herman & Lent [1976] and at the time of writing has become
accepted as a de-facto standard for Tomographic PIV. Petra et al. discuss a
maximum-sparsity formulation which yields potentially more accurate re-
sults than the MART process, but encounter problems related to the degree
of sparsity in the solution for practical particle densities.
Mishra et al. [1999] propose a modification to the MART algorithm - Simulta-
neous Mart (SMART) - which allows computation of pixel-voxel weightings
simultaneously. There is a strong argument in favour of adopting SMART
above MART; either approach leads to a high reconstruction time but the
algorithmic structure of SMART lends itself to parallelisation (across mul-
tiple computational cores), whereas implementations of MART are difficult
to parallelise efficiently (except when running explicitly parallel processes,
12 Literature Review
which results in high memory consumption). Run-times are typically large
for MART reconstruction - especially for large reconstruction volumes such
as those of Hain et al. [2008] and the computing industry is increasingly utilis-
ing multi-core architectures to allow speed improvements, so parallelisation
is an immediate requirement for those aiming to process large datasets.
Worth & Nickels [2008] compute the complexity of the MART algorithm as
O(4ktNcamsNvox[8L + 4]), where kt is the number of iterations, Ncams is the
number of cameras, Nvox is the average number of voxels above the threshold
(across multiple iterations) and L is the average length of the line of sight of
a pixel. The Multiplicative First Guess (MFG) approach proposed by Worth
& Nickels in the same paper precomputes a first guess for the intensity dis-
tribution. This allows iterations of the MART algorithm to take advantage of
sparsity (i.e. the Nvox term is substantially reduced), accelerating the solution
(especially on the first iteration, where a MART computation with uniform
first guess is not sparse at all). The computational complexity (for the MFG
step) is O(Nvox[8Ncams + 1]) and its advantage compared to an uninitialised
MART calculation is shown in figure 2.3. Note the decreasing benefit at in-
creased particle densities (due to decreasing sparsity).
FIGURE 2.3 Performance of MFG algorithm compared to MART (left) with
number of iterations and (right) with particle density [Worth &
Nickels, 2008].
Atkinson & Soria [2009] propose an alternative first guess technique - the
Multiplicative Line Of Sight (MLOS) algorithm. The effect is equivalent to the
MFG approach discussed above (i.e. weighted pixel intensities are multiplied
to form an initial guess). Despite a higher cyclomatic complexity than the
2.2 Tomographic PIV in Turbulent Boundary Layers 13
MFG, Atkinson & Soria appear to create a more computationally efficient
formulation by projection from voxel to pixel space, then interpolation of
the pixel intensity on that line of sight. The authors use the Simultaneous
MART algorithm rather than standard MART for their subsequent iterations.
For a similar quality of reconstruction (MLOS + 40 iterations SMART, vs. 5
iterations MART), the MLOS-SMART algorithm is reported to be up to 4.4×
faster than a direct MART reconstruction. Higher improvement factors are
reported for various configurations (trading off accuracy).
Atkinson & Soria [2009] claim an improvement with respect to the MFG, due
to there being no requirement for assembly of a weightings matrix. How-
ever, for both MFG and MART the weightings matrix need not be assembled
directly, but can be computed on the fly from partially-precomputed coef-
ficients (as is the case with most modern MART codes) which results in a
more efficient implementation. It appears from the times reported in Atkin-
son & Soria, table 5 that the reference is taken using a MART algorithm with-
out such a modification; in which case MFG and MLOS could be considered
broadly equivalent. On the ’zeroth’ iteration (i.e. making comparison with-
out the confusion of different ’downstream’ algorithms) both approaches of-
fer a speed-up of around ×8 (see figure 2.3 and Atkinson & Soria, figure 7b)
- subsequent variation in performance is likely due to the relative implemen-
tation of the iterative algorithm (be it MART or SMART).
2.2 TOMOGRAPHIC PIV IN TURBULENT BOUNDARY LAYERS
2.2.1 Coherent structure in the boundary layer
In the last four decades, a great deal of work on turbulent boundary layers
has surrounded the analysis of so-called ’coherent structures’ - turbulent ed-
dies which retain some self-coherence over both space and time. For bound-
ary layers, these structures have become known as ’hairpin’ vortices.
The concept of a hairpin (or arch) vortex was introduced by Theodorsen
[1952] who developed the idea into a model: ’Vortical ”tornadoes” form astride
near-wall regions of low-velocity fluid and grow outward with heads inclined at 45◦,
and with spanwise dimensions proportional to distance from the wall’
It was pointed out by Offen & Kline [1975] and Hinze [1975] that the dy-
14 Literature Review
FIGURE 2.4 Visualisations of structural information in the boundary layer using
DNS (top) hairpin vortices shown during a numerical study of
structural evolution, [Adrian, 2007, courtesy of K.Kim] (bottom) a
conditionally averaged hairpin structure surrounding ’ejection’ (Q2)
events in a boundary layer from Adrian & Liu [2002]
namics of hairpin structures could be related to the production of turbulence
in the near-wall region. Discussion centred on a ’bursting’ process which
(even to date) has not been precisely defined. Hinze referred to the condition
where low speed, near-wall fluid is lifted between the vortex legs and vio-
lently breaks down (’bursts’) into a patch of high turbulence intensity. This
idea has been lent weight by various DNS studies [Guezennec et al., 1989;
Robinson, 1990; Jeong et al., 1997] with the conditional averaging technique
employed by Adrian & Liu [2002] (see figure 2.4) representing a step forward
in identification of structural form.
Over time, the formation of hairpin vortices has been extensively studied
and modelled [Offen & Kline, 1975; Zhou et al., 1996; Adrian et al., 2000]. An
emergent idea is that of ’autogeneration’, in which passage of a single hair-
pin tends to generate aligned trains (’packets’) of hairpin vortices in its wake.
The zone between the ’legs’ of the hairpins has a characteristically low speed
- leading to the concept of hairpin packets straddling low speed zones. This
2.2 Tomographic PIV in Turbulent Boundary Layers 15
is highlighted in figure 2.6, in which green hairpin vortices (identified using
a ’swirling’ criterion) sit atop a low-speed region of flow (denoted by a blue
isosurface in the same figure. The relationship between large-scale hairpins
(i.e. order of the boundary layer thickness) and so-called ’Very Large Scale
Motions’ (VLSMs) much longer than the boundary layer thickness is still un-
der investigation [Kim & Adrian, 1999; Marusic, 2001; Ganapathisubramani
et al., 2003] but is outside the scope of discussion here.
Study of coherent vortices to date has been based on the canonical case of
a zero-pressure gradient boundary layer at a smooth wall. In general, the
Reynolds numbers studied have been low (Reθ < 5000). Some variations
from this simple case have been made: for a review the user is referred to
Robinson [1991], although the case of wall roughness is considered later. Re-
cent computational studies have reached high enough Reynolds Numbers
for the logarithmic layer of the boundary layer to develop fully [Schlatter &
Orlu, 2010] and good agreement between DNS and experimental data have
been found [Schlatter et al., 2009].
Low-order models are also used to describe the process. The first predic-
tive model, an ’Attached Eddy Model’ (AEM) was proposed by Townsend
[1976], who suggested a random distribution of eddies, growing in a self-
similar way from the wall (and thus ’attached’). This was greatly extended
by Perry & Chong [1982], who assumed an eddy shape and demonstrated
that the Attached Eddy Model could be used to replicate mean flow proper-
ties, broadband turbulence intensity distributions and turbulence spectra.
Fundamentally, the attached eddy model of Perry & Chong consists of an
array of geometrically similar ’Λ’ vortices (figure 2.5). The representative
structures are distributed at different scales according to an inverse power-
law probability density function (i.e. many small-scale vortices, few at a
large scale), and positioned randomly in the streamwise and cross-stream
directions. The Biot-Savart law is used to compute the velocity field, whose
properties can then be compared with experiment. Comparison with experi-
ment in High Reynolds Number flow facilities [Marusic & Perry, 1995] have
shown remarkable agreement, leading toward increasingly robust predictive
measurements [Marusic et al., 2010a; Mathis et al., 2011].
Since an understanding of structural behaviour can demonstrably result in
16 Literature Review
FIGURE 2.5 ’Lambda’ vortex structure, reproduced from Perry & Chong [1982]
low-order predictive models for otherwise intractable high Reynolds Num-
ber flows, as in the case above, it follows that study of turbulent flow topol-
ogy and structure can of great value to the scientific and engineering commu-
nities. Although the canonical case of zero-pressure gradient, high Reynolds
Number boundary layer flow is well modelled by the AEM, the subject of
boundary layer theory - as well as the study of a wide range of turbulent
flows - can still benefit from an improved understanding of topological be-
haviour. Rough-walled boundary layer flow is one of many examples [Maru-
sic et al., 2010b].
Despite being somewhat limited by accuracy concerns (see the error values
reported in section 2.1.1 and 2.1.2), the strength of Tomographic PIV lies in
the direct measurement of flow structures (which were previously inferred
from 2D measurements, dye/smoke visualisations). Computational simula-
tions are able to reach increasingly high Reynolds Numbers, but are still lim-
ited by the processing power of modern computer systems and the boundary
conditions relating to particular scenarios.
2.2 Tomographic PIV in Turbulent Boundary Layers 17
2.2.2 Tomographic PIV Experiments
The number of experiments performed with Tomographic PIV is continu-
ally increasing. To allow comparison between different approaches and to
explore the range of volumes and other experimental parameters for which
Tomographic PIV has been successfully used, parameters extracted from a
range of published studies have been tabulated at the end of this chapter. Of
the experiments shown in the tables [Elsinga et al., 2006, 2007; Schro¨der et al.,
2008, 2011; Elsinga et al., 2010; Worth et al., 2010; Ku¨hn et al., 2011], the most
relevant to this work is that of Schro¨der et al. [2011] whose free-stream ve-
locity and Reθ are in a similar range to those expected in the C.U.E.D. water
tunnel where the present work will be conducted (see section 3.2).
Table 2.1 considers the properties of turbulent boundary layers in several
investigations, highlighting extents of the reconstruction volumes and the
Reynolds numbers investigated.
a. b. c.
Free-stream velocity (m/s) U∞ 9.9 0.53 [Mach2]
Boundary layer thickness (m) δ99 0.024 0.038 0.020
Wall friction velocity (m/s) uτ 0.45 0.0219 19.5
Wall scale (m) + 2.94× 10−5 5.414× 10−5 2.8× 10−6
Wall scale Reynolds Number δ+ 710 800 7080
Momentum Thickness (m) θ 2.9× 10−3 5.5× 10−3 1.2× 10−3
Momentum Thickness Reynolds Number Reθ 1900 2460 34000
Wall-normal measurement range (z+) 34 : 13 : 408 13 : 15 : 328 1062 : 177 : 3363
TABLE 2.1 Boundary layer parameters for existing TPIV experiment data. a.
Elsinga et al. [2007]. b. Schro¨der et al. [2011]. c. Elsinga et al. [2010].
From table 2.1 it is clear that although boundary layers have been investi-
gated ’close to the wall’ as in the case of experiments a. [Elsinga et al., 2007]
and b. [Schro¨der et al., 2011] - the latter resolving down to z+ = 13 with
spacing of 15+ - the investigations have been performed for low Reynolds
Number based on wall scale (δ+). Experiments at higher Reynolds Number
may be beneficial - particularly since structures become less well organised
at higher Re and since models such as the Attached Eddy Model (see 2.2.1)
make the assumption of a large Reynolds Number.
In contrast, experiment c. [Elsinga et al., 2010] successfully resolves a large
proportion of the boundary layer for the purpose of large-scale structural
visualisation - the authors are able to make extensive commentary on the
18 Literature Review
FIGURE 2.6 A conceptual sketch of boundary layer structure at high Re and
experimental evidence from TPIV to support it [Elsinga et al., 2010].
structural composition in the boundary layer, verifying conceptual models
with TPIV data (see figure 2.6). Due to the limited depth of field available
with Tomographic PIV [Elsinga et al., 2006], a very thin boundary layer is
used in order to measure across the range of z/δ. The measurement range
therefore comes at the expense of resolution (in terms of wall units), with the
measurement zone starting outside the logarithmic region.
Experiments have not been performed in a facility large enough to allow
users to visualise the logarithmic region with a good resolution and at a
higher Reynolds Number than δ+ = 800. A good resolution (i.e. fine enough
to visualise the dominant structures in the region of interest) is achieved in
experiments a. and b. by matching the boundary layer size to the capability
of the measurement technique; to increase the Reynolds Number for a given
spatial resolution requires an improvement in dynamic range of TPIV.
The results of Schro¨der et al. [2011] (who use the TPIV data to investigate
Lagrangian particle trajectories and 2-point correlations) are also used by
2.3 Post-processing for Tomographic PIV results 19
Elsinga & Marusic [2010] in a study of flow topology and strain rates in the
boundary layer. Together with the work of Elsinga et al. [2010] which utilises
filtering techniques to observe influence of structures at different scales, these
works demonstrate the wide variety of analyses which can be performed us-
ing Tomographic PIV. The present work will extend such analyses to a higher
Reynolds Number.
2.3 POST-PROCESSING FOR TOMOGRAPHIC PIV RESULTS
The term ’post-processing’ (e.g. of a measured field) can describe a wide
range of operations - such as application of derivative operators, filtering,
extraction of features, point measurements and visualisation techniques.
Post-processing is applied to velocity field results. Here, we limit the defi-
nition of ’post-processing’ to exclude application specific operations and the
determination of derivatives, in order to focus on operations which condition
the field itself in a way which has general application (such as filtering).
Being a relatively new measurement technique, very few post-processing op-
erations have been established specifically for Tomographic PIV. We turn to
the literature for methods of processing images (i.e. 2D scalar fields), aiming
to extend methods to 3D scalar fields then to multiple velocity components
(thus forming post-processing methods for 3D vector fields).
This section aims to apply techniques used in 2D image processing to 3D3C velocity
field results following cross-correlation, in order to improve the final results. This is
not to be confused with pre-processing of the raw particle images discussed in 4.3.
We focus on three general classes of image processing:
Denoising Reduction in noise within the image
Band Limiting Filtering to specifically exclude or retain features of a cer-
tain scale or frequency component
Restoration Repair of images in local regions containing a large error
20 Literature Review
Since TPIV data suffers from a high degree of noise, false vectors occasion-
ally arise - post-processing techniques will ideally minimise the effect of false
vectors on the surrounding (good) data.
In addition to extraction of features of different scales from the measured
flow (coherent vortices), we refer to the discussion on spatio-temporal reso-
lution (see chapter 1 and section 3.2.2). Since it may not be possible to mea-
sure all scales for a desired Reynolds number (even if the accuracy of TPIV
can be improved) data should be handled correctly (i.e. enforcing a scale or
band limit appropriate to the available resolution) in such an event.
2.3.1 Denoising
The practice of image/signal denoising can be broadly split into two cate-
gories: The use of filtering and the use of fit functions. In the former, it is
assumed that the characteristic length scale of the noise is different to that of
the underlying signal. This is dealt with by band-limiting (see below) or local
averaging. In the latter, an underlying functional form is assumed and fitted
to the noisy date(usually by a least-squares approximation).
In using fit functions, care must be taken that the underlying function repre-
sents the nature of the data - for the case of turbulent fluid flow, such gen-
eralised functions do not exist. However, Lourenco & Krothapalli [1995] use
local polynomials (of varying truncation error) to fit components of PIV data.
The authors demonstrate successfully that smooth derivatives can be pro-
duced even with noisy data. The approach is used to form a derivative field
directly, but could be used to determine a de-noised velocity field.
A similar system is utilised by Elsinga et al. [2010], who implement a 3-
dimensional Savitzky-Golay filter. This is a multi-dimensional windowed
regression fit of a polynomial form to the gridded data. At each grid point,
the 5× 5× 5 surrounding region is used to fit a second order 3D polynomial
function from which a smoothed velocity value and field derivatives can be
taken. The technique is shown to be very effective in denoising, but does
band limit the solution somewhat3. Although robust to noise, the sensitivity
3The band limit is noted by Elsinga et al. to be similar in length-scale to the effect of
windowing the cross correlation. Use of the filter is thus a practical option for many purposes
since the spatial resolution of the final measurement may not be significantly diminished.
2.3 Post-processing for Tomographic PIV results 21
of the filter to outlier (false) vectors can be high due to the polynomial nature
of the fit; especially in cases where false vectors cluster close together.
Denoising of scalar fields is more usually done by utilising the frequency or
wavelet (spectral) domain representations of the field, and band-limiting that
representation. Amongst many sources Butterworth [1930] discusses the use
of band limiting filters (in that case, for audio processing purposes).
The band limiting process is done in the frequency domain - more recently,
examples have arisen which utilise the wavelet domain. Wavelet based de-
noising is used for scale-based denoising (wavelet transforms retaining spa-
tial information more effectively than fourier decompositions), although the
principle is very similar to band-limiting. The use of wavelets in denoising is
exposed by Mallat [1999] and Daubechies [1992]. Again, there is a large body
of literature covering many signal and image processing applications.
Wavelets have recently been used in the field of turbulence [Farge, 1992;
Farge & Schneider, 2001; Griebel & Koster, 2000] for simulation purposes,
but this application is not to be confused with denoising: wavelet transforms
can also be used to express a field with greater numerical sparsity than can
be achieved in the Real domain - this characteristic is useful in general, but
not relevant to the present discussion.
Although available and immediately generalisable to 2D and 3D, wavelet-
based denoising techniques are not commonly used in post-processing of PIV
data despite their well-characterised behaviour for similar applications and
large bodies of supporting literature.
However, denoising of velocity fields through Gaussian smoothing is fre-
quently used for PIV applications. In a scalar image, Gaussian smoothing
(blurring) is the use of a localised averaging technique to weight data ele-
ments with an average of their surrounding elements. A 3× 3 Gaussian ker-
nel is frequently used for this purpose - many PIV software programs include
Gaussian smoothing as a standard option during post-processing stages. The
practice is also recommended by Ganapathisubramani et al. [2007] and adop-
ted for visualisation purposes by Schro¨der et al. [2006]. Worth et al. [2010]
note that smoothing of velocity fields in this way can improve the accuracy
of Tomographic PIV due to the reduction in noise. To improve computational
performance, blurring is often applied using spectral techniques, but is not
22 Literature Review
to be confused with band-limiting.
Selection of the blurring kernel is essentially arbitrary. In the case of Worth
et al. [2010] the effect was documented for low Reynolds Numbers (i.e. scale
of flow features large relative to the grid spacing), and no comparison was
drawn with other forms of denoising. [REF BARATH] performed an evalu-
ation of the effect of blurring on flow features whose scale is of the order of
the grid size. It is pointed out by Worth et al. that selecting kernel size below
the cutoff frequency limits the effect on scales which can be captured at the
experimental resolution.
2.3.2 Band and Scale Limiting
Other than the application of band/scale limiting for denoising purposes,
the practice of band limiting images/signals can be used effectively for sep-
arating aspects of the flow field. Coming from signal processing, application
of a band limit allows specific spatial or temporal bands to be excluded or
retained. Translating this principle to the interrogation of turbulent fields
is possible, as we are able to decompose the field into different components
comprising small and large scale features - an example of this is the sepa-
ration of large motions from individual structures in Elsinga et al. [2010].
This is known as scale decomposition4 and is the principle underlying the
work of Farge [1992] for computational applications. Scale decomposition
has not been used (in any universally appropriate manner) for PIV record-
ings to date, although the practice of filtering results is not uncommon for
particular applications - for example, Elsinga et al. [2010] use spatial band
limiting filters
The work of Schram et al. [2004] uses a wavelet decomposition to identify
vortices and educe information about eddy structure in a turbulent flow
measured with PIV (performing the transform on the flow’s enstrophy field).
The technique is extremely similar to scale decomposition; their discrimina-
tor being simply energy rather than scale. However, the approach was not
used in an attempt to post-process the field (in the present meaning of ’post-
processing’) but to identify vortices.
4’Scale Decomposition’ has a closely analogous but more formal definition for use within
wavelet theory - see Daubechies [1992]
2.3 Post-processing for Tomographic PIV results 23
2.3.3 Temporal Filtering
The above discussion refers to spatial filtering; temporal filtering is also pos-
sible and is commonly used in video processing (e.g. in MPEG encoding/de-
coding) as well as in obvious 1D time series filtering applications. In video
processing, temporal filters typically perform a local average of values (at a
point in the image) as the image changes over time. The principle can be
applied to time-resolved PIV fields; locally averaging components in time
rather than space.
The use of temporal filtering is discussed in 4.1.1 and is taken advantage of
by Fore et al. [2005]. In the latter, the authors use a three different nonlinear
filters to detect outliers - including a median filter.
Ve´tel et al. [2011] use temporal filtering for reduction of noise. In addition
to the methods highlighted above (i.e. local averaging and wavelet-based
transforms) an optimal Wiener filter is shown to be effective in removing
noise. The authors also point out that since noise is autocorrelated in the spa-
tial domain, temporal filters are more effective than their spatial analogues
for the purpose of noise removal in PIV recordings.
The most common application of a temporal filter in high speed PIV is the
use of a high-pass filter in order to perform a Reynolds decomposition on the
flow [Baur & Koengeter, 2000].
2.3.4 Restoration
Image restoration is the subject of ’fixing’ images which contain local cor-
ruption; examples include speckled film, torn photographs and compression
artefacts in digital image files. ’Corruption’ can also include non-local ran-
dom noise. In a typical scenario, the unknown (corrupt) region of data is
restored using properties of the image such as frequency content. The sub-
ject of restoration can be generalised to include denoising and band-limiting
processes if required.
First attempts at image restoration were made using an algorithm for band-
limited extrapolation, by Gerchberg [1974] and Papoulis [1975]. Now known
as the ’Gerchberg-Papoulis’ algorithm, it has a variety of applications in sig-
nal processing. The approach is linear - it was noted by Youla [1978] that in
24 Literature Review
a general case, it may exhibit pathologically slow convergence and be poorly
conditioned in the presence of noise.
Youla & Webb [1982] propose an alternative to the Gerchberg-Papoulis algo-
rithm utilising the Method of Convex Projections or Projection Onto Convex
Sets (POCS), which was formulated for faster convergence and robustness
against noise. In a second part to the same piece of work, Sezan & Stark
[1982] utilise the theory of Youla & Webb in the restoration of an image whose
data had been modified (in the spectral domain). Results demonstrated a sig-
nificant improvement of the convex projections method with respect to the
Gerchberg-Papoulis algorithm. The authors did not add noise to the image
so the effect of noise is unknown.
The method was demonstrated by Simard & Mailloux [1988] for use with
a 2D vector field, rather than a scalar image. Most projections (imposition
of a smoothness constraint, say) are trivially extensible to vector spaces by
application to each component of the vector field in turn - the series of pro-
jections proposed by Sezan & Stark can all be applied in this way. In contrast,
Simard & Mailloux utilised both components (u, v) in a projection operator
formulated to reduce divergence. In overview, the nearest field to the in-
put (in a least squares sense) which satisfied incompressibility is determined.
Detailed review of the methodology is made in section 6.2.2 so will not be
exposed here. This nontrivial use of vector information paves the way for
a variety of possible operators which are applicable only to vector spaces.
Simard & Mailloux [1990] investigate the behaviour of the method of convex
projections for a range of 2D vector fields and confirm the performance of
their divergence reduction operator in 2D. To date, the method has not been
applied in measured 3D fields.
In all types of PIV (and especially where results are noisy as in the case of
TPIV), results contain ’false’ vectors caused by spurious peaks in the cross-
correlation. In cases where false vectors are not removed, they can alter the
appearance of local flow topology. Since POCS is frequently used to restore
corrupted (e.g. torn) images, it is suggested that false vectors in the field
could be treated within the POCS framework, although there is no literature
dealing with this application.
Where false vectors are removed, (having been identified by criteria such as
2.3 Post-processing for Tomographic PIV results 25
the Normalised Median Test), the missing data must be filled in (usually by
interpolation). Simard & Mailloux [1990] evaluate the performance of the
divergence reduction operator for restoration of missing vectors, although
this approach has not yet been applied specifically for PIV data.
2.3.5 Summary
The current ’standard’ Tomographic PIV is limited in the spatio-temporal res-
olution (or similarly the measurable dynamic range) of the technique. This
limitation is as a result of noise (to be investigated further in and restricts the
Reynolds Numbers for which the fluid flow can be fully spatially resolved.
Bias error in Tomographic PIV causes a systematic error to be present when
measuring flows containing shear (velocity gradients).
The level of noise in Tomographic PIV is large - up to 3× that of conven-
tional PIV, which is the planar analogy to TPIV (see 2.1.2). This problem is
exacerbated when attempting to take derivatives of the field, as differencing
increases the relative significance of noise.
False vectors arise in PIV results as a result of spurious peaks in the cross
correlation plane (through noise or particle mis-matches). Where these are
successfully identified and removed, accurate methods of filling in the empty
data sites are desirable. Where not successfully identified, it is desirable that
post-processing operators are insensitive to their presence.
At high Reynolds Numbers, even a substantially improved TPIV technique
(or any form of PIV - see Adrian [2011]) may not have sufficient dynamic
range to capture all scales existing within the flow.
This work aims to address these problems in three ways:
- New methods will be developed to improve the noise and bias errors
(and therefore the measurable dynamic range).
- Experiments and numerical studies will be carried out to quantify the
effectiveness of the new methods.
26 Literature Review
- Post-processing techniques for Tomographic PIV results will be for-
mulated, based on those used in other fields such as image process-
ing. Post-processing will be aimed at reducing divergence in measured
fields, denoising results, reducing the impact of false vectors on appar-
ent flow topology and correctly handling cases where spatio-temporal
resolution is insufficient for the flow being measured.
2.3 Post-processing for Tomographic PIV results 27
2.3.6 Parameters of previous experiments
X 40mm (730vox)
Y 40mm (730vox)
Z 10mm (184vox)
δt 35µs
Particle displacement 0.18mm (3.2vox)
Window Size (overlap) [41× 41× 21] (75%)
Deforming Windows? Yes
Type of flow Cylinder Wake, ReD=2700
Seeding Smoke - 1µm droplets
Particle Density 0.05ppp
Laser Pulse Energy 400mJ
Number of cameras 4 [1280× 1024, 12 bit]
Algorithm MART, 5 iterations, µ = 1
TABLE 2.2 Experimental parameters from Elsinga et al. [2006]
X 33mm, 34mm, 42mm
Y 1-12mm, 8-20mm, 0-30mm
Z 26mm, 26mm, 10mm
δt 100µs
Particle displacement 26vox (free stream) - 27vox/mm
Window Size (overlap) [423] (75%)
Deforming Windows? Yes
Type of flow Turbulent boundary layer
Seeding Smoke - 1µm droplets
Particle Density 0.05
Laser Pulse Energy 400mJ
Number of cameras 4 [1376× 1040, 12 bit]
Algorithm MART, 5 iterations, µ = 1
TABLE 2.3 Experimental parameters from Elsinga et al. [2007]
X 32mm
Y 18mm
Z 29mm
δt 200µs
Particle displacement unstated
Window Size (overlap) [483] (75%), 24vox/mm
Deforming Windows? Yes
Type of flow Turbulent spot and tripped boundary layer. U∞
= 7m/s
Seeding Olive oil droplets
Particle Density Not stated
Laser Pulse Energy 21mJ with light ’amplification’ mirrors
Number of cameras 4 [1024× 1024 10 bit]
Algorithm MART, unstated iter no. and µ
TABLE 2.4 Experimental parameters from Schro¨der et al. [2008]
28 Literature Review
X 63mm (734vox) 1380+
Y 15mm (176vox) 328+
Z 68mm (793vox) 1490+
δt 1ms, 333µs
Particle displacement not stated
Window Size (overlap) [323] (75%), 11.6vox/mm
Deforming Windows? Yes
Type of flow Turbulent boundary layer
Seeding 56µm polyamide
Particle Density 0.05ppp
Laser Pulse Energy 25mJ
Number of cameras 4 [1024× 1024 10 bit]
Algorithm MART, unstated iter no. and µ
TABLE 2.5 Experimental parameters from Schro¨der et al. [2011]
X 70mm
Y 3-9.5mm (0.15-0.47 y/δ)
Z 35mm
δt 2.4µs
Particle displacement 20vox
Window Size (overlap) [403] (75%), 23vox/mm
Deforming Windows? Yes
Type of flow Mach 2 Turbulent boundary layer
Seeding 240nm TiO2
Particle Density 0.05ppp
Laser Pulse Energy 400mJ
Number of cameras 4 [2048× 2048, 14 bit]
Algorithm MART, 5 iterations, µ = 1
TABLE 2.6 Experimental parameters from Elsinga et al. [2010]
X 70mm
Y 70mm
Z 11mm
δt not stated
Particle displacement not stated
Window Size (overlap) [323] (75%), 15vox/mm
Deforming Windows? No
Type of flow Isotropic homogeneous turbulence tank
Seeding Dantec dynamics 10µm S-HGS
Particle Density 0.01ppp
Laser Pulse Energy 120mJ
Number of cameras 4 [1024× 1024, 10 bit]
Algorithm MFG-MART
TABLE 2.7 Experimental parameters from Worth et al. [2010]
2.3 Post-processing for Tomographic PIV results 29
X 0.69m
Y 0.42m
Z 0.24m
δt not stated
Particle displacement not stated
Window Size (overlap) [64× 64× 32] (75%)
Deforming Windows? not stated
Type of flow Forced convection in a convection cell 1m in size
Seeding Helium filled soap bubbles 0.2-0.3mm
Particle Density not stated
Laser Pulse Energy n/a (Strobe illumination)
Number of cameras 3 [1280× 1024] + 1 [1376× 1040]
Algorithm MART, 5 iterations, µ = 1
TABLE 2.8 Experimental parameters from Ku¨hn et al. [2011]

3 Experimental Methodology
3.1 DESIGN OF EXPERIMENTS
Experimental work was conducted to form a basis for development, to assess
the effectiveness of and to demonstrate new Tomographic PIV interrogation
and post-processing techniques, developed in later chapters (4 and 6). More
specifically, there were three key aims for the experimental work:
1. To validate the accuracy of Tomographic PIV using independent and
simultaneous measurements (3.1.1).
2. To verify and validate enhancements to the ’standard’ Tomographic
PIV procedure, which consists of a ’Multiplicative Algebraic Recon-
struction Technique’ (MART) reconstruction followed by a ’VOlume
Deforming Iterative Method’ (VODIM) windowed 3D cross correlation.
This general procedure has been established principally by the work of
Elsinga et al. [2006] and Scarano [2002] (3.1.2).
3. To investigate the behaviour of coherent structures in a high Reynolds
Number, fully turbulent boundary layer flow (3.1.3).
Consideration of these items in more detail (see below) led to the following
experimental requirements:
- Simultaneously triggered 2D PIV and Tomographic PIV measurements
- 2D PIV plane must be oriented in plane with the ’depth’ (Z) direction
of the Tomographic PIV setup
- Particle density within suitable bounds for conventional TPIV
- Non-time-resolved data captured for statistical purposes
- Time-resolved data captured to follow evolution of coherent structures
- Flow field used should be a high Reynolds Number boundary layer
- Well-resolved logarithmic region of the boundary layer (at the expense
of not visualising other regions if necessary).
31
32 Experimental Methodology
3.1.1 Validation of Tomographic PIV
The Tomographic PIV method has been validated in several studies (see sec-
tion 2.1). Typically, such validation work falls into three brackets; numerical
evaluation of reconstruction accuracy, numeric evaluation of whole-process
accuracy and comparison with experiment. The use of numerical techniques
to evaluate accuracy of the process will be discussed in more detail in section
5.2. Section 2.2.2 highlights a deficiency of experimental work in the valida-
tion of Tomographic PIV - experiments have typically been performed at a
low Reynolds Number. Worse, comparison between TPIV data and the ’true’
solution is based on previous experiments (such as hot wire measurements)
taken at a different point in time to the TPIV test, opening possibility that
experimental conditions (e.g. temperature, performance of equipment etc)
have changed.
The experiments here aim to rectify the lack of verifiable data by using two
independent measurements, conducted simultaneously in the same flow. In
particular, two dimensional PIV measurements (whose accuracy and robust-
ness are well established for a wide range of flows) are used in conjunction
with Tomographic PIV.
The accuracy of Tomographic PIV is typically reduced along the camera line
of sight, in the Z (’camera normal’) plane as opposed to the other two (X and
Y) directions. This is due to the camera separation angle being limited by the
required depth of field. The effect is of a ’smeared’ reconstruction in the Z di-
rection (Figure 3.1) which affects the correlation peak location [Worth et al.,
2010]. Since the accuracy of TPIV is at its worst in the Z direction, the 2D PIV
system was oriented to the side (looking at the X-Z plane passing through
the centre of the tomographic measurement volume). That arrangement al-
lowed two comparisons to be made: X components (2D PIV vs. ’best-case’
TPIV measurements) and Z components (2D PIV data vs. ’worst-case’ TPIV
measurements).
3.1.2 Verification of enhancements to Tomographic PIV
As discussed in later chapters (see chapter 4), several improvements have
been made to the processing algorithms for Tomographic PIV. Experimental
3.1 Design of Experiments 33
FIGURE 3.1 Elongation (or ‘smearing’) of reconstructed particles in the Z direction
results are used in addition to the numerical studies contained in chapter 5
in order to demonstrate the validity of these improvements.
Before undertaking the experiments, the most important parameter to se-
lect a-priori was the particle density. The Correlation Tracking Enhancement
(CTE) technique (see 4.1) gives rise to a higher signal to noise ratio in the
correlation plane than conventional cross correlation - and a lower bias error.
A higher particle density can be afforded than can be robustly used with a
traditional MART-VODIM algorithm (approximately 0.05 particles per pixel,
as found by Elsinga et al. [2006]). Consideration was given to a second set of
tests with higher seeding density to demonstrate this. However, increasing
the particle density prevents direct comparison between velocity fields pro-
duced by CTE and VODIM cross correlation algorithms, since the VODIM
algorithm no longer works correctly. A particle density of 0.05 (suitable for
application of VODIM) was therefore used throughout all experiments.
34 Experimental Methodology
3.1.3 Investigation of Coherent Structures
Considerable effort has been made in the literature (see section 2.2.1) to in-
vestigate the concept of a ’hairpin’ vortex - a canonical structure which is
replicated at different scales throughout the turbulent boundary layer. Nu-
merical and experimental studies have to date been limited to low Reynolds
Number (see section 2.2.2). Discussion is ongoing with regard to the exis-
tence and significance of hairpin vortices at higher Reynolds Numbers, i.e.
in boundary layers with well established logarithmic regions. Measurement
and simulation of these structures at high Reynolds Numbers is imperative
for improvement and validation of theoretical results and empirical models.
To help clarify this issue, the current measurements also aim to observe vor-
tical structure within a boundary layer at sufficiently high Reynolds Number
for the logarithmic region to be well established. In the event that hairpin
vortices are prevalent (or not) in the results, some discussion could be laid to
rest. It is therefore a requirement that velocity fields produced should be well
resolved spatially in order to visualise the key structures. In section 3.2.1 it
is shown that resolving the logarithmic region comes at the expense of the
ability to visualise the entire boundary layer. Attention was therefore paid to
resolving the logarithmic region, which is of most interest in this case.
At higher Reynolds Numbers, it has been shown (by conditional averaging
of velocity fields) that the hairpin structure concept is statistically useful, re-
gardless of the existence of hairpin structures in an instantaneous sense. To
allow further investigation (especially in the event that instantaneous struc-
tures are not observable), non-time-resolved (’statistical’) data is required in
order to determine well-converged statistical quantities (using time-resolved
data for statistical analyses results in samples which are not independent).
The literature review (see 2.2.1) comments on various proposed mechanisms
governing the formation and life-time dynamics of hairpins. A corollary of
investigating only the near-wall and buffer regions using this technique is
that the formation process of structures in the boundary layer can be ob-
served (due to the high spatial resolution in the region where formation oc-
curs). Since formation is a dynamic process, it was also a requirement that
time-resolved data was captured. The triggering system was therefore de-
signed to facilitate both time-resolved and statistical measurements.
3.2 Experimental Facility and Boundary Layer Characteristics 35
3.2 EXPERIMENTAL FACILITY AND BOUNDARY LAYER
CHARACTERISTICS
Note: The values for free stream speed, friction velocity uτ and spatial resolutions given here
are preliminary parameters used to set up up the experiment. Actual values for these tests
differ. They are computed from the experimental results and are reported in table 5.1
The facility used was the Boundary Layer Water Tunnel at Cambridge Uni-
versity Engineering Department. The tunnel is a recirculating design with an
open channel; i.e. there is zero streamwise pressure gradient. It is specifically
designed for high Reynolds Number experiments, having an 8m long work-
ing section and cross section of 900mm (wide) x 500mm (deep), as illustrated
in Figure 3.2.
FIGURE 3.2 C.U.E.D. Boundary Layer Water Tunnel
The boundary layer is tripped to turbulence at entry to the working section
and can reach a thickness of δ99 = 150mm by the end. The maximum free-
stream speed is 900mm/s.
The arrangement of the experimental setup is shown in figure 3.3.
36 Experimental Methodology
FIGURE 3.3 General arrangement of the Tomographic PIV apparatus
3.2.1 Spatial Resolution
To enable initial specification of the experiment, boundary layer character-
istics were determined using data from 2D PIV measurements previously
taken in the same tunnel by Tee & Nickels [2008]. The measurement position
used for that study was 4m from the trip (halfway along the working section)
and key parameters of the boundary layer are reported in Table 3.1.
The literature review (see table 2.1) together with Raffel et al. [2007, p.239,
Table 7.1] gather data from a range of Tomographic PIV experiments predat-
ing this work [Scarano et al., 2006; Humble et al., 2007; Michaelis et al., 2006;
Schro¨der et al., 2006]. The spatial resolution of those experiments (assuming
for direct comparability that interrogation windows of size 32 voxels3 and
75% overlap were used) varies between 0.44mm and 0.26mm. The experi-
ments had an average measurement volume of 38x38x12mm.
3.2 Experimental Facility and Boundary Layer Characteristics 37
Free-stream velocity U∞ = 0.51 m/s
Boundary layer thickness δ99 = 0.122 m
Wall friction velocity uτ = 0.021 m/s
Wall scale 1+ = 0.049× 10−3 m
Wall scale Reynolds Number δ+ = 2240
Momentum Thickness θ = 0.0073 m
Momentum Thickness Reynolds Number Reθ = 3.7× 103
TABLE 3.1 Boundary layer parameters measured in the C.U.E.D. tunnel [Tee &
Nickels, 2008]
Using this information from previous Tomographic PIV experiments with
the boundary layer parameters in Table 3.1, an initial estimate of the achiev-
able resolution was made: For Tomographic PIV in the CUED Water Tunnel,
vector spacing was estimated to be 5-10 wall units. The wall normal extent
of the volume was estimated to occupy 2.5 ≤ Z+ ≤ 244 with streamwise (X)
and spanwise (Y) extents of 38mm (or 775+).
The estimated reconstruction volume (shown in figure 3.4) therefore spanned
the majority of the logarithmic region (which ranges between Z+ = 40 and
Z+ = 370) with a spatial resolution of 5-10 wall units.
FIGURE 3.4 Estimated reconstruction volume (expressed in wall units) with
boundary layer profile and key layers highlighted
38 Experimental Methodology
3.2.2 Temporal Resolution
The value of δt used between image pairs for PIV must be chosen to ensure
that all components of velocity are accurately measured. A conventional ’rule
of thumb’ for PIV is that the particle displacement in the field is at least 3 pix-
els (or voxels), and less than 1/3 of the window size (i.e. the dynamic range
is limited by window size and vice-versa). The nominal accuracy of PIV has
been studied at length; claimed values of error (in particle displacement) typ-
ically vary between 0.01 and 0.1 pixel (the former for noise-free, simulated
data). Equivalent values for Tomographic PIV are less numerous, but three
studies (Worth & Nickels [2008]; Worth et al. [2010]; Wieneke & Taylor [2006])
have determined accuracy to be within 0.2 voxels (X and Y directions) and 0.3
voxels (Z direction).
Where instantaneous measurement of coherent structures is made, operators
such as the swirling, or ’Q’ criterion, [Cucitore et al., 1999] are typically ap-
plied to identify the existence of a structure independent of the mean flow.
Thus, error must not only be small relative to the velocity magnitude (which,
in the case of a shear layer, is high) but small relative to the magnitude of the
turbulent fluctuations.
For statistical purposes, the higher degree of random error associated with
the Z component simply means that more samples are required to converge
on mean quantities derived from uz as compared to those derived from ux
and uy. However, Root Mean Square values of fluctuating quantities may
require special treatment since error in RMS values cannot be decreased by a
straightforward increase in the number of samples.
The measurement volume encompasses the regions (logarithmic and buffer
layers) in which turbulence production and Reynolds Stresses are highest.
Discussion of Reynolds normal stress distributions for very high Reynolds
Numbers in Schlichting & Gersten [2003, section 17.1.2] suggests that the
streamwise turbulent intensity is highest, with a peak of u′2X/u
2
τ = 7.5 at
around Z+ = 20. This is shown in Figure 3.5. Conversely, the minimum
turbulent intensity exists in the wall normal velocity components. The wall-
normal stresses vary from 0 at the wall, asymptoting towards u′2Z /u
2
τ ≈ 1 in
the outer regions of the boundary layer.
The δt value was selected based on these values. The proposed goal is 10%
3.2 Experimental Facility and Boundary Layer Characteristics 39
FIGURE 3.5 Variation of Normalised Reynolds Stresses with wall distance
(reproduced and modified from Schlichting & Gersten [2003], Fig.
17.6).
error but the problem as-posed is singular: resolving the wall-normal com-
ponent of turbulent fluctuations at the wall within 10% requires zero measure-
ment error. A more realistic target uses the turbulent fluctuations at the mid-
point of the measurement region (Z+ = 122). Taking the averaged square of
the turbulent fluctuations (the Reynolds Stress) as a representative velocity:
u′2Z /u
2
τ ≈ 0.8
|u′Z| ≈ 1.88× 10−2 ms−1 (3.1)
For a target error of 10%, using displacement error (in voxels) eZ = 0.3 and
voxel size of 0.0653mm (determined through calibration and setup of a pre-
liminary test):
δtmin =
6.53× 10−5 · eZ
0.1× |u′Z|
s
= 10.4× 10−3 s (3.2)
Assuming a 323 window size, a similar argument is made to determine the
maximum value of δt from the axial component of turbulent intensity:
40 Experimental Methodology
u′2X/u
2
τ ≈ 8.0
|u′X| ≈ 5.9× 10−2 ms−1
δtmax =
(32/3) · 6.53× 10−5
|u′X|
s
δtmax = 11.7× 10−3 s (3.3)
These arguments suggest that the preferred value of δt is around 10ms. How-
ever, during initial testing stages the MART+VODIM algorithm was applied
for a range of values 1ms <= δt <= 12ms, with a window size of 323 voxels.
Results for tests with δt above 4ms were found to deliver a large proportion of
false vectors. The test at 4ms was somewhat improved, while the test at 2ms
delivered meaningful results (i.e. a plausible vector field). Above, we used
the Reynolds Stress to determine a representative magnitude for the velocity
components which describe a coherent structure. Although well-reasoned
in terms of achieving a required accuracy, this approach clearly does not ac-
count for the sensitivity of the robustness of PIV to the selection of δt.
Note that a further constraint applies to the value of δt based on the displace-
ment gradient, which should remain at less than 0.1 pixels/pixel.
FIGURE 3.6 Power density spectrum reproduced from Davidson et al. [2006].
These values are for a boundary layer at Reθ = 12600, somewhat
higher than the value of 3700 used here. Lowering the Reynolds
number decreases the range of scales over which turbulent
fluctuations occur.
For a more formal consideration of PIV ’robustness’ as a function of δt, we
turn to the work of Nickels et al. [2005], Davidson et al. [2006] and Hinze
[1975, Fig. 7-30]. Amongst many others, these sources demonstrate that un-
steady motion in a turbulent boundary layer exists across a range of scales.
Graphs such as that reproduced in figure 3.6 describe the power spectrum of
3.2 Experimental Facility and Boundary Layer Characteristics 41
the boundary layer in terms of the scale of turbulence. Where the power den-
sity spectrum is nonzero at wavenumbers lower than the spatial resolution of the
PIV system, then turbulent motions must exist within an interrogation win-
dow. Considering the impact of this observation on the PIV process, there
are three relevant scale ranges:
- Scales larger than the the window size are captured by the PIV pro-
cess (although the Nyquist Criterion must also be satisfied in order to
extract spectral information from the result)
- Scales of the same order as the window size substantially affect the
particle position within the window, despite application of a VODIM
cross correlation algorithm
- Scales smaller than the order of the window size have some effect
on the placement of individual particles within a window but do not
significantly change their relative position
In the first case, there is no significant impact on either the accuracy or the
robustness of PIV.
Scales of the same order of the window size critically affect the robustness
of the PIV process by altering relative particle positions within the window.
This has the effect of diminishing the ’true’ peak in the cross correlation
plane/volume and introducing ’false’ peaks1.
Scales significantly smaller than the window size affect the accuracy of the
PIV (typically by shifting and widening the correlation peak) but do not have
a significant effect on robustness since multiple peaks are unlikely to be in-
troduced.
Using Taylor’s Frozen Eddy Hypothesis to relate eddy size (or wavenumber)
with fluctuation frequency, the above argument can be formulated in the time
domain rather than a spatial one - thereby imposing a limit on the value of
δt used. Using values derived from Tee & Nickels [2008] (Table 3.1) at height
1These terms are intentionally placed in quotation marks. In fact, the ’false’ peaks are
true - since they represent actual motion, while the ’true’ peak is true only in the sense of an
ensemble average of the motion within the window (i.e. is based only on the on fluid motion
at larger scales than the window size).
42 Experimental Methodology
Z+ = 122 for eddies at the order of the window size (32/3 <= rvoxels <= 32),
where r is the scale length of the turbulence (in m or voxel units):
uX = 0.63 ·U∞
rm = rvoxels × 6.53× 10−5
δt|max = rm/uX (3.4)
δt|max =
6.50ms where rvoxels >= 322.17ms where rvoxels >= 32/3 (3.5)
Since fully developed turbulent boundary layer flow at these Re consists of
turbulent motions of considerably smaller spatial scales than the window
size, the lower boundary in eq. 3.5 is appropriate. This is consistent with the
experimental finding that 2ms is the largest value of δt which can be used to
produce a robust solution.
All tests were run with δt = 2ms. Convergence in the mean boundary layer
profile was ensured by taking a large number of independent samples (1500)
for the statistical work (see section 5.1). For time-resolved data it was antic-
ipated that an improvement in random error could be brought about using
alternative processing techniques (see chapter 4).
One corollary of this analysis is that velocity fields measured using Tomo-
graphic PIV can be safely low-pass filtered without loss of (or even with an
improvement in) accuracy, providing that the cut-off point of the applied spa-
tial filter corresponds to the temporal resolution effected by the selection of
δt. This is reflected in the work of Worth & Nickels [2008]; Worth et al. [2010]
who noticed in passing that application of a Gaussian smoothing filter im-
proves the accuracy of the results. However, a Gaussian kernel has a similar
but not directly equivalent effect to that of a low pass filter so the problem
merits further investigation and more rigorous treatment.
3.3 EXPERIMENTAL SETUP
3.3.1 PIV Hardware
All laser and camera hardware was triggered centrally using a National In-
struments (NI) digital I/O controller and LabView software.
3.3 Experimental Setup 43
For 2D PIV, a New-Wave Gemini double-pulsed low speed laser was used.
This laser is capable of operation at repetition rates up to 15Hz. The pulse
energy was 200mJ (100mJ/pulse) at 532nm wavelength.
For Tomographic PIV, the laser selected was a New-Wave Pegasus double-
pulsed Nd:Yag (527nm) high speed laser. The laser is capable of repetition
rates up to 5kHz, although the highest used in this testing was 1kHz. The
laser delivers 10mJ/head/pulse. Both heads were triggered simultaneously,
in ’tied’ operation (20mJ/pulse, single head) throughout this work to max-
imise particle intensity in the images.
The 2D PIV setup required a single camera. For this purpose, a Photron Fast-
cam APX was used. This camera is capable of 10 bit image acquisition, has a
resolution of [1024x1024] and capture rate of 2000fps (at full resolution). The
2Gb memory buffer allowed storage of 2048 images before download to disc
was required. The CMOS type sensor has a 17µm pixel size. The camera was
operated at full resolution and triggered (in ’Random’ mode with a sync sig-
nal) using TTL signals with the NI controller mentioned above. At download
time, the middle 8 bits of the 10 bit image were stored in Tagged Image File
Format (TIFF).
The model of camera used for Tomographic PIV was the Photron Fastcam
SA1.1. This camera consists of a 12 bit CMOS sensor with 20µm pixel size,
making it ideal for Tomographic PIV where light intensity is critical. The
maximum operating speed is 5400fps and the memory buffer stores 5047
images. The triggering was performed identically to the Fastcam APX. At
download, images were converted to 16 bit TIFF files to prevent loss of data.
It was later found that the tomographic reconstruction was not sensitive to
the difference between 8 and 16 bit images. Since a large amount of hard
drive space is occupied by these images, this observation can make a signifi-
cant time and cost saving for future users.
3.3.2 Particle Selection
The particles used were Dantec 10µm silver plated hollow glass spheres.
These particles are widely used in water PIV applications and have been
standardised upon for use in the CUED water tunnel. Being highly reflec-
tive, they typically outperform glass or polymer beads of much larger size.
44 Experimental Methodology
This is due to the reduced attenuation of the light path when compared to
larger particles (for a constant particle density, which must be maintained in
order to satisfy the requirements of the Tomographic PIV).
The particle density in the water tunnel was tuned to a value close to 0.05ppp,
by iteratively capturing images and using the nppp() function in the Tomo-
graphic PIV Toolbox (see 3.4.2) to evaluate particle density before making ad-
justment. The particle density typically varies throughout a test but is main-
tained within an error of 0.004ppp for all images in all tests. For a window
size of 323 voxels, this corresponds to approximately 3 particles per window.
3.3.3 Laser Optics
A customised beam expander was developed to illuminate the test volume.
The volume was placed halfway across the tunnel (see figure 3.3), so the
beam path length (through water) was 450mm.
The beam expander consisted of kinematic mirrors for fine alignment, then a
beam expander comprising spherical lenses (to create a fat beam) followed by
a collimating cylindrical lens pair to expand the beam to a fat sheet parallel
to the floor of the tunnel. A knife-edge mask is used to convert the elliptical
beam profile to a rectangular one. This is shown in diagrammatic form in
figure 3.7.
Figure 3.8 shows the arrangement of these expansion optics. Note the align-
ment equipment used at the calibration plate position and at the entry point
of the laser to the water tunnel - these masks serve as ’sights’ when orienting
the laser beam and fixing beam thickness. The laser beam for tomographic
volume illumination entered from the side of the tunnel, preventing any need
for disturbances in the flow (such as a mirror) and minimising the beam path
length to reduce attenuation.
The Gemini laser (for 2D PIV) was fitted with a standard cylindrical lens
arrangement utilising a variable focal length lens. The focal length was set
to 1250mm (the length of the beam path), providing a thin (1mm) laser sheet
at the measurement position. The sheet thickness for 2D PIV was therefore
greater than a single 323 interrogation window for Tomographic PIV (see 3.2).
The sheet width at the measurement position was set to 180mm - larger than
3.3 Experimental Setup 45
Both the spherical and cylindrical beam expanders 
are adjustable to allow control of beam collimation.
A pair of adjustable ('kinematic') 
mirrors allow the beam angle 
and centreline location to be 
continually adjusted - preventing 
the need for fine alignment of 
the entire laser body.
Beam expanded and 
collimated (to approx. 
12mm thickness)
Collimation is achieved using a 'target block'
with markings aligned with the knife-edge 
mask. Beam angle, centreline location and 
collimation are adjusted (in that order) until 
only the marked region of the target block 
is illuminated
Beam expanded and collimated
through cylindrical lens pair to 
illuminate whole volume
Knife-edge mask
Target Block
FIGURE 3.7 Alignment and general arrangement of laser beam optics.
the 100mm wide field of view in order to prevent low-light edge effects.
3.3.4 Camera Optics
The camera and lens combination used for Tomographic PIV must have a
field of view approximately 50mm x 50mm in size. This limit is imposed
by the particle size: imaging particles less than one pixel in size results in
an extremely poorly conditioned tomographic reconstruction. To achieve
such a small field of view, Sigma 180mm macro lenses were used. The line
of sight from reconstruction volume to camera CMOS sensor was approxi-
mately 750mm in length.
Having selected a lens to satisfy the constraint on field of view, the depth of
field available is considered. For successful imaging of the entire volume,
this must exceed the depth of the reconstruction volume. Once lens selection
had been made, the f-stop of the lenses was adjusted to meet the depth of
field requirement whilst maximising the light available to the camera. Typi-
cally, the f-stop value was 11, leading to very dark images - gain and filtering
operations were later applied in order to see the particles.
46 Experimental Methodology
FIGURE 3.8 Customised laser beam expansion optics. The orange target block is
temporarily positioned in the reconstruction volume for alignment
purposes. Note that the laser shown, a Quantronix Darwin Duo, was
not ultimately used due to reliability concerns.
For 2D PIV, the lens (a Sigma 105mm macro) was selected to provide a field of
view approximately 120mm square. This allowed the 2D PIV to capture the
entire boundary layer thickness. Despite larger field of view, the improved
resolution of 2D PIV compared to Tomographic PIV allowed the spatial res-
olution of the two measurement techniques to be comparable.
Although the floor of the tunnel is nominally clear, abrasion due to particle
movement over time has left the perspex surface cloudy. It was decided that
superior optical access was achievable by placing cameras above the water
tunnel, looking through the free surface. The entire system was mounted on
a Thor Labs XT95 rail system, with Manfrotto 110 tripod heads allowing fine
adjustment of camera angle (figure 3.9).
The camera separation angle used was 30◦ from the vertical in both stream-
wise and cross-stream directions. This angle was selected based on a study
3.3 Experimental Setup 47
FIGURE 3.9 Mounting detail for SA1.1 cameras.
by Elsinga et al. [2006, Fig. 6] concluding that the error in reconstruction (due
to shallow angle effects) is minimised at 30◦.
Scheimpflug adapters were used on all cameras in order to orient the focal
plane with the floor of the tunnel, maximising the useful depth of field.
The test volume cannot be viewed directly through the free surface, as motion
of the surface distorts the image. To prevent this effect, a stationary perspex
tray (or ’bird bath’) was placed into the water just touching the free surface.
Additionally to the ’bird bath’, prisms containing water were (at first) used
to prevent the optical aberration caused by looking through the surface at
an angle. This led to a difficulty: adjusting the Scheimpflug angle and mak-
ing small adjustments to the cameras during alignment caused the required
prism angle to change. To prevent the need for an extensive iterative process
to align prism, camera and Scheimpflug angle the bird bath was filled with
stationary water. ’Goggles’ were constructed which attached to the lenses
then partly submerged into the bird bath. The flat bottom of the goggles
48 Experimental Methodology
forms a prism, which is always normal to the camera line of sight - thus a
degree of freedom is removed from the alignment problem. See Figure 3.10.
FIGURE 3.10 Setup of bird bath and goggles to view through the free surface
3.3.5 Calibration
Calibration was performed for 2D and Tomographic PIV apparatus using
the respective calibration routines of the DaVis software (3.4.1) and the To-
mographic PIV Toolbox (3.4.2).
Since the 2D and tomographic measurements were to be directly compared,
the calibration plates were mounted on a plastic template to relate the two
frames of reference.
The calibration for Tomographic PIV was completed using a perspex cali-
bration plate. These plates are transparent, with known dot locations on the
front and back faces. This is inexpensive and considerably easier than trans-
lating a single plane calibration plate through the flow. However, the draw-
3.3 Experimental Setup 49
FIGURE 3.11 Wide-grooved aluminium calibration plate.
back is that the refractive index of the perspex must be corrected for when
making calibration maps. The correction is handled by the Tomographic PIV
Toolbox (see 3.4.2), but does introduce a source of error in cases where the
refractive index is not precisely known.
For the later experiments a solid aluminium calibration plate was designed,
with optimally wide grooves to maximise the depth between planes whilst
maintaining a high camera separation angle (Figure 3.11).
A second step in calibration for the tomographic experiments was to per-
form a particle self calibration. A particle fitting algorithm (which can also be
used for particle tracking velocimetry) was embedded into the Tomographic
PIV Toolbox, following the procedure of Mann et al. [1999]. The disparity-
based recalibration process described by Wieneke [2008] was implemented
and used to reduce calibration errors.
In order for the particle fitting algorithm to be successful, image preprocess-
ing is required: background (’dark’) image subtraction is carried out, then
images are band pass filtered and thresholded. This improves the robustness
of the 2D particle finder, which applies a 3× 3 Gaussian fit to local maxima
in the images. TomoPIV Toolbox pseudocode containing the exact processing
options is presented in Figure 3.12 along with zoomed regions of the input
and output images.
50 Experimental Methodology
%% TomoPIV Toolbox image processing script
% Pre−processes images for input to self−calibration routine.
% SET UP IMAGE PROCESSING PARAMETERS
% Load the dark image set for cameras 1−4
[darkImgCell] = tomoPIV loaddarks(darkImageLocation, [1 2 3 4]);
% Apply the gain settings for each camera
gain = [150 50 50 50];
% Use nominal blur radius of 3 pixels, adjusted by the pixel to
% voxel ratio to give the same 'real space' blur for each camera)
blurRadius = [3 3 3 3]./[0.9273 1.0566 0.8651 1.1614];
% Order of processing operations to be performed
order = { 'dark im corr' ;... % Subtract background images
'bypass' ;... % Band pass filter
'gain and saturate' ;... % Apply gain
'gaussian blur' }; % Apply final blur
% LOOP FOR EACH SET OF IMAGES IN THE TIME SERIES
for iFrame = 0:2047
% Load the raw image set
[imgCell] = tomoPIV loadset(imageLocation, iFrame, cameras);
% Do the processing for each camera
for iCam = 1:4
% Define the image processing options:
imOpts = defineImageOptions(order,...
'darkImage', darkImgCell{iCam}, ...
'lowThreshold', 0.3, ...
'hiThreshold', 3, ...
'gain', gain(iCam), ...
'saturation', 1.0,...
'blurSigma', blurRadius(iCam), ...
'nBlur', 3);
% Call the image processing routine
processedImgs{iCam,1} = imageProcessor(imOpts,...
imgCell{iCam});
end % end for iCam = 1:4
% Save the set of images for visual reference later
tomoPIV saveprocessed(processedImgs, ...
processedImageLocation, iFrame);
end % End for iFrame
FIGURE 3.12 TomoPIV Toolbox script: Image pre-processing for self calibration
3.3 Experimental Setup 51
−0.1 −0.05 0 0.05 0.1
0
5
10
15
20
25
Disparity (R Direction) in pixels
−0.4 −0.2 0 0.2 0.4
0
5
10
15
20
25
30
35
Disparity (C Direction) in pixels
−2 −1 0 1
0
10
20
30
40
50
60
70
−0.5 0 0.5 1
0
10
20
30
40
50
60
−0.5 0 0.5 1
0
5
10
15
20
25
30
35
40
45
50
−0.4 −0.2 0 0.2 0.4
0
5
10
15
20
25
30
35
40
45
−0.4 −0.2 0 0.2 0.4
0
10
20
30
40
50
60
70
−0.4 −0.2 0 0.2 0.4
0
5
10
15
20
25
30
FIGURE 3.13 Histograms of calibration error, using a 7× 7× 4 grid of disparity
bins following self-calibration. The outlier in camera 2 is caused by
poor laser illumination in one corner leading to no true particle
matches within that disparity bin (corner excluded from
reconstruction volume in later testing).
52 Experimental Methodology
−40
−20
0
20
40
−50
0
50
0
2
4
6
8
10
X (mm)
Final Disparity Field (camera 1) Median Disparity = 0.092519
Y (mm)
Z 
(m
m)
−40
−20
0
20
40
−50
0
50
0
2
4
6
8
10
X (mm)
Final Disparity Field (camera 2) Median Disparity = 0.091057
Y (mm)
Z 
(m
m)
−40
−20
0
20
40
−50
0
50
0
2
4
6
8
10
X (mm)
Final Disparity Field (camera 3) Median Disparity = 0.090298
Y (mm)
Z 
(m
m)
−40
−20
0
20
40
−50
0
50
0
2
4
6
8
10
X (mm)
Final Disparity Field (camera 4) Median Disparity = 0.094159
Y (mm)
Z 
(m
m)
FIGURE 3.14 Disparity maps showing final calibration error within the
tomographic reconstruction volume.
The error in the initial calibration (’disparity’) was found to be relatively
large; up to 5.6 pixels for one camera. Such a large error completely dis-
rupts the tomographic reconstruction since the error is larger than particle
size itself: the reconstructed volume comprises entirely ghost particles. The
process of self calibration is therefore viewed as an essential part of the To-
mographic PIV process.
The median calibration error (corrected for by the self calibration process)
was [4.3, 5.6, 3.0, 4.4] pixels for cameras 1-4 respectively. Figures 3.13 and
3.14 show the magnitude and distribution of disparity error within the vol-
ume. The remaining error after 12 iterations of the self calibration process
was [0.093, 0.091, 0.090, 0.094] pixels (median error) - almost an order of
magnitude improvement compared to the initial calibration and within the
acceptable bound of calibration error as specified by Wieneke [2008].
3.4 Processing Software 53
3.4 PROCESSING SOFTWARE
3.4.1 Two-Dimensional P.I.V.
The software used to process 2D PIV data was LaVision’s DaVis v.7.2.2, which
has become a de facto standard software for PIV analysis. DaVis has multi-
pass window deformation PIV algorithms and uses the Cardinal / Whittaker
approach for window deformation and correlation peak location.
MATLAB (with LaVision’s proprietary file format reader as a plugin) was
used as a post-processing tool for further analysis of the velocity fields output
from DaVis (section 5.1.2).
3.4.2 Tomographic P.I.V.
Tomographic PIV data was processed using the Tomographic PIV Toolbox
for MATLAB, an extensive collection of MATLAB based functions written
intentionally for this work.
The toolbox consists of the following features:
- Image processing wrappers for handling and preprocessing of raw To-
mographic PIV images
- 3D Calibration functions, including a Graphical User Interface
- Self calibration capability based on the methods of Mann et al. [1999]
and Wieneke [2008]
- Setup functions for determination of ’case files’ detailing geometric and
precomputed reconstruction parameters
- Tomographic reconstruction routines written in FORTRAN (interfaced
directly to MATLAB) for high performance processing
- A variety of reconstruction algorithms including MART, MFG [Worth
& Nickels, 2008] and WRS (see 4.2)
- Cross correlation functions implemented in MATLAB, FORTRAN and
CUDA for high performance processing
54 Experimental Methodology
- Multiple pass window deformation, using The´venaz-prefiltered Cubic
Spline Interpolation [Ruijters & The´venaz, 2010], discrete window off-
set and cardinal / Whittaker reconstruction
- Gaussian-fit correlation peak location
- Correlation Tracking Enhancement capability (see 4.1)
- Multiple post-processing operators including deformation tensor com-
putation and control volume based computation of circulation.
- Implementation of UNIRELAX algorithm and various operators for
POCS-based post-processing of vector fields [developed in chapter 6]
- Flexible structure and code harnesses for optimisation (see 4.3)
The Tomographic PIV Toolbox is designed to be modular in use; different
stages each comprise an input file (or files) and an output file (or data struc-
ture) with a predefined format. Different stages in the process can therefore
be looped and scripted as the user requires. Particular advantages of this
include the ease with which alternative algorithms can be implemented and
the ability to nest processing functions into MATLAB’s inbuilt optimisers, for
automated tuning of processing parameters (as described in 4.3). The general
procedure followed using the Toolbox is outlined in Figure 3.15.
The Tomographic PIV Toolbox is state of the art in all respects except the
correlation peak location routine. Lourenco & Krothapalli [1995] found that
using a Gaussian-based correlation peak location method introduced a phe-
nomenon known as ’peak-locking’, where the returned location is biased
towards integer locations in the cross correlation volume. The Whittaker
Reconstruction was proposed as a preferable alternative, as it mitigates the
problem by sampling at sub-pixel locations using a band-limited ’Cardinal’
interpolant [Stearns & Hush, 1990]. However, at the time of that publication
the use of window deformation in PIV was not standard. Introduction of iter-
ative window deformation by Scarano [2002] caused the location of the cross
correlation peak to converge on the integer location at the centre of the corre-
lation plane. Hence where window deformation is used (as in the TomoPIV
Toolbox), the peak locking error from a Gaussian fit becomes negligible.
3.4 Processing Software 55
Calibration Images
Sets of Particle Images
Calibration (GUI)
- Carry out interactive dot fitting and an 
automated 3rd order polynomial fit to 
derive a calibration mapping for each 
camera
Calibration Map File
Calibration Plate Details
Image pairs / time series
Self Calibration
- Use a 3D particle fitting code and 
disparity map calculation to improve 
accuracy of calibration maps
Image Processing
- Band pass filter and threshold the 
images to improve performance of the 
particle fitting algorithm
Calibration Map File
Setup File
Case Setup (GUI)
- Determine voxel size and extent of 
reconstruction domain
- Precompute pixel line of sight 
coefficients
Tomographic Reconstruction
- Using MART, MFG or WRS 
Cross Correlation
- Using first or second order accuracy.
- Choice of MATLAB native, CUDA 
(GPU based) or FORTRAN-MEX 
algorithms
Velocity Field Files
Image Processing
- Apply background subtraction, gain 
correction, band pass filters and/or 
nonlinear/custom filters to images, to 
improve the reconstruction quality and 
robustness of the cross correlation
Reconstruction 
Files (temporary)
FIGURE 3.15 Workflow chart for the Tomographic PIV Toolbox for MATLAB

4 Enhancement of Tomographic PIV
In the literature survey (see section 2.1), discussion revealed that dynamic
range can be limited by the degree of error present in the measurement tech-
nique. The survey also revealed that the computational cost of Tomographic
PIV can be high.
Here, a novel approach called the Correlation Tracking Enhancement (CTE) is
described, for improving the performance of Tomographic PIV with respect
to random and bias errors. The aim of the Correlation Tracking Enhancement
is thus to improve the dynamic range achievable with Tomographic PIV by
reducing the error.
A new technique for improving the computational cost of Tomographic PIV
is presented in section 4.2.
In section 4.3, a closer look is taken at the image pre-processing stage, and
an optimisation-based approach is formulated in an attempt to improve the
error and robustness of the Tomographic PIV process.
4.1 CORRELATION TRACKING ENHANCEMENT
4.1.1 Motivation
As discussed in section 3.2.2, the timestep δt required for robust PIV in high
Reynolds Number experiments (such as the boundary layer flow described
in 3.2) is too small to ensure an accurate measurement of turbulent intensities within
the flow.
Since increasing the timestep may not be a viable option for improving mea-
surement accuracy, an alternative technique is sought to allow measurement
within an acceptable range of uncertainty. Although an increase in compu-
tational load must be expected, an important caveat is that the process must
still be feasible (i.e. well within an order of magnitude of the current MART-
VODIM analysis technique).
The requirement of low (or at least reasonable) computational load eliminates
the opportunity to use the Motion Tracking Enhancement (MTE) algorithm
57
58 Enhancement of Tomographic PIV
discussed in section 2.1.1, which increases computational work dramatically.
It is, however, possible to filter results in time. This practise has been adopted
for some datasets [Fore et al., 2005], and is a trivial extension of the spatial
filtering approaches discussed in Chapter 6. Noise apparent in the output
measurements can undoubtedly be reduced. However, as considered in more
detail below (section 4.1.2), filtering cannot account for nonlinearities in the
cross correlation.
The VOlume Deformation Iterative Multigrid (VODIM) method was intro-
duced by Scarano & Riethmuller [2000] and is a standard technique for PIV
analysis. Performance of the algorithm is detailed by Scarano [2002] and
shown to be particularly effective in shear and vortical flows. Essentially, the
principle of deforming the input images improves the ability of the cross-
correlation to detect matches in particle patterns by reducing the relative
motion of several particles within the same window. The use of an itera-
tive technique allows an improvement in the trade-off between spatial reso-
lution and dynamic range of the method1. Originally developed for 2D PIV,
the VODIM technique was immediately extended to 3D (trivial) when the
Tomographic PIV method was introduced by Elsinga et al. [2006], although
for reasons of computational speed some TPIV implementations such as that
used by Worth [2010] retained a simpler approach of window offset rather
than deformation.
It is worth noting that several accuracy studies [Worth et al., 2010; Atkinson
et al., 2011] make use of a discrete window offset algorithm rather than using
deforming windows. Atkinson et al. [2011], who also used a deforming win-
dow algorithm report that a reduction in random error of up to 0.2 pixels is
possible when using a window deformation technique - this is a significant
fraction of the overall error reported, and it is noted that further improve-
ment could be introduced using a higher order image interpolation method
than the linear approach used. An accuracy study (numerical or experimen-
tal) utilising state-of-the-art codes (i.e. incorporating window deformation,
high order sub-pixel interpolation and direct cross correlation) has not been
published at the time of writing.
1Small interrogation windows limit the range of velocities measurable for a given flow,
but give good spatial resolution. An iterative approach with a reducing window size allows
capture of a wider range of window sizes whilst retaining fine spatial resolution
4.1 Correlation Tracking Enhancement 59
Fundamentally, four approaches for improving the measurement accuracy
can be taken:
- Improving the quality of the reconstruction requires either a different
hardware setup (e.g. more cameras for better determination of the vol-
ume, more precise calibration, greater light intensity leading to clearer
images, higher quality optics with less aberration etc) or use of a post-
processing technique, such as application of MTE (tending to reduce
the instance of ghost particles) or band pass filtering (ghost particles
typically having a smaller characteristic size than real particles).
- Improved accuracy of the PIV interrogation algorithm is achieved by
using deformed windows with improved sub-pixel interpolation algo-
rithms [Lourenco & Krothapalli, 1995; Roesgen, 2003; Nobach, 2004]
and direct (rather than the more usual FFT-based) cross correlation [Pi-
irto et al., 2005], bias errors and robustness of the cross correlation pro-
cess can be improved considerably.
- Use of higher-order techniques High order polynomial regression has
been shown [Lourenco & Krothapalli, 1995] to improve accuracy in
scenarios such as evaluating derivatives of the resulting velocity field.
With the exception of work on spatial derivatives by Scarano [2004], lit-
tle attempt has been made in moving to a high order scheme within the
cross correlation itself.
- Increase the seeding density Assuming that the reconstruction can be
made to a sufficient quality, increasing the number of particles in a win-
dow of any given size will improve the definition of the cross correla-
tion peak, leading to improved accuracy.
Addition of cameras and increasing laser power is considered prohibitively
expensive. Every effort is already made to remove optical aberrations from
the experimental setup and self-calibration techniques are used to minimise
calibration error (see 3.3). As discussed above, the MTE technique for en-
hancing reconstruction quality is considered too computationally intensive
to be of practical use at the time of writing. The first option in the list above
(improving reconstruction quality) is therefore discounted for use here.
60 Enhancement of Tomographic PIV
The remaining option to be investigated is a move to a higher order scheme
within the cross correlation. Scarano [2004] performs a direct cross correla-
tion for spatial derivatives; the work shows considerable improvement in the
computation of spatial gradients in the flow, with slight detriment to noise (in
a 2D measurement). For TomoPIV the method proposed by Scarano could be
extremely valuable (especially since window sizes are often relatively large
compared to flow features being measured, thereby affecting spatial veloc-
ity gradients) but the expected increase in noise error may detract from the
present goal, especially since noise is already a principal source of error in
TomoPIV measurements (see section 2.1.2).
Here, we follow the alternative (although not mutually exclusive) path and
move toward a higher-order method in time, in an effort to reduce random
error. This higher order approach forms the basis of the Correlation Tracking
Enhancement (CTE) described below.
4.1.2 The CTE Process
The CTE process derives from an attempt to make the cross-correlation sec-
ond order accurate. The steps in the PIV process are linear, except for one: Re-
construction, window extraction/deformation and discrete normalised cross
correlation are all examples of linear operations; the step of particle peak lo-
cation is not, since the maximum value of the correlation plane/volume is
selected - equation 4.1 highlights the nonlinearity, where a and b are correla-
tion planes.
max(a) + max(b) 6=max(a + b) (4.1)
As discussed in section 4.1.1, the application of a linear filter to the results is a
possibility. For example, an expression for the velocity field evaluated using
second order accuracy is as follows:
uo2=
uA⊗B + 2uB⊗C + uC⊗D
4
(4.2)
Where u denotes a 2 or 3 component velocity field, subscript o2 denotes a
second order solution, A, B, C and D denote four particle fields spaced δt
apart in time and ⊗ is the cross correlation operator.
4.1 Correlation Tracking Enhancement 61
Critically, the nonlinearity in the cross correlation results in a sensitivity of
the measurement to noise and false peaks in the correlation plane. This is the
principle reason why VODIM is generally well conditioned - the deformation
reduces false peaks in the correlation plane, and conditions the shape of the
peak to be more regular - thereby suppressing nonlinear behaviour in the
peak-finding process.
Another key aspect for ensuring high accuracy is the robustness of the tech-
nique, especially where iterative interrogation is used (since a poor result in
an intermediate pass can lead to a corrupted field in subsequent passes). In
this sense, robustness cannot be extricated from accuracy. The attempt to fil-
ter the final time-series of velocity fields by application of equation 4.2 cannot
account for nonlinearity within the cross correlation and therefore has little
effect on the robustness of the VODIM process. Robustness is particularly
important in 3D compared to 2D, due to (in general) a reduced number of
particles within the window and an additional dimension - the former caus-
ing false peaks to be stronger in magnitude, the latter resulting in a greater
number of potential particle matches (and therefore false peaks) in the corre-
lation volume.
From these considerations of sensitivity to noise and of robustness, it is clear
that application of a higher order technique would be most advantageous
when applied before the peak location process rather than after, in order to
stabilise nonlinearity in the peak finder.
In Correlation Tracking Enhancement, a high order scheme is put in place
before the application of the peak finding process, in order to improve the
signal to noise ratio (and shape of the peak) in the correlation volume/plane.
Instead of relying on two volumes, CTE requires three or four successive
reconstructions to make the cross correlation. Here, only operations with
four reconstructions are considered.
The fundamental approach (shown diagramatically in figure 4.1) is to apply
a VODIM algorithm. However, rather than extract deformed windows from
two fields, window locations are convected an additional step forward and
backward in time, to be extracted from four fields - here denoted A (taken
at time t− δt), B (at time t), C (at time t + δt) and D (at time t + 2δt). Note
that to achieve full temporal resolution, 3δt << T where T is the minimum
62 Enhancement of Tomographic PIV
Ghost particles
True particles
Window B Window DWindow CWindow A
(AxB) (BxC) (CxD)
(AxB) . (BxC) . (CxD)
Time t - δt
Time t
Time t + δt
Time t + 2δt
FIGURE 4.1 Windowing pattern for CTE compared to conventional VODIM
algorithm (windows for VODIM outlined in red).
characteristic time period in the flow.
Cross correlation is made between pairs of fields A⊗B, B⊗C and C⊗D. At
this point, three cross correlation volumes are available. Each window pair
consists of the same particle group, extracted and deformed from the particle
field at a different point in time. The three resultant correlation volumes
therefore describe the temporal evolution of the flow observed from a frame
of reference moving with the flow at three points in time (t− 0.5δt, t + 0.5δt,
t + 1.5δt).
With reference to equation 4.3 (in which subscript n refers to the iteration
number within a multiple-pass PIV interrogation), the CTE technique allows
4.1 Correlation Tracking Enhancement 63
a high order evaluation of the time dependent term ∂x/∂t of the material
derivative.
un+1=
Dxn+1
Dt
≡∂xn+1
∂t
+ un∇˙xn+1 (4.3)
Spatial terms are still evaluated to a first order accuracy (the linear predic-
tor methods exposed by Scarano & Riethmuller [2000] are used to convect
window locations). However, since the spatial resolution increases with suc-
cessive iterations of a multi-pass grid, error in spatial derivative terms de-
creases with successive passes in both the VODIM and CTE techniques. In
some sense, then, a multigrid approach is analogous to an interrogation with
high order spatial terms. However, multigrid approaches currently utilise a
fixed time-step. It is therefore beneficial to incorporate higher-order terms in
time by using the CTE technique.
The use of additional particle recordings in the higher order scheme increases
the number of particles imaged for a given measurement2 resulting in a greater
signal to noise ratio (and therefore improved measurement noise and robust-
ness) in addition to the benefits of a higher order scheme.
Bias error resulting from correlated ghost particles is also improved by the
CTE. Ghost particles exist over shorter timescales than true particles [Schro¨der
et al., 2011], so the use of several correlations spanning a longer time period
increases the ratio of correlated-true to correlated-ghost particles.
Elsinga et al. [2011] note that bias error decreases with δt as it is more likely
that ghost particles will become de-correlated. Without strictly increasing the
δt used in each correlation, the CTE has the same effect, since the combination
of correlation planes effectively suppresses peaks arising from de-correlated
particles, whilst retaining peaks arising from true particles (see figure 4.1).
A further improvement in accuracy can be leveraged through use of CTE.
Section 3.2.2 discusses an upper limit on δt which is as a result of robustness;
improved robustness in the CTE technique can therefore facilitate use of a
greater δt without breakdown of the algorithm; improving both bias error
(through the increase of δt as in Elsinga et al. [2011]) as well as an increase in
dynamic velocity range resulting from decreased measurement noise.
2In effect, since the same particles are imaged in a series of independent measurements
64 Enhancement of Tomographic PIV
The effect of CTE is similar to that of the Motion Tracking Enhancement pro-
posed by Novara et al. [2010], which works by removal of ghost particles
that are not moving with the local flow (i.e. are de-correlated). By iterative
deformation of reconstruction volumes, Novara et al. use a volume derived
at time t (here denoted B), convect it according to the velocity field to time
t+ δt and use it to correct the existing reconstruction at t+ δt (C). This effec-
tively improves reconstruction quality by suppressing ghost particles.
The CTE technique can be used to improve spatial resolution (which is cur-
rently limited by particle density as discussed in 3.2.1) in two ways. Firstly,
use of twice the number of windows allows use of a reduced number of par-
ticles per window, thereby reducing the minimum window size3. Secondly,
having a reduced sensitivity to ghost particles, the CTE technique can be
used with denser particle fields (and therefore smaller window sizes) than a
typical MART-VODIM setup for a given acceptable level of bias error.
The additional benefit made by the CTE technique comes in computational
performance. At worst, i.e. in the case where a single correlation is made at a
time, use of the CTE technique doubles the processing requirement. At best
(where a time series of data is being processed, so the series of reconstruc-
tions is made anyway) the CTE technique can be applied to the final passes
of a cross correlation process only, minimising computational impact4. This
is a substantial improvement with respect to MTE, which requires succes-
sive applications of VODIM as well as convection operations applied to the
reconstructed particle fields.
Discussion should be made regarding the application of the higher-order
scheme. Attentive readers will have noticed that the final operation indicated
in figure 4.1 is a dot product, inconsistent with implications of the discussion
above. At the point where three cross correlation volumes have been com-
puted, two different operations may be applied:
3although care should be taken if utilising this practise, since an additional source of error
is introduced as the centroid of the imaged particles is less likely to coincide with the window
centre
4Application of CTE rather than VODIM for a cross correlation pass slightly more than
doubles the computational load for that pass - double the number of window fetches and
triple the number of correlations are made
4.1 Correlation Tracking Enhancement 65
⊗ABCD=A⊗ B + 2B⊗ C + C⊗D4 (4.4)
⊗ABCD=A⊗ B . B⊗ C . C⊗D (4.5)
The former expression (equation 4.4) is strictly a second order operation.
However, the maximum theoretical improvement in signal to noise ratio is
×3 for that additive process. In the presence of noise and where ghost parti-
cles are prevalent, the latter (equation 4.5) may be used since multiplication
provides greatly improved suppression of uncorrelated ghost particles and
noise - the improvement in signal to noise ratio for real, experimental data is
typically an order of magnitude (see figure 4.2).
2 4 6 8 10 12 14 16 18 20
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
Signal to Noise Ratio (Correlation Plane)
PD
F
 
 
VODIM, 413 final pass size
VODIM, 323 final pass size
CTE, 413 final pass size
CTE, 323 final pass size
FIGURE 4.2 PDFs showing the signal to noise ratio of correlation volumes in a
VODIM algorithm (top) and in the same field with coincident
window locations using the CTE method (bottom). PDFs are shown
for four applications of the PIV algorithm, having final-pass window
sizes of 41, and 32 voxels3.
Novara et al. [2010] include an investigation to compare results between mul-
tiplicative and additive MTE approaches (remembering that MTE is applied
in the real domain rather than the correlation plane). It is reported that the
multiplicative approach performs well at the first iteration but leads to a
poorer solution over several iterations due to slight errors in the convection
scheme causing cancellation of real particle intensities. For CTE, the same
problem does not occur: all three pairs of windows are offset based on the
same velocity field. Error in the window deformation process is consistent
between the three correlations and cancellation of the peaks does not occur.
Any error is then compensated for by subsequent interrogation passes in the
normal way.
66 Enhancement of Tomographic PIV
It is left to the user to select whether the algorithm is used in additive or
multiplicative form. For this work, the multiplicative form was used since,
unlike the multiplicative MTE, there is no apparent disadvantage.
An outstanding issue having impact on accuracy of the MTE technique and
potentially any cross correlation for Tomographic PIV is that smoothing of
the Tomographic reconstruction volume (between successive iterations of
MART) has been found to result in improvement in the accuracy [Michaelis,
2011; Discetti & Astarita, 2011]. It is speculated here that this is due to the
filter acting to extract particles (true particles being Gaussian in nature, simi-
lar to a smoothing kernel). Due to there being a nonzero error in the velocity
field used for convection, the additive nature of the MTE technique results in
a substantial smoothing effect (using the particle field itself as a kernel!). The
influence of smoothing should be investigated independently from the appli-
cation of MTE or other technique, since potential improvements are available
for little (almost no) computational cost.
4.1.3 Performance of the CTE
Figure 4.3 shows correlation volumes (from experimental data rather than
simulated) with the shape of the correlation peaks highlighted as isosurfaces
of magnitude (normalised relative to the maximum magnitude in the corre-
lation volume). The sub-pixel accurate peak location process is affected by
the shape of the peak, which is clearly better defined for the CTE cases - this
is also highlighted by the misshapen peak shown (for the 41 voxels3 case)
in figure 4.4. Figure 5.11 and its commentary show that the CTE technique is
particularly powerful in reducing random error in the uz component of veloc-
ity (i.e. normal to the ’plane’ of measurement). Figure 4.3 gives good insight
into why this is - in the correlations computed using VODIM, the shape of
the peak is elongated in the vertical direction (even for the larger window
size) leading to less precise determination of that component. This is largely
mitigated with CTE.
In general, the CTE process improves accuracy by allowing a better defined
peak shape in the correlation domain - this is due to an improved signal
to noise ratio (SNR) in the correlation plane and leads to a more accurate
subpixel peak location process. Here, SNR is defined as the ratio between
4.1 Correlation Tracking Enhancement 67
FIGURE 4.3 Correlation volumes at the final pass in a multigrid cross correlation.
Volumes a. and c. are for VODIM, with interrogation windows 413
and 323 voxels in size respectively. Volumes b. and d. are for CTE,
using the same interrogation windows of 413 and 323 voxels
respectively. Windows are taken at the same point in the same field
for both CTE and VODIM. Isosurfaces of magnitude (normalised by
peak magnitude) are plotted with different opacities, showing the
misshapen peak for application of VODIM to smaller window sizes.
68 Enhancement of Tomographic PIV
FIGURE 4.4 Surface plot of correlation peak, sliced through the max location
within correlation volumes 41 voxels3 in size (as in figure 4.3). The
colour scale is selected to highlight the deformed shape at the base of
the peak which affects the subpixel location accuracy.
peak and mean value in the correlation volume. Figure 4.2 shows the clear
improvement in signal to noise ratio for CTE compared to VODIM using a
probability density function of the SNR.
4.2 ACCELERATION OF TOMOGRAPHIC PIV BY WEIGHTING
REDUCTION
4.2.1 Background to the Weighting Reduction Scheme
Since the accuracy of Tomographic PIV can be improved by schemes such
as the CTE technique (section 4.1), it is desirable to do so. However, such
schemes are typically associated with an increase in computational cost com-
pared to the standard MART-VODIM approach. In some cases, this is mit-
igated by the particular experimental requirements. For example, the time-
series experiments described in 3.1 are imaged at a constant period, so ap-
plication of the CTE technique requires no additional reconstructions. The
non-time-resolved experiments require groups of 4 (rather than 2) image sets;
double the number of reconstructions in order to use CTE.
4.2 Acceleration of Tomographic PIV by Weighting Reduction 69
Even ignoring the potential for improvement using enhanced methods like
CTE, application of the MART algorithm is highly computationally intensive.
Effective parallelisation of the algorithm is achievable (either by performing
batches of several MART reconstructions simultaneously, or by adopting the
SMART algorithm of Mishra et al. [1999]). Nevertheless, reconstruction of a
single volume typically takes several minutes on modern workstations5.
Even for the shortest of experiments, the number of reconstructions ranges in
the 100s. For larger experiments (e.g. attempting to achieve converged flow
statistics or observing structures over a longer time period), the computation
time for Tomographic PIV becomes prohibitive - a user may have to wait
weeks or even months to complete the processing.
The literature review (2.1.3) discusses various methods for improving com-
putational cost. The Multiplicative First Guess of Worth & Nickels [2008] and
the similar MLOS algorithm of Atkinson & Soria [2009] improve the speed of
computation by a factor of ×8. Despite this advantage, there is still a large
processing overhead and it is desirable to find an alternative to the MFG/M-
LOS + MART/SMART approach which is more computationally efficient.
As described in the next section (4.2.2), the underlying cause of this prob-
lem is the size and storage format of the weightings matrix W, which is pro-
hibitively large and cannot be stored in computational memory in a form
conducive to direct solution. It is also observed that in the case of reconstruc-
tion for PIV, the problem is characterised by ’compound sparsity’ - current
solution approaches takes partial but not full advantage of the high degree
of sparsity in the problem.
Here, a reduction method is described (tantamount to preconditioning W)
which allows the programmer to take full advantage of the high degree of
sparsity in the system. The weightings matrix W may now be retained in
memory in its entirety - in a form suitable for use with the wide range of
computationally efficient matrix solvers already available in software such
as MATLAB, Intel’s Math Kernel Library and NVIDIA’s CUSPARSE library
for GPU computation.
5The computer used for this work was equipped with 8 cores (dual Intel XEON Harper-
town processors at a clock speed of 2.83GHz), 12Gb DDR2 ram (clocked at 667MHz) and a
SAS hard drive for rapid disc access.
70 Enhancement of Tomographic PIV
4.2.2 Sparsity of the weightings matrix
As described by Elsinga et al. [2006], the tomographic reconstruction is a so-
lution of the equation system Wx = p where x is a vector describing the
unknown voxel intensities, p is a vector containing the intensity distribution
recorded by the cameras (the ’pixels’ array) and W is a weightings coefficient
array which records the contribution of intensity in each pixel to the intensity
in each voxel (i.e. a weighting of 1 indicates that the pixel looks directly at
the voxel, a weighting of 0 indicates that the pixel cannot see the voxel).
For Tomographic PIV, the problem is highly underdetermined - in the present
experiments, 4×1MP cameras are used to observe a reconstruction volume
700 × 600 × 160 voxels in size. Thus the degree of determination is 4 ×
106/67.2× 106 = 0.059. This is for a volume of typical depth - in some cases
the degree of determination is even worse.
The Wij matrix (where subscript i refers to a pixel, and j refers to a voxel)
therefore consists of 22.6× 1014 elements, requiring a total of 1× 106Gb of
storage in full form. However, the matrix is highly sparse. Pixels (projected
into physical space) are approximately the same size as voxels (for optimal
reconstruction performance - see Elsinga et al. [2006]). The line-of-sight (LoS)
for each pixel is computed and the pixel weighted with those voxels sur-
rounding the LoS through the volume. Most TPIV codes (including the To-
moPIV Toolbox) take a 3× 3 region of voxels surrounding the intersection of
the LoS at each plane in the voxels array. For the case of the present experi-
ments, each pixel is therefore weighted with 160× 9 voxels. There are thus
4× 106 × 160× 9 = 5.76× 109 nonzero coefficients in W.
The sparsity in W is therefore very high - it is 99.9978% sparse. Still, the
coefficients require approximately 21Gb of memory (at single-precision). The
indexing array required to store the matrix in sparse form is in at least 32 bit
integer precision for modern computer systems, so takes another 21Gb.
Sparse problems whose weightings matrices are too large to store are well
studied (CFD problems are examples of these) and the MART approach is
one method for solving the problem. For general problems, there is no way
around the issue and weightings have to be computed at run-time within
the solver (or stored on disc and retrieved element-by-element with the large
time-lag associated with optical/magnetic hardware).
4.2 Acceleration of Tomographic PIV by Weighting Reduction 71
As currently adopted, the MART/SMART approaches are suitable for a gen-
eral intensity distribution in the pixels array. For example, if every pixel
were non-zero, the MART approach would still work (an example of this
is in scalar field reconstruction).
For the specific case of PIV, the inputs themselves are sparse. A typical To-
mographic PIV image (shown in figure 4.5) is 58% sparse.
FIGURE 4.5 Typical image (left) used as input to Tomographic PIV, shown with
zoomed central region (right). Pixels array from this image is 58%
sparse.
Consider the case that we do not compute the voxels which are weighted
with a zero-level pixel, assuming the intensity in that voxel also to be 0. The
size of the weightings array is reduced (in this case) by almost 60%.
In a typical implementation of the MART algorithm, an if-statement is in-
cluded to take advantage of pixel array sparsity by skipping calculation of
the MART equation for zero level pixels. However, weightings must still be
computed for all zero level pixels since the voxels which they affect must be
set to zero6.
A significant computational benefit can therefore be achieved by utilising a
method which does not require weightings of zero-level pixels to be com-
puted at all.
6A further disadvantage is that if-statements inside loops tend to adversely affect the
speed of execution, since the branch itself requires clock cycles and typically prevents some
compiler optimisations.
72 Enhancement of Tomographic PIV
Sparsity in the pixels array has a further implication. Consider the case of
a pixel which causes a zero-level to occur in a voxel. That voxel may be
weighted with several other pixels (e.g. neighbouring pixels, or ones from
other cameras). In the case of a MART algorithm, computation must be made
for the other pixels despite the knowledge that being weighted with a zero level
pixel, the intensity of this voxel must tend toward zero.7 In order to ascertain the
ultimate value of the voxel, determination of weightings between it and any
other pixels are not necessary - neither is the actual computation. Here, these
unnecessary weightings are known as redundant coefficients.
Voxels array (idealised to 2d)
Camera 2 - pixels array (idealised to 1d)
C
a
m
e
ra
 1
 -
 p
ix
e
ls
 a
rr
a
y
 (
id
e
a
li
se
d
 t
o
 1
d
)
Nonzero voxels (both cameras)
Voxels nonzero in Camera 1, zeroed by Camera 2
Voxels nonzero in Camera 2, zeroed by Camera 1
Nonzero pixels (both cameras)
FIGURE 4.6 Voxel and pixel arrays for an idealised tomographic reconstruction, to
demonstrate the Wrs principle. Here, dimensionality is reduced to 1D
images and 2D voxel domain for simplicity. Weightings take values
[0, 1] in this example, but in general can be any non-negative value in
the range 0 <= w <= 1.
By elimination of these redundant coefficients from the weightings matrix,
the already high sparsity of W can be substantially increased, to the point
where the entire matrix can reasonably be held in computer memory. An
example (using 1D pixel arrays and a 2D voxel array, with binary weighting
values) is made in figure 4.6 to demonstrate the principle.
7By inclusion of an if-statement to determine voxel intensity in the inner loop of the MART
algorithm, this computation can be somewhat mitigated. Algorithmically, however, this con-
sists of a retrieval of the value from the voxels array ( 2-300Mb memory retrieval) and a branch
conditional statement within the inner loop, both of which are computationally intensive in
their own right.
4.2 Acceleration of Tomographic PIV by Weighting Reduction 73
4.2.3 Weighting Reduction
We define a reduced weightings array, still denoted W but utilising sub-
indices r and s. In order to reduce Wij and sort the kept coefficients into
the smaller matrix Wrs, two logical masks are required:
δir The same size as the pixels array. Contains logical true where pixels are
nonzero, false otherwise.
δjs The same size as the voxels array. Contains logical true where all non-
zero weighting coefficients for a voxel relate to pixels whose intensity
is nonzero, false otherwise.
In theory, the reduction process consists of application of the two logical
masks in succession:
Wis = Wij ∀δjs = true (4.6)
Wrs = Wis ∀δir = true (4.7)
In practise, Wrs is computed directly using the logical masks, since Wij is too
large to construct.
Voxels array (reduced using camera 1)
Camera 2 - pixels array (idealised to 1d)
C
a
m
e
ra
 1
 -
 p
ix
e
ls
 a
rr
a
y
 (
n
o
 lo
n
g
e
r 
sp
a
rs
e
)
Voxels array (reduced using camera 2)
Camera 2 - pixels array (no longer sparse)
C
a
m
e
ra
 1
 -
 p
ix
e
ls
 a
rr
a
y
 (
id
e
a
li
se
d
 t
o
 1
d
)
FIGURE 4.7 Elimination of redundant coefficients, using camera 1 (left) and
camera 2 (right).
74 Enhancement of Tomographic PIV
Voxels array (reduced by both cameras)
Camera 2 - pixels array (no longer sparse)
C
a
m
e
ra
 1
 -
 p
ix
e
ls
 a
rr
a
y
 (
n
o
 lo
n
g
e
r 
sp
a
rs
e
)
Voxels and pixels remapped to full, 
dense storage in computer memory
FIGURE 4.8 Following elimination of coefficients from the weightings matrix,
retained pixel and voxel arrays are full and nonzero. The reduced
weightings matrix Wrs remains highly sparse (figure 4.9), although is
in general less sparse than the input Wij.
Figures 4.6 to 4.8 illustrate the use of the pixel and voxel domains in con-
structing the logical masks. To show the effect of masking on the coefficients
matrix, figure 4.9 shows the weightings arrays for this simplified example.
Computationally, the δjs mask is constructed by projecting pixel lines of sight
through the voxels domain. The logical operation can be made either by pro-
jecting only the zero-level voxels, or by only the nonzero voxels. The process
is closely equivalent to that detailed in Worth & Nickels [2008] for making a
Multiplicative First Guess (MFG). Since only a subset of the pixels must be
projected, determination of the mask takes at most 50% of the computational
effort required to perform the MFG.
Once the masks are computed, the size of the Wrs matrix is known (nnz(δir)×
nnz(δjs)). The matrix is initialised in sparse form then lines of sight for the
nonzero pixels are computed from calibration coefficients. These are used to
determine weighting coefficients and thus populate the Wrs array.
4.2 Acceleration of Tomographic PIV by Weighting Reduction 75
Original Wij coefficients
Coefficients retained in W
rs
Coefficients eliminated due to zero pixels in C1
Coefficients eliminated due to zero pixels in C2
Coefficients eliminated due to zero pixels in both C1 and C2
0 50 100 150 200 250
0
20
512  nonzero coefficients, 93.75% sparse
0 50 100 150 200 250
0
20
304  coefficients eliminated
 
 
0 20 40 60 80 100 120 140 160
0
10
20
262  nonzero coefficients, 93.75% sparse
FIGURE 4.9 (top) Weightings matrix Wij corresponding to the reconstruction
problem in figure 4.6. (middle) Weightings matrix coloured to show
redundant and kept coefficients for the present example. (bottom)
Resultant Wrs matrix.
4.2.4 Degree of determination in the Wrs matrix
The assembly of Wrs is similar in computational effort to the process of mak-
ing a Multiplicative First Guess. In the case of the TomoPIV Toolbox, coeffi-
cients for lines of sight are precomputed to further improve efficiency of the
process (i.e. polynomial evaluations of calibration functions are not required
within the computation). In the absence of a high overhead for weighting
computation, the MFG approach is used as a basis (instead of MLOS) due to
its lower cyclomatic complexity [Atkinson & Soria, 2009].
This reduction technique has a curiosity seen in the matrix of figure 4.10. The
matrix has more rows than columns; i.e. it is overdetermined. To explain this
physically, consider that each particle is ≈3 voxels in diameter (in 3D space)
76 Enhancement of Tomographic PIV
FIGURE 4.10 Sparsity pattern of Wrs for a small (3003 domain size) experiment,
with 0.03 particles per pixel and a characteristic particle size of 4
pixels. Although the matrix has no classical structure, the three
horizontal bands correspond to the three cameras
4.2 Acceleration of Tomographic PIV by Weighting Reduction 77
and≈3 pixels in diameter (assuming that it occupies 9 pixels) when projected
into 2D. Using formula for volume of a sphere, the particle occupies 14 or 15
nonzero voxels; compared to Ncams × 9 pixels (27 in this case), where Ncams is
the number of cameras used. Thus more pixels ’see’ the particle than there
are voxels in it. The problem must be overdetermined and this is verified by
the dimensions of the above matrix.
In contrast, cases containing a higher particle density (or larger particles) can
result in an underdetermined Wrs matrix as seen later in figure 4.13.
In general, the arrangement of particles within the volume also contributes to
the degree of determinacy. For systems like these with a sufficient number of
randomly arranged particles, the degree of determinacy in the whole system
is insensitive to this concern.
An overdetermined problem has implications for the solution method, since
the structure of the required solver will be very different and for accuracy of
the reconstruction, since overdetermination theoretically allows perfect re-
construction (assuming no noise in the measurement and that no noise or
artefacts are introduced by discretisation of the problem). It is possible to take
advantage of overdetermination by increasing the voxel size (relative to pixel
size) somewhat, by application of a sharpening filter to the input images (i.e.
making particles smaller) or by increasing the threshold below which pixel
intensities are considered zero (i.e. by artificially increasing sparsity of the
pixels array, sacrificing definition at the edges of reconstructed particles).
4.2.5 Solution of Wrs
Detailed discussion of matrix solution methods is outside the scope of this
work. The key requirement is that a non-negative constraint is applied to the
solution - since it is nonphysical for light intensity to take a negative value.
A variety of methods which allow such a constraint have been used to solve
Wrs and results have been found to be consistent, although methods vary in
convergence rate, run-time and memory overhead:
NNLS The non-negative least squares algorithm described by Law-
son & Hanson [1995] and implemented in the MATLAB inter-
nal function lsqnonneg().
78 Enhancement of Tomographic PIV
CG An extremely simple Conjugate Gradient (CG) solver imple-
mented in MATLAB. Code structure is given in figure 4.11.
PCG A preconditioned conjugate gradient solver (variation of CG).
SMART The SMART (or MART) techniques can be trivially implemented
to utilise Wrs instead of Wij. Ostensibly, the implementation
is the same as for the non-reduced problem. However, since
weightings are now stored in a quickly accessible manner and
since accessing a column of the matrix directly gives the indices
of weighted pixels used within SMART, the numerical imple-
mentation of SMART becomes much more computationally ef-
ficient. Weightings are no longer computed at all (following the
first guess) and the iterations are performed as sparse-matrix-
multiply operations, for which most modern Mathematics Li-
braries (e.g. Intel’s Math Kernel Library) contain efficiently
vectorised and massively parallelised routines. Moreover, opti-
mising compliers such as the Intel Fortran Complier are able to
prefetch coefficients for many matrix and scalar multiply oper-
ations, which streamlines memory bandwidth demands com-
pared to existing reconstruction codes.
The last three solvers are suitable for both under-and-over determined cases.
For overdetermined cases only, the first solver provides excellent ’out-of-the-
box’ performance since MATLAB’s routines are parallelised.
Since the solver implementations here are made for the purpose of demon-
strating validity of the technique - and since the goal of this work is not to
produce commercial standard code - no aggressive optimisations have been
made. There is scope for a wide variety of improvements on the prototype
code presented here, many of which are able to be implemented in a straight-
forward way:
- Implementation on a massively parallel architecture, such as a GPU
using the CUSPARSE library.
- Implementation within FORTRAN utilising the Intel Math Kernel Li-
brary for handling sparse matrices - even with the simple CG technique
this is likely to result in considerable speed improvement.
4.2 Acceleration of Tomographic PIV by Weighting Reduction 79
function [x] = nnlslisting(Wrs, nzPixels, nIterations)
%NNLSLISTING Example solver for the Wrs*nzVoxels = nzPixels problem
% Transpose to multiply
wsr = wrs';
% Initialise x as x0>=0
x=zeros(size(wrs,2),1);
% Initialise the x old previous step solution (arbitrarily high ...
value)
x old = ones(size(x)).*100000;
% Initialise a vector of residuals
residual = zeros(nIterations,1);
% Majorise curvatures
curvatures = wsr*(sum(wrs,2));
% Iterate solution
for i=1:nIterations
% Gradient of curvatures
gradient = wsr*(wrs*x−b);
x = x−gradient./curvatures;
% Apply non−negative constraint
x(x<0)=0;
% Compute residual. Normalise whilst excluding NaNs (/0 errors
% otherwise).
mask = x old ˜= 0;
residual(i) = sum(abs(x(mask) − x old(mask))./abs(x old(mask)));
% Save solution for next iteration
x old = x;
end % End iteration loop
end % End main function
FIGURE 4.11 TomoPIV Toolbox Function: Simple conjugate gradient (CG)
least-squares solver with non-negative constraint for sparse
matrices, implemented in MATLAB for testing purposes.
80 Enhancement of Tomographic PIV
function [x] = smart(wrs, nzPixels, nIterations)
%SMARTLISTING Example solver for the Wrs*nzVoxels = nzPixels problem
% Relaxation factor (0 < mu < 2)
mu = 1;
% Start with unity intensity distribution
xOld = ones(size(wrs,2),1);
% Solution P i
P = nzPixels;
% Do iterations of the SMART algorithm
for smartIterCtr = 1:nIterations
% Using wrs instead of wij
A = wrs*xOld;
% Projection term at the current iterate
B = P./A;
% Loop SMART iteration − parallelisable loop!
for voxCtr = 1:size(wrs,2)
% Figure out which pixels see this voxel
mcjMask = wrs(:,voxCtr)>0;
% Get the number of pixels weighted w/ current voxel
mcj = nnz(mcjMask);
% Current correction terms
correction = B(mcjMask).ˆ(mu*wrs(mcjMask,voxCtr));
% Product term
product = prod(correction).ˆ(1/mcj);
% Update the intensity
x(voxCtr) = x(voxCtr)*product;
end
disp([' smart.m: Completed SMART iteration ' ...
num2str(smartIterCtr) ' of ' num2str(nIterations)])
% end SMART iteration
end
FIGURE 4.12 TomoPIV Toolbox Function: The SMART algorithm is trivially
implemented with parallelised sparse matrix libraries for solution of
the Wrs matrix. Here, a solver is shown in prototype MATLAB code.
4.2 Acceleration of Tomographic PIV by Weighting Reduction 81
- Investigation of more appropriate preconditioners and direct solution
techniques for the problem.
- Implementation of sparsity maximisation solvers may improve quality
of the output.
- Current implementations rely on MATLAB’s native sparse array ob-
ject, which stores values in double precision. Since the precision is not
required to be greater than single, a customised sparse object could re-
duce storage demand by up to a factor of 2, reducing demand for mem-
ory bandwidth from the solver.
4.2.6 Performance of Wrs with a CG solver
It is worth tabulating results for the Wrs approach on real data. We use two
different cases to compare performance with MART:
Case 1 A 4-camera setup (used to produce experimental results in section
5.2) with relatively sparse seeding (0.02 particles per pixel). Particles
are of characteristic size 3.5 pixels. This leads to an overdetermined
weightings array.
Case 2 The same 4-camera setup with dense seeding (0.05 particles per
pixel). Particles are also of characteristic size 3.5 pixels. This leads
to an underdetermined weightings array.
Since the effectiveness of the Wrs technique is limited by sparsity of the solu-
tion, the speed of computation (and the ultimate ’capacity’ in terms of parti-
cle density) can be controlled using a thresholding value on the input images.
In the cases used here, a threshold of 12 counts (12/255 for normalised im-
ages) was applied - assuming that particles have a Gaussian intensity profile,
this corresponds to a cut-off at approximately 5% of the peak intensity.
The Wrs matrices for the two cases are shown in figure 4.13. Computational
times are quoted for single-core execution for comparability between systems
- the first guess technique is easily parallelisable8. Memory storage require-
8The first guess is a simpler version of the MFG requiring two steps. In the first, binary
masks are computed and in the second, weightings are retrieved to assign the contents of
82 Enhancement of Tomographic PIV
ments are quoted for single precision weighting values and 32-bit sparse ar-
ray indexing.
FIGURE 4.13 Sparsity patterns of Wrs for a four-camera setup. (top)
Overdetermined matrix for a 0.02ppp case, occupying 53.1 Mb in
memory and computed in 23.6s. Reconstruction is 98.97% sparse in
the voxels domain and 58.28% sparse in the pixels domain. (bottom)
Underdetermined matrix for 0.05ppp case occupies 895.3 Mb in
memory. Its reconstruction is 82.32% sparse in the voxels domain
and 27.15% sparse in the pixels domain.
Note an additional benefit that increasing the number of cameras increases
the sparsity of the voxels domain (for a given particle density). Not only
is the quality of the reconstruction improved, but the speed is actually im-
proved too, as the Wrs matrix reduces in size.
the Wrs matrix. Both steps are order independent since the MFG is not actually made, but a
binarised equivalent is used. Matrix computation can thus be implemented in an efficiently
parallelised way.
4.2 Acceleration of Tomographic PIV by Weighting Reduction 83
Case 1 Case 2
Number of Rows in Wrs 383735 6584601
Number of Columns in Wrs 473867 827493
Size of Wrs (Mb) 53.1 895.3
Sparsity of Wrs 99.9963% 99.9957
Pixel Array Sparsity (%) 58.3 27.2
Seeding Density (ppp) 0.02 0.05
Wrs Computation time (serial) 23.6 s 55.5 s
TABLE 4.1 Parameters for the two cases used to assess performance of the Wrs
technique.
We consider the performance of the technique for the case of lowest particle
density9. Figure 4.14 shows the variation of Q-factor (defined by Elsinga et al.
[2006] as a metric of reconstruction quality) and Normalised Median Resid-
ual (taken over the nonzero elements of the voxels array) in the solution with
Time. Despite a sub-optimal implementation, the Wrs algorithm is quicker
to converge toward a reasonable value of Q than the MART implementation.
The peak Q-factor achieved is lower than that which can be achieved with a
MART algorithm by 2.1% for this case. The Q factor following the 5th itera-
tion of MART is 0.771, which is the same as that reached by the CG solver -
i.e. for practical purposes there is no difference in Q factor between a MART
implementation and a CG solution.
Figure 4.15 shows a subdomain of the reconstruction, made with MART and
with Wrs. Isosurfaces of intensity are plotted at E = 0.02, where E is the
unnormalised intensity field resulting from reconstruction. The domains are
broadly similar, but the Wrs technique appears somewhat more susceptible to
ghost particles despite almost identical Q-factors between these comparisons.
In fact, the intensity retained by the CG solver is somewhat higher than by
MART - hence additional artefacts are apparent when isosurfaces are plotted
at the same intensity.
Taking into consideration the implementation of solvers - where in this case
the CG solver is wholly unoptimised and built from prototype code in MAT-
9Similar performance is achievable for both low and high density test cases; however
only the lowest density case was used to evaluate Q factors of the reconstruction for practical
reasons relating to the forward-projection algorithm required to determine the input images.
84 Enhancement of Tomographic PIV
LAB, where the MART algorithm is a finely tuned function written in FOR-
TRAN, it is difficult to provide a reliable estimate for the ultimate relative
performance of the two. Reported speed-ups for re-implementing MATLAB
code in FORTRAN range between 2 and 10,000 from a wide variety of sources
- it is reasonable to suppose that a full implementation of the Wrs technique
can be substantially quicker than MART, although the speed up is likely to
be at the lower end of this range10.
Since all techniques can be effectively parallelised, serial timings are pre-
sented throughout for comparability.
0 50 100 150 200
0.5
0.55
0.6
0.65
0.7
0.75
0.8
0.85
0.9
Time (s)
Q
 
 
QMART
QWRS
0 50 100 150 200
0
0.2
0.4
0.6
0.8
1
1.2
1.4
Time (s)
N
or
m
al
is
ed
 M
ed
ia
n 
Re
sid
ua
l η
 
 
ηMART
ηWRS
FIGURE 4.14 Quality factor and Normalised Median Residual plotted against
Time for a Conjugate Gradient (Wrs) solver - one of many possible
solvers which can be implemented based on the Weighting
Reduction Scheme - and for a standard MART solution. The CG
solver is implemented in unoptimised form in MATLAB compared
to a highly-tuned MART code based in FORTRAN.
One remaining concern with the Wrs scheme is that sparsity in the image ar-
ray is a key factor in the viability of the reduction process, as dense image
arrays result in very little reduction in the number of coefficients. This is
somewhat irreconcilable with the general aim within the TPIV community
to move toward denser particle fields and therefore improved spatial res-
olution. Increasing camera resolution whilst decreasing particle size is the
only way of reconciling these aims with the use of the weighting reduction
scheme.
10MATLAB’s sparse matrix libraries are already well-optimised, and the code above is
partly vectorised to improve the performance of the MATLAB-based function.
4.3 Image Pre-processing for Tomographic PIV 85
FIGURE 4.15 Central 32 voxels3 of a domain reconstructed using MART (5
iterations, µ = 0.99) (left), and Wrs with CG at 250 iterations, at the
same Q factor as iterations of MART (right). Blue spheres indicate
’true’ particles.
4.3 IMAGE PRE-PROCESSING FOR TOMOGRAPHIC PIV
4.3.1 Background
Before using recorded images in a tomographic reconstruction, an image pre-
processing stage is performed. Image-preprocessing is a common technique
for 2D PIV but is of greater importance in Tomographic PIV, since it is pos-
sible to significantly alter the characteristics of the reconstructed field using
image pre-processing.
There is no literature survey exhaustively detailing the effects of image pre-
processing on the reconstruction. However, is is noted in Elsinga et al. [2006]
that ’The reconstruction process is improved by means of image preprocess-
ing with background intensity removal, particle intensity equalisation and a
Gaussian smooth (3×3 kernel size)’. These image processing techniques have
become accepted as necessary for a well-conditioned reconstruction without
substantial or systematic investigation regarding their effect - with the excep-
tion of background subtraction, whose performance was
Currently, the image processing operations applied are arbitrary, in the sense
that there is no consistent methodology used to ascertain the optimal set of
86 Enhancement of Tomographic PIV
processing parameters. Background subtraction can be applied in a consis-
tent way (subtracting the average of a time series from each individual image
in the series). However, determination of the radius of a smoothing kernel or
scale of a band pass limit (for example) is left to the whim of the user - these
are typically determined by trial and error.
Here, a consistent methodology is formulated for determination of optimal
image processing parameters.
4.3.2 Performance of the Reconstruction
There are three fundamental problems associated with the reconstruction
process which pre-processing may mitigate:
- The particles to be reconstructed are small relative to the voxel size
(typically 3 voxels in diameter). Their shape and centre of intensity
therefore cannot be well-defined following voxelisation.
- Once projected into the reconstruction volume, the line of sight of a
pixel is weighted to affect the voxels through which it passes and their
neighbours within the reconstruction. Lines of sight corresponding to
zero (or low) intensity pixels form the boundary - a cage - around each
reconstructed particle. Since lines of sight are weighted to affect neigh-
bouring voxels, the reconstructed particle is ’sharpened’ (inversely di-
lated/blurred) by the application of MART11.
- The existence of ’ghost’ particles affects the quality of the cross correla-
tion. Since the number of ghost particles increases with seeding density,
their presence limits the ultimate spatial resolution of Tomographic PIV.
Figure 4.16 highlights these problems. Using an experimental calibration
(which had been self-calibrated as discussed in section 3.3.5) to encompass
the likely degree of error associated in the reconstruction problem, a single
particle was reconstructed from artificial images with no noise present. No
image processing or thresholding is applied. It is clear that for particles of
the size considered in most PIV experiments (second from top in figure 4.16),
11This is in addtition to the inherent sharpening effect of MART, which arises through
repeated multiplication of a gaussian distribution with itself.
4.3 Image Pre-processing for Tomographic PIV 87
FIGURE 4.16 Reconstruction of a single particle at a range of characteristic sizes.
(left) Zoomed particle images in 4 cameras, (right) Reconstructed
particle shown green, with an ideal particle corresponding to an
idealised particle in blue whose characteristic size corresponds to
that in the generated images. (top to bottom) Particle characteristic
diameter 2, 3, 4, 5, pixels (in first camera - size corrected for pixel to
voxel ratio in the other camera images to maintain consistent
reconstructed particle size). Axis ticks denote 1 voxel distance.
88 Enhancement of Tomographic PIV
the shape can be significantly deformed by the discretisation of the problem
and the particle decreased in size. This is a large source of the measurement
noise associated with Tomographic PIV (discussed in 3.2.1), particularly at
high resolutions (small interrogation window sizes) where there are very few
particles per window.
The application of band-pass filters, gaussian smoothing and background
subtraction of the input images tends to mitigate all three of the above con-
cerns. Smoothing dilates particles in the image plane - they become larger
relative to the voxel size, addressing the first two problems above. Back-
ground subtraction and band pass filtering eliminate low-intensity particles,
noise and DC levels in the input images, increasing the sparsity of the pixels
array (as discussed in 4.2) with a corresponding reduction in the number of
ghost particles present.12
4.3.3 Selection of Image Processing Parameters
There are three methods which can be used to ascertain appropriate parame-
ters for use in image processing:
- Trial and error (the status quo).
- Parameter Searches - vary parameters systematically. Perform image
processing, reconstruction and cross correlation for each unique set
of parameters and evaluate the result based on a performance metric.
Choose the set of parameters with the best metric.
- Use an optimisation algorithm to vary parameters in an automated
way. At each iteration of the algorithm, perform the image processing,
reconstruction and cross correlation and evaluate the result based on a
performance metric. The algorithm adjusts parameters to minimise (or
maximise) the value of the metric.
There is any number of potential parameters which could be used (corre-
sponding with an infinitum of possible filters). For simplicity, this study
ignores the use of highly nonlinear filters such as local normalisation and
12Although a trade-off exists: the process of smoothing / blurring the images decreases
image sparsity, increasing the number of ghost particles.
4.3 Image Pre-processing for Tomographic PIV 89
background subtraction - these filters have well defined purposes (brightness
equalisation and DC level cancellation respectively) and their application is
not dependent on (or sensitive to, at least) ’tuning’ factors. This study, then,
is limited to the following key parameters: gain, low-pass lengthscale, high-
pass lengthscale, blur radius, intensity threshold and saturation level. Each
of these can be applied separately for different cameras - for a 4 camera setup,
this requires 6× 4 = 24 variables.
Since reconstruction and cross correlation are computationally expensive,
this list must be parsed down - a design space search of 24 independent vari-
ables could take months.
Firstly, it is observed that lengthscales can be linked. The Pixel to Voxel Ratio
(PVR) is a quantity used in any Tomographic PIV code to relate the pixel size
to the voxel size, when objects are projected from the image plane into the
reconstruction domain. Elsinga et al. [2006] suggest that the optimal PVR is
approximately unity; within real experiments the PVR of each camera typ-
ically varies as 1 ± 0.1. Defining length scales (blur radius and band-pass
scales) in voxel distance units, they are related to image (pixel) units for each
camera through that camera’s PVR. Thus the Blur, Low-pass and High-pass
parameters reduce to a single variable for all cameras.
Secondly, the low-pass lengthscale is held constant since the nature of the
bypass filter used in the Tomographic PIV Toolbox requires an odd integer
value >1. Since particles in the raw images were approximately 3 pixels, a
value of 3 was chosen.
The saturation level (maximum intensity allowable) is dependent on the gain
- decreasing the saturation point is equivalent to increasing the gain. Since all
images imported into the TomoPIV Toolbox are normalised by the bit count
to vary in the range [0, 1] a constant (default) saturation value of 1 was used.
Except where very high thresholds are used, variation of the intensity thresh-
old (minimum intensity allowable - values below this are set = 0) has little ef-
fect on the shape of reconstructed particles (although as discussed above, in-
creasing threshold typically reduces size and number of ghost particles found
in the reconstruction).
Sensible ranges for parameters were defined based on trial and error - ranges
are intended to be extremes (i.e. obviously incorrect values at both ends of
90 Enhancement of Tomographic PIV
the range, with appropriate values between). Value ranges are listed for re-
maining parameters in table 4.2. Note that for different experiments, ranges
can vary substantially - for example, the uppermost values of gain are high
in this case since 12 bit cameras were used (requiring 16 bit filetype), and
output .tif files saved from the cameras using the lower bit registers - upon
normalisation (by 65,535 for a 16 bit image) all intensities appear low.
Gain (Cam 1) 1 ≤ g1 ≤ 1000
Gain (Cam 2) 1 ≤ g2 ≤ 1000
Gain (Cam 3) 1 ≤ g3 ≤ 1000
Gain (Cam 4) 1 ≤ g4 ≤ 1000
Low Pass Length-scale 0.1 ≤ slp ≤ 1
Blur Radius 0 ≤ σ ≤ 3
TABLE 4.2 Value ranges for image processing parameters.
This reduction of parameters is sufficient for automated optimisation pur-
poses. However, to facilitate full design space searches, it is possible to match
intensities between cameras - i.e. take the median intensity in Camera 1 and
scale the gain of other cameras to reach the same median intensity. A his-
togram matching algorithm can equally be used.
4.3.4 Metric of performance
In order to perform an automated or design-space search for ’best perfor-
mance’, a performance metric must be defined. Clearly, it is desirable to im-
prove the accuracy of the ’output’ field, so an RMS value of the difference
between measured (including image processing) and true fields would be
appropriate.
However, the intention is to develop a practical technique for improvement
of experimental results - for which the true fields are not known a priori. An
alternative metric must be chosen. In this case, a stringent false vector crite-
rion was applied to the output field (the normalised median test developed
by Westerweel & Scarano [2005]) and percentage of false vectors contained in
the field used as a score. While suitable for demonstration of the technique,
this metric has some pitfalls, discussed in 4.3.5.
4.3 Image Pre-processing for Tomographic PIV 91
4.3.5 Design Space and Simplex Search Methods
σ
Lo
w
 P
as
s 
Le
ng
th
sc
al
e
260
0.25
0.26
0.27
0.28
0.29
0.3
0.31
0.32
0.33
0.34
0.35
 (Gaussian Blur)
Gain (Camera 1)
1
1.5
2
2.5
160
180
200
220
240
50% false
60% false
70% false
Lo
w
 P
as
s 
Le
ng
th
sc
al
e
FIGURE 4.17 Isosurfaces of performance in parameter space
An initial attempt was made to investigate the effect of altering parameters
using a simple design space search. The additional simplification (mentioned
above) of matching camera intensities was taken, and the parameter space
searched systematically. Figure 4.17 shows isosurfaces of performance index
for different combinations of image processing parameters.
Experimental data was taken from the tests described in 3.3; i.e. particle im-
ages taken in the logarithmic region of a high Reynolds Number turbulent
boundary layer, using a high-speed laser and Photron SA1.1 1MP cameras.
The isosurfaces plotted actually show very poor performance of the TPIV
process - 50%, 60% and 70% false vectors identified using a normalised me-
dian test with standard parameters (recommended by Westerweel & Scarano
[2005]). However, the ’peak’ performance is not captured in the region of de-
sign space investigated. A natural step is to extend design space to the point
where isosurfaces of performance become complete hulls, then locate a peak.
It was found, however, that the search was poorly-conditioned. output image
varied widely in dynamic range i.e. two cameras had dim, smoothed images,
92 Enhancement of Tomographic PIV
the other two brighter and sharper13. This caused the histogram matching
algorithms (both an iterative algorithm from the MATLAB Image Processing
Toolbox and a straightforward scaling correction were tried) to over-saturate
particles in some images. In that case, an increasing low-pass length scale
gave rise to large, apparently smooth particles and very low sparsity in the
images. The result was that a large number of ghost particles were present in
the reconstruction, having profound effect on the accuracy and robustness of
the cross correlation. The design space search was halted.
The next step required a design space search incorporating the individual
values of gain for each camera - a 6 parameter search. Since each variable
should be evaluated at (say) 10 points, a design space search (using a hexa-
hedral grid of points) requires 106 evaluations; prohibitive to any computer
system other than a very large cluster. In order to proceed, an automated
optimisation algorithm was used.
In selecting an optimisation approach, a number of requirements were con-
sidered:
Non-Locality Optimisations which descend to the closest local minimum
in a design space may not be appropriate - figure 4.17 does
not show evidence of the design space having local min-
ima, but the characteristics of the space may alter with the
addition of 3 further dimensions.
Efficiency The algorithm must converge quickly, since each iteration
consists of a tomographic reconstruction and cross correla-
tion.
No derivatives Derivatives of the scoring function (the score at each point
in parameter space) are unreliable due to the discrete na-
ture of aspects of the problem. Moreover, calculation of
a derivative function typically requires many evaluations,
affecting the requirement for efficiency.
13This is not as a result of the mie-scattering effect as input images are normalised. Which
two cameras are dimmed in the output of the optimisation depends on the starting guess for
image processing parameters
4.3 Image Pre-processing for Tomographic PIV 93
Nonlinearity While highly nonlinear behaviour is not expected, the pro-
cess of cross correlation is nonlinear so it must be assumed
that the scoring function is also nonlinear.
Constrained The ability to apply constraints on variable values (see ta-
ble 4.2) is desirable for some cases - however for the ranges
defined here, constraints should not be required since the
optimisation reaching the range boundaries is a sign of the
scoring function being improperly bounded.
Two classes of optimiser were considered which had this capability; a Nelder-
Mead Simplex Search [Nelder & Mead, 1965] and a Genetic algorithm. Using
MATLAB’s Optimisation Toolbox, a first attempt found that the Genetic algo-
rithm required far too many function evaluations (to be effective, the initial
population size had to be > 100 so even a few generations takes consider-
able time). Worse, a poor convergence rate is a well-known aspect of Genetic
algorithm performance, so this technique was discounted due to computa-
tional effort required. A simplex method was implemented which proved
suitable - although not strictly global, simplex methods are non-local enough
that highly local minima (usually caused by discontinuities in the score func-
tion) are handled by the routine, which is sufficient in this case.
The Tomographic PIV Toolbox was extended to include a function harness
for optimisation in this way; see figure 4.18 for details. A prerequisite of the
algorithm to be implemented was that an initial guess is supplied. Here, the
initial guess for values of gain was chosen such that the images appeared
similarly bright without saturating; the low-pass filter size was taken to be
0.3 pixels and the blur radius set at unity.
Figure 4.19 shows convergence of the simplex search in terms of the variation
of parameters over 104 iterations. In this case, no bounds were supplied -
i.e. the converged values reflect an optimal solution, rather than artificially
imposed boundaries. The proportion of false vectors was improved from
approximately 25% at the start down to <5% at convergence. Figure 4.20
shows the vector field before and after the image optimisation process.
The algorithm is clearly successful in reducing the number of false vectors,
the qualitative ’condition’ of the vector field and in converging on a solution
without manipulation of the bounds or constraints. However, the objective
94 Enhancement of Tomographic PIV
function imageOptimisation
% −−> Initialise − load required image and calibration files
% Initialise arrays for output and storage of each step
maxIterations = 300;
nParameters = 6;
% Define the processing filters to be applied (in order)
order = {'dark im corr' ;...
'bypass' ;...
'gain and saturate';...
'gaussian blur' };
% −−> Define options used in tomographic reconstruction
% −−> Specify settings used in cross correlation (VODIM algorithm).
% Define control options for optimisation algorithm
optimOpts = optimset( 'Display', 'iter',...
'FunValCheck', 'on',...
'MaxFunEvals', maxIterations,...
'TolX', 0.01);
% Set the first guess and run optimisation
x0 = [100 100 100 100 1 0.3];
[x,fval,exitflag,output] = fminsearch(@evaluationFcn,x0,optimOpts);
% −−> Save results and plot iterations / convergence
% NESTED EVALUATION FUNCTION
function [score] = evaluationFcn(input pars)
% Sort the input parameters vector and scale lengths by PVR
gain(1:nCams) = input pars(1:4);
blur(1:nCams) = input pars(5)*pvr(1:nCams)/pvr(1);
lp(1:nCams) = input pars(6)*pvr(1:nCams)/pvr(1);
% −−> Perform Image Processing with the current processing options
% −−> Perform the reconstruction and cross correlation
% Evaluate score (proportion of false vectors)
score = nnz(velocity(end).nanMask(:))/numel(velocity(end).nanMask);
end % END EVALUATION FUNCTION
end % END MAIN FUNCTION
FIGURE 4.18 Pseudocode function prototype for image optimisation in the
TomoPIV Toolbox for MATLAB. A full function would include error
handling statements, storage of intermediate results and progress
indicators.
4.3 Image Pre-processing for Tomographic PIV 95
0 20 40 60 80 100 120
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
Iteration #
Optimisation of image processing parameters
by derivative−free unconstrained non−linear Nelder−Mead Simplex Method
 
 
Band Pass Noise Lengthscale (Pixels)  (Cameras 1 and 2)
0.01*Gain  (Cameras 1 and 2)
0.1*Gaussian blur radius (σ in Pixels) (Cameras 1 and 2)
Band Pass Noise Lengthscale (Pixels) (Cameras 3 and 4)
0.01*Gain (Cameras 3 and 4)
0.1* Gaussian Blur radius (σ in Pixels) (Cameras 3 and 4)
FIGURE 4.19 Convergence of a Nelder-Mead simplex search for optimal
performance. Note in this case the gain of cameras 1 and 2 were
varied independently of cameras 3 and 4 (at different angles relative
to the laser beam, the first two cameras are illuminated more
brightly than the last two) - gain values may alternatively be linked
using a histogram matching tool.
96 Enhancement of Tomographic PIV
FIGURE 4.20 Vector field slices showing condition of a raw vector field (i.e. no
false vector selection applied) with (top) images processed according
to the input parameter set and with (bottom) images processed using
the parameter set output from the optimisation routine.
4.3 Image Pre-processing for Tomographic PIV 97
function of the search algorithm remains a concern: the optimisation toward
a minimum number of false vectors tends to increase values of gain and blur
beyond those deemed acceptable - it was noticed that many particles satu-
rated before the final application of blur, losing their definition and therefore
reducing the accuracy of particle locations. In fact, the objective chosen re-
sults in less sparse images (larger particles) than is necessary.
Equation 5.4 discussed in Chapter 5 shows that increasing particle size geo-
metrically increases the number of ghost particles and hence bias error in the
cross correlations of reconstructed domains. As a result it seems unwise to
utilise the image optimisation procedure presented here ’in-anger’ unless an
alternative objective function is formulated to reduce both bias and random
errors in the reconstruction.
A possible alternative objective function is the signal to noise ratio in the
correlation plane, which may take into account the effect of increasing ghost
particle density as well as the effect of particle shape and size. A further
alternative is the use of sparsity maximisation in the reconstruction.

5 Effect of the Correlation Tracking
Enhancement (CTE)
5.1 EXPERIMENTAL RESULTS
The work of Elsinga et al. [2011] discussed the bias error introduced into To-
mographic PIV by the coherent motion of ’ghost’ particles in the reconstruc-
tion. Their work is reinforced by Atkinson et al. [2011] who investigate the
error in a specific flow - a boundary layer. These studies use simulated par-
ticle fields to demonstrate the error, with supporting experimental evidence.
Here, experimental results are used verify the conclusions of these studies
and investigate the effect of the CTE on bias error.
Tomographic PIV experiments were performed to measure the buffer and
logarithmic regions of a boundary layer as described in chapter 3.
In Elsinga et al. it is noted that despite the presence of bias error, topological
behaviour of the flow may still be investigated. To that end, the time-series
results produced from the series of experiments are utilised in chapter 6 in
a discussion on denoising and visualisation. Here, we are concerned with
the non-time-resolved results from which boundary layer statistics can be
determined and bias errors measured.
5.1.1 Measurements and Convergence of Statistics
The CTE technique was used to process results from experiments described
in chapter 3 due to the improved robustness and spatial resolution achiev-
able with respect to the VODIM technique. This allowed interrogation of the
domain using windows of size 32 voxels3. Using 75% window overlap led
to 16 vectors in the wall normal direction, spaced at intervals of 9.25 wall (+)
units (0.52 mm) in the range 23.13 ≤ z+ ≤ 161.91.
Interrogation of the field was completed in five passes with window size
(and overlap) as follows: 643 (0%), 41 3 (0%), 323 (0%), 323 (75%) and 323
(75%). The predictor-corrector for volume deformation was stabilised using
a Gaussian smoothing filter [Schrijer & Scarano, 2008] .
99
100 Effect of the Correlation Tracking Enhancement (CTE)
Straightforward 2D PIV images were captured in sync with (i.e. immediately
preceding) the tomographic PIV results and processed using 16 pixels2 inter-
rogation windows and 75% overlap - spatial resolution was comparable to
(slightly higher than) the Tomographic PIV results and a larger field of view
was captured.
A wide variety of quantities can be determined and used for analysis - espe-
cially in the case of the 3D3C Tomographic PIV results. The mean flow, tur-
bulent products, flow derivatives (all terms of the deformation tensor) and
invariant quantities are considered here - no attempt is made to converge on
second derivatives of the flow field.
Images were captured for a total of 3400 velocity fields, sampled at time in-
tervals ∆t = 1s. The noise within Tomographic PIV being higher than 2D PIV,
it was anticipated that convergence for TPIV would be slower. However, the
rate was found to match (or slightly exceed) the 2D PIV data - averaging field
quantities in the X (streamwise) and Y (spanwise) directions allows many
samples to be taken from a single field - so accelerates the rate of conver-
gence. Due to computational constraints, 1500 of the 3400 TPIV fields were
used (full 3400 field dataset was used for 2D PIV statistics).
Convergence plots for computation of the mean flow profile are displayed
in figure 5.1, using the normalised residual as an indicator for convergence.
Equation 5.1 shows the normalised residual for a mean flow field, where NR
is the normalised residual and u1:N represents a velocity field averaged over
N independent samples. On this metric, mean flow velocity and turbulent
intensities are converged to within 1% in 1500 samples.
NR =
u1:N − u1:N−1
u1:N−1
(5.1)
Each sample ui is separated in time by ∆t = 1s. The large scale eddy turnover
time for the boundary layer is δ99/U∞ = 0.25. Since ∆t ≥ 0.25s the samples
can be assumed to be independent.
Convergence plots (also using the normalised residual to indicate conver-
gence) for all nine entries of the deformation tensor are shown in figure 5.2.
An unexplained behaviour is that terms in the cross-stream direction (e.g. uy,
∂uy/∂x and ∂uy/∂z) show slower convergence than terms in the wall-normal
and streamwise direction, despite a similar level of noise in the data. The af-
5.1 Experimental Results 101
0 500 1000 1500 2000 2500 3000 3500
10−6
10−5
10−4
10−3
10−2
10−1
100
101
Number of iterations
N
or
m
al
is
ed
 R
es
id
ua
l
 
 
u x (2d PIV)
u z (2d PIV)
u x (TPIV)
u y (TPIV)
u z (TPIV)
1% residual
0 500 1000 1500 2000 2500 3000 3500
10−5
10−4
10−3
10−2
10−1
100
101
Number of iterations
N
or
m
al
is
ed
 R
es
id
ua
l
 
 
u ′2
x
(2d PIV)
u ′2
z
(2d PIV)
u ′
x
u ′
z
(2d PIV)
u ′2
x
(TPIV)
u ′2
y
(TPIV)
u ′2
z
(TPIV)
u ′
x
u ′
y
(TPIV)
u ′
x
u ′
z
(TPIV)
u ′
y
u ′
z
(TPIV)
1% residual
FIGURE 5.1 Convergence (normalised residuals) of mean flow velocity
components and turbulent intensities measured by TPIV and 2D PIV.
fected terms in the deformation tensor are converged to within 10% of their
long-term mean value; unaffected terms within 1%. Topological invariants
are converged to within 5% with the exception of R, which retains a 20%
confidence interval despite the high number of samples.
5.1.2 Boundary Layer Characteristics
In order to provide a basis for evaluation of Tomographic PIV, the 2D PIV re-
sults were used to characterise the boundary layer with an analytical profile.
The profile selected (van Driest’s) is also used for numerical studies and is
discussed in 5.2.1. Figure 5.3 shows a Clauser plot of the boundary layer con-
102 Effect of the Correlation Tracking Enhancement (CTE)
0 500 1000 1500 2000 2500 3000 3500
10−5
10−4
10−3
10−2
10−1
100
101
102
103
Number of iterations
N
or
m
al
is
ed
 R
es
id
ua
l
 
 
∂ u x/∂x (2d PIV)
∂ u x/∂ z (2d PIV)
∂ u z/∂ x (2d PIV)
∂ u z/∂ z (2d PIV)
∂ u x/∂x (TPIV)
∂ u x/∂ z (TPIV)
∂ u z/∂ x (TPIV)
∂ u z/∂ z (TPIV)
0 500 1000 1500
10−4
10−3
10−2
10−1
100
101
102
Number of iterations
N
or
m
al
is
ed
 R
es
id
ua
l
 
 
∂ u x/∂ y (TPIV)
∂ u y/∂x (TPIV)
∂ u y/∂ y (TPIV)
∂ u y/∂ z (TPIV)
∂ u z/∂ y (TPIV)
0 500 1000 1500
10−5
10−4
10−3
10−2
10−1
100
101
102
Number of iterations
N
or
m
al
is
ed
 re
sid
ua
l
 
 
P
Q
R
λ2
λ ci
5% residual
FIGURE 5.2 Convergence (normalised residuals) of velocity field derivatives and
topological invariants measured by TPIV and 2D PIV.
5.1 Experimental Results 103
taining the 2D PIV data - the logarithmic region is clearly visible. Using this
figure, the value of skin friction velocity uτ was determined and the value
of the log-law constant (C) verified to match that of Coles [1956]. The mean
2D flow profile was used to ascertain boundary layer thickness δ99 and free-
stream speed U∞ (these latter measurements made possible since the field of
view for 2D PIV extended out of the boundary layer).
102 103 104 105
0.4
0.5
0.6
0.7
0.8
0.9
1
1.1
yU
∞
/ν
u
(y)
/U
∞
 
 
Universal law of the wall
2D PIV
FIGURE 5.3 Clauser plot for determination of skin friction velocity from 2D PIV
data.
Having ascertained a fit for the logarithmic region of the flow, the constant
A+ associated with the van Driest profile (in the buffer layer) was adjusted to
maintain the required continuous distribution, as shown in figure 5.4. Table
5.1 shows the resultant values.
The geometry of the setup is such that the laser sheet for 2D PIV is incident
on the wall. Being a boundary layer flow, particles (small enough to appear
as dust) accumulate at the wall since the turbulence is not intense enough (in
the laminar subregion) to sweep the wall clean. A considerable part of the
laser beam is reflected and 2D PIV images become highly saturated close to
the wall, despite efforts to compensate by adjusting the lens aperture. This
results in erroneous measurements close to the wall, as seen in figure 5.3
where the 2D measurements differ largely from the expected profile - which
104 Effect of the Correlation Tracking Enhancement (CTE)
Free-stream velocity U∞ = 0.45 m/s
Boundary layer thickness δ99 = 0.115 m
Wall friction velocity uτ = 0.019 m/s
Wall scale + = 6.1× 10−5 m
Wall scale Reynolds Number δ+ = 1890
Momentum Thickness θ = 9.6× 10−3 m
Momentum Thickness Reynolds Number Reθ = 3657
TABLE 5.1 Boundary layer parameters measured in the CUED tunnel using 2D
PIV synchronised with TPIV acquisition, in good agreement with
earlier data from the same tunnel (table 3.1).
should diminish to zero at the wall.
5.1.3 Bias errors (mean profiles)
Where a flow is locally sheared, the effect of windowing for cross correlation
acts as a local averaging filter, smoothing out strong gradients in the mea-
sured flow field. This is considered in great detail in Elsinga et al. [2011].
Since the effect causes measurements to deviate from an expected analytical
profile, we wish to represent the magnitude of this effect in order to differen-
tiate between bias arising through the effect of ghost particles and bias inher-
ent in a windowed measurement. The analytical velocity profile described
above is filtered using a top-hat filter the same size as an interrogation win-
dow. That is, the analytical profile is locally averaged to mimic the effect of
windowing (32 voxel3 windows) - the Tomographic PIV profiles should be
expected to fall onto the locally averaged, rather than analytic, profile. The
analytical and filtered profiles are shown in figure 5.4 for comparison with
experimental data.
Figure 5.4 shows the mean streamwise velocity profile determined by Tomo-
graphic PIV, plotted over the analytical and 2D PIV profiles. CTE analysis
using windows 323 voxels in size is shown compared to a VODIM analysis
with windows 483 voxels in size. Different window sizes are deliberately
used make the measurements comparable: the number of true particles in
the cross correlation is similar, whereas a CTE correlation with 483 voxel win-
dows would capture a larger number of particles. This effect would unfairly
advantage the CTE analysis in a comparison with VODIM.
5.1 Experimental Results 105
101 102
4
6
8
10
12
14
16
18
20
22
z+
u
+
 
 
Analytic
Filtered Analytic
2D PIV
VODIM: 483 windows, 75% overlap
CTE: 323 windows, 75% overlap
1
2
4
3
5
FIGURE 5.4 Performance of CTE and VODIM techniques in measuring the mean
profile of a boundary layer, including an analytical velocity profile
filtered to mimic the effect of windowed measurement. Numbered
features are explained in section 5.1.3.
There are several salient features of the data, numbered in figure 5.4.
1 Averaging the filtered analytical profile in the range of z+ covered by
the TPIV-CTE domain gives a value of u+ = 15.38. This is the value at
which the bias error must tend toward zero, since bias typically skews
measurements toward the mean velocity in the volume. The measured
(CTE) velocity profile is therefore expected to cross the filtered profile.
For both CTE and VODIM techniques, the crossing point is close to the
centre of the measurement volume. A more detailed discussion is made
below.
2 Bias of the velocity profile toward the domain mean value is clear, even
utilising the CTE technique to minimise the effect of ghost particles.
However, bias error is substantially reduced with the CTE technique at
the same time as improving measurement resolution.
106 Effect of the Correlation Tracking Enhancement (CTE)
3 Elsinga et al. [2011] note that bias errors decrease where local displace-
ment is more than 1 particle diameter different to the mean displace-
ment (expressing this in more formal terms, the peak in the correla-
tion plane caused by coherent motion of the ghost particles becomes
segregated from the true peak - and therefore does not bias the peak
location). At the interrogation position close to the wall, particle dis-
placements are approximately 3 voxels lower than at the half-height of
the volume. The return of the TPIV profile toward the analytical curve
at this location supports the arguments of Elsinga et al., suggesting a
reduction in the degree of bias at this location.
4 The laser beam profile reduces in intensity at the higher values of z+
(over 120); resulting in a reduced particle intensity (and ultimately no
reconstructed particles above this height). Ghost particles therefore be-
come more significant and the profile strongly tends toward the mean
flow value. It is routine to remove these extreme data points for any
investigation of the actual flow - they are shown here to highlight the
effect of a varying laser illumination intensity: that is, a nonuniform
laser beam can alter the magnitude of an existing bias error.
5 The 2D PIV measurement close to the wall (below z+ = 50) is corrupted
by saturation of the recorded images, caused by reflection of the laser
sheet incident on the wall. Hence the 2D PIV profile deviates from the
analytical solution.
As expected for this flow, the mean cross-stream and wall-normal velocities
are equal to zero within the 1% error determined by their convergence limit.
Figure 5.5 shows the crossover between analytic (filtered as above to exclude
the effect of windowing error) and measured velocity profiles, for CTE and
VODIM techniques. The velocity profile from both VODIM and CTE analysis
crosses within 0.3% of the expected location, which is within the confidence
interval (1%) associated with the statistical convergence of these results. Note
that the expected crossover point is determined by integrating the filtered
analytical profile over the range covered by the laser-illuminated volume, from
Z = 0mm to Z = 12mm.
A side-note should be made concerning the effect of ’peak-locking’ in PIV.
Piirto et al. [2005] demonstrate (using image interpolation and peak location
5.1 Experimental Results 107
102
9
10
11
12
13
14
15
16
17
18
Z+
u
+
 
 
Filtered Analytic (323 voxel windows)
Filtered Analytic (483 voxel windows)
VODIM: 483 windows, 75% overlap
CTE: 323 windows, 75% overlap
Mean of filtered analytic within laser volume
Cross − over (zero bias) point 
of CTE and VODIM profiles.
Expected cross−over is at u+=15.38. 
Actual cross−over of CTE and VODIM
profiles is within 0.3% of the expected 
value.
FIGURE 5.5 Crossover of measured and analytic velocity profiles occurs at the
mean velocity within the measurement region, since bias error tends
to zero at this point.
routines used in the present analyses) that there is a so-called ’bias error’
inherent in the peak location process. This is not to be confused with the
bias error discussed here; to prevent ambiguity it will be referred to as peak
locking error. Peak locking error is expected to be below 0.05 voxels (for the
image interpolation and peak location methods used within the Tomographic
PIV Toolbox - see 3.4.2). Being considerably smaller than the random and bias
errors which are addressed within this work (see section 2.1), peak locking is
not considered here.
Figures 5.6 and 5.7 show the bias error expressed in voxels, allowing compar-
ison with other Tomographic PIV experiments and flows. In each case, bias
error in the mean flow is ascertained by subtraction of the filtered analytical
velocity profile F(z) from the measured average flow. The subtracted analyt-
ical profile is filtered according to the window size of each case (32 and 48
voxels for CTE and VODIM tests respectively), to remove dependence of the
bias error on window size: i.e. the difference between curves for CTE and
108 Effect of the Correlation Tracking Enhancement (CTE)
20 30 40 50 60 70 80 90 100 110
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
Z (voxels)
δ
x
(z
)−
F
(z
)
(v
o
x
e
ls
)
 
 
CTE
VODIM
FIGURE 5.6 Bias error (expressed as displacement in voxels) as a function of
wall-normal distance.
−3 −2.5 −2 −1.5 −1 −0.5 0 0.5 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
δx(z )− 1Z
∫ Z
0
F (z )dz (voxe ls)
δ
x
(z
)−
F
(z
)
(v
o
x
e
ls
)
 
 
CTE
VODIM
FIGURE 5.7 Bias error (expressed as displacement in voxels) as a function of
particle displacement away from the volumetric mean displacement.
5.1 Experimental Results 109
VODIM is solely as a result of the characteristics of the technique and not a
function of window size or particle density.
Use of the CTE technique substantially reduces bias error compared to the
VODIM approach as shown in table 5.2.
Displacement CTE Bias (vox) VODIM Bias (vox) Improvement(%)
0.66 0.26 0.43 40
0.48 0.18 0.31 41
0.29 0.11 0.18 39
0.06 0.03 0.05 37
-0.20 0.06 0.06 9
-0.58 0.19 0.25 24
-1.18 0.38 0.63 40
TABLE 5.2 Percentage improvement of CTE bias compared to VODIM bias for
varying displacements
The domains measured have a ’crossover’ velocity (as shown in figure 5.5)
at which bias error is equal to zero. In figure 5.7, this crossover point is used
to determine the displacement of particles (at a given wall-normal location)
relative to the mean flow in the measurement volume. Since (on average)
ghost particles convect with the mean flow velocity (as depicted in figure 2.2),
this demonstrates the effect (item 3, figure 5.4) of ghost particles in distorting
the position of a ’true’ correlation peak.
Elsinga et al. [2011] suggest that ’This phenomenon [bias error] is alleviated when
the difference between particles displacement along the volume depth is increased be-
yond a particle image diameter’. Here, a slightly different hypothesis is made:
Assume a perfect correlation of ghost particles (i.e. the pattern and relative
displacement of ghost particles is constant between two windows) and as-
sume the same perfect correlation of true particles. At least two peaks will
be present in the correlation volume, one relating to the ghost particle dis-
placement and the other to the true particle displacement. Having assumed
perfect matching, the peaks are single spikes. In a discretised volume, these
spikes (being usually at a noninteger location) are described by a region of
nonzero elements having size 2× 2× 2 within the correlation volume. Po-
sitions in the correlation volume further than 1 element away from the peak
location have zero intensity. Thus if the ’ghost peak’ is more than 2 voxels
110 Effect of the Correlation Tracking Enhancement (CTE)
away from the true peak, the shape of the true peak remains unaffected and
vice versa; otherwise the peaks merge, causing bias error when determining
the peak location.
Since the extents of both ghost and true peaks are increased by imperfection
in particle matching, it is expected that in a true case, a distance of somewhat
more than 2 voxels is required; but that the relative displacement need not
be as large as the particle diameter. The evidence presented in figure 5.7 sup-
ports this hypothesis, showing a substantial drop in bias error1 for relative
displacement greater than 2 voxels (in this case 2.5 voxels), despite a particle
image diameter of approximately 4 voxels.
An important conclusion from these observations is that bias error is affected
by (although not caused by) discretisation in the correlation plane. The range
of displacements over which bias error is present can therefore be reduced
by upsampling during the correlation process. This operation can be per-
formed trivially: Some PIV software utilises the (computationally efficient)
FFT-based correlation for the majority of the correlation plane/volume but
performs a (computationally expensive but more accurate) direct cross corre-
lation to ascertain the central array surrounding the correlation peak (which
in the final passes of a VODIM/CTE interrogation process lies at the centre of
the volume). The latter direct correlation stage could be used to perform an
efficient local upsampling of the central position without increasing the com-
putational requirement by upsampling the entire correlation volume. This is
equivalent to a multigrid technique with refinement in the central region.
It is noted that the flow field used here is close to being a worst-case for predic-
tion of bias error. This is due to the high shear present in the velocity field.
Moreover, the region of highest shear is at the bottom of the volume close
to the wall, so all smoothing techniques (including the predictor-corrector
utilised within multi-pass PIV algorithm itself) are expected to exacerbate
the problem.
5.1 Experimental Results 111
0 50 100 150 200
0
2
4
6
8
10
y+
 
 
u ′2x/u
2
τ (2d PIV)
u ′2z /u
2
τ (2d PIV)
u ′xu ′z/u
2
τ (2d PIV)
u ′2x/u
2
τ (TPIV)
u ′2z /u
2
τ (TPIV)
u ′xu ′z/u
2
τ (TPIV)
0 50 100 150 200
−0.1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
y+
 
 
u ′2y /u
2
τ (TPIV)
u ′xu ′y/u
2
τ (TPIV)
u ′yu ′z/u
2
τ (TPIV)
0 50 100 150 200
2
4
6
8
10
12
14
16
y+
 
 
k+ (TPIV)
FIGURE 5.8 Turbulent intensities. (Top) Comparison of TPIV using CTE with 2D
PIV (Middle) CTE-TPIV Values unobtainable with 2D PIV (Bottom)
Summated turbulent intensity from CTE-TPIV results.
112 Effect of the Correlation Tracking Enhancement (CTE)
5.1.4 Bias errors (turbulent intensities)
Although the mean profile is biased, it is worth investigating whether turbu-
lent intensities recorded by TPIV are also affected. The 2D PIV results can
only be used to ascertain intensities in u′ and w′, while TPIV can measure
products of fluctuations in all directions.
Figure 5.8 shows the variation of turbulent intensity terms with wall-normal
distance. The consistency between TPIV and 2D PIV results is poor - 2D PIV
produces curves with a similar shape and magnitude to that expected (as
discussed in section 3.2.2, figure 3.5). Tomographic PIV measurements do not
agree well - exhibiting more than 50% error between the TPIV measurement
and the PIV measurement at most heights above the wall for all quantities
which can be compared (i.e. terms in x and z). Whilst shapes of curves are
similar, the magnitudes vary wildly from those measured with 2D PIV.
Consider first the u′2z /u2τ term. The CTE-TPIV measurement of this term sig-
nificantly under-predicts the 2D PIV measurement (by up to a factor of 4).
In these experiments, bias error is most obvious for its effect in the stream-
wise direction but it also affects velocity measurement in the wall-normal
and cross stream directions; it has the effect of reducing the turbulent inten-
sities based on uy and uz terms. It is therefore this bias effect which causes
under-prediction of u′2z /u2τ.
The u′2x /u2τ term shows different behaviour to u′2z /u2τ. In the streamwise di-
rection, intensities are over-predicted compared to the 2D case. The main
contributory factor is again the bias error. Unlike the wall normal direction
(with uz = 0), there is a mean flow; ux 6= 0. The bias error here strongly
affects the mean profile (see figure 5.4). Error is propagated from the mean
flow to the turbulent intensity via the Reynolds decomposition u′x = ux − ux.
The deviation between true and biased mean profiles is recorded as turbulent
intensity; leading to over-prediction of u′2x /u2τ. This is especially significant
at the top and bottom of the volume where bias error is worst.
Using solely the mean flow bias error explanation for over-prediction of u′2x /u2τ,
agreement between 2D PIV and TPIV could be expected around z+ = 80
at the centre of the volume, where bias error in the mean flow decreases to
1For CTE results only - the spatial resolution of VODIM was too low to show this effect,
all measurement points having displacement less than 1.5 voxels from the mean
5.2 Numerical Accuracy Study 113
zero (see figure 5.6). However, bias error also affects the turbulent intensi-
ties directly. In a mean sense, the mean bias error is zero at the centre point.
Measured axial velocity components at this mid-height location are biased
by particles at the top and bottom of the volume. In a mean sense, these two
sources of bias cancel (hence the ’cross-over’ point of zero mean bias). In
an instantaneous sense, however, the momentary arrangement of particles
causes bias to occur one way or another - away from the true mean value.
The turbulent intensity is therefore increased by the bias error in this case.
In light of the point that bias error can, in some cases, increase the turbulent
intensity, it is worth revisiting the discussion for u′2z /u2τ. It is plausible that
large scale turbulence in the flow can cause some instantaneous measure-
ments to be biased away from the mean, contributing to an increased turbu-
lent intensity. However, since the mean component uz is uniform throughout
the flow the effect of biasing toward the mean is much more significant, lead-
ing to the under-prediction of u′2z /u2τ which is seen.
There are additional sources of error present. The noise inherent in PIV (es-
pecially TPIV) measurements will act to increase turbulent intensities. For
all components, the effect of using larger interrogation windows than for 2D
PIV tends to reduce turbulent intensity of the TPIV measurement compared
to 2D PIV since the windowing effect is similar to a local averaging filter.
A side note is that the magnitude of k+ (the total turbulent intensity) is close
to the values which might be expected from figure 3.5. However, since inten-
sities in the y and z directions were under-predicted while turbulent intensity
in the x direction was over-predicted, it is believed that this result is due to a
coincidental cancellation of error.
5.2 NUMERICAL ACCURACY STUDY
A numerical study was carried out in order to investigate the performance of
Tomographic PIV with the new CTE procedure in a shear flow, since flows
with high velocity gradients are most affected by bias error. The method used
is that of Atkinson et al. [2011]: artificial images are generated using the mean
profile of a turbulent boundary layer then the Tomographic PIV process is
applied to the images. The output velocity distribution is compared to the
analytical solution in order to compare accuracy.
114 Effect of the Correlation Tracking Enhancement (CTE)
Since Atkinson et al. consider the accuracy of the VODIM algorithm using
this technique, the aim of this analysis is not to replicate that work but to
use the same method to investigate the accuracy of the CTE technique when
compared to a VODIM algorithm.
5.2.1 Analytical Mean Velocity Profile
The distribution used is consistent with that used by Atkinson et al., compris-
ing a linear region in the viscous sublayer, a van Driest profile [van Driest,
1956] in the overlap layer and a logarithmic region conforming to the univer-
sal law of the wall.
u+ =

z+ z+ < 5∫ z+
0
2
1+
√
1+ 4κ2z+2(1− e−z+/A+)2
dz+ 5 ≤ z+ < 50
1
κ
ln(z+) + C 50 ≤ z+ < 200
(5.2)
Constants κ = 0.42, A+ = 25.77 and C = 5.20 are taken from a fit to 2D PIV
data for the experiments reported in section 5.1.2.
5.2.2 Simulated Displacement Field
In order to convert the analytical boundary layer profile above to a displace-
ment field for simulation of the TPIV process, a scaling (of wall units to
voxels) equal to the experimental value was applied, thereby replicating the
same mean velocity profile in the reconstruction domain as used in the ex-
periments.
A timestep δt = 2ms is used in order to convect particles between consecutive
fields. This study therefore uses the same δt and makes a close approximation
to the mean flow of section 5.1; values of accuracy reported here may be
applied to the experimental data.2
2Estimations of accuracy are subject to the limitation that there is no fine-scale turbulent
component in the analytical flow field. Since tomographic PIV results are known to be biased
in the presence of strong shear, it is anticipated that the dominant source of error is encapsu-
lated by the present analysis.
5.2 Numerical Accuracy Study 115
5.2.3 Implementation
A minor shortcoming in the work of Atkinson et al. is that numerical studies
are made with no error in the calibration of the cameras. The forward pro-
jection described by Atkinson & Soria [2009] is used to create images directly
from the weightings array and hence a ’perfect’ image conforming to that
weightings array is determined.
Here, an alternative approach is taken, consistent with the work of Worth
et al. [2010], which utilises a real-world calibration in order to encapsulate the
reconstruction error associated with slightly misaligned cameras. A volume
calibration was produced from experimental data utilising a self-calibration
process, with an error of less than 0.5 pixels throughout the entire volume
(median absolute error less than 0.095 pixels in all cameras) as described in
3.3.5. The inverse calibration is determined from the mapping polynomi-
als and particle locations projected from real to image space. Particles are
assumed to have gaussian intensity distribution and images are produced
according to equation 5.3.
Ir,c =
N
∑
p=1
exp
((
r− rp
)2
+
(
c− cp
)2
σ2
)
(5.3)
Where N is the number of particles, [r, c] are the row and column indices
(in image space) of a pixel and [rp, cp] give the location of each particle in
image space. The standard deviation σ is set at a nominal value of 0.75 voxels
(giving a particle diameter of approximately 3 voxels) and scaled by the pixel
to voxel ratio of each camera to give rise to the standard deviation of particles
in image space. Pixel to voxel ratios for all cameras are within 10% of unity
so this factor has a noticeable (but not dominant) effect on the regularity of
reconstructed particles.
The domain is set at 10mm in depth. In voxels, the domain size lx,y,z is [500×
500× 149]. The number of particles per pixel (ppp) is varied. The number of
cameras Nc = 4. Using equation 5.4 (c.f. Elsinga et al. [2011, eqn. 3]), the ratio
between Np (number of true particles) and Ng (number of ghost particles) 3
3The number of ghost particles (e.g. with 0.05ppp) can be higher than the number of
true particles. However, as noted by Elsinga et al. this number does not account for the size
or intensity of the ghost particles, which are characteristically smaller and have significantly
lower intensity than true particles - allowing meaningful cross correlation.
116 Effect of the Correlation Tracking Enhancement (CTE)
is tabulated for this domain and a range of ppp values in table 5.3. Accuracy
studies presented here use this range of ppp values. The particle area Ap is
taken to be 3.0 pixels - as noted by Elsinga et al., the effective area is less
than the area based on diameter of the diffraction spot (which in this case is
7.1pixels2).
Np
Ng
=
1
pppNc−1.ANcp .lz
(5.4)
ppp 0.001 0.01 0.02 0.04 0.05 0.06 0.08 0.10
Np/Ng 82800 82.8 10.3 1.29 0.663 0.384 0.162 0.0828
TABLE 5.3 Variation of true / ghost particle ratio with number of particles per
pixel
Figures 5.9 to 5.11 characterise the differences between CTE and VODIM
methods for three scenarios:
a. Component of flow ’in-plane’ with the volume, where the field includes
mean shear leading to bias error (X components).
b. Component of flow ’in-plane’ with the volume but with no mean shear,
i.e. no bias error (Y components).
c. Component of flow normal to the ’plane’ of the volume, having no
mean shear (Z components).
5.2.4 Numerical Results
The present results are reported for instantaneous fields (averaged along the
streamwise and cross-stream directions) rather than averaging results from
an entire series of data as in section 5.1, to reduce computational load. This
results in some noise in figures 5.9 to 5.11 (seen as occasional jagged features
in the colormaps) but trends and typical values of the data are clear.
The three cases listed above are valuable in comparing the VODIM and CTE
techniques. Let us first discuss the suitability of this flow field for evaluating
error. For evaluating bias error, the use of a shear or boundary layer is the
most extreme case possible.
5.2 Numerical Accuracy Study 117
σ = 0 . 750
Z
(
v
o
x
)
 
 
0.02 0.04 0.06 0.08 0.1
60
80
100
120
−0.4
−0.2
0
0.2
0.4
0.6
σ = 0 . 875
 
 
0.02 0.04 0.06 0.08 0.1
60
80
100
120
σ = 1 . 000
Z
(
v
o
x
)
N ppp
 
 
0.02 0.04 0.06 0.08 0.1
60
80
100
120
σ = 1 . 250
N ppp
 
 
0.02 0.04 0.06 0.08 0.1
60
80
100
120
C
T
E
σ = 0 . 750
Z
(
v
o
x
)
 
 
0.02 0.04 0.06 0.08 0.1
60
80
100
120
−0.4
−0.2
0
0.2
0.4
0.6
σ = 0 . 875
 
 
0.02 0.04 0.06 0.08 0.1
60
80
100
120
σ = 1 . 000
Z
(
v
o
x
)
N ppp
 
 
0.02 0.04 0.06 0.08 0.1
60
80
100
120
σ = 1 . 250
N ppp
 
 
0.02 0.04 0.06 0.08 0.1
60
80
100
120
V
O
D
IM
V
O
D
IM
 a
lg
o
ri
th
m
 u
n
st
a
b
le
(n
o
 r
e
su
lt
s)
 i
n
 b
la
n
k
 r
e
g
io
n
FIGURE 5.9 Variation of streamwise (X-direction) error (in voxels, including the
biasing effect) with wall-normal distance and particle density. (Upper
four) CTE interrogation with 483 windows, (Lower four) VODIM
interrogation with 483 windows, each for four characteristic particle
diameters (diameter Dp in the image plane ≈ 4σ).
118 Effect of the Correlation Tracking Enhancement (CTE)
σ = 0 . 750
Z
(
v
o
x
)
 
 
0.02 0.04 0.06 0.08 0.1
60
80
100
120
−0.02
0
0.02
0.04
σ = 0 . 875
 
 
0.02 0.04 0.06 0.08 0.1
60
80
100
120
σ = 1 . 000
Z
(
v
o
x
)
N ppp
 
 
0.02 0.04 0.06 0.08 0.1
60
80
100
120
σ = 1 . 250
N ppp
 
 
0.02 0.04 0.06 0.08 0.1
60
80
100
120
C
T
E
σ = 0 . 750
Z
(
v
o
x
)
 
 
0.02 0.04 0.06 0.08 0.1
60
80
100
120
−0.02
0
0.02
0.04
σ = 0 . 875
 
 
0.02 0.04 0.06 0.08 0.1
60
80
100
120
σ = 1 . 000
Z
(
v
o
x
)
N ppp
 
 
0.02 0.04 0.06 0.08 0.1
60
80
100
120
σ = 1 . 250
N ppp
 
 
0.02 0.04 0.06 0.08 0.1
60
80
100
120
V
O
D
IM
V
O
D
IM
 a
lg
o
ri
th
m
 u
n
st
a
b
le
(n
o
 r
e
su
lt
s)
 i
n
 b
la
n
k
 r
e
g
io
n
FIGURE 5.10 Variation of cross-stream (Y-direction) error (in voxels) with
wall-normal distance and particle density. (Upper 4) CTE
interrogation with 483 windows, (Lower four) VODIM interrogation
with 483 windows, each for four characteristic particle diameters
(diameter Dp in the image plane ≈ 4σ).
5.2 Numerical Accuracy Study 119
σ = 0 . 750
Z
(
v
o
x
)
 
 
0.02 0.04 0.06 0.08 0.1
60
80
100
120
0
0.02
0.04
0.06
0.08
0.1
σ = 0 . 875
 
 
0.02 0.04 0.06 0.08 0.1
60
80
100
120
σ = 1 . 000
Z
(
v
o
x
)
N ppp
 
 
0.02 0.04 0.06 0.08 0.1
60
80
100
120
σ = 1 . 250
N ppp
 
 
0.02 0.04 0.06 0.08 0.1
60
80
100
120
C
T
E
σ = 0 . 750
Z
(
v
o
x
)
 
 
0.02 0.04 0.06 0.08 0.1
60
80
100
120
0
0.02
0.04
0.06
0.08
0.1
σ = 0 . 875
 
 
0.02 0.04 0.06 0.08 0.1
60
80
100
120
σ = 1 . 000
Z
(
v
o
x
)
N ppp
 
 
0.02 0.04 0.06 0.08 0.1
60
80
100
120
σ = 1 . 250
N ppp
 
 
0.02 0.04 0.06 0.08 0.1
60
80
100
120
Capped
colour limit
(max = 0.24)
V
O
D
IM
V
O
D
IM
 a
lg
o
ri
th
m
 u
n
st
a
b
le
(n
o
 r
e
su
lt
s)
 i
n
 b
la
n
k
 r
e
g
io
n
FIGURE 5.11 Variation of wall-normal (Z-direction) error (in voxels) with
wall-normal distance and particle density. (Upper 4) CTE
interrogation with 483 windows, (Lower four) VODIM interrogation
with 483 windows, each for four characteristic particle diameters
(diameter Dp in the image plane ≈ 4σ).
120 Effect of the Correlation Tracking Enhancement (CTE)
Consider scenario a. A group of particles, lying at some height H from the
wall result in a ghost particle at some other height H′. In this artificial sim-
ulation, there is no representation of turbulence in the field so the group of
particles will be convected by the same distance for every step δt in time.
That is to say, identical4 ghost particles appear for both (VODIM) or all four
(CTE) fields. Since mean shear is present, bias error is expected. Since ghost
particles are identical over all time steps, the CTE technique is expected to
result in the same bias as the VODIM technique. This is borne out by figure
5.9, although CTE does show a slight reduction in bias - due to the consis-
tent discretisation error in true particles but inconsistent discretisation error
in ghost particle reconstruction causing the magnitude of the ’ghost peak’ to
diminish compared to the magnitude of the true peak.
Consider scenario b. With no bias due to mean shear, the magnitude of er-
ror is expected to be comprised of image distortion, calibration error, recon-
struction/discretisation error and PIV error. We here classify these sources
of error as ’noise’ (i.e. as having no mean biasing effect). Having a zero Y-
velocity component in the analytical field results in somewhat reduced noise
error compared to an arbitrary flow5 and is susceptible to the increased parti-
cle density effect discussed below. Figure 5.10 shows that the CTE technique
decreases the level of noise by a factor of at least 1.5 (¿2 in 90% of cases) over
all combinations of particle size and density tested compared to the VODIM
results. Moreover, the CTE technique is seen to stabilise the processing algo-
rithm: blacked out areas in the colormaps show regions where the VODIM
technique produced completely erroneous data but CTE continued to work.
Scenario c. is closely similar to scenario b., with the same limitations. How-
ever, as discussed in section 3.2.2, the accuracy of Tomographic PIV is less-
ened in the direction which is nominally normal to the line of sight (the
’depth’ direction in the volume). This is due to reconstruction effects, where
particles are elongated in that direction (despite a camera separation angle
4’Identical’ in the sense that the ghost particle is placed analytically in the same position.
discretisation of the images and volume mean that the ghost particle does not appear identical
between time steps, especially where the displacement of the ’true’ particle group is not close
to an integer number of voxels
5Having zero velocity component in the direction of interest results in no error due to
different discretisation effects in the reconstruction (all particles are at the same Y location
relative to the voxel grid for all time steps). Any peak locking present in the PIV interrogation
is also minimised by this choice of field.
5.2 Numerical Accuracy Study 121
which in this case exceeded 30◦ for all cameras as recommended by Elsinga
et al. [2006]). Comparing figure 5.11 with 5.10 highlights that reduction in
accuracy, while comparing CTE (top four colormaps) with VODIM (bottom
four colormaps) in figure 5.11 shows a large reduction in error - by a factor
between 2 and 7. This is to be expected, since the error source is as a result of
inaccuracy in reconstructing individual particles. Using CTE, the same par-
ticles are captured at multiple locations - effectively averaging out the error
in each instantaneous reconstruction.
Note that the presence of ghost particles (which are ’identically’ placed be-
tween time steps as discussed above) actually increases the perceived accu-
racy of displacement measurement in the Y and Z directions - since ghost
particles have the same velocity component as true particles, their effect is
the same as increasing the particle density. Account must be made for this
if using the results presented here to indicate noise levels within generalised
Tomographic PIV data (considering noise as experimental error excluding
bias), but in this case results suffering the same effects are compared (be-
tween CTE and VODIM) so no such correction is necessary. Account could
be made by using table 5.3 to determine the density of true + ghost particles
in the present analyses, and using that value to compare against the density
of true particles in the case for which an estimate of noise error is required.
To summarise these scenarios:
- The numerical simulation over-predicts bias error compared to the ex-
perimental case, due to artificially strong correlation of ghost particles.
- The numerical simulation is likely to under-predict the level of noise in
a Tomographic PIV measurement of an arbitrary flow unless the parti-
cle density effect is accounted for.
A useful additional conclusion is that turbulence in a flow reduces bias error
in the mean flow statistics, although introduces local bias effects to instanta-
neous fields which may not be present once averaging has taken place. This
relates to the discussion on effects of biasing on turbulent intensities in sec-
tion 5.1.4, and it is important to recognise that the presence of turbulence has
a nonlinear effect on the accuracy of TPIV.
122 Effect of the Correlation Tracking Enhancement (CTE)
The non-turbulent nature of the analytical field used prevents reduction of
bias error when using the CTE technique - as opposed to the experimental
work above where reduction in bias is observed, due to the effect of turbu-
lence in the field allowing the CTE method to work as intended. However,
the CTE technique is shown to substantially reduce measurement noise, es-
pecially for velocity components in the ’depth’ direction of the field.
5.3 REVIEW OF RESULTS
The CTE technique, which was formulated and described in chapter 4 has
been tested both experimentally and numerically.
Results in section 5.2 show a strong reduction in random error compared
to existing Tomgraphic PIV analysis types. In the numerical study (which
disadvantaged the CTE technique compared to VODIM), random error was
reduced by a factor of at least 2 for all particle densities considered.6 In cases
of highest particle density, CTE works where VODIM algorithms fail due to
extremely noisy correlation volumes.
Since the CTE technique is more robust to small-scale movements within
windows (see description of the problem in section 3.2.2), it allows use of
longer δt values than otherwise acceptable in high Reynolds Number flow.
Although all the results presented within this work are interrogated with
the same value of δt for comparability, investigation revealed that δt could
be doubled (compared to that achievable with VODIM) in the flow consid-
ered within this work. This is to the detriment of temporal resolution, but
compounds the benefit of improved accuracy in the measurement by further
doubling the dynamic range: a trade-off worth making in many cases.
In summary, the dynamic range of the measurement can be improved com-
pared to conventional TPIV by a factor of at least four7, subject to the bias
error being within allowable limits for a given flow.
Bias error associated with ghost particles in the reconstruction was reduced
6Improvement of random error is assessed on the basis of improvement in eY and eZ,
since estimation of eX is strongly affected by bias error
7Factor of 4 comprises a factor of > 2 improvement in accuracy at the ’bottom end’ of the
velocity range compounded with a factor of 2 in the maximum velocity robustly measurable
(due to doubling of the allowable δt value)
5.3 Review of Results 123
by approximately 40% throughout the volume (in a region of highly sheared
fluid flow); decreasing to a lesser reduction in error at the ’mid-point’ where
bias error decreases to zero anyway. This reduction is consistent in both ex-
perimental and numerical analyses.
The CTE technique increases computational demand compared to existing
TPIV by less than or equal to a factor of 2; i.e. the approach is substan-
tially less computationally intensive than alternative methods of improving
the SNR in TomoPIV such as the MTE-MART approach [Novara et al., 2010]
which requires up to 10 iterations of the MART-VODIM process and with
which the CTE approach shows some similarities.
Since random error has been substantially improved, a recommendation for
future work is the task of integrating the CTE with the spatial derivative
correlation technique of Scarano [2004]. It is anticipated that combination
of the two approaches could be powerful in improving measured velocity
gradients whilst controlling the random measurement noise.

6 Constrained Restoration of
Tomographic PIV Data
Tomographic PIV allows the measurement of a fully 3D-3C velocity field,
but a significant degree of noise is introduced to the data. In 2D PIV and
hot-wire measurements, noise is often removed by smoothing and filtering
respectively. However, having knowledge of the fully 3D flow field allows
more advanced techniques to be used in post-processing of the data, which
utilise a-priori knowledge of the flow field.
Following post-processing, results must remain meaningful. Say we know
that a ’true’ field is likely to be smooth (to some degree), but the correspond-
ing measured field contains noise. Smoothness of the measured field is triv-
ially achieved by setting all velocity components to zero; this is a valid post-
processing operation but clearly not meaningful. Thus, we wish to find a
smooth field which best represents the true data. Since we don’t have true
data, but we know that the measured data is close to it, we satisfy ourselves
by finding a smooth field which best represents the measured data.
Put more formally, the ideal post-processing operation removes all experi-
mental error without changing the ’true’ data. This is not possible in a gen-
eral case, so we require operations to be applied (to reduce experimental er-
ror) whilst minimising the change in true data. In fact, since the true field is
not known, minimising the change in true data cannot be achieved; but such
a procedure is closely approximated by minimising the change in the measured
field.
A requirement for minimisation implies that an optimisation is carried out.
Here, such optimisation has many degrees of freedom (large number of vec-
tors in the flow field, each with 3 components) and is constrained by the
nature of each operation applied.
The principle of Projection Onto Convex Sets (POCS) is used in order to apply
an optimisation (to minimise the magnitude of correction to the input field)
whilst applying constraints relating to a-priori knowledge of the input field.
The principle can either be used to reduce error in, or to impose required
125
126 Constrained Restoration of Tomographic PIV Data
constraints on, the measured field (or both). That is, the action of the post-
processing operator need not strictly reduce error.
Here, we utilise POCS in order to apply two constraints of particular interest:
Incompressible The fluid flow measurements considered elsewhere in this
work are for an incompressible fluid. However, evaluating
divergence of the results (section 5.1) reveals that measured
fields are not divergence free, due to experimental noise.
We thus apply the incompressibility constraint to enforce
this known property.
Scale-limiting Turbulent boundary layer measurements (section 5.1) are
at a higher Reynolds Number than can be fully spatially
resolved using Tomographic PIV, despite improvements in
spatial resolution brought about by Correlation Tracking
(see chapters 4 and 5). Scales in the flow smaller than the
spatial resolution appear as noise in the results and should
be handled appropriately to allow correct visualisation and
interpretation of scales which are measurable.
Making the assumption that noise is not spatially correlated with itself, the
process of scale-limiting can also be considered as an effective de-noising
technique. The latter constraint (above) thus addresses discussion in the lit-
erature survey (see 2.3) relating to both denoising and band/scale limiting.
The process of POCS is also applied in section 6.5.3 to perform restoration
(see 2.3) of a locally ’corrupt’ input - in this case a field in which some vectors
have been identified as invalid and require ’filling-in’.
6.1 PROJECTION ONTO CONVEX SETS
The technique of Projection Onto Convex Sets (POCS) is a constrained op-
timisation technique allowing an effectively infinite number of dimensions.
In particular, the method allows application of one or more known a priori
constraints, which may be nonlinear in nature. Thus it is suitable for use
on datasets such as 3D-3C velocity fields from Tomographic PIV, where the
6.1 Projection Onto Convex Sets 127
number of vectors in the field is high and where the measurement technique
does not account for a number of valid physical constraints.
Starting Point
(Initial measurement / guess)
Convex Set C1
Convex Set C2
Othogonal Projection P1 (onto C1)
Valid Solutions (belonging to both C1 and C2)
Converged solution
FIGURE 6.1 Projection onto convex sets in Hilbert Space, where the convex regions
C1 and C2 represent constraints applied to the input.
Rather than minimising an objective function, the optimisation is achieved
by projecting an initial guess onto a ’set’, which represents the group of all
possible solutions obeying some pre-defined constraint. This projection is
achieved through the use of a ’projection operator’, a linear or (more usu-
ally) nonlinear map onto the set. Any number of constraints may be applied
through successive application of different operators. Iteration of the pro-
cess converges the guess to a solution which is a member of all sets, thereby
meeting all constraints (Figure 6.1).
In formal terms, a set is defined as a Closed Linear Manifold (CLM) in Hilbert
Space. Sets must be closed, convex (a tangent line at any point on the set
boundary never passes through the set) and at least one set must be finite to
ensure convergence of the algorithm. Figure 6.2 shows how the use of non-
convex (concave) sets can lead to non-convergence of the algorithm. The
requirement for closedness follows from convexity since a hypersurface is
generally not convex1.
The principle of POCS has been most frequently applied in image process-
1A hyperplane is the only example of a convex hypersurface, however hyperplanes are
not strictly convex as defined in Youla & Webb [1982, Figure 2], which is a requirement for
rapid convergence
128 Constrained Restoration of Tomographic PIV Data
Starting Point
(Initial measurement / guess)
Concave Set C1
Concave Set C2
Othogonal Projection P1 (onto C1)
Valid Solutions (belonging to both C1 and C2)
Iteration never converges to a valid solution
FIGURE 6.2 Effect of set concavity on the convergence of a POCS algorithm
ing techniques (application to scalar fields). The papers reviewed in 2.3.4
introduce the application of the technique for image restoration; in particu-
lar, Youla & Webb [1982] give a thorough exposition of the associated theory.
The work of Simard & Mailloux [1988, 1990] extends the work to operators
in a vector space (2D). The reader is referred to those sources for detailed
proofs but a review is given here, adopting the notation of Youla & Webb
except where stated.
While the literature refers to the process of POCS as ’restoration’, the number
of degrees of freedom in a general problem is substantially higher than the
number of constraints applied. Here, the term ’reconstruction’ is used since
it is more appropriate to the solution of an underdetermined problem.
6.1.1 Review of Theory
Consider an original vector g. It is known a priori that g belongs to a linear
subspaceR within Hilbert SpaceH. The vector g represents a ’true’ solution
but is not known to or measurable by an observer.
Here, we define g as comprisingM elements which in turn represent a scalar
or vector field of dimension N , where 1 <= M,N ,<= ∞ and M >= N
always.
An observer measures a vector f , which also lies withinH. SinceH contains
both f and g, it can be stated that f represents an orthogonal projection of
g onto some other linear subspace in H. In general the projection is both
6.1 Projection Onto Convex Sets 129
unknown and nonlinear, thus no inverse projection operator is available and
f consists of partial data from g.
Given f , application of a sufficient number of constraints allows unique de-
termination (restoration) of g. Each constraint corresponds to a closed convex
set in H, so g occupies the intersection of the m sets C1, C2, . . . , Cm (eqn. 6.1,
from Youla & Webb eqn. 1).2
g ∈ C0 =
m⋂
i=1
C1 (6.1)
In the case that the sets intersect at a singular point, g can be uniquely re-
covered; otherwise g occupies the intersection of the sets C0 and although
not unique, a closer approximation to g than f may be obtained (assuming
f /∈ C0 already).
In summary, the problem of POCS reconstruction is as follows: Given f ,
apply a sufficient number of physically appropriate constraints in order to
uniquely recover (or more closely approximate) g.
6.1.2 Projection operators in a 2D vector space
To project f onto a linear subspace Ci ⊂ H, a corresponding operator Pi
must be defined. The derivation of Pi is dependant on the type of constraint
to be applied - a variety of common operators are described in the literature,
especially for scalar fields.
For the purposes of this work, interest lies in the application of POCS to vec-
tor spaces (i.e. velocity fields). There are three immediately relevant oper-
ators, each defined here for a function s, which is analogous to f above but
existing in a 2 dimensional vector space.
2Note that notation f , g is reversed compared to the convention used in Youla & Webb to
retain consistency with later works such as Simard & Mailloux [1988].
130 Constrained Restoration of Tomographic PIV Data
Operator P1 Applies a frequency band limit3 to an input function s. P1 is
the projection onto subspace C1 ⊂ H which is the set of all
functions whose Fourier transform is zero outside a prescribed
region δ in the Fourier plane having coordinates (u, v).
P1s =
∣∣∣∣ S(u, v) (u, v) ∈ δ0 (u, v) /∈ δ (6.2)
Operator P2 Applies a prescribed value to elements of the input function.
Formally, C2 ⊂ H is the set of all s ∈ H whose components
(sx, sy) projected onto arbitrary subsets Ix and Iy ⊂ H assume
the prescribed u(x, y) and v(x, y) respectively.
P2s =
∣∣∣∣∣
sx(x, y) = u(x, y) (x, y) ∈ Ix
sy(x, y) = v(x, y) (x, y) ∈ Iy
s otherwise
(6.3)
Operator P3 Applies a divergence-free condition to an input function, i.e.
C3 is the subset of functions in H for which ∇ · s = 0. Deriva-
tion of the operator for N = 2 and extension to N = 3 is
addressed in section 6.2.2.
P3s = s−∇p (6.4)
Where p denotes potential of a ’corrector’ velocity field deter-
mined by solving the poisson equation (see section 6.2.2).
The above operators are reproduced from Simard & Mailloux [1988, paras.
1, 2, 3 and eqns. 5, 6 and 8 respectively] for vector fields in 2D space with
components u and v (N = 2). Notation s is used here to prevent conflicting
terminology between the use of f above and by Simard & Mailloux.
3Direct application of the P1 operator (which is a ’top-hat’ function) in a numerical envi-
ronment will result in artefacts (ringing) due to edge effects of the Discrete Fourier Transform.
In reality, a filter must be constructed to correctly apply the bandwidth limit.
6.2 A Divergence Reduction Operator in 3 Dimensions 131
6.2 A DIVERGENCE REDUCTION OPERATOR IN 3 DIMENSIONS
6.2.1 A Note on Hilbert and Sobolev Spaces
Consider a three-component vector field, whose component values are finite
and nonsingular, existing within a finite region Ω in three dimensional, real
Euclidean space R3. A function f which describes the values (three compo-
nents of velocity ui for i = 1, 2, 3) of the field at any point x within Ω has the
property of being quadratically integrable; i.e.
∫
Ω
| f (x)|2dx<∞ (6.5)
By definition of a quadratically integrable function, there exists a complete
inner product space. Thus the function space (consisting of all possible func-
tions) is a Hilbert Space H, denoted L2(Ω). In simplistic terms, each point
in H corresponds to a unique function and therefore a unique vector field
within Ω.
We dispense with the confusing notation (u, v) used by Simard & Mailloux
[eqn. 2], which utilises the two-dimensionality of the problem to express the
inner product as a dot product. The generalised inner product and norm of
H are given by equations 6.6 and 6.7 respectively.
〈 f , g〉=
∫
Ω
n
∑
i=1
fi(x)gi(x)dx ∀ f , g ∈ H (6.6)
|| f ||=〈 f , f 〉 (6.7)
where n = 3 in this three component example.
Functions are continuous for this case, but continuity is not a requirement
for eqn. 6.5 to be satisfied. The number of unique and valid functions char-
acterising fluid flow within Ω is infinite, thus the space H used here is also
infinite.
Before continuing with derivation of the P3 operator, it is worth noting that
continuity and finiteness (both of which are assumed here for the case of
fluid flow) of a function f ∈ H together imply that well-behaved derivative
fields of f can be obtained. Here, functions are required to be smooth to
at least the second derivative in order that the equations of incompressible
132 Constrained Restoration of Tomographic PIV Data
fluid flow are satisfied. By observation, the first derivative field D f is a finite
continuous vector field which also lies within L2(Ω). We thus introduce the
notion of a Sobolev SpaceW ; a vector space consisting of the set of functions
in Lp(Ω) whose derivatives (up to a given order) also reside in Lp(Ω). This
is expressed formally in equation 6.8.
W s,p(Ω)={ f ∈ Lp(Ω) : ∀|α| ≤ s, ∂αx( f ) ∈ Lp(Ω)} (6.8)
where in d dimensions α = (α1, . . . , αd), |α| = α1 + · · ·+ αd, and pertains to
the degree of differentiability d of the function in a given dimension. The
function derivatives ∂αx = ∂1x, . . . , ∂dx are taken in a weak sense (i.e. the func-
tion is a distribution and a solution for the derivative in the classical sense
may not exist).
The value of p denotes that a function is p− integrable. In the case that p =
2, W s,2 is a Hilbert Space Hs(Ω). For construction of the P3 operator, we
will see below that the first derivative of f is required (s = 1). A further
observation can be made that taking derivative of a field is a linear operation
in an entirely real domain. Thus, D f can be further restricted to lie within
H−1(Ω), the dual vector space ofH1(Ω) (eqn. 6.9).
f ∈ H−1(Ω) (6.9)
6.2.2 Derivation of the P3 operator
We aim to derive a projection operator P3 equivalent to that in 6.1.2, but ap-
plicable in three dimensions. Section 6.2.1 constrains f and D f to define an
input vector field and its first derivative respectively in an n dimensional re-
gion Ω. Following the general methodology described in Simard & Mailloux
[1988], we here review the derivation whilst making extension to 3D.
The first condition is that the divergence-free vector space onto which we
wish to project is convex. Convexity of such spaces is shown in Stark [1987]
and discussed further in Suter [1994] for an arbitrary number of dimensions.
The work is not repeated here since no extension (e.g. to 3D) is required.
Projection of an input field f is made onto the subspace C3(Ω) (eqn. 6.10).
C3(Ω)={g ∈ L2(Ω), ∇.g} (6.10)
6.2 A Divergence Reduction Operator in 3 Dimensions 133
Velocity fields may be decomposed linearly into components consisting of
hyperbolic, rotational and divergent fields [Jeong & Hussain, 1995]. A cor-
rector field f c is therefore formulated (eqn. 6.11) in order to remove the di-
vergent component of f .
g= f − f c, ∇. f = ∇. f c (6.11)
Arbitrary point in L2
Non orthogonal projection
Equipotentials of 'distance' in Hilbert Space from the 
starting point to the set boundary (distance defined
 using the norm of an n-dimensional displacement vector)
Boundary of set
Orthogonal projection (minimal norm  
required to project onto the set)
FIGURE 6.3 Requirement for orthogonality of projections in Hilbert space
Any arbitrary velocity field whose divergence is equal to the divergence of
f satisfies this requirement on f c presented here. This projects f onto C3.
However, it is a requirement that the projection is orthogonal. Geometrically,
this corresponds to projection to the closest part of the set (figure 6.3). Alge-
braically, this is a question of minimising the inner product (eqn. 6.6) between
the input f and the projection g:
min
g∈C3
|| f − g|| (6.12)
134 Constrained Restoration of Tomographic PIV Data
Simard & Mailloux minimise 6.12 by re-expression as a saddle point problem
and utilisation of a Lagrange multiplication. Replicating their equation 11:
min
g∈L2(Ω)
max
p∈H1
1/2 ‖ f − g ‖2 +
∫
Ω
p∇.g (6.13)
which gives rise to the following optimal conditions...
∇.g=0
fc=∇p (6.14)
The Lagrange multiplier p (and hence fc) is determined by solution of the
Poisson equation (6.15, from Simard & Mailloux, eqn. 9) with a homogeneous
Dirichlet (constant p value) boundary condition (eqn. 6.16).
∇2 p = ∇. f ∈ H−1(Ω), x ∈ Ω (6.15)
p = 0, x|γ (6.16)
with f c = ∇p and where γ denotes the boundary of Ω. The selection of this
boundary condition is intuitive for a 3D measurement domain in a homo-
geneous fluid, since selection of a Neumann-type or nonuniform boundary
condition requires use of the input measurement (subject to noise and inac-
curacy) in order to derive the boundary. Solution of the correction is thus
invalidated by dependency on the input field.
The domain boundary γ surrounds the region of measurement without in-
tersecting it (i.e. passing through any measurement points - see figure 6.4).
This is a requirement of the equation system above since imposing 6.15 at the
boundaries 6.16 results in multiple constraints leading to numerical instabil-
ity. In practical terms, this is implemented through solution of a standard
poisson matrix on a 1-padded domain.
The exposition above has closely followed Simard & Mailloux [1988] whilst
adding an explanation based on geometric means. Through the geometric
arguments made it can be seen that since all the above steps are made in
Hilbert space, no limitation is placed on dimensionality of the problem. Thus
the procedure above applies for a three dimensional volume as well as the
two dimensional area considered in the original application of this technique.
6.3 Wavelet-Based Band Limiting: The P4 Operator 135
TPIV Measurement Domain
1-padded boundary grid, p = 0
Measured scalar field p on TPIV 
grid (window centres)
FIGURE 6.4 Application of homogeneous Dirichlet-type boundary conditions on a
1-padded grid for solution of the Poisson equation (simplified to 2D)
6.2.3 Uniqueness of solution
The application of three dimensionality in similar problems can cause diffi-
culty. For example, determination of the third component of velocity from a
3-Dimensional, 2-Component volume is indeterminate. Here also, the correc-
tor field fc as presented in eqn 6.11 is indeterminate (i.e. not unique). How-
ever, in contrast to the 3D2C example, the technique presented here is con-
strained by application of the minimisation requirement (eqn 6.12). Results
are therefore meaningful, but the caveat must be stated that this technique
determines the closest field exhibiting the desired property (incompressibil-
ity), rather than the exact field corresponding to the noisy input.
6.3 WAVELET-BASED BAND LIMITING: THE P4 OPERATOR
6.3.1 Application of the P1 operator to Tomographic PIV data
The P1 operator defined in section 6.1.2 extends trivially to three dimensions.
The operator was implemented and applied to Tomographic PIV data using
a basic ’top-hat’ filter, chosen to eliminate scale (’frequency’) content above
the Nyquist criterion corresponding to the grid spacing.
Figure 6.5 shows the effect that this projection has on an input velocity field.
While the results are smooth, the boundaries of the domain vary widely from
136 Constrained Restoration of Tomographic PIV Data
the input field due to the periodicity of the operator. Critically, the presence
of erroneous vectors in the input field (of which several are obvious) causes
substantial change in topology in the output field - that is to say that noise/er-
ror (at high frequency) is spread to the surrounding region by this operator
in order to meet the band limit, instead of effectively cut out.
−20 −10 0 10 20
−20
−15
−10
−5
0
5
10
15
20
X (mm)
Y 
(m
m)
−20 −10 0 10 20
−20
−15
−10
−5
0
5
10
15
20
X (mm)
Y 
(m
m)
FIGURE 6.5 Application of the P1 operator to a velocity field. (Left) Slice of a
3D3C input field. (Right) Slice at same location following application
of P1. Note the smearing of underlying flow topology close to
erroneous vectors and boundaries.
It is possible to design a filter which applies a band limit in the spectral do-
main whilst suppressing the ringing effect. One such example is the ’Butter-
worth’ filter [Butterworth, 1930]
H(jω) =
1
1+ (−(jω)
2
ωc
)n
(6.17)
where ωc is the cutoff frequency of the filter, and n is the order of the filter.
A wide variety of alternative filter designs are proposed by the literature.
For the purpose of band-limiting Tomographic PIV data, selection of a filter
is somewhat arbitrary since the filter characteristics need not bear any re-
lation to the physics of the problem. Results can vary widely according to
the selection of filter type and tuning parameters (such as the order of the
Butterworth filter n).
6.3 Wavelet-Based Band Limiting: The P4 Operator 137
One further problem is that of periodicity. Spectral decompositions are in-
herently periodic in nature. The small domain, and the (sometimes) large
difference between values of velocity components across the field can result
in large errors where periodic ’wrap-around’ is made in filter construction.
6.3.2 Derivation of an alternative band-limiting (’P4’) operator
To address the problems described in 6.3.1, a wavelet-based filter was imple-
mented. Wavelet Transforms are analogous to Fourier Transforms, utilising
generalised instead of sinusoidal basis functions. Rather than give a thor-
ough exposition here, the reader is referred to a number of key texts; see
Daubechies [1992] and Farge [1992].
The work of Griebel & Koster [2000] and Farge & Schneider [2001] is of partic-
ular relevance here, since a wavelet transform is used in order to decompose
a three dimensional fluid flow into components - unlike a Fourier Transform,
information regarding both spatial positioning and scale of features is re-
tained in a wavelet transform (made possible due to the locality of the basis
functions selected). The fluid flow can be filtered by either energy content
of features in the flow (such as vortical structures) or by scale-length of the
features (figure 6.6).
FIGURE 6.6 Filtering a 2D scalar field. (Left) Input field. (Right) Band-limited field
using wavelet transform thresholded by scale.
Retention of spatial information and the ease with which non-periodic in-
tervals can be applied (see figure 6.7 compared to figure 6.6) make wavelet-
based filters far preferable to spectral filters for this application. While the
138 Constrained Restoration of Tomographic PIV Data
Ringing artefacts
FIGURE 6.7 (left) Spectral equivalent of the scale-filtered field in figure 6.6. A
top-hat filter with cutoff at the Nyquist frequency was used. (right)
Zoomed region showing the top-right corner, highlighting ringing
effects and periodic artefacts close to the boundary.
range of possible filters remains large, selection of the precise filter type can
be made using a greater degree of physical insight than a spectral filter (see
section 6.3.3).
The work of Farge, Schneider and colleagues is mostly aimed toward im-
proving solution techniques for CFD by increasing the sparsity of the prob-
lem. This is irrelevant here, since the velocity fields are extremely small in
comparison with those required for general CFD problems. However, the
principle that the velocity field can be filtered in the wavelet domain trans-
lates well.
False vectors and noise arising from application of PIV to a velocity field can
be large in magnitude. Translated into the wavelet domain, this results in en-
ergetic components within the field. However, coherent structures may also
be highly energetic (depending on the nature of the flow field) so filtering on
the basis of energy content could be unwise. However, it is known that noise
is typically random thus has a characteristically short length-scale.
As a result, the top-hat P1 operator is replaced here with an operator to per-
form scale-based filtering in the wavelet domain:
6.3 Wavelet-Based Band Limiting: The P4 Operator 139
Operator P4 Limits the range of scales (feature sizes) retained in an output
field. P4 is the projection onto subspace C4 ⊂ H which is the
set of all functions s whose wavelet transform W(s) is zero in-
side a prescribed region δ in the wavelet domain.
W(P1s)⇔
∣∣∣∣ W(s) x /∈ δ0 x ∈ δ (6.18)
For a cuboid region Ω with cartesian coordinate system x =
(xi, xj, xk), limitation of features (by scale) to half the measure-
ment resolution uses a δ mask defined by:
δ =
∣∣∣∣ 1 ceil(Ni,j,k2 ) <= i, j, k <= Ni,j,k0 elsewhere (6.19)
where N denotes number of grid points in the i, j, k directions.
6.3.3 Selection of the wavelet filter
One criticism of the P1 operator was the arbitrariness of filter selection. Sim-
ilar criticism can be made of this P4 operator since it is possible to select
different wavelet types in order to perform the transform: the shape of the
wavelet (which is effectively a filter) will affect the content of the signal con-
tent which is filtered out. However, it is possible to select wavelet types on
the basis of more physical reasoning than the selection of tuning parameters
for a spectral filter.
To compare with precedent, a variety of wavelets have been used in similar
problems, including quintic-spline [Farge et al., 1999, 2002], univariate in-
terpolets [Griebel & Koster, 2000] and mexican-hat wavelets [Schram et al.,
2004] (all for fluid simulation or coherent vortex extraction). The signal pro-
cessing community frequently uses the Daubechies family of wavelets [Mal-
lat, 1999] for denoising purposes.
The current problem is effectively that of signal denoising, rather than flow
simulation. Higher order filters such as the quintic splines tend to retain
finer detail and are therefore less effective for denoising. In contrast, the
underlying signal is a fluid velocity field - very low order filters such as the
Daubechies 4th order wavelet [Daubechies, 1992] do not reflect the physics of
a fluid continuum. A second order, biorthogonal spline is used in order to:
140 Constrained Restoration of Tomographic PIV Data
- Maintain the appropriate level of continuity in the reconstructed signal
- Effectively minimise noise in the filtered signal
- Allow invertible wavelet transforms (biorthogonality)
The wavelets discussed here were investigated by their application to the
streamwise component of an instantaneous flow field (taken from the bound-
ary layer experiments described in section 3.1). Figure 6.8 shows the origi-
nal data together with the removed component and final filtered field for
Daubechies 4th order wavelets, biorthogonal quintic spline wavelets and the
chosen biorthogonal second order spline wavelets. The biorthogonal second
order spline wavelets clearly give the best signal retention whilst effectively
denoising - in formal terms, the data removed appears entirely incoherent.
FIGURE 6.8 Denoising a flow field variable using wavelet filtering. Top to bottom:
Daubechies 4th Order, Biorthogonal Quintic Spline, Biorthogonal
Second Order (2:2) Spline. Left to right: Original Field, Filtered Result,
Removed Components. Intensity represents magnitude of streamwise
velocity component in a slice plane parallel to the wall.
The literature review (see 2.3.5) discusses the requirement for false vectors to
be more effectively handled during post-processing stages. In contrast with
the spectral filter of figure 6.5, figure 6.9 (which is applied to the same input
6.4 The RELAX Algorithm and Variants 141
−20 −10 0 10 20
−20
−15
−10
−5
0
5
10
15
20
X (mm)
Y 
(m
m)
−20 −10 0 10 20
−20
−15
−10
−5
0
5
10
15
20
X (mm)
Y 
(m
m)
FIGURE 6.9 Application of the P4 operator to a velocity field. (Left) Slice of a
3D3C input field. (Right) Slice at same location following application
of P4. Note the improved retention of underlying flow topology and
lack of error due to periodic conditions at the boundaries.
field) shows the effective removal of erroneous data without the underly-
ing flow topology being substantially altered. Treatment of data close to the
boundaries is also improved (the wavelet transform used is adjusted for the
interval to prevent the application of periodic conditions at the boundaries).
To ’survive’ the filter, false vectors must be in clusters which are at least 2×
2× 2 grid points in size - i.e. the cluster must have a length scale greater than
the cut-off scale in all three dimensions, which is rare - hence the technique is
more effective than a spectral filter at removing clusters of false vectors.
6.4 THE RELAX ALGORITHM AND VARIANTS
6.4.1 Implementation Framework
In terms of numerical implementation, this work is concerned with the ma-
nipulation of 3 component, 3 dimensional vector fields. Such data is stored
readily as a set of three 3D arrays and the parsing of data between different
subroutines and functions is trivial. Moreover, the velocity fields considered
are not large in size so no special considerations are required to handle the
amount of data. Therefore, a simple function framework was devised to per-
142 Constrained Restoration of Tomographic PIV Data
form the iterative POCS procedure.
Projection operators as described in 6.1.2 are each implemented as separate
functions within, but independent of, the Tomographic PIV Toolbox (3.4.2).
A template for the projection function interface, which allows addition of an
arbitrary number of operators, is shown in figure 6.10.
Projection operators are called by an overarching algorithm whose properties
are discussed in section 6.4.2.
6.4.2 Data flow
A variety of algorithms are available to implement restoration in a simi-
lar way to POCS. One such algorithm is the Gerchberg-Papoulis (GP) tech-
nique [Gerchberg, 1974; Papoulis, 1975]. However, the work of Sezan & Stark
[1982] indicated that the *RELAX family of algorithms allowed a significantly
quicker convergence rate than the GP technique. Consequently the *RELAX
algorithm was implemented for this work.
6.4 The RELAX Algorithm and Variants 143
function [p3sx p3sy p3sz] = P#(sx, sy, sz, delta, options)
%P# Projection operator example function interface
% Detailed operator description
%
% Syntax:
% [p#sx p#sy p#sz] = P#(sx, sy, sz, delta, options)
%
% Inputs:
% sx [ny x nx x nz] 3D array containing x−components
% of the input vector field s
% sy [ny x nx x nz] 3D array containing y−components
% of the input vector field s
% sz [ny x nx x nz] 3D array containing z−components
% of the input vector field s
% delta [1 x 1] Spacing (arbitrary units) of the
% monotonic, regular grid on
% which sx, sy, sz lie.
% options structure Structure containing run control
% parameters for *RELAX algorithms
% and all operators.
%
% Note:
% nx is the number of grid points in the x direction,
% ny the number in the y direction, etc. Thus the
% x direction increases across the columns of the arrays
% (i.e. their second dimension) and the y direction
% increases across the rows (the first dimension) of the
% array.
%
% Outputs:
% p#fx, p#fy, p#fz
% [ny x nx x nz] x, y, and z components
% (respectively) of the field
% resulting from the projection
%
% PERFORM PROJECTION ON INPUT FIELDS
% ASSIGN TO OUTPUT
FIGURE 6.10 Function prototype for a projection operator. Use of this interface
allows straightforward implementation of additional operators
144 Constrained Restoration of Tomographic PIV Data
‘RELAX’ algorithm:   fn+1 = (1-λ)fn + λP(fn)
Input 
Converged?
Output 
D
ivergence = 0
B
and Pass
Lim
it
Length-scale 
Filter
Etc... Etc...
FIGURE 6.11 Flow of data within the *RELAX family of POCS algorithms
Figure 6.11 shows the structure of the *RELAX algorithm, of which there are
three variants (see 6.4.3 below). The algorithm is extremely simple, iterating
between (any number of) projections toward convergence. The algorithm
is guaranteed to converge provided the conditions on convexity, closedness
and intersection described in section 6.1.1 are met. The rate of convergence
(figure 6.12) is dependent on the strength of convexity as well as tangency
between the two (or more) sets, which is not known a priori. Convergence
rate can be controlled using relaxation parameters.
Starting Point
(Initial measurement / guess)
Slow convergence to final solution
FIGURE 6.12 Geometric explanation of the decrease in convergence rate with
decreasing tangency angle between two sets. Each green line
represents a projection.
6.4 The RELAX Algorithm and Variants 145
6.4.3 The *RELAX algorithm family
The family of RELAX algorithms consists of three closely related approaches.
The differences are marked by the relaxation factor λ which varies between
the three:
UNIRELAX The simplest of the family. In this case, there is no relaxation
of the algorithm for any operators (i.e. in figure 6.11, λ = 1).
RELAX A generalisation of UNIRELAX in that the relaxation can
vary from unity and is not necessarily the same for all op-
erators. As stated in Sezan & Stark [1982], using 0 < λi < 2
(where subscript i indicates the corresponding number of
the operator) can accelerate the convergence rate (see sec-
tion 6.4.4).
MULTIRELAX A full generalisation of RELAX developed for this work, in
which every member of the function set s = sxyz has an
independent relaxation value. Thus single values of λi be-
come scalar fields λixyz of the same size as the input func-
tion s. As per the RELAX algorithm, values vary as 0 <=
λixyz < 2. MULTIRELAX allows progress of the algorithm
to be adjusted for different regions of the input fields. For ex-
ample, a lower confidence (due to some known and quan-
tifiable source of error) in some elements of s can be ac-
counted for by increasing their relaxation factor accordingly
- thus weighting the restoration toward retaining elements
whose confidence score is high. A minor additional benefit
of MULTIRELAX is that the P2 operator described in section
6.1.2 can be trivially implemented by setting λ2xyz = 0 for all
(x, y, z) ∈ (Ix, Iy, Iz).
6.4.4 Selection of the relaxation parameter
It is suggested by Sezan & Stark that values of λi must be chosen based on
the tangency of two sets, since a low tangency angle slows the rate of conver-
gence (figure 6.12). The selection of λi is not trivial; at first sight the selection
146 Constrained Restoration of Tomographic PIV Data
of a value close to 2 seems ideal (assuming stability is achievable). However,
the benefit of fine-tuning relaxation parameters for this work was deemed
negligible, and possible implications (re. instability) of using λ > 1 were a
concern. The present work is therefore carried out with a value of λi = 1.0
for all i except where otherwise stated (e.g. section 6.5.1).
Section 6.5.3 provides an insight into the selection of λ distributions in order
to use the MULTIRELAX approach.
6.5 APPLICATION OF POCS TO TOMOGRAPHIC PIV DATA
6.5.1 Visualisation of data: Application of RELAX34
Statistically representative structural information can be obtained using tech-
niques such as conditional/ensemble averaging. This has been effectively
demonstrated for Tomographic PIV results by Schro¨der et al. [2011]. De-
pending on the process taking place, measurement noise may be diminished
through averaging - leaving a smooth field for visual interrogation as well as
application of quantitative processing techniques.
However, where TomoPIV is used for time-resolved visualisation of struc-
tural evolution (interrogating instantaneous fields) or where ensemble aver-
aging uses field variables other than mean velocities, the effect of measure-
ment noise propagates error into the results.
Denoising techniques can be a powerful method for diminishing the effect
of noise on both visualised fields (making evolution of structure with time
clearer to the human eye) and quantitative interrogation processes (making,
for example, determination of Lagrangian particle trajectories more robust).
Here operators P3 and P4 are applied using the RELAX algorithm frame-
work. The resulting technique is referred to as RELAX34. Results shown here
are therefore band-limited and divergence-free. Fields were iterated for 100
iterations of RELAX34, enough to consistently achieve normalised residuals
below 1× 10−4.
The RELAX algorithm allows different relaxation factors to be applied to dif-
ferent projections. λ3 was set to unity for the divergence-free projection oper-
ator, but the value of λ4 was set to 0.01 (divided by the number of iterations)
6.5 Application of POCS to Tomographic PIV Data 147
in order to scale the relative magnitudes of successive projections. This im-
proved the convergence rate after it was found that the P4 operator tended
to suppress the required iterative steps of P3.
The first check is that the divergence of the field actually reduces; this check
is made in figure 6.13, which shows the probability density functions for the
normalised divergence ratio ξ in fields before and after the application of
RELAX34. Regions of high gradient magnitude can produce relatively large
absolute divergence contributions; but they could be small relative to the lo-
cal gradient in the flow. Due to the large change in shear across the flow field
being assessed, the divergence is locally normalised as per equation 6.20,
such that reduction is represented for the entire field rather than in zones
of high shear. The convergence properties of the algorithm are such that it is
computationally prohibitive to eliminate all divergence; it is nonetheless sub-
stantially reduced as seen in the figure - the algorithm is particularly effective
at suppressing outlying regions of high divergence (unsurprising given the
least-squares nature of the problem).
ξ =
(∂u/∂x + ∂v/∂y + ∂w/∂z)2
(∂u/∂x)2 + (∂v/∂y)2 + (∂w/∂z)2
(6.20)
Observe the fields in figure 6.14. These represent slices of the velocity and
vorticity fields at the same position (parallel to the wall, z+=42) in the turbu-
lent boundary layer flow described in section 3.2. The velocity fields are very
similar; the post-POCS field is clearly smoother but retains all the key struc-
tures which a human observer picks out of the pre-POCS field. The vorticity
fields, in contrast, are very different. The differential nature of the vorticity
field is such that measurement noise is exacerbated. The POCS techniques
removes this noise, allowing clear visualisation of the structural information.
Erroneous ’vortices’ caused by noise and occasional false vectors (intention-
ally included in this example field for visualisation purposes) are substan-
tially diminished (except in extreme cases such as for groups of false vectors
or systematic error, which cannot be effectively removed by the scale filter or
divergence reduction operator).
Figure 6.15 illustrates structural information using isosurfaces of the second
invariant Λ2. Note that Λ is used here as opposed to the lower case version
used by Jeong & Hussain [1995] who review a variety of criteria, to avoid con-
148 Constrained Restoration of Tomographic PIV Data
0 0.5 1 1.5 2 2.5 3
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
Normalised Divergence Ratio ξ
PD
F
 
 
Before POCS
After POCS (Relax34)
FIGURE 6.13 Probability Density Functions for normalised divergence values
within an experimentally determined velocity field, showing
substantial reduction in divergence following application of
RELAX34
flict with the notation for Lagrange multipliers in the present text. Although
the trend is similar (e.g. areas of low and high turbulence remain clear be-
fore and after the POCS application), the decreased noise in the post-POCS
field not only diminishes ’false vortices’ resulting from noise (seen as many
small artefacts in the isosurface) - but also improves the apparent coherence
of structures which seem to be present in the input, but whose continuous
nature is broken by noise and artefacts.
For comparison, the input field is twice-smoothed before computation of the
invariants, using a 3 × 3 × 3 kernel with σ = 1.0. Figure 6.16 shows the
effect of this smoothing in 2D; 3D comparisons (using isosurfaces of Λ2) are
shown in figure 6.17. Small scale structures are merged together, cancelled
and diminished. It is clear that the RELAX34 method allows a useful degree
of noise cancellation/smoothing, whilst retaining coherency in the flow for
6.5 Application of POCS to Tomographic PIV Data 149
−20 −15 −10 −5 0 5 10 15 20
−20
−15
−10
−5
0
5
10
15
20
X (mm)
Y 
(m
m)
−20 −15 −10 −5 0 5 10 15 20
−20
−15
−10
−5
0
5
10
15
20
X (mm)
Y 
(m
m)
X (mm)
Y 
(m
m)
−20 −15 −10 −5 0 5 10 15 20
−20
−15
−10
−5
0
5
10
15
20
X (mm)
Y 
(m
m)
 
 
−20 −15 −10 −5 0 5 10 15 20
−20
−15
−10
−5
0
5
10
15
20
−80
−60
−40
−20
0
20
40
60
80
X (mm)
Y 
(m
m)
−20 −15 −10 −5 0 5 10 15 20
−20
−15
−10
−5
0
5
10
15
20
X (mm)
Y 
(m
m)
 
 
−20 −15 −10 −5 0 5 10 15 20
−20
−15
−10
−5
0
5
10
15
20
10
20
30
40
50
60
70
80
90
100
FIGURE 6.14 Velocity (mean removed, 1 in 4 vectors shown), vorticity magnitude,
and wall-normal vorticity component data for a turbulent boundary
layer presented as a wall-parallel slice at z+=42, pre (left) and post
(right) application of RELAX34. Note the consistency in overall
topology seen in the velocity field, but drastic reduction in noise
within the vorticity field allowing clearer visualisation of structural
information.
150 Constrained Restoration of Tomographic PIV Data
all measurable scales4 as opposed to just those dictated by the size of the
smoothing kernel.
Some dominant structures can be recognised in the pre-POCS visualisations
(of a derivative-based field, Λ2), but are discontinuous due to the effect of
noise and incorrectly resolved small-scale flow features in the flow. Follow-
ing application of POCS, a significant amount of noise and small-scale infor-
mation5 - the human eye is now able to interrogate and evaluate the flow
field, observing the dynamics and evolution of structures.
Although band passed (appearing smoother), the general topology of the
flow is not substantially altered by the application of POCS. This has sub-
stantial benefit (indeed is a prerequisite for the usefulness of this technique).
Two uses for future study are immediately apparent:
- The POCS technique can be used as a pre-processor for structural in-
terrogation techniques such as the Lagrangian trajectories computed
by Schro¨der et al. [2011], allowing a much more stable computation of
trajectories with little or no need for smoothing of the input.
- Instantaneous streamlines can be improved. Figure 6.18 shows the im-
provement in output of a streamline algorithm, in this case applied
some way from the wall (z+ = 42) using the same field visualised in
figure 6.15. In Chapter 7, this approach is used to close to the wall in
boundary layer data to provide skin-friction lines (which are stream-
lines, plotted in a plane close to the wall in the laminar sublayer); being
orthogonal to the vorticity vector, skin friction lines provide insight into
the flow topology close to the wall. Visualisation (and stable computa-
tion) of the lines and their corresponding topologies from experimen-
tal measurement can lead to an improved understanding of turbulent
wall-shear.
4Since the Reynolds Number in the present flow is high enough for turbulent eddies to
exist at scales below our grid spacing, then vortical motions at around the scale of the grid
should be present, as is the case in both the input and POCS-conditioned fields).
5Data is only removed at scales below the nyquist criterion; i.e. small scale flow features
are discarded if not measured at a sufficient resolution to be meaningful
6.5 Application of POCS to Tomographic PIV Data 151
FIGURE 6.15 Isosurface plots (3D, viewed from above) of the 2nd invariant
(Λ2 = −200) from raw results (top) and following application of
RELAX34 (bottom). Results are superimposed on a slice of the
velocity field (mean removed, 1 in 4 vectors plotted) from the input
field
152 Constrained Restoration of Tomographic PIV Data
FIGURE 6.16 Isosurface plot (3D, viewed from above) of the 2nd invariant
(Λ2 = −200), following application of a 3× 3× 3 gaussian
smoothing filter (σ = 1.0). Note the loss of fine-scale coherent
structures when compared to the output of RELAX34 in figure 6.15.
6.5 Application of POCS to Tomographic PIV Data 153
FIGURE 6.17 Isosurface plots of the 2nd invariant (Λ2 = −200) from raw results
(red) and following application of RELAX34 (green). Raw results are
shown together with similar isosurface taken from a smoothed field
(blue), indicating the dominant regions of swirl.
154 Constrained Restoration of Tomographic PIV Data
FIGURE 6.18 Instantaneous stream lines (shown in a plane at z+ = 42, z = 2.6mm)
are useful for observing flow topology. The use of POCS
reconstruction prevents noise from affecting algorithms which use
derivatives and tracking techniques to ascertain topological
behaviour. Top: Extents of flow field. Bottom: Zoomed-in region.
6.5 Application of POCS to Tomographic PIV Data 155
6.5.2 Divergence reduction: Using the P3 Operator
Where a requirement exists to reduce the divergence of the data without the
presence of smoothing, the P3 Operator can be applied on its own. In such a
case, the change to the field is not easily visible, since vectors are altered by a
small amount due to the differential nature of the correction. The change in
vector components is substantially smaller than the error associated with the
measurements as demonstrated by figure 6.19.
Error in measurement (based on 1/10 Voxel rule-of-thumb)
N
um
be
r o
f c
ou
nt
s
% Change in velocity (all components) due to divergence reduction (blue) and smoothing (red)
0
3
2.5
2
1.5
1
0.5
0
x104
-15% -10% 5% 0% 5% 10% 15%
FIGURE 6.19 Change in velocity components resulting from application of
UNIRELAX3 (blue). Also plotted are the bounds of error
measurement. For comparison, change in field resulting from
application of a 3× 3× 3 gaussian smoothing kernel is also shown
(red).
In itself, reduction of divergence may not be considered particularly impor-
tant. However, an improved topological structure is a direct corollary of the
divergence reduction technique. Worth [2010] used the RELAX3 algorithm
presented in section 6.2 on TPIV data whose resolution was approaching the
Kolmogorov lengthscale of a mixing-tank flow. The divergence reduction
156 Constrained Restoration of Tomographic PIV Data
process had a significant effect on the appearance of the QR plane, seen in fig-
ure 6.20: The now ’classical’ teardrop shape dicussed in Chong et al. [1990]
sharpens and appears closer to that of analytical or DNS based solutions.
Noise and random error in the results, as shown by Worth, causes a spread
of contours in the QR plane toward a more rounded shape - particularly af-
fecting the sharpness of the teardrop.
FIGURE 6.20 Improved topological classification in the Q-R plane, adapted from
Worth [2010]. (left) Before POCS. (right) After application of POCS
(RELAX3)
6.5.3 False vector fill-in: Using UNIRELAX with the P2 Operator
PIV in general (and Tomographic PIV in particular) produces false vectors
in the input field, identified by criteria such as the normalised median test
[Westerweel & Scarano, 2005]. False vectors are often identified at corners
and edges of the domain due to reduced light intensity (a characteristic of
many optical setups). Once identified and removed, these false vectors must
be replaced.
Typically, vector ’fill-in’ is performed using a form of interpolation such as
trilinear or tricubic. While straightforward, these methods can be poorly
conditioned where a neighbouring vector is false (but unidentified), there
are groups/clusters of false vectors and especially at the faces, edges and
corners of the domain where the interpolation must extrapolate (in 1, 2 and
3 directions respectively).
The P2 operator defined in section 6.1.2 acts as a mask, preventing some data
6.5 Application of POCS to Tomographic PIV Data 157
from being modified. A mask is defined for P2 which allows only those vec-
tors identified as false to be modified upon application of RELAX.
Application of UNIRELAX with P2 and P3 (denoted UNIRELAX23) allows
fill-in of the vector field using the requirement that the resulting local field
be divergence-free. The technique works robustly in the case that false vectors
are singletons. Where false vectors are grouped together, the equation system
becomes indeterminate in a manner similar to that discussed in section 6.2.3.
The implementation of UNIRELAX is such that a first guess is required. De-
spite the indeterminacy of the system, if the first guess is close to the ultimate
solution (e.g. obtained by linear interpolation or local averaging) the solution
is able to converge. However, in such a case the solution is highly dependent
on the first guess and therefore violates the core principle of POCS restoration
that data be restored/corrected using physical principles alone.
Use of UNIRELAX23 results in instability in two cases:
- where the first guess is far from an appropriate solution and there are
neighbouring false vectors to be filled.
- where there are large groupings of false vectors (regardless of the first
guess).
The latter may be somewhat suppressed by using the RELAX algorithm with
reduced values of λ but the approach is far from universal.
To address the non-physical nature of the restoration and to improve conver-
gence characteristics the P4 operator was included to form a new approach,
UNIRELAX234. The effect is that the solution at the sites of false vectors is
constrained to a band-limited version of the input field, which has a strong
stabilising effect.
The multi-scale nature of the wavelet-based P4 operator is ideal for this sce-
nario where clusters of vectors are missing - coefficients of the wavelet trans-
form for larger scales than the cluster size still provide band-limited infill
(effectively using a downsampled grid).
It should be noted that the system is still strictly indeterminate for clusters of
vectors but for most purposes stable restoration is possible, relying solely on
physical reasoning. Restoration is insensitive to the initial guess.
158 Constrained Restoration of Tomographic PIV Data
Figure 6.22 shows a slice of the result from application of UNIRELAX234. In
this case, a random selection of 17% of vectors in the input field were labelled
as false. In addition to those, 3% were highlighted as false using the nor-
malised median test (with a threshold of 2, see Westerweel & Scarano [2005])
on the raw experimental data. Thus a total of 20% of the original data was
discarded (velocity components set to 0) before restoration.
−100 −80 −60 −40 −20 0 20 40 60 80 100
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
Error (mm/s) between reconstructed and original velocity components
PD
F
 
 
εX
Error in original field (0.2 voxel)
−100 −80 −60 −40 −20 0 20 40 60 80 100
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
Error (mm/s) between reconstructed and original velocity components
PD
F
 
 
εX
Error in original field (0.2 vox)
−100 −80 −60 −40 −20 0 20 40 60 80 100
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
Error (mm/s) between reconstructed and original velocity components
PD
F
 
 
εY
Error in original field (0.2 voxel)
−100 −80 −60 −40 −20 0 20 40 60 80 100
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
Error (mm/s) between reconstructed and original velocity components
PD
F
 
 
εY/uY
Error of original field (0.2 vox)
FIGURE 6.21 Histograms showing the magnitude of difference between original
and filled-in vectors, showing the error margin associated with the
original vectors (assuming 0.2 voxel error). (Top) Error in X
direction. (Bottom) Error in Y direction. Higher error in X
components is assumed to result from the effect of ghost particles on
the noise in the original (input) data (see section 5.2.4).
Of the discarded (but valid) 17% of original vectors, each was restored with
a median error in velocity component of 5.8mms−1, where error was deter-
mined by comparing the input field (subject to experimental noise) to the
output of UNIRELAX234. The experimental error in the input field is slightly
larger than the average difference between input and reconstruction (exper-
imental error assumed to be 0.2 voxels displacement as discussed in section
6.5 Application of POCS to Tomographic PIV Data 159
3.2.2, which corresponds to a value of 6.5mms−1 in the present units). The
distribution of this ’reconstruction error’ is shown in figure 6.21.
To investigate the sensitivity of this approach to a high level of indeterminacy
in the equation system, restorations were made discarding different propor-
tions of the data from 3% (the actual identified false vectors only) to 100%.
Figure 6.24 shows the variation in median error with percentage of missing
data. Despite up to 70% of the data being missing, the field can be recon-
structed within twice the error associated with the input field. For higher
percentages, computations consistently diverge due to the almost complete
indeterminacy in the system.
Figure 6.23 is similar to figure 6.22 but with a higher percentage (50%) of
discarded vectors; showing the power of the technique in reconstructing a
large portion of the field. The same technique is useful for upsampling data
(obviously subject to the same band limit as the original data).
6.5.4 Weighted application of POCS: Using the MULTIRELAX algorithm
As noted in section 6.4.3, the MULTIRELAX algorithm can be used to weight
different regions of the field according to an ’accuracy function’ determined
a-priori.
For TomoPIV, the signal to noise ratio in the correlation volume provides a
metric of confidence in the accuracy (or at least robustness) of an individual
vector: False vectors are more likely to originate from correlation volumes
with low signal to noise ratio.
The signal to noise ratio for each vector is used to construct a ’weighting’
field as follows:
λixyz = 1− snrmax snr (6.21)
A slice through the input weighting field and results of applying the al-
gorithm (with P3 and P4 operators) to convergence (normalised maximum
residuals within 0.0001) is shown in figure 6.25.
A similar approach is made using the normalised divergence as a metric of
’how inaccurate’ an individual vector is. The results of this application are
shown in figure 6.26.
160 Constrained Restoration of Tomographic PIV Data
−20 −15 −10 −5 0 5 10 15 20
−15
−10
−5
0
5
10
15
X (mm)
Y 
(m
m)
 
 
Original Field (retained)
Original Field (discarded)
−20 −15 −10 −5 0 5 10 15 20
−15
−10
−5
0
5
10
15
X (mm)
Y 
(m
m)
 
 
Original Field (retained)
Reconstructed Field
FIGURE 6.22 Restoration of field with 20% missing data. Red vectors were zeroed
before application of UNIRELAX234
6.5 Application of POCS to Tomographic PIV Data 161
−20 −15 −10 −5 0 5 10 15 20
−15
−10
−5
0
5
10
15
X (mm)
Y 
(m
m)
 
 
Original Field (retained)
Original Field (discarded)
−20 −15 −10 −5 0 5 10 15 20
−15
−10
−5
0
5
10
15
X (mm)
Y 
(m
m)
 
 
Original Field (retained)
Reconstructed Field
FIGURE 6.23 Restoration of field with 50% missing data. Red vectors were zeroed
before application of UNIRELAX234
162 Constrained Restoration of Tomographic PIV Data
20 30 40 50 60 70 80 90 100
0
5
10
15
20
25
30
% vectors removed then reconstructed
M
ed
ia
n 
er
ro
r (
mm
/s)
 be
tw
ee
n r
ec
on
str
uc
ted
 
a
n
d 
or
ig
in
al
 v
el
oc
ity
 c
om
po
ne
nt
s
 
 
εX
εY
εZ
Input Error (≈ 0.2 voxel displacement)
FIGURE 6.24 Median error in restoration of a field with varying amount of
missing data.
X (mm)
Y (
m
m
)
-20 -15 -10 -5 0 5 10 15 20
-15
-10
-5
0
5
10
15
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Original Field
Local Relaxation Factor
X (mm)
Y (
m
m
)
-20 -15 -10 -5 0 5 10 15 20
-15
-10
-5
0
5
10
15
Restored Field (Iteration 100)
FIGURE 6.25 Restoration of data using MULTIRELAX. Highly weighted areas
indicate poor signal to noise ratio.
6.5 Application of POCS to Tomographic PIV Data 163
X (mm)
Y (
m
m
)
-20 -15 -10 -5 0 5 10 15 20
-15
-10
-5
0
5
10
15
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Original Field
Local Relaxation Factor
X (mm)
Y (
m
m
)
-20 -15 -10 -5 0 5 10 15 20
-15
-10
-5
0
5
10
15
Restored Field (Iteration 500)
FIGURE 6.26 Restoration of data using MULTIRELAX. Highly weighted areas
indicate high divergence in the input field.
There are two aspects to ’inaccuracy’ in an input field: false vectors not re-
moved in earlier processing of the field and the noise inherent in the measure-
ments. Expecting the POCS technique to remove false vectors is equivalent to
stating that the weighting function is a good metric for identification of false
vectors, which may not be true. The various false vector identification crite-
ria [Raffel et al., 2007, p.185] have not been tried as input weightings since
their determination is not necessarily physical in nature.
Other than not removing obvious false vectors in the field, it is apparent that
the output fields following application of MULTIRELAX are substantially
more noisy than the inputs. One of two conclusions can be drawn:
- Neither the signal to noise ratio nor the divergence present in the ini-
tial field are reliable indicators of which vectors are most affected by
measurement noise
or
- The application of MULTIRELAX violates convexity in a fashion to that
noted in 6.5.3 and therefore does not converge to the correct result.
164 Constrained Restoration of Tomographic PIV Data
0 10 20 30 40
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
Iteration Number
R
el
at
iv
e 
R
es
id
ua
l
Convergence not guaranteed
 in a poorly−conditioned case
FIGURE 6.27 Unsteady convergence of MULTIRELAX.
Although no applications of MULTIRELAX were divergent in the sense of a
continually increasing residual from the computation, the convergence pro-
cess was seen to be less steady (figure 6.27) than in ’normal’ applications,
which exhibit continuous convergence. Moreover, since divergence in the
input vector field must occur as a result of measurement noise, the latter con-
clusion is deemed to be the valid one.
6.6 Review of Operators 165
6.6 REVIEW OF OPERATORS
The operators considered in previous sections are reviewed here in order to
provide an overview of the capabilities established within this work.
Operators should be viewed as independent constraints and can be applied
separately or in conjunction within the RELAX algorithm framework.
P1 Uses a spectral decomposition to apply a spatial band-limit to input
data. Although effective at denoising data, the technique is sensitive
to the presence of false vectors. Periodic artefacts are to be expected
(without special handling of boundary points in the input data).
P2 Allows fixed values to be applied to the input. This is used in conjunc-
tion with other operators to restore data elements known to be missing
or corrupt by keeping surrounding (valid) data elements fixed.
P3 Reduces divergence in the input field, to meet the incompressibility
condition.
P4 Uses a wavelet decomposition to apply a scale limit. This removes
noise which is not spatially correlated, together with components of
the flow field existing at scales smaller than the spatial resolution.
Of these four operators, the latter two (P3 and P4) are of greatest significance,
as they can be applied to general Tomographic PIV configurations without
special treatment (assuming incompressibility in the fluid).
Operator P4 is a preferable alternative to the P1 operator since:
- Spurious (false) vectors are effectively removed with minimal effect on
the underlying flow topology
- The cut-off scale is precisely defined at the Nyquist criterion without
the ringing effects associated with top-hat filters
- Spurious/false vectors are effectively removed (even in clusters, pro-
viding clusters have limited extent in at least one dimension) without
substantial effect on the underlying flow topology.
166 Constrained Restoration of Tomographic PIV Data
Additional operators are be the subject of further investigation. In particular,
a Savitzky-Golay filter (whose application by Elsinga et al. [2010] is discussed
in 2.3.1) could be considered as an alternative toP4 for applications where the
requirement for a strict band limit (as applied by P4) is of less importance
than smooth derivative fields.
7 Coherent Structures in High
Reynolds Number Boundary Layers
The strength of Tomographic PIV lies in the ability to access and visualise
topological data. Despite the presence of bias error (see chapter 5), topologi-
cal information from the measured flow fields is still highly valuable.
Here, we use results from the experiments described in chapter 3 to highlight
the combined power of TPIV techniques introduced elsewhere in this work.
7.1 REVIEW OF BOUNDARY LAYER AND EXPERIMENTAL PARAMETERS
Tomographic PIV measurements were conducted in a high Reynolds Num-
ber turbulent boundary layer. The experiments and setup are described in
detail elsewhere (see chapter 3) but it is worth reviewing the setup geometry
and boundary layer properties here.
7.1.1 Boundary Layer and Measurement Volume
To form comparison with existing work on boundary layers which utilises
Tomographic PIV, the present experimental parameters are tabulated along-
side those from Schro¨der et al. [2011] (whose setup and flow parameters,
amongst various sources in the literature1, are most similar to the present
tests) and from Elsinga et al. [2010] (whose boundary layer experiment is the
only work to date performed at a significantly higher Re than considered
here). Figure 7.1 shows the present measurement volume (expressed in wall
units) in the context of these two previous experiments: although the inter-
rogation volume size is geometrically similar between all three experiments,
the boundary layer thickness in the present facility has allowed more detailed
visualisation of the near-wall flow than in the previous experiments.
1See section 2.3.6 for a detailed review of existing work.
167
168 Coherent Structures in High Reynolds Number Boundary Layers
Present work Schro¨der et al. Elsinga et al.
U∞ 0.45 0.53 Mach 2 m/s
δ99 0.115 0.038 0.020 m
uτ 0.019 0.0219 19.5 m/s
+ 6.1× 10−5 5.4× 10−5 2.8× 10−6 m
δ+ 1890 800 7080
θ 9.6× 10−3 5.5× 10−3 1.2× 10−3 m
Reθ 3657 2460 34000
δGRID 0.52 (8.6+) 0.69 (15+) 0.49 (177+) mm
LVOLUME 40.8 (669+) 65 (1380+) 70 (25000+) mm
DVOLUME 7.8 (128+) 15 (328+) 6.5 (2301+) mm
TABLE 7.1 Comparison of experimental parameters between the present work and
the TPIV-based boundary layer studies of Schro¨der et al. [2011] and
Elsinga et al. [2010]. δGRID is the spatial resolution (here taken to be the
spacing between adjacent grid points). All experiments utilise a 75%
window overlap. LVOLUME and DVOLUME are the length and depth of
the reconstruction volume respectively.
7.1.2 Basis of Comparison
In chapters 4 to 6, a variety of enhancements to the Tomographic PIV process
are presented, including:
- The Correlation Tracking Enhancement for improved accuracy
- The Weighting Reduction Scheme (WRS) allowing efficient solution of
the tomographic reconstruction problem
- An image optimisation methodology for improved conditioning of the
reconstruction process
- A framework (’RELAX’) for vector field post-processing based on Pro-
jection Onto Convex Sets
- A divergence reduction operator for use within the POCS framework
- A wavelet-based scale limiting filter operator for use within the POCS
framework
7.1 Review of Boundary Layer and Experimental Parameters 169
Equivalent size (using wall scaling) of present volume
FIGURE 7.1 Measurement volume (in wall units) of the present Tomographic PIV
experiments superimposed onto figures from Schro¨der et al. [2011]
(top) and Elsinga et al. [2010] (bottom).
170 Coherent Structures in High Reynolds Number Boundary Layers
To evaluate net effect of improvements to the process, a ’baseline’ standard
approach was used which consisted of image processing, MART reconstruc-
tion and VODIM-based cross correlation. For comparison (here referred to
as the ’enhanced’ approach) we present results utilising the same image pro-
cessing and reconstruction technique, but comprising a CTE-based cross cor-
relation (section 4.1) and application of the RELAX34 divergence & wavelet
noise reduction algorithm (section 6.4.3).
The image optimisation methodology (section 4.3) was found to be ill-posed
in its current form - leading to excessive bias error in final vector fields. In
light of this concern the developed technique was not ultimately used and
more typical image preprocessing practices were applied.
As discussed in section 4.2.5, solvers for the reduced reconstruction were im-
plemented in prototype form (aspects of computational science being outside
the scope of this work) so for the present comparisons, experimental data (see
chapter 3) was processed using a MART reconstruction algorithm.
7.1.3 Summary of Experimental Method
Image preprocessing operations (as used in the majority of TPIV experiments
in the literature) are applied, comprising: background subtraction; bypass
filtering; gain (to equalise images) and gaussian blurring.
The MART reconstruction algorithm was applied in 5 iterations with relax-
ation parameter µ = 0.9. No smoothing was applied to the reconstruction
volume between iterations.
Window sizes for the comparison between standard and enhanced techniques
are selected based on the minimum size which is robust and which con-
tains at least three particles - this corresponds with windows 323 in size for
the setup used. When using the CTE process, results can be achieved with
smaller windows than this: for windows 243 voxels in size, approximately
one true (i.e. non-ghost) particle is contained in each window (see figure 7.2).
Despite apparently robust cross correlation (structures are consistent between
the two independent time-steps shown in figure 7.2), the errors due to bias,
intensity variation and position of particles within windows discussed by
Nobach & Bodenschatz [2009] become excessive at this window size.
7.1 Review of Boundary Layer and Experimental Parameters 171
FIGURE 7.2 Field showing isosurfaces of λci = −200, following use of CTE
technique with final pass window size of 243. Green and blue surfaces
are taken at times t = 0 and t = 3δt respectively; i.e. are computed
from independent sets of images.
The use of windows 323 voxels in size with the standard technique is not
wholly robust for this case but is nevertheless possible. For both cases, a five-
pass interrogation was used to allow the final window size to be reached
stably: Passes used window sizes of [643, 413, 323, 323, 323] and overlaps of
[0%, 0%, 0%, 75%, 75%] respectively.
It was found that progressing through to the final window size using zero
overlap is an effective method for reducing clustered groups of false vectors
in the final result. Groups can otherwise result from a spurious feature be-
ing present in several adjacent windows and thus not identified during false
vector identification stages.
Both enhanced and standard techniques used identical false vector selection
criteria. The normalised median test was applied with noise parameter of
0.15 and threshold of 2. In addition, absolute values were used to identify
false vectors; limits on the maximum absolute value (for each velocity com-
ponent) were applied outside the range of plausible values and vectors dis-
carded if the limit was exceeded by any component.
The value of δt used was 2ms for the ’standard’ analysis; since this was found
172 Coherent Structures in High Reynolds Number Boundary Layers
by trial and error to be the largest value for which computation of the vec-
tor field was robust using VODIM cross correlation (see the more extensive
discussion on this aspect in section 3.2.2).
Selection of δt for the ’enhanced’ analysis is less trivial. It is tempting when
showcasing results, to select an increased δt of 4ms - the improved robust-
ness of the CTE technique allows use of this value without destabilising the
algorithm. This is a considerable advantage; doubling the dynamic range of
the technique even before the improved accuracy of peak location is consid-
ered. However, for direct comparability with the standard analysis, δt was
selected to be 2ms.
For application of the RELAX34 postprocessing algorithm, the only tuneable
parameters are the relaxation factors for theP3 and P4 operators which were
selected as λ3 = 1.0 and λ4 = 0.01. Selection of unity relaxation factor for
P3 follows directly from the UNIRELAX algorithm. Since 100 iterations of
the algorithm were used (sufficient to reduce divergence to within the finite
sampling and truncation error of the differencing algorithm), the selection
of λ4 = 1/100 weights the application of the wavelet filter (intended as a
single-application filter) across the iterations.
7.2 COMPARISON OF STANDARD AND ENHANCED TECHNIQUES
7.2.1 Reduction of Divergence
Figure 7.3 shows isosurfaces of divergence. The standard approach clearly
yields nonzero divergence in the results, despite incompressibility of the flow.
The maximum divergence in the standard field is 827.0s−1 compared to a
maximum of 6.1s−1 in the enhanced field; an improvement of×134. The reg-
ular pattern shown for the enhanced field in figure 7.3 oscillates about every
other grid point. This oscillation is caused by the remaining finite sampling
and truncation error: Central differencing is used to take velocity deriva-
tives in the implementation used here; if necessary, higher order schemes
for computation of derivatives or interrogation methods which directly yield
derivatives (such as those discussed by Scarano 2004) can be used to reduce
the residual divergence further.
7.2 Comparison of Standard and Enhanced Techniques 173
FIGURE 7.3 Comparison of field divergence between standard (left) and enhanced
(right) TPIV analyses. Isosurfaces are plotted at 413.5 and 3.06
respectively; an improvement of ×134 in 100 iterations of RELAX34.
7.2.2 Velocity Field Comparison
Figures 7.4 and 7.5 show velocity fields for the standard (VODIM) and en-
hanced (CTE and RELAX34) cases. The images, pre-processing, reconstruc-
tion and PIV options (such as false vector selection thresholds) and visuali-
sation methods are otherwise identical.
Two fields are produced from the same images with identical processing and
false vector selection options; one using the standard analysis technique and
one using the enhanced technique. The velocity field, Q-criterion and skin
friction lines are shown to compare the abilities of the two techniques.
In figure 7.4, the vector plots (slice positioned at the closest plane to the wall
at z+ = 21.4) show a marked difference between the standard and enhanced
cases; with the standard case appearing considerably less continuous and
having substantially greater noise.
An isosurface of the Q-criterion is used to identify swirling motion in figure
7.5. The same isovalue is used for both fields. The enhanced field shows
coherent flow structures consistent with current theory (section 2.2.1)such
as low speed regions and vortex packets - flow features visualised with the
enhanced technique are discussed in section 7.3. In contrast, little or no co-
herency is seen in the standard field: since the Q-criterion is based on field
derivatives, noise associated with the technique dwarfs meaningful features
which may be present in the data.
174 Coherent Structures in High Reynolds Number Boundary Layers
−300 −200 −100 0 100 200 300
−300
−200
−100
0
100
200
300
X+
Y+
−300 −200 −100 0 100 200 300
−300
−200
−100
0
100
200
300
X+
Y+
FIGURE 7.4 Vector slice plots at Z+ = 21.4, with mean velocity removed for
standard (top) and enhanced (bottom) TPIV analyses. Scaling is
arbitrary between the two figures due to the differing bias error (i.e.
different apparent mean flow). 1 in 4 vectors shown for clarity.
7.2 Comparison of Standard and Enhanced Techniques 175
FIGURE 7.5 Isosurfaces of Q at 5×median(Qenhanced) = 380 for standard (top) and
enhanced (bottom) TPIV analyses. Both fields are produced from the
same images with identical processing and false vector selection
options. Field is the same as that shown in figure 7.4.
176 Coherent Structures in High Reynolds Number Boundary Layers
Red paths indicate vortex packets.Skin friction line 
paterns are consistent with the efect of packets 
close to the wal.
FIGURE 7.6 Skin Friction lines for standard (top) and enhanced (bottom) analyses.
Isosurfaces of Q are shown (with transparency) for comparison of
topology.
7.3 Observation of Structural Elements in the Flow 177
The skin friction lines shown in figure 7.6 are computed using the 2D instan-
taneous velocity field in a plane parallel to (and as close as possible to) the
wall at Z+ = 21.4. A lagrangian particle tracking routine is applied to the
velocity field to determine the lines, which can also can be thought of as in-
stantaneous streaklines parallel and close to the wall. A frame of reference
with the time-mean flow removed is used.
Computed in this way, the lines help to highlight the topology of the near-
wall flow and accentuate flow features. Similarly to figures 7.4 and 7.5, there
is significant improvement between the standard and enhanced methods. In
the lower (enhanced) frame, the Q field corresponds well with expected flow
topology seen in the skin friction lines - this is not the case in the standard
technique (upper frame) where regions of swirl identified by the Q-criterion
frequently do not correspond to the topology expected from the skin friction
lines. This agreement for the enhanced technique underlines the findings
in section 6.5.1 in which particle tracking type techniques were found to be
more robust following application of RELAX34.
7.3 OBSERVATION OF STRUCTURAL ELEMENTS IN THE FLOW
The boundary layer measurements presented here highlight a number of key
features consistent with the observations of Theodorsen [1952] who discusses
the concept of a vortex ’arch’, Offen & Kline [1975] and Hinze [1975] who dis-
cuss ’sweeps’ and ’bursts’ in the logarithmic region of a turbulent boundary
layer, and a variety of other sources in the literature (see section 2.2.1) which
discuss the emergent idea of vortex ’packets’.
7.3.1 Hairpin packets and near-wall streaks
Figure 7.7 shows structures in a flow with low-speed regions. The green
isosurfaces effectively highlight packets of hairpin vortices. The action of
hairpin vortices in the near wall region is inextricably linked with the charac-
teristic skin friction velocity uτ in a turbulent boundary layer. uτ is thus used
to set the isosurfaces of velocity (at 5× uτ) indicate low speed fluid which
has been lifted [Offen & Kline, 1975] from the viscous layer away from the
wall - this is a low speed near-wall streak.
178 Coherent Structures in High Reynolds Number Boundary Layers
FIGURE 7.7 An isosurface of Q (green) is shown with an isosurface of velocity
(blue). (top) Field at time t = 0s. (bottom) Field at time t = 20s
7.3 Observation of Structural Elements in the Flow 179
The two visualisations in figure 7.7 belong to the same dataset but are taken
20s apart in time; intermediate fields have a consistent low speed region. The
streak with its corresponding series of arched hairpins prevails for more than
80s until the end of the data capture, although is distorted in some frames by
the effect of larger scale motions passing above (as seen in figure 7.8 taken
at t = 60s). Simple application of Taylor’s ’Frozen Eddy’ hypothesis over
80s would suggest that the structure were longer than the water tunnel itself;
indicating that an autogeneration process is perpetuating the streak.
The subject of extremely long, auto-generating trains of hairpins (Very Large
Scale Motions or VLSMs) is frequently discussed within the literature (see
section 2.2.1) but commonly refers to larger scale motion than could not be
visualised in this volume. Near the wall, Offen & Kline [1974] do discuss a
self-perpetuating process for smaller scale (inner region) structures. Using
high speed PIV, Adrian et al. [2000] report the occurrence of regular vortex
packets in the near-wall region of the flow spanning streamwise lengths up
to 2δ. However, the evidence of a long, low-speed streak shown here with
its associated long train of vortices suggests that that motions with scale far
greater than the length of the tunnel itself could be possible even at these
small scales.
Although the motions are distorted by larger scale flow features, these very
long, small motions do not appear to meander significantly (i.e. change loca-
tion in the cross-stream direction over time). This observation would require
much longer test times to confirm - since structures are close to the wall,
meandering is likely to be linked to the wall (rather than the outer) velocity
scales and therefore occur over greater timescales than observed here.
In all visualisations, the flow was sharply segregated into patches of high
turbulent intensity or quiescent flow, consistent with the observations of Of-
fen & Kline [1974] and Ganapathisubramani et al. [2003]. See for example
figure 7.8 in which the turbulent intensity is significantly lower (few visible
structures at or above the chosen value of Q) outside the vicinity of the train.
Offen & Kline [1974] utilise dye-flow and hydrogen bubble experiments to
visualise a ’bursting’ process in the turbulent boundary layer. ’Sweeps’ are
highlighted, where high speed fluid is drawn downward to the wall. Regions
where low speed fluid is lifted are also identified (seen as blue isosurfaces in
180 Coherent Structures in High Reynolds Number Boundary Layers
FIGURE 7.8 An isosurface of Q (green) is shown with an isosurface of velocity
(blue). This field is taken at time t = 20s after those in figure 7.7. Note
distortion of the train by the presence of larger structures (it is slanted
to the side) and quiescent regions away from the central, highly
turbulent, train.
figures 7.7 and 7.8). The ’outer’ dye stream in their experiments is positioned
at location z = 100+ to effectively highlight the sweeping motions observed.
In figure 7.9 a similar study is made using streamlines in a plane perpendic-
ular to the wall. Being instantaneous (i.e. for a snapshot of the flow), the
streamlines are not equivalent to the dye lines in the work of Offen & Kline
which are Lagrangian paths. Nevertheless, the streamlines indicate the direc-
tion of flow normal to the wall so the sweep to which Offen & Kline referred
is highlighted. The packet of vortices causes a small series of sweeping mo-
tions which corresponds to the ’sawtooth’ of sweeps described by Offen &
Kline. A uniform velocity equivalent to the mean flow at height z = 100+ is
removed before calculation of the streamlines to highlight the topology and
the vertical component of fluid motion.
7.3 Observation of Structural Elements in the Flow 181
FIGURE 7.9 A field is visualised using an isosurface of Q (green) and isosurface of
velocity (blue) to highlight vortices and low speed near-wall streaks
respectively. The red box indicates a small packet of hairpin vortices
which are isolated in the lower image. Instantaneous streamlines are
used to highlight the sweep and ejection regions.
182 Coherent Structures in High Reynolds Number Boundary Layers
7.3.2 Near-wall topology
A useful technique for visualisation of topology - mentioned in section 6.5.1 -
is the use of streaklines close the the wall (or ’skin friction lines’). Figures 7.10
and 7.11 show skin friction lines (for an entire field and for a hairpin train,
respectively) at the wall, together with the now-familiar isosurfaces of Q and
ux. Note that strong foci (regions of intense Q) are frequently reflected in the
skin friction lines (although Q is galilean-invariant, whereas the streaklines
are not - in this case the frame of reference chosen is moving with the time-
average flow speed at Z+ = 21.4).
Both figures 7.10 and 7.11 show part of the low-speed streak discussed in
the previous section. The topology and arrangement of vortex arches high-
lights how regular the inner flow can be, despite the high Reynolds Number
- Schlatter & Orlu [2010] remark on the chaotic appearance of the outer flow
at a Reynolds Number (Reθ = 2500) lower than that of the present measure-
ment (Reθ = 3650).
A key implication of regularity in the inner flow is that creation of turbu-
lence occurs in a structured way close to the wall, even as Reynolds Number
increases (from the lower values for which computational visualisations are
currently available). This observation builds increased confidence in predic-
tive models for high Reynolds Number flow which rely on structural infor-
mation, such as those of Marusic et al. [2010b] and Mathis et al. [2011].
An instance of a small hairpin group standing alone from a larger train was
found - this is visualised in the same way in figure 7.12, showing the detailed
topology surrounding hairpins close to the wall.
7.3 Observation of Structural Elements in the Flow 183
FIGURE 7.10 Skin friction lines at Z+ = 21.4 for an entire field, shown with
isosurfaces of Q and ux to highlight key motions.
184 Coherent Structures in High Reynolds Number Boundary Layers
FIGURE 7.11 Skin friction lines at Z+ = 21.4 beneath a train of hairpins.
7.3 Observation of Structural Elements in the Flow 185
FIGURE 7.12 Skin friction lines at Z+ = 21.4 beneath a small group of otherwise
isolated hairpins.
186 Coherent Structures in High Reynolds Number Boundary Layers
7.3.3 Appearance of closed loops of vorticity
FIGURE 7.13 The enhanced flow field of figure 7.5, with isosurface Q = 380.
Isosurface is coloured with the wall-normal height Z. Note the
vortex ring structures whose axis is in the streamwise direction
An interesting feature is highlighted by figure 7.13. In this field, several
’trains’ of vortices are present. The isosurface threshold and colour is se-
lected so that two strong trains are highlighted in the image (away from the
wall, toward the red end of the colour scale).
Consistent with current theory, it appears that vortex arches are entrained
in the streamwise direction. However, arches (according to current theory)
have extended legs - hence the name ’hairpin’ was coined to describe the
structures. In this visualisation, the ’arches’ appear to be full vortex rings.
The vortex rings occur frequently - a series of fields was visualised and at
least 60% contained one or more rings. Rings typically occur in regions of
high turbulent intensity rather than the quiescent Offen & Kline [1974] as
’quiescent regions’ of little activity (descrirains of rings are less prevalent
than single rings and usually consist of 4 or less rings (only one train of more
7.3 Observation of Structural Elements in the Flow 187
than 4 rings, pictured in figure 7.13, was seen).
Axes of the vortex rings are aligned in the streamwise direction - the flow is
similar in nature to a jet. The characteristic diameters and axis heights of the
ring in figure 7.14 are 40+ and 90+ respectively. These values can vary by 10+
for general rings seen in the flow.
FIGURE 7.14 An isosurface of Q shows an isolated vortex ring in the field.
Streaklines are plotted at constant height, in a frame of reference
adjusted to be local to the ring’s velocity, showing the typical
cross-sectional visualisation of a vortex ring. The red quiver plot
shows discretisatoin of the measured field.
An individual ring and a train of rings are visualised in figures 7.14 and 7.15
respectively, using an isosurface of the Q-criterion. Streaklines (of constant
Z, in a frame of reference local to the ring) are shown in figure 7.14 to high-
light the flow field. However, it is important to note the quiver plot (yellow
arrows) also contained in figure 7.14, highlighting the spacing of vectors at
75% window overlap: The spatial resolution is not sufficient to visualise the rings
using independent measurement points. This calls into question the validity of
the measurement in this case; further experiments whose spatial resolution
is at least a factor of 2 smaller must be made (either in a larger boundary
layer, or using a different optical setup) before conclusion can be drawn on
whether rings are truly present.
188 Coherent Structures in High Reynolds Number Boundary Layers
FIGURE 7.15 Isosurfaces of the Q-criterion (Q = 3000) are used to highlight a train
of vortices, extracted from figure 7.13.
If the ring structures exist, it is speculated that they are the tops of hairpin
vortices whose legs have undergone vortex reconnection. In contrast, struc-
tures whose characteristic spanwise width is less than their distance from the
wall do not form vortex rings and appear as arches: Packets of these vortex
arches can be seen in later figures (7.7 and fig:streak2) but almost always have
a characteristic height less than or equal to their width. Trains of vortex rings
are more prevalent (or at least as prevalent) in the data than large structures
(i.e. hairpins too large to wholly visualise) whose legs extend to the floor in
the manner expected.
A non-physical explanation for the rings lies with a cluster of false vectors2
pointing in the upstream direction. Although typically false vectors are as-
sociated with shear rather than rotational motion, it is possible for the Q-
criterion to identify such a cluster as a vortex ring.
An alternative (less likely) explanation is that the rings could relate to bias er-
ror affecting the bottom of arches in such a way that they appear to reconnect
- however, bias error (see section 5.1.3) tends to have an effect over longer
lengthscales than the 1 to 3 grid points over which the topology is chang-
ing the the bottom of the rings. It is therefore difficult to imagine bias error
affecting structural form in this way.
2The source must be a cluster of several vectors - at least two in each direction - since
single vectors are either identified during false vector selection or filtered out in the RELAX34
procedure
7.4 Preliminary Conclusions on Flow Structure 189
7.4 PRELIMINARY CONCLUSIONS ON FLOW STRUCTURE
In most cases, improvements made to a measurement technique are con-
firmed and verified by additional work on the new technique across a broader
scientific community. Since the main thrust of this work has been toward im-
provement the measurement technique, it seems unwise (and is outside the
present scope) to present detailed conclusions on the actual fluid flow within
the same piece of work.
However, the work has produced visualisation of high Reynolds Number
flow topology in the near-wall region of the boundary layer, whose spatial
resolution clearly exceeds that of existing experimental work. From the ob-
servations made above, some preliminary conclusions can be drawn which
may be useful in further discussion of this type of flow field:
- Very Long Near-wall Streaks are present, consisting of trains of hairpin
vortices within the buffer/overlap layer. These long, small structures
are highly regular, exist over lengthscales longer than the tunnel (i.e.
are perpetuated by an autogeneration mechanism), do not appear to
meander significantly and co-exist with much larger structures within
the flow.
- Hairpins and Packets have been visualised at small scale close to the
wall. Typically, the width of the arches (in the span-wise direction)
is around 60 wall units. The height of the arches in the wall-normal
direction is almost always less than the width.
- Vortex Rings have been seen in the flow, entrained with axes oriented
in the streamwise direction. It is hypothesised that these are formed
due to vortex reconnection of hairpin legs. The rings are not visualised
with independent measurements so may arise due to clusters of false
vectors; further investigation at a higher spatial resolution is required
to conclude on their existence.

8 Conclusions
8.1 ON SOFTWARE
A state-of-the-art software toolbox (for MATLAB) has been developed to per-
form Tomographic PIV, including functions incorporating the techniques of
chapters 4 and 6. The toolbox has been rigorously validated as part of the
present work and is available free to the scientific community. It includes:
- A wide range of image processing capabilities (independent of, but
benefiting from, MATLAB’s own Image Processing Toolbox).
- GUI-based calibration and setup routines.
- Parallelised reconstruction functions utilising MART [Herman & Lent,
1976], the MFG [Worth & Nickels, 2008], and the Wrs matrix reduction
technique described in 4.2.
- Parallelised cross correlation functions for VODIM [Scarano & Rieth-
muller, 2000] with multiple pass algorithms and Whittaker/Cardinal
window deformation.
- Post-processing routines including band-pass, Savitzky-Golay and lo-
cal polynomial regression based filtering1 together with visualisation
tools based on MATLAB’s inherent plotting capabilities.
The novel developments presented within this work are also available within
the toolbox:
- The matrix reduction (WRS) technique described in 4.2.
- The Correlation Tracking Enhancement of section 4.1.
- The Projection Onto Convex Sets code framework and operators de-
scribed in chapter 6)
1Filter functions, with the exception of wavelet based filtering, are cast to N dimensions
allowing a fourth dimension - time - to be incorporated into the filtering process
191
192 Conclusions
8.2 ON ENHANCEMENTS TO TOMOGRAPHIC PIV
8.2.1 The Correlation Tracking Enhancement
The CTE technique has been formulated, described and tested both experi-
mentally and numerically. Random error was reduced by a factor of at be-
tween 2 and 7 for the majority of particle densities considered. In cases of
highest particle density, CTE works where VODIM algorithms fail due to ex-
tremely noisy correlation volumes.
An improvement in spatial resolution has been demonstrated, with the tech-
nique able to resolve flow to acceptable accuracy using a smaller window
size than is possible with a VODIM algorithm (subject to the limit that parti-
cle density is sufficient to allow a decrease in window size).
Dynamic range of the measurement technique is improved by a factor of be-
tween 2 and 7, due to decreased random error in the measurements. Using
CTE allows robust increase of timestep δt for PIV in a highly turbulent flow
(see section 3.2.2). Taking this into account, the dynamic range is improved
by a factor of at least 4 overall.
Bias error associated with ghost particles in the reconstruction was reduced
by approximately 40% in the majority of the volume.
The CTE technique increases computational demand compared to existing
TPIV by less than or equal to a factor of 2.
8.2.2 The Weightings Reduction Scheme
A Weighting Reduction Scheme (WRS) has been developed to remove re-
dundant coefficients from the tomographic reconstruction weightings ma-
trix. This facilitates extremely fast implementation of tomographic recon-
struction on multi-and-many-core hardware, utilising existing sparse matrix
libraries and techniques.
Solver implementations of the Simultaneous MART and Conjugate Gradient
(CG) solution techniques (for the reduced problem) are described.
Using a sub-optimal Conjugate Gradient solver, the solution time and qual-
ity using a WRS scheme is shown to be broadly comparable with a MART
8.3 On the Accuracy of Tomographic PIV 193
computation. Significant speed-up relative to both MART and SMART tech-
niques is expected in an optimal implementation.
The Weighting Reduction Scheme is shown to be useful at a typical particle
density of 0.05ppp but becomes impractical at higher densities.
8.2.3 Image pre-processing
An attempt is made to formulate a consistent approach to image preprocess-
ing for Tomographic PIV, using a Nelder-Mead optimisation. The attempt
appears successful, yielding robust cross correlations and exhibiting well-
bounded, convergent optimisation.
However, discussion reveals a concern that these benefits are achieved to the
detriment of bias error in the solution. The technique is therefore not recom-
mended for general application until a well-posed alternative objective func-
tion is found (possibly utilising sparsity maximisation for a given allowable
percentage of false vectors).
8.3 ON THE ACCURACY OF TOMOGRAPHIC PIV
The Tomographic PIV technique, including the CTE modification has been
validated experimentally by direct comparison with measurements made us-
ing a 2D PIV apparatus (synchronised with the TPIV image acquisition) and
numerically through reconstruction of artificially generated data.
8.3.1 Experimental Findings
Bias is found to be a significant source of error for this case, due to the high
mean shear in the flow field. From experimental results, the peak value of
bias error is 0.5 voxels very close to the wall for the CTE technique. Bias
error is worse for VODIM analysis. However, for the majority of the volume
(where shear is lessened) bias error is less than 0.3 voxels (VODIM-based
techniques) or 0.2 voxels (CTE-based techniques).
The CTE technique is shown experimentally to improve the bias error com-
pared to VODIM analysis by 40% (away the ’mean flow position’).
194 Conclusions
The CTE technique was found to be robust at higher values of timestep δt
than VODIM, allowing increase in dynamic range of the technique by a factor
of 2 on top of increases through improved accuracy (see below).
8.3.2 Numerical Findings
An analytical flow field was used to artificially generate images in order to
numerically predict the accuracy of TPIV.
Section 5.2.4 reveals that bias error is over-predicted by the numerical study
compared to experiments due to the effect of turbulence (in the experimental
flow) de-correlating ghost particles. A corollary is that presence of turbulence
in a flow reduces bias error in the mean flow statistics. However, turbulence
must introduce local bias effects in instantaneous fields which have not been
quantified here.
Random error was evaluated within the numerical study. The study was
found to under-predict the level of noise in a Tomographic PIV measurement
of an arbitrary flow, due to the presence of strongly correlated ghost particles.
Results are expected to have the same level of under-prediction, so relative
performance between different techniques can be evaluated.
In all tests, CTE substantially outperformed VODIM - random error was re-
duced by a factor varying between 2 and 7 depending on the particle sizes
and densities used (figures 5.10 and 5.11). Dynamic range increases by a cor-
responding factor (in addition to improvement in dynamic range available
through increased timestep when using CTE).
8.4 On Constrained Restoration of Tomographic PIV Data 195
8.4 ON CONSTRAINED RESTORATION OF TOMOGRAPHIC PIV DATA
The principle of Projection Onto Convex Sets was applied in order to post-
process tomographic PIV data. The POCS process involves successive appli-
cation of different ’Operators’ (constraint functions) in order to apply con-
straints to the vector field. A wide variety of constraints can be formu-
lated (providing a least-squares-minimising correction formulation can be
derived) and several operators were considered within this work:
Operator P1 Applies a spatial ’frequency’ band limit to the input field in the
Fourier domain.
Operator P2 Applies a prescribed value to elements of the input field.
Operator P3 Applies a divergence-free condition to an input field.
Operator P4 Applies wavelet-based scale filtering to structures in the input
field.
Most constraints here are simply filters trivially extended from 1D or 2D data
processing, however the P3 operator for divergence reduction is not trivial.
The operator is re-derived (based on the approach used in the literature for
2D fields). Dirichlet-type boundary conditions are used and shown to be ap-
propriate for this problem (which demands solution of the Poisson problem
on a cuboid domain in 3 dimensions).
Potential criticism that the solution is not unique for a 3 dimensional case
is addressed with a logic-based argument that, given the fully determined
3D3C flow field as an input (as opposed to a 3-dimension 2-component fluid
flow), the corrector field is uniquely determined given the caveat that the so-
lution of the equation system is formulated as a minimisation: that is, the
’corrector field’ does not merely exhibit the required property (of having di-
vergence consistent with the input), but that it is the closest field to the input
which exhibits that property (rather than the ’true’ solution, which remains
unknown).
A variety of different operators have been used in combination and sepa-
rately under the framework of the ’RELAX’ algorithm. The procedure has
been shown to effectively de-noise, to filter coherent structures according to
196 Conclusions
their scale and to removes divergence in a measured input field. Velocity vec-
tor components are altered within their experimental error in the majority of
cases.
Having demonstrated the ability to de-noise, limit scale and reduce diver-
gence; postprocessing algorithms which use derivatives and tracking tech-
niques to ascertain topological behaviour are used to demonstrate the effec-
tiveness of the POCS technique in isolation (chapter 6) and in combination
with other enhancements (chapter 7).
The POCS method is used to restore missing or inaccurate data (subject to
a scale limit on coherent structures contained within the flow) for velocity
fields which have up to 70% missing data, with an accuracy of the order
of the original measurement accuracy. Beyond a certain threshold, the non-
uniqueness of the reconstruction problem causes divergence in the POCS al-
gorithm, resulting in extremely high error in the restoration.
8.5 ON HIGH REYNOLDS NUMBER TURBULENT BOUNDARY LAYERS
Tomographic PIV has been used to perform statistical flow analyses for a
high Reynolds Number turbulent boundary layer (at δ+ = 1890, Reθ = 3660).
High speed time-resolved experiments have been performed for the same
flow. The region of interest was very close to the wall (within the overlap
and logarithmic regions of the flow). Grid points were spaced2 at 8.57+ (wall
units) or 0.52 mm, in the range 21.43 ≤ z+ ≤ 150.04. Table 7.1 gives further
details of the boundary layer parameters and measurement range used.
Both time-series and statistically converged results have a higher spatial res-
olution3 and are at a higher Reynolds Number δ+ than previously published
data visualising the near-wall region of the boundary layer. Results are avail-
able for use by the scientific community.
The following strictly preliminary conclusions have been drawn, subject to
wider examination and acceptance of the measurement technique used.
2Using 75% window overlap
3High spatial resolution is due predominantly to the physical characteristics of the bound-
ary layer but nevertheless benefits from the improved accuracy of the CTE technique, which
facilitates use of a smaller window size than otherwise achievable
8.6 Outlook 197
- Very Long Near-wall Streaks are present, consisting of trains of hairpin
vortices within the buffer/overlap layer. These long, small structures
are highly regular, exist over lengthscales longer than the tunnel (i.e.
are perpetuated by an autogeneration mechanism), do not appear to
meander significantly and co-exist with much larger structures within
the flow.
- Hairpins and Packets have been visualised at small scale close to the
wall. Typically, the width of the arches (in the span-wise direction)
is around 60 wall units. The height of the arches in the wall-normal
direction is almost always less than the width.
- Vortex Rings have been seen in the flow, entrained with axes oriented
in the streamwise direction. It is hypothesised that these are formed
due to vortex reconnection of hairpin legs. The rings are not visualised
with independent measurements so may arise due to clusters of false
vectors; further investigation at a higher spatial resolution is required
to conclude on their existence.
8.6 OUTLOOK
For Tomographic PIV, still further improvements are required before data be-
comes as powerful as computational results for investigations of flow struc-
ture. The development of multi-pulse systems such as the Correlation Track-
ing Enhancement shown here and similar methods involving specialist hard-
ware are presently being discussed [Adrian, 2011]. Such approaches, as indi-
cated within this work, can substantially improve the dynamic range of the
technique.
For 3D PIV in general, Tomographic PIV is very hardware-intensive, requires
careful calibration and extensive setups. Depth of field is the most limiting
factor (both the optical limitations and the bias error affect the ability to mea-
sure a deep volume). The bias error, whilst improved by application of the
CTE, is still present. Alternative systems are under development such as
high-speed scanning PIV systems which may yield a better depth of field
with less hardware.
198 Conclusions
With a decreased experimental noise available using multi-pulse systems, To-
mographic PIV would benefit from the application of derivatives correlation
[Scarano, 2004] in which the field derivatives were directly measured. Al-
though the method resulted in higher noise (exacerbating noise associated
with conventional TPIV techniques) the technique is far better able to cap-
ture gradients in the flow - a strong advantage when using TPIV for it’s core
strength - investigating structures, which are primarily identified using field
derivatives.
From section 7.3, there are three clear avenues for future research in turbulent
boundary layers:
- Very long (greater than the water tunnel length) near-wall streaks with
overarching vortex packets have been found. Predictive models based
on the Attached Eddy hypothesis of Perry & Chong [1982] do not yet
take into account these very long structures, assuming stochastically
distributed hairpins or smaller packets. Attempting to take into ac-
count the very long near wall motions (and the regions of quiescent
flow away from the near-wall streaks shown in this work) may yield
interesting results.
- The regularity of the near-wall region shown in figures 7.10 and 7.11, to-
gether with the streamwise extent of the structure in figures fig:quxlongtime
and fig:distortedtrain begs questions related to the interaction of in-
ner and outer flows. In particular, the subject of Very Large Scale Mo-
tions (VLSMs) in the outer layer is presently under discussion [Kim &
Adrian, 1999; Adrian et al., 2000; Marusic, 2001; Ganapathisubramani
et al., 2003]. The concept of autogeneration is central to the explana-
tion of VLSMs, but with increasingly chaotic flow at the larger scales
as Reynolds Number increases it is difficult to imagine a VLSM aris-
ing without meandering significantly or becoming intermittent. Hav-
ing shown a very long packet of hairpins generating turbulence near
the wall, it is worth investigating whether a VLSM at the outer scales
could be spatially correlated with or even caused by such a near-wall
structure.
- Vortex ring type structures were visualised in the flow. Both physical
and non-physical explanations for their appearance in the date were
8.6 Outlook 199
made, but the spatial resolution of the experiments was insufficient to
rule out the possibility that the rings are non-physical. Experiments
with a higher spatial resolution (easily achievable using Tomographic
PIV with greater zoom on the optics and a smaller depth of field) should
be undertaken, in the region z = 40+ to z = 120+, to confirm presence
of these structures.

References
Adrian R. (2011). Conceptual directions for the next generation of piv. In
PSFVIP-8: The 8th Pacific Symposium on Flow Visualization and Image Pro-
cessing, Moscow, Russia. 25, 197
Adrian R. & Liu Z. (2002). Observation of vortex packets in direct numerical
simulation of fully turbulent channel flow. J. Vis, 5(1):9–19, Also in: 6th
Asian Symposium on Visualization (ASV6), Pusan, South Korea, May 27-
31, 2001. xi, 14
Adrian R., Christensen K. & Liu Z. (2000). Analysis and interpretation of in-
stantaneous turbulent velocity fields. Exp. Fluids, 29(3):275–290. 14, 179,
198
Adrian R.J. (2007). Hairpin vortex organization in wall turbulence. Phys. Flu-
ids, 19(4). xi, 14
Atkinson C. & Soria J. (2009). An efficient simultaneous reconstruction tech-
nique for tomographic particle image velocimetry. Exp. Fluids, 47(4-5):553–
568. 12, 13, 69, 75, 115
Atkinson C., Coudert S., Foucaut J.M., Stanislas M. & Soria J. (2011). The
accuracy of tomographic particle image velocimetry for measurements of
a turbulent boundary layer. Exp. Fluids, 50(4, SI):1031–1056, Also in: 8th
International Symposium on Particle Image Velocimetry (PIV 09), Monash
Univ, Melbourne, Australia, Aug 25-28, 2009. 4, 9, 10, 58, 99, 113, 114, 115
Baur T. & Koengeter J. (2000). High-speed piv and the post-processing of
time-series results. In EUROMECH 411, Rouen. 23
Butterworth S. (1930). On the theory of filter amplifiers. Experimental Wireless
and the Wireless Engineer, 7:536–541. 21, 136
Chong M., Perry A. & Cantwell B. (1990). A general classification of three
dimensional flow fields. Phys. Fluids A., 2(5):765–777. 156
Coles D. (1956). The law of the wake in the turbulent boundary layer. J. Fluid
Mech., 1(02):191–226. xvii, 103
201
202 References
Cucitore R., Quadrio M. & Baron A. (1999). On the effectiveness and limita-
tions of local criteria for the identification of a vortex. Eur. J. Mech. B-Fluids,
18(2):261–282. xvii, 38
Daubechies I. (1992). Ten Lectures on Wavelets. SIAM: CBMS-NSF regional con-
ference series in applied mathematics; 61, 1st ed. 21, 22, 137, 139
Davidson P., Nickels T. & Krogstad P. (2006). The logarithmic structure func-
tion law in wall-layer turbulence. J. Fluid Mech., 550:51–60. xi, 40
Discetti S. & Astarita T. (2011). Spatial filtering improved tomographic piv. In
9th International Symposium of Particle Image Velocimetry, Kobe, Japan. 66
Elsinga G. & Marusic I. (2010). Evolution and lifetimes of flow topology in a
turbulent boundary layer. Phys. Fluids, 22(1). 19
Elsinga G., Scarano F., Wieneke B. & van Oudheusden B. (2006). Tomographic
particle image velocimetry. Exp. Fluids, 41(6):933–947. xi, xv, 2, 7, 8, 11, 17,
18, 27, 31, 33, 47, 58, 70, 83, 85, 89, 121
Elsinga G., Kuik D., van Oudheusden B. & Scarano F. (2007). Investigation
of the three-dimensional coherent structures in a turbulent boundary layer
with tomographic-piv. In 45th AIAA Aerospace Sciences Meeting and Exhibit,
Reno, Nevada. xv, 17, 27
Elsinga G., Adrian R., van Oudheusden B. & Scarano F. (2010). Three-
dimensional vortex organization in a high-reynolds-number supersonic
turbulent boundary layer. J. Fluid Mech., 644:35–60. xi, xv, 2, 3, 17, 18, 19,
20, 22, 28, 166, 167, 168, 169
Elsinga G.E., Westerweel J., Scarano F. & Novara M. (2011). On the velocity
of ghost particles and the bias errors in tomographic-piv. Exp. Fluids, 50(4,
Sp. Iss. SI):825–838, also in: 8th International Symposium on Particle Image
Velocimetry (PIV 09), Melbourne, Australia, Aug 25-28, 2009. 8, 63, 99, 104,
105, 106, 109, 115, 116
Farge M. (1992). Wavelet transforms and their applications to turbulence.
Annu. Rev. Fluid Mech., 24:395–457. 21, 22, 137
Farge M. & Schneider K. (2001). Coherent vortex simulation (cvs), a semi-
deterministic turbulence model using wavelets. Flow Turbul. Combust.,
66(4):393–426. 21, 137
References 203
Farge M., Schneider K. & Kevlahan N. (1999). Non-gaussianity and coherent
vortex simulation for two-dimensional turbulence using an adaptive or-
thogonal wavelet basis. Phys. Fluids, 11(8):2187–2201, also in: International
Conference on Turbulence, Los Alamos Nat’l Lab., Los Alamos, New Mex-
ico, Mar, 1998. 139
Farge M., Schneider K. & Abry P. (2002). Analysing and compressing turbu-
lent fields with wavelets. Tech. Rep. 20, ISSN 1626-8334, Institut Pierre Simon
Laplace. 139
Fore L.B., Tung A.B., Buchanan J.R. & Welch J.W. (2005). Nonlinear tempo-
ral filtering of time-resolved digital particle image velocimetry data. Exp.
Fluids, 39(1):22–31. 23, 58
Ganapathisubramani B., Longmire E. & Marusic I. (2003). Characteristics of
vortex packets in turbulent boundary layers. J. Fluid Mech., 478:35–46. 15,
179, 198
Ganapathisubramani B., Lakshminarasimhan K. & Clemens N. (2007). Deter-
mination of complete velocity gradient tensor by using cinematographic
stereoscopic piv in a turbulent jet. Exp. Fluids, 42:923–939, 10.1007/s00348-
007-0303-5. 21
Gerchberg R. (1974). Super-resolution through error reduction. Opt. Acta,
21:709–720. 23, 142
Griebel M. & Koster F. (2000). Adaptive wavelet solvers for the unsteady in-
compressible navier-stokes equations. In J. Malek, J. Necas & M. Rokyta
(eds.) Advances in Mathematical Fluid Mechanics, pp. 67–118, Springer-Verlag
Berlin, ISBN 3-540-67786-0, also in: 6th International School on Mathemat-
ical Theory in Fluid Mechanics, Paseky, Czech Republic, Sep 19-26, 1999.
21, 137, 139
Guezennec Y., Piomelli U. & Kim J. (1989). On the shape and dynamics of
wall structures in turbulent channel flow. Phys. Fluids A, 1(4):764–766. 14
Hain R., Ka¨hler C. & Michaelis D. (2008). Tomographic and time resolved
piv measurements on a finite cylinder mounted on a flat plate. Exp. Fluids,
45:715–724, 10.1007/s00348-008-0553-x. 12
204 References
Herman G.T. & Lent A. (1976). Iterative reconstruction algorithms. Comp. Biol.
Med., 6:273–294. 11, 191
Hinze J. (1975). Turbulence. McGraw-Hill, 2nd ed. 13, 14, 40, 177
Humble R., Scarano F. & van Oudheusden B. (2007). Instantaneous 3d flow
organization of a shock wave/turbulent boundary layer interaction using
tomo-piv. In Proc. 37th AIAA Fluid Dynamics Conference and Exhibit. Miami
(USA). 36
Jeong J. & Hussain F. (1995). On the identification of a vortex. J. Fluid Mech.,
285:69–94. 133, 147
Jeong J., Hussain F., Schoppa W. & Kim J. (1997). Coherent structures near the
wall in a turbulent channel flow. J. Fluid Mech., 332:185–214. 14
Kim K. & Adrian R. (1999). Very large-scale motion in the outer layer. Phys.
Fluids, 11(2):417–422. 15, 198
Ku¨hn M., Ehrenfried K., Bosbach J. & Wagner C. (2011). Large-scale tomo-
graphic particle image velocimetry using helium-filled soap bubbles. Exp.
Fluids, 50(4, SI):929–948, also in: 8th International Symposium on Particle
Image Velocimetry (PIV 09), Monash Univ, Melbourne, Australia, Aug 25-
28, 2009. xv, 17, 29
Lawson C. & Hanson R. (1995). Solving Least Squares Problems. SIAM: Classics
in Applied Mathematics, 1st ed. 77
Lourenco L. & Krothapalli A. (1995). On the accuracy of velocity and vorticity
measurements with piv. Exp. Fluids, 18(6):421–428. 8, 20, 54, 59
Mallat S. (1999). A Wavelet Tour of Signal Processing. Academic Press, 2nd ed.
21, 139
Mann J., Ott S. & Andersen J.S. (1999). Experimental study of relative, turbu-
lent diffusion. Tech. Rep. R-1036(EN), Riso National Laboratory, Roskilde,
Denmark. 49, 53
Marusic I. (2001). On the role of large-scale structures in wall turbulence.
Phys. Fluids, 13(3):735–743. 15, 198
References 205
Marusic I. & Perry A. (1995). A wall-wake model for the turbulence structure
of boundary-layers 2. further experimental support. J. Fluid Mech, 298:389–
407. 15
Marusic I., Mathis R. & Hutchins N. (2010a). Predictive model for wall-
bounded turbulent flow. Science, 329(5988):193–196. 15
Marusic I., McKeon B.J., Monkewitz P.A., Nagib H.M., Smits A.J. & Sreeni-
vasan K.R. (2010b). Wall-bounded turbulent flows at high reynolds num-
bers: Recent advances and key issues. Phys. Fluids, 22(6). 1, 16, 182
Mathis R., Hutchins N. & Marusic I. (2011). A predictive inner-outer model
for streamwise turbulence statistics in wall-bounded flows. J. Fluid Mech.,
681:537–566. 15, 182
Michaelis D. (2011). personal communication. 66
Michaelis D., Poelma C., Scarano F., Westerweel J. & Wieneke B. (2006). A
3d-time resolved cylinder wake survey by tomographic piv. In Proc. 12th
Int. Symposium on Flow Visualisation, ISFV12, Go¨ttingen, Germany. 36
Mishra D., Muralidhar K. & Munshi P. (1999). A robust mart algorithm for to-
mographic applications. Numer. Heat Transf. B-Fundam., 35(4):485–506. 11,
69
Nelder J.A. & Mead R. (1965). A simplex method for function minimization.
The Computer Journal, 7(4):308–313. 93
Nickels T., Marusic I., Hafez S. & Chong M. (2005). Evidence of the k(1)(-1)
law in a high-reynolds-number turbulent boundary layer. Phys. Rev. Lett.,
95(7). 40
Nobach H. (2004). Accuracy of sub-pixel interpolation in PIV and PTV image
processing. Tech. Rep. 001/2004, Fachbereich Maschinenbau, Fachgebiet
Stro¨mungslehre und Aerodynamik, Technische Universita¨t Darmstadt. 8,
59
Nobach H. & Bodenschatz E. (2009). Limitations of accuracy in piv due to
individual variations of particle image intensities. Exp. Fluids, 47:27–38. 10,
170
206 References
Novara M., Batenburg K.J. & Scarano F. (2010). Motion tracking-enhanced
mart for tomographic piv. Meas. Sci. Technol., 21(3). 9, 10, 64, 65, 123
Offen G. & Kline S. (1974). Combined dye-streak and hydrogen-bubble visual
observations of a turbulent boundary layer. J. Fluid Mech., 62(Jan):223–&.
179, 180, 186
Offen G. & Kline S. (1975). Proposed model of bursting process in turbulent
boundary-layers. J. Fluid Mech., 70(Jul):209–228. 13, 14, 177
Papoulis A. (1975). A new algorithm in spectral analysis and band-limited
extrapolation. IEEE Trans. Circuits Syst., CAS-22:735–742. 23, 142
Perry A. & Chong M. (1982). On the mechanism of wall turbulence. J. Fluid
Mech., 119(Jun):173–217. xi, xvii, 15, 16, 198
Petra S., Schro¨der A. & Schno¨rr C. (2009). 3D Tomography from few projec-
tions in experimental fluid dynamics. In W. Nitsche & C. Dobriloff (eds.)
Imaging Measurement Methods for Flow Analysis, vol. 106 of Notes on Numer-
ical Fluid Mechanics and Multidisciplinary Design, pp. 63–72, Springer Berlin
/ Heidelberg, ISBN 978-3-642-01105-4. 11
Piirto M., Eloranta H., Saarenrinne P. & Karvinen R. (2005). A comparative
study of five different piv interrogation algorithms. Exp. Fluids, 39(3):571–
588. 8, 59, 106
Raffel M., Willert C., Wereley S. & Kompenhans J. (2007). Particle Image Ve-
locimetry A Practical Guide (Second Edition). Springer, 2nd ed. 10, 36, 163
Robinson S. (1990). The kinematics of turbulent boundary layer structure. Ph.D.
thesis, National Aeronautics and Space Administration. 14
Robinson S. (1991). Coherent motions in the turbulent boundary-layer. Annu.
Rev. Fluid Mech., 23:601–639. 1, 15
Roesgen T. (2003). Optimal subpixel interpolation in particle image velocime-
try. Exp. FLuids, 35(3):252–256. 8, 59
Ruijters D. & The´venaz P. (2010). Gpu prefilter for accurate cubic b-spline
interpolation. The Computer Journal. 54
References 207
Scarano F. (2002). Iterative image deformation methods in piv. Meas. Sci. Tech-
nol., 13(1):1–19. 8, 31, 54, 58
Scarano F. (2004). A super-resolution particle image velocimetry interroga-
tion approach by means of velocity second derivatives correlation. Meas.
Sci. Technol., 15(2):475–486. 59, 60, 123, 172, 198
Scarano F. & Riethmuller M. (2000). Advances in iterative multigrid piv im-
age processing. Exp. Fluids, 29(S):S51–S60, Also in: 3rd International Work-
shop on Particle Image Velocimetry, Santa Barbara, CA, Sep 16-18, 1999. 8,
58, 63, 191
Scarano F., Elsinga G., Bocci E. & van Oudheusden B. (2006). Three-
dimensional turbulent cylinder wake structure investigation with tomo-
piv. In Proc. 13th Intl. Symp. on Applications of Laser Techniques to Fluid Me-
chanics, Lisbon (Portugal). 36
Schlatter P. & Orlu R. (2010). Assessment of direct numerical simulation data
of turbulent boundary layers. J. Fluid Mech., 659:116–126. 15, 182
Schlatter P., Orlu R., Li Q., Brethouwer G., Fransson J., Johansson A., Al-
fredsson P. & Henningson D. (2009). Turbulent boundary layers up to
re(theta)=2500 studied through simulation and experiment. Phys. Fluids,
21(5). 15
Schlichting H. & Gersten K. (2003). Boundary Layer Theory. Springer, 8th re-
vised and enlarged ed. xi, 38, 39
Schram C., Rambaud P. & Riethmuller M. (2004). Wavelet based eddy struc-
ture eduction from a backward facing step flow investigated using particle
image velocimetry. Exp. Fluids, 36(2):233–245. 22, 139
Schrijer F.F.J. & Scarano F. (2008). Effect of predictor-corrector filtering on the
stability and spatial resolution of iterative piv interrogation. Exp. Fluids,
45(5):927–941. 99
Schro¨der A., Geisler R., Elsinga G., Scarano F. & Dierksheide U. (2006). Inves-
tigation of a turbulent spot using time-resolved tomographic piv. In Proc.
13th Intl. Symp. on Applications of Laser Techniques to Fluid Mechanics, Lisbon
(Portugal). 21, 36
208 References
Schro¨der A., Geisler R., Elsinga G., Scarano F. & Dierksheide U. (2008). Inves-
tigation of a turbulent spot and a tripped turbulent boundary layer flow
using time-resolved tomographic piv. Exp. Fluids, 44(2):305–316. xv, 17, 27
Schro¨der A., Geisler R., Staack K., Elsinga G.E., Scarano F., Wieneke B., Hen-
ning A., Poelma C. & Westerweel J. (2011). Eulerian and lagrangian views
of a turbulent boundary layer flow using time-resolved tomographic piv.
Exp. Fluids, 50(4, Sp. Iss. SI):1071–1091, also in: 8th International Sympo-
sium on Particle Image Velocimetry (PIV 09), Melbourne, Australia, Aug
25-28, 2009. xi, xv, 2, 3, 17, 18, 28, 63, 146, 150, 167, 168, 169
Sezan M.I. & Stark H. (1982). Image restoration by the method of convex
projections: part 2 applications and numerical results. IEEE Trans. Med.
Imaging, 1(2):95–101. 24, 142, 145
Simard P. & Mailloux G. (1988). A projection operator for the restoration
of divergence-free vector fields. IEEE Trans. Pattern Anal. Mach. Intell.,
10(2):248–256. 24, 128, 129, 130, 131, 132, 134
Simard P. & Mailloux G. (1990). Vector field restoration by the method of con-
vex projections. Computer Vision, Graphics and Image Processing, 52(3):360–
385. 24, 25, 128
Stark H. (ed.) (1987). Image Recovery: Theory and Application, pp. 29–77. Aca-
demic Press, New York, 1st ed. 132
Stearns S. & Hush D. (1990). Digital Signal Analysis. Prentice-Hall, 2nd ed. 54
Suter D. (1994). Divergence-free wavelets made easy. Tech. Rep. MECSE 1994-2,
Dept. of Electrical and Computer Systems Engineering, Monash Univer-
sity. 132
Tee B. & Nickels T. (2008). Experimental investigation of zero pressure-
gradient turbulent boundary layers using particle image velocimetry. APS
Division of Fluid Dynamics Meeting Abstracts, pp. G9+. xv, 36, 37, 41
Theodorsen T. (1952). Mechanism of turbulence. In 2nd Midwestern Conf. on
Fluid Mech., pp. 1–18, Ohio State University. xi, 1, 13, 177
Townsend A. (1976). The Structure of Turbulent Shear Flows. Cambridge Uni-
versity Press, 2nd ed. 15
References 209
van Driest E. (1956). On turbulent flow near a wall. J. Aero. Sci., 23(11):1007–
1011,1036. 114
Ve´tel J., Garon A. & Pelletier D. (2011). Denoising methods for time-resolved
piv measurements. Exp. Fluids, 51(4):813–916. 23
Westerweel J. & Scarano F. (2005). Universal outlier detection for piv data.
Exp. Fluids, 39(6):1096–1100. 8, 90, 91, 156, 158
Wieneke B. (2005). Stereo-piv using self-calibration on particle images. Exp.
Fluids, 39(2):267–280. 8
Wieneke B. (2008). Volume self-calibration for 3d particle image velocimetry.
Exp. Fluids, 45(4):549–556, also in: 7th International Symposium on Particle
Image Velocimetry, Univ Roma, Rome, Italy, Sep 11-14, 2007. 8, 49, 52, 53
Wieneke B. & Taylor S. (2006). Fat-sheet piv with computation of full 3d-
strain tensor using tomographic reconstruction. In Proc. 13th Intl. Symp. on
Applications of Laser Techniques to Fluid Mechanics, Lisbon (Portugal). 10, 38
Worth N.A. (2010). Tomographic-PIV Measurement of Coherent Dissipation Scale
Structures. Ph.D. thesis, University Of Cambridge. 58, 155, 156
Worth N.A. & Nickels T.B. (2007). A computational study of tomographic re-
construction accuracy and the effects of particle blocking. ASME Conference
Proceedings, 2007(42886):691–700. 4, 10
Worth N.A. & Nickels T.B. (2008). Acceleration of tomo-piv by estimating the
initial volume intensity distribution. Exp. Fluids, 45(5):847–856. xi, 12, 38,
42, 53, 69, 74, 191
Worth N.A., Nickels T.B. & Swaminathan N. (2010). A tomographic piv res-
olution study based on homogeneous isotropic turbulence dns data. Exp.
Fluids, 49(3):637–656. xv, 4, 17, 21, 22, 28, 32, 38, 42, 58, 115
Youla D. (1978). Generalised image restoration by the method of alternating
projections. IEEE Trans. Circuits Syst., CAS-25. 23
Youla D.C. & Webb H. (1982). Image restoration by the method of convex
projections: Part 1 Theory. IEEE Trans. Med. Imaging, 1(2):81–94. 24, 127,
128, 129
210 References
Zhou J., Adrian R. & Balachandar S. (1996). Autogeneration of near-wall vor-
tical structures in channel flow. Phys. Fluids, 8(1):288–290. 14