A MACHINE LEARNING APPROACH TO FOLLOWING
TRACE PARTICLES IN TURBULENT FLOW
By
RANDEL ALLAN CROWE
A DISSERTATION PRESENTED TO THE GRADUATE COUNCIL OF
THE UNIVERSITY OF FLORIDA
IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE
DEGREE OF DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
1977
UNIVERSITY OF FLORIDA
3 1262 08552 3156
To Patricia, A Most Special Person
ACKNOWLEDGEMENTS
Deserving much of the credit for the completion of this work is,
of course, my advisor, Professor Gale E. Nevill. Often in dark times,
he gave me needed encouragement and guidance and has my most sincere
thanks. Others on my committee should be noted for their patience as
well as valuable suggestions. Dr. Kurzweg reviewed the fluid mechanics
aspects and Dr. Boykin the dynamics and modeling. Little could have
been accomplished in my graduate work without the very early support
and friendship of Dr. Hemp. Enduring the final push, Dr. Shaffer
has my deepest appreciation. Yet many others should also be included.
Surely, Professor E. Rune Lindgren, and his associates, Drs.
Elkins and Jackman whose project provided this problem, together with
Professor Roland C. Anderson who aided in the optical analysis, should
be recognized for providing unending assistance and counseling. Quietly
in the background are the members of the "MNDC,' a local social group,
too numerous to mention individually, but all very supportive and there
when needed. Unusual was the serendipity and productivity of my con
versations with Mickey Rogers. All words fail me when I try to express
my appreciation and admiration of my wife, Pat, who gave me incessant
drive and inspiration. To my parents as well as those mentioned and
unmentioned, I would like to express my everlasting gratitude and
indebtedness.
TABLE OF CONTENTS
ACKNOWLEDGEMENTS.... .............................................
LIST OF FIGURES...................................................
LIST OF TABLES....................................................
ABSTRACT..... .................................. ..................
CHAPTER
I INTRODUCTION.............................................
Scope of Work................................ ...... ....
Pertinent Background for the Particle Following
Machine.. ..............................................
II THE TRACE PARTICLE TECHNIQUE AND DATA ACQUISITION SYSTEM..
The Use of Trace Particles for Flow Visualization......
Earlier Particle Following Approaches .................
Qualitative Observations of Particle Motion............
III A PARTICLE FOLLOWING MACHINE .............................
Overview ..............................................
PFM: Input and Output.................................
Image Path Attributes.................................
Following Image Paths: The Decision Process...........
The Image Feature Vector...............................
Modeling the Residual.................................
Control of the PFM ....................................
Measuring Performance.................................
Summary...............................................
Page
iii
vi
viii
ix
CHAPTER Page
IV PFM: INITIALIZATION AND IMPLEMENTATION................... 58
Initialization........................................ 58
The Data ............................................... 66
The Initialization Program ............................. 69
The Particle Following Program......................... 69
Program LEARN: Determination of the Joint
Probabilities......................................... 74
V RESULTS OF PERFORMANCE TESTS.............................. 76
Initialization........................................ 76
Particle Following.................................... 85
Performance of the PFP................................. 90
VI DISCUSSION AND CONCLUSIONS............................... 95
APPENDIX A ANALYSIS OF THE DATA ACQUISITION SYSTEM ............... 98
The Data Acquisition System........................ 98
Approximate Geometric Ray Trace.................... 104
Qualitative Observations of Images................. 107
Geometric Ray Trace................................ 108
Spot Diagrams.................................... 128
Digitization of Film Data.......................... 133
Path Characteristics Using Orthogonal Projections.. 136
APPENDIX B APPROXIMATE GEOMETRIC RAY TRACE PRISM EQUATIONS........ 144
APPENDIX C TABLET DIGITIZER ERROR................................ 155
APPENDIX D FORTRAN PROGRAM LISTINGS.............................. 164
BIBLIOGRAPHY.............................. ......................... 183
BIOGRAPHICAL SKETCH................................... ............ 187
LIST OF FIGURES
FIGURE Page
11 DATA ACQUISITION SYSTEM BLOCK DIAGRAM ..................... 3
21 THE DATA ACQUISITION SYSTEM ............................... 11
31 THE PARTICLE FOLLOWING MACHINE ............................ 22
32 DIGITIZED IMAGE PLANE COORDINATES ......................... 22
33 IMAGE CANDIDATES AND ERRORS............................... 35
34 LINK CLASSES, DECISION STATES AND HYPOTHESES.............. 37
35 EXAMPLE OF JOINT TRANSITION PROBABILITY................... 49
36 BLOCK DIAGRAM OF THE PARTICLE FOLLOWING MACHINE............ 57
41 FRAMEBYFRAME DATA....................................... 67
42 INVERTIBLE PATH DATA...................................... 67
43 SIMPLIFIED FLOWCHART OF "INIT"............................ 70
44 SIMPLIFIED FLOWCHART OF "PFP"............................. 71
45 SIMPLIFIED FLOWCHART OF "LEARN" ........................... 75
51 TANGENTIAL AND NORMAL ACCELERATION COSTS.................. 77
52 COST FUNCTION CONTOURS................................... 79
53 INITIAL IMAGE PATH SEGMENTS: TOP HALF OF IMAGE PLANE..... 80
54 INITIAL IMAGE PATH SEGMENTS: BOTTOM HALF OF IMAGE PLANE.. 81
55 RESIDUAL JOINT TRANSITION PROBABILITY, s = 1.0............ 86
56 RESIDUAL JOINT TRANSITION PROBABILITY, s = 1.2............ 87
57 RESIDUAL JOINT TRANSITION PROBABILITY, s = 1.5............ 88
58 RESIDUAL JOINT TRANSITION PROBABILITY, s = 2.0............ 89
Ai DEFINITION OF COORDINATE AXES............................. 99
A2 APPROXIMATE MERIDIONAL RAY TRACE .......................... 106
FIGURE Page
A3 GRID OF OBJECT POINTS.................................... 109
A4 MERIDIONAL RAY TRACE: LARGE PRISM ........................ 110
A5 OFFAXIS RAY TRACE: LARGE PRISM.......................... 110
A6 RAY TRACE TOP VIEW....................................... 111
A7 PRINCIPAL RAYS FROM OBJECT GRID INTERSECTING IMAGE PLANE.. 111
A8 MERIDIONAL RAY TRACE: SMALL PRISM ........................ 113
A9 LARGE PRISM RAY LENGTHS FOR FULL LIGHTING.................. 115
A10 UNNECESSARY CONFUSION FOR LARGE PRISM WITH FULL LIGHTING.. 115
A11 SMALL PRISM RAY LENGTHS FOR FULL LIGHTING................. 116
A12 UNNECESSARY CONFUSION FOR SMALL PRISM WITH FULL LIGHTING.. 116
A13 LIGHTING SCHEMES......................................... 117
A14 UNNECESSARY CONFUSION FOR HALF AND CIRCLE ILLUMINATION.... 118
A15 VOLUME SENSITIVITIES (cm3/mm2)........................... 123
A16 PROBABILITIES FOR CONFUSION AND OVERLAP.................... 124
A17 HEXAPOLAR ARRAY ........................................... 130
A18 OBJECT POINT LOCATIONS.... ............................... 130
A19 SENSITIVITY OF FOCUSING POINT A............................ 131
A20 SPOT DIAGRAMS FOR POINT A AND A', BEST FOCUS.............. 132
A21 SPOT DIAGRAMS FOR POINT C AND C', BEST FOCUS SET FOR
POINT A.................................................... 134
A22 SPOT DIAGRAMS FOR POINTS B AND B', BEST FOCUS ON A........ 135
A23 ORTHOGONAL STEREOSCOPIC PROJECTION OF A PATH............... 137
B1 RAY TRACE GEOMETRY.................................... .... 146
C1 DISTRIBUTION OF PEN LOCATION ABOUT DESIRED IMAGE
LOCATION.................................................. 156
C2 TWO DIMENSIONAL REGIONS................................... 156
LIST OF TABLES
TABLE Page
51 KEY TO FIGURES 5.3 AND 5.4................................ 82
52 INITIALIZATION PERFORMANCE ON 207 FRAMEONE IMAGES........ 84
53 PFP PERFORMANCE, s = 2.0.................................. 92
54 ISOLATION OF DECISIONS................................... 94
Ala APPROXIMATE PARTICLE MOVEMENT (MM) IN SHUTTER TIME......... 103
Alb APPROXIMATION USING MAX VELOCITY= UBARO.791............... 103
A2 PROBABILITY OF TWO PARTICLES OCCURRING IN VI .............. 119
A3 AVERAGE CONFUSION AND OVERLAP PROBABILITIES FOR SMALL
PRISM..................................................... 126
A4 AVERAGE CONFUSION AND OVERLAP PROBABILITIES FOR LARGE
PRISM..................................................... 127
C1 TWO DIMENSIONAL TABLET ERROR.............................. 161
C2 PROBABILITY OF BIT ERROR FOR GIVEN E....................... 162
viii
Abstract of Dissertation Presented to the Graduate Council
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy
A MACHINE LEARNING APPROACH TO FOLLOWING
TRACE PARTICLES IN TURBULENT FLOW
By
Randel Allan Crowe
August 1977
Chairman: Gale E. Nevill, Jr.
Major Department: Engineering Sciences
Machine learning is used in a statistical decision theoretic scheme
to follow trace particles in turbulent flow. Data used to demonstrate
performance came from previous experiments studying turbulent flow in a
pipe (Reynolds number less than 6500), and consisted of digitized stereo
scopic trace particle image locations. Accuracy of the data acquisition
technique is examined with emphasis on optical errors. Particle
location in the pipe is determined from its two stereo images. How
ever, high trace particle count prevents immediate pairing of particle's
stereo images. Therefore, identification of image paths (sequentially
sampled particle locations) is carried out separately in each view.
The particle following scheme is initialized by short image path segments
found using search procedures based on finding fiveimage sequences with
minimum velocity variation. Image dynamics are modeled using a constant
acceleration assumption. An image's dynamic state is estimated using
a Kalman filter with finite fading memory. The next image belonging
to the path is sought in the neighborhood of a prediction. A minimum
error rate Bayes decision rule is applied to choose the next image
from the candidates in the neighborhood. The feature vector is
composed of the two most recent filter residuals. Probabilities of
the residual transitions are learned by analyzing previously completed
and verified sequences. It is found that this scheme gives better
performance than earlier attempts, which utilized a constant velocity
assumption, including improved handling of image overlap, confusion,
and measurement noise. It is concluded that a semiautomated system
is necessary where the particle following procedure's output is.
verified by an operator who can change or correct decisions made.
CHAPTER I
INTRODUCTION
Scope of Work
Turbulence has evaded detailed and complete understanding by scores
of scientists over many generations. Recent attempts at its description
suffer from the difficulty of obtaining good experimental verification
which is needed in any modeling effort. In particular, two separate
experiments that have been reported utilized trace particles to facilitate
the detailed quantitative description of the structure of turbulence in
pipes. Johnson (1974) and Jackman (1976) used neutrally buoyant particles
in water and a novel prismcinematographic data recording system to
study transition and low Reynolds number flows (3500 Re 6500).
Breton (1975) used small glass spheres in trichloroethelyne (TCE)
and a unique "distributed camera" system* to study a flow with Reynolds
number of the order of 100,000. These workers were interested in a
Lagrangian description of the turbulence. That is, a turbulent model
written in terms of the individual traces of fluid elements. These
experiments generated stereoscopic film records of temporally sampled
particle paths, and both suffered from the inability to automatically
reduce these data to the individual path sequences.
*This consisted of two sets of twenty lenses placed along a test section
of the pipe.
2
Quantitative flow visualization requires taking data by an optical
system that records two different views (stereo views) of the particles
in the flow. Each view is a projection of the threedimensional par
ticle paths onto a twodimensional film plane. By measuring the image
locations in the film plane, the two views of a particle can be used to
determine its threedimensional position. This is what is termed photo
grammetry. Breton (1975) used a sufficiently small number of particles
to allow straightforward determination of stereo pairs. Johnson (1975)
and Jackman (1976), however, had high particle counts making the immedi
ate determination of stereo image pairs impossible. Instead, they deter
mined the twodimensional image paths ( a process called particle following)
of all images and then compared the paths to determine which were stereo
pairs (conjugate paths). A high particle count yields a better and more
reliable visualization. But the use of more particles increases the
burden on the data reduction procedures. Breton (1975) and Jackman (1976)
used automated procedures as much as possible but were not successful at
attaining full automation. Automatic data analysis is difficult and will
not become a reality without a thorough understanding of the data
acquisition system, its errors, the particle path characteristics and
the concepts of particle following. The work presented here considers
Jackman's experimental arrangement and the errors in his data acquisition
system. Considered in detail is the particle following problem.
Data reduction consists of several identifiable subproblems which,
together, form a linear data acquisition and reduction system (See
Figure 11). Once the experimental system has been defined, the first
problem is to acquire quality stereo images of the trace particles.
DATA ACQUISITION SYSTEM BLOCK DIAGRAM
FIGURE 11
This is an acute optics problem in the case of Johnson's (1975) prism
system. The present work includes an in depth analysis of the optical
properties of the pipe and prism system. Appendix A presents details of
the analysis. Second, the images must be recorded on film and their lo
cations digitized for use by the analysis programs. Jackman's (1976)
use of a graphic tablet system for data entry is treated in Appendix C.
The primary concern of this work is the third step; converting the
digitized film data into image paths. The approach taken here involves
three procedures: prediction, decision making, and learning. The
combination of these processes to perform the image following task is
termed the Particle Following Machine (PFM). The PFM takes the digitized
image locations and determines which images represent the same particle
from frame to frame (sample to sample). It outputs image paths which can
be compared to find the conjugate paths. Then, a transformation is used
to determine the threedimensional particle paths. These last two aspects
of the data analysis are not considered in the current work but they are
discussed along with the actual evaluation of the results in terms of
fluid turbulence theory by Johnson (1974), Breton (1975), and Jackman (1976).
The principal contribution of this work is the development of a
theoretical basis for the PFM. This approach to following images of
trace particles makes use of a combination of estimation, decision and
pattern recognition theory. The PFM is presented in a general format
and could be applied to other problems where a group of identical objects
are being tracked and the motion of an object is fairly consistent over
time. Performance tests of a specific implementation of the PFM were
made using a computer code written in FORTRAN. Some modifications would
be necessary to make this implementation an economic production program.
A secondary contribution, the optical analysis of the prism system
and digitization error of a data tablet, is important to the continued
vitality of the experimental procedure. The results of the analysis
performed gives insight into the data acquisition system and the charac
ter of the input to the PFM.
Pertinent Background for the Particle Following Machine
The particle following task is simply stated as determining which
images in the digitized film data belong to a given particle. An image
path is found by determining which image in each frame of data belongs
to the specified particle. Therefore, for each frame, a decision must
be made as to which image belongs to the particle being followed. This
leads to the task being considered a pattern recognition problem. The
image paths are the desired patterns to be identified. The images of
a frame of data must be classified as belonging to the given particle
or not. To make these classifications, a feature vector is necessary.
Here, the entire feature vector is used in the decision process. The
class probabilities for a given feature vector are learned by using
previously identified image paths. The use of the entire feature vector
and a set of given classifications is referred to as machine learning,
i.e. parallel pattern recognition with supervised learning (Hunt, 1975).
Synonymously, Hunt uses the term adaptive pattern recognition while
Duda and Hart (1972) refer to it as simply adaptive learning. Further,
the particle following task is a compound decision problem since more
than one decision is being made at each step.
The PFM uses residuals of an adaptive tracking filter as components
of the feature vector. Finite memory filters can be used as tracking
filters as shown by Jazwinski (1970). An alternate technique, used in
the present work is the finite fading memory filter of Tarn and
Zaborszky (1970). These filters are successful at tracking because
they eliminate possible divergence of the Kalman filter. The divergence
of the Kalman filter is considered by Price (1968). Adaptive radar
tracking systems following manned maneuvering targets use similar track
ing filters. Singer (1970) uses a Kalman filter to track targets with
known maneuvering capabilities. Singer and Behnke (1970) present a
comparison of five tracking filters concluding that the Kalman filter
requires the most computation but provides measures of tracking error
statistics. McAulay and Denlinger (1973) present a tracking filter
that tracks by changing its dynamic model of the target when a maneuver
is detected. This system requires knowing a target's allowed maneuvers
beforehand. Adaptive Kalman filters are treated in Gelb (1974). Back
ground on learning and adaptive systems can be found in Tsypkin (1971
and 1973). The general area of adaptive systems is scanned by Lainiotis
(1976).
Since the feature vector is constructed from filter residuals, it
is a time series. It is assumed here to be a Markov process with the
next state only dependent on the immediately previous state. This
use of the history of the feature vector makes the PFM an example of
making decisions in the context of the current situation (i.e. the way
the system got to its current state is important). Modeling the fea
ture vector is then partially involved with time series analysis (Box
and Jenkins, 1970, and Robinson, 1967) and partially estimation theory
(Jazwinski, 1970).
The most interesting aspect of the PFM is its need to first
decide which image in a frame belongs to the particle being followed.
Once this decision has been made, normal estimation theory applies, i.e.
the last estimate can be updated with information from the next measure
ment. The particle following task is therefore a combination of appli
cations of decision theory and estimation theory. No directly applicable
literature was found that discussed the underlying theory of such a
system. This unique problem, then, requires as basic a development
as possible to provide a foundation on which to build a viable production
system. This is evident from the lack of success of other attempts
at automatic particle following. Breton (1975) avoided the problem
entirely by resorting to manual entry of conjugate images and having a
very low particle/image density. Jackman (1976) tried to use a deter
ministic approach that assumed constant velocity particles and no
measurement noise. This effectively reduced the amount of decision
making, but when a choice had to be made, it was arbitrary. Its per
formance was fairly poor and the images were eventually followed by
manual effort. Since the desired image paths are known to be subsets
of the data (collections of images), one approach might be to test all
image combinations with a heuristic cost function. It would then
be assumed that paths with minimum cost would be the desired image
paths. With the image density that Jackman used, the number of
possible combinations is very large; at ten times the density, the
problem would be enormous. Nilsson (1971), and Newell and Simon (1972)
discuss heuristic problem solving. For this approach, the search of
possible combinations is heuristically guided by rules that "seem"
appropriate. These rules are adhoc but possibly simpler than the
decision theoretic approach taken here. It is difficult to base such
rules on theory and, typically, they provide neither insight to their
accuracy nor guidance for their improvement. The heuristic approach
is useful in solving the initialization problem. (See Chapter IV.)
Initially, there is no information available as to the dynamic state
of an image. By severely limiting the search and accepting some'error,
a heuristic approach is taken to identify short image paths. These
initial path segments are then sufficient to start the particle following
procedure. Application of this initialization procedure to the particle
following task would ignore much information that is readily available
such as image state and measurement noise.
After a brief description of the data acquisition system and a
discussion of some aspects of the analysis of the optical system in
Chapter II, the PFM is described in Chapter III. Chapter IV presents
the details of initialization and programming a specific implementation
of the particle following machine while results and conclusions are
presented in Chapters V and VI respectively.
CHAPTER II
THE TRACE PARTICLE TECHNIQUE AND DATA ACQUISITION SYSTEM
Presented in this chapter is the background of the trace particle
technique and data acquisition system used by Johnson (1974) and
Jackman (1976). In addition, results of qualitative and quantitative
analyses of the system are discussed. These results provide a foundation
for the development of the particle following machine in Chapter III.
Details of the analysis are given in Appendix A.
The Use of Trace Particles for Flow Visualization
There exist numerous techniques for flow visualization. Their
purpose is to reveal the motion of a fluid element over time in com
plicated flows. Additives such as dye, wax and rosin spheres, dust,
hydrogen bubbles, aluminum powder and mica flakes have all been tried.
Ideally, the trace material used follows the fluid motion precisely
without altering the flow structure by its presence. Examples of flow
visualization may be found in Reynolds (1883), Prandtl (1904), Fage and
Townend (1932), Prandtl and Tietjens (1934), Lindgren (19541969), Coles
(1965), Kline et al (1967), Corino and Brodkey (1969), Kim et al.(1971),
Johnson (1974), Breton (1975), Jackman (1976), Johnson et al.(1976), and
Elkins et al.(1977).
The use of pliolite trace particles is due to Nychas, Hershey, and
Brodkey (1973). These particles are desirable because they are opaque
(reflect light well), of selectable size, and almost neutrally buoyant.
Johnson (1974) gives a lengthy argument as to their abilities to follow
the fluid motion faithfully.
Johnson (1974) introduced the use of a prism to generate the stereo
images required for quantitative measurements. This technique requires
only one camera, simplifying somewhat the required equipment compared
to the special distributed camera of Breton (1975).
Quantitative descriptions are typically not obtainable from flow
visualization techniques due to the complicated equipment and analysis
required. Breton (1975) traced only a few particles compared to Johnson
(1974). Jackman (1976) was able to attain one or two particles per cubic
centimeter which approaches a reasonable particle count (henceforth re
ferred to as density). A goal of the current work of Lindgren and his
associates is to increase the density by an order of magnitude (Lindgren,
1977). As the density is increased, the work required in data reduction
increases tremendously and fully or semiautomatic procedures are vital.
Johnson performed all particle following by hand. Jackman was able to
develop semiautomated data entry, and path matching (finding conjugate
paths) procedures but resorted to following the images manually. Breton
reported failure of his automatic system to follow image:path pairs and
essentially used a machineassisted manual entry procedure.
The present work is a study of the particle following problem.
In order to understand the basis for the development of the particle
follower, the experimental system and data acquisition equipment need
to be understood. The remainder of Chapter II describes the apparatus
and discusses the results from analysis of the optical system.
Earlier Particle Following Approaches
Of primary concern to the current work is the data acquisition
system used by Jackman. Mounted at the observation station of his
experimental apparatus is a 90 prism as shown in Figure 21. Light
LLi
CJ LLJ
\U <
LJL LL
0
LU
o)
00
C
AJ
I
Lii
I
LU
Lid
1.
Lii Lii
a
LUJ
0)
Ct
0
cr 
L
LL
3:
projected down the pipe reflects off the trace particles making them
highly visible. An observer looking into the pipe through the prism
sees two views of the particles in the pipe, one view in each of the
two prism faces. A particle is invertible if it has an image in both
views. Its location inside the pipe can then be calculated. The
region for which a particle's images are invertible depends on the
prism parameters, the pipe parameters, and the observer's position.
To enhance the capability of the system, two prisms were used. The
smaller is appropriate for studies of motion close to the wall; the
larger to make observations deep into the pipe. A Bell and Howell
35mm movie camera was used to record the particle images. Its film
rate was adjusted to suit the flow rate (faster flowsmore particle
motionfaster film rate required). Two small light sources mounted
on the prism were used as reference points for frame alignment during
subsequent data conditioning and analysis. The display of a Hewlett
Packard digital counter was also recorded on each frame to allow
accurate determination of the sample rate.
Jackman used a graphic tablet digitizer to encode the image lo
cations. The graphic tablet system consisted of a Scriptographics
Tablet Digitizer connected online to a Tektronix graphic terminal
interfaced to a central IBM 370/165 computer facility. Using a photo
graphic enlarger, the film of a turbulent flow (Re=3500) was projected
onto the eleveninch square tablet surface. The particle images were
then entered onebyone using the tablet's locator pen. The images,
once encoded, were translated and rotated to a reference position by
use of the frame's reference points. His data consisted of 55 frames
and roughly 13,000 points whose entry into the computer through the
tablet was straightforward but tedious and highly susceptible to
operator error. Graphic output was used to verify data and quickly
showed most needed corrections. Once the image paths were identified
(manually), FORTRAN programs were used to find conjugate paths (particle
paths seen through both prism faces). More processing yielded the three
dimensional particle paths. (Jackman found approximately 150 pairs of
invertible image paths. Breton's system made all images invertible,
but he used only 1020 particles.) Use of the computer greatly reduced
the time used for data analysis; from the many months that Johnson*
used to a matter of hours. The procedures established by Jackman were
faster, but required a great deal of humancomputer interaction for
data entry and manual particle following.
Jackman's attempted but abandoned automatic particle following
procedure estimated the location of the next image from the last two
by assuming the particle's velocity constant and the measurement per
fectly accurate. A window was constructed around the estimate and
the next image of a path was sought in this region. This was a straight
forward, deterministic approach, yet it failed to follow particle paths
for a significant distance and made many errors in assigning an image
*Johnson mounted the film as slides and using a standard projector,
enlarged film records of a Re=6500 turbulent flow onto graph paper.
He then recorded by hand the particle images' locations, and manually
determined the motion in time of each image path.
to its correct path. As is evident from this (and has been known to
researchers in computer vision for a long time) programming a machine
to do even as simple a human task as merging points into lines is not
easy. If frames are superimposed, the particle paths are immediately
evident to a human observer, yet to program a machine to find these
correct paths is not elementary.
Qualitative Observations of Particle Motion
There are a number of useful qualitative observations pertinent
to the discussion of the data acquisition system and image path at
tributes. Direct observation of turbulence for Reynolds numbers from
3500 to 6500 was performed to obtain the following information. The
system was designed to study particle motion using a Lagrangian or
referential description of fluid flow. That is, motion of the fluid
elements is being observed which is strikingly different from more
common fluid flow descriptions using Eulerian or spatial descriptions.
For flows in the pipe, a strong axial velocity exists so that particles
neither back up nor traverse circles. The observer is given a feeling
of circular fluid motion, but if any specific particle is followed, none
is seen. Furthermore, particles are not observed (by the eye) to have
any random "jumpy" motion as is the case in Brownian motion and their
paths are seen to be smooth trajectories, typically of small curvature.
(Meant here is the curvature in the projection; the threedimensional
curvature is nearly impossible to observe since the particle's stereo
image and hence location in the pipe is not known.)
Particles are observed throughout the field of view to have widely
different speeds with faster particles giving the effect of being deeper
in the pipe. Bright, slowlymoving particles can sometimes be paired
to their stereo image. Observation of these shows that the motions
of the two stereo images are not particularly similar.
At low Reynolds numbers, large periods of laminar flow are seen.
Here, as expected, particles travel in paths parallel to the pipe
axis; slow particles near the wall, faster particles in the center.
As the Reynolds number is increased through transition, regions of
disturbed motion become more frequent. This turbulence occurs for
isolated periods and is referred to as a slug flow. As a turbulent
slug approaches the viewing station, observed particles begin to deviate
from their laminar trajectories with increasing violence of motion
until the slug has passed, when a sudden calming occurs. Observation
of particle paths at a higher Reynolds number (Re=4500) shows that the
projected paths vary in direction over time, smoothly, not erratically,
with at most 10 to 15 degrees of slope. This can be expected to in
crease for higher Reynolds numbers but it does indicate how dominant
the mean velocity is.
These observations imply that following images should be straight
forward since the expected paths are smooth and slowly changing. This
is not the case however. With Jackman's data, one major problem was the
large measurement noise incurred with the use of the tablet for data
entry and too few references for frame alignment. As the optical analy
sis shows, there are other significant errors (See Appendix A).
There are errors due to the pipeprism refractive properties,
perspective, camera imperfections, and particle movement. It is shown
in Appendix A, that during the exposure time, a particle may move on
the order of one millimeter in the pipe. This causes elongation of the
16
particle images and a loss of accuracy in the position of their center.
This error has been eliminated from future experiments by the use of
a strobed lighting system synchronized to the camera (Lindgren, 1977).
In addition, coma, caused by use of too low an f number on the camera
lens can be corrected by increasing the f number. This requires more
illumination which has been increased to some extent. An increased
f number also reduces the effect of astigmatism due to the pipe and
prism. The spot diagrams in Appendix A demonstrates the focusing
problems typical of astigmatism. A larger f number has a greater depth
of field bringing more of the pipe into focus and keeping the image size
small. Other errors are more significant.
Determining the location of a particle from the location of its
two stereo images requires transforming the image locations into inter
secting rays in the pipe. Both Johnson and Jackman wrote analytic
approximations that were adequate in a region near the meridional
plane. A full threedimensional analysis requires an iterative tech
nique to determine a ray's direction and is therefore much more compli
cated than the straightforward "prism equations" given in Appendix B.
It reveals that considerable error can occur in a particle's position
when using these approximations. This error increases as the camera
is moved closer to the pipe. The camera cannot be too far from the
pipe because the illumination is not adequate to sufficiently expose
the film at the film rate required and the lenses available. Thus,
full threedimensional analysis is necessary. The close proximity
of the camera and the pipe causes this perspective problem. It also
causes another error. When projecting a line parallel to the pipe
axis into the image plane, the two stereo images are not scaled the
same, i.e. the foreshortening of the two images is different. This
causes sampled image locations of the two image paths not to be above
one another, increasing the difficulty of finding stereo image pairs.
Perspective errors are obscured in Jackman's data because of the large
measurement noise. As the system is improved and measurement noise
reduced, they will become more important and if not considered, will
reduce the capabilities of the data acquisition system.
There are other problems intrinsic to systems using projections.
The process of projecting the threedimensional images onto an image
plane causes a loss of information. Two such projections are required
just to recover the threedimensional position of the particle. When
many particles are involved, images overlap and particles shadow
other particles. Stereo pairs of images cannot be readily identified.
A particle following machine must handle these phenomena to achieve a
reasonable performance level. The statistical approach taken in
Chapter III requires knowledge of the probability of occurrence of
these problems. In particular, two separate situations need to be
considered: overlap and confusion.
Overlap occurs when a particle comes between the camera and anoth
er particle deeper in the pipe. Since the optical system has a finite
resolution, there is a region, something like a shadow, that exists for
each particle in which a second particle will be overlapped. This region
varies in size and shape depending on the particle's location. Any
particle in the shadow of another particle cannot be resolved by the
data acquisition system. If the two particles' images overlap only
partially, the center of the resultant image does not represent the
true location of either image and an error will be incurred when the
location of this image is digitized. This problem will always occur
for a system with finite resolution and must be considered in any error
analysis. In qualitative observations, the overlap of particles was
observed infrequently. As the particle density is increased, overlaps
can be expected to increase in occurrence. It can be concluded that
this "instrumentation" noise will always exist, perturbing the true
locations of the images.
Confusion is the term applied to the occurrence of images existing
near each other in the film plane even though their respective particles
may not be close to each other in the pipe. This is a result of the
projection of threedimensional locations into two dimensions. There
are two types of confusion: necessary and unnecessary. Necessary
confusion occurs when invertible paths (particle paths with two image
paths) have images at a given time located in the same region. This
cannot be avoided. Unnecessary confusion occurs when paths that are
not invertible (paths with only one image path and hence not resolvable
into their threedimensional positions) have images which occur in
the neighborhood of an invertible image. This problem can be corrected
to some extent by limiting the illumination in the pipe to a volume
that can be resolved by the prism (only invertible images are produced).
Of course, images occurring in the same region but at different times
cause no problems. Confusion cannot be observed directly since the
individual projections give no information as to particle location;
any image could possibly belong to an invertible particle.
In summary, the possibility of overlap means that image locations
cannot be considered as perfectly accurate, some measurement noise will
always exist. The fact that confusion will always occur to some extent
means that image paths may not necessarily be considered isolated. There
fore, it cannot be assumed that a small window drawn around an image
that is larger than the system resolution will contain only one image.
This has implications about any image following scheme. If an estimate
is made of an image location, no matter how small a variance this esti
mate may have, it is possible that other than the correct image is
arbitrarily close to the estimate. Hence, no estimator can "isolate"
a path, i.e. a decision as to which image is correct will always have
to be made. Finally, to make matters worse, both of these problems
increase in severity as particle density increases. The analysis in
Appendix A shows that data collected by Jackman can be considered
relatively sparse. He had a particle density on the order of one
particle per cubic centimeter. The probabilities of overlap and con
fusion for this case are very small and allow the data to be reduced
even with the large measurement noise. Higher densities would have
less chance of being reduced correctly and hence would require a smaller
particle density or a smaller illuminated region. The sparse nature
of the data explains why it can be reduced at all.
The data acquisition system described has been shown to have
a number of errors present; some correctable, some not. Optical error
sources were identified which can typically be reduced in severity or
eliminated by better design. Basic phenomena of stereo projections were
discussed and quantified for use in the particle follower theory.
20
A simplified projection system is considered in Appendix A.
This shows that the characteristics of the two projections, while
directly related to the particle path, are not particularly similar.
This means that there is very little information available to identify
a particle's stereo images. This has impact on particle following in
the sense that it limits the possible techniques available. The pro
cedures developed in this work are limited to identification of image
paths as opposed to a more global approach that incorporates some
threedimensional information.
CHAPTER III
A PARTICLE FOLLOWING MACHINE
Overview
Presentation of the theory of a Particle Following Machine (PFM)
begins with the formal definition of its input and output. Then image
path attributes are presented and models developed. Next, the decision
process used by the PFM to follow a path is given. This process requires
use of a Kalman filter and, from its residual, an image feature vector
is defined. After the residuals have been modeled, the details of the
decision process are stated. The use of learning to find some of the
probabilities needed for the decision process is then discussed. Finally,
control of the PFM and measurement of its performance are presented.
Figure 31 is an abstract description of the PFM. The image data
array Sik. (described in the next section) forms the input to the PFM.
1j k
Using its "concept of paths," the PFM sorts images into groups which are
composed of all images of one particle. These image sequences are denoted
as paths P p = 1, 2, ...N, where N is the total number of paths found.
The PFM is then a decision maker; identifying patterns (paths) described
only by their "concept." This "concept" is formulated in terms of pro
vided rules (models) and learned probabilities. The remainder of this
chapter formally defines the PFM and develops procedures to accomplish its
task. A summary of the process is given in the last section. Figure 36
shows the details of the PFM described in the summary but may be used to
INPUT
THE PARTICLE FOLLOWING MACHINE
FIGURE 31
REGION FOR IMAGES
THROUGH TOP OF PRISM
imox
~ I 
im ax `*.
THROUGH BOTTOM OF PRISM
IMAGE PLANE COORDINATES
FIGURE 32
OPTICAL AXIS
(Uaa)
Jmin
j Sijk
DATA
FRAME
BOUNDARY
flow
DIGITIZED
23
aid understanding throughout the chapter as the specific formulation is
developed.
PFM: Input and Output
The input data to the PFM is assumed to be a digitized and "reduced"
version of cinematographic records. These records consist of pictures
(copies of the front image plane) taken every T seconds. They could be
produced by an analog or digital video source as well as a film camera.
Reduction of the picture data consists of estimating the center of the
finite sized images. (Even with known errors in shape and location of
the image, this is currently the only reasonable way of determining the
image location.) Unavoidable errors exist in these locations (measure
ments) which are due to the phenomena discussed in Chapter II (overlap,
optical abberation and particle motion). In the case of Jackman's data,
the error of manual entry of the image centers and reference points was
also present.
Continuous coordinate axes, (u,v), are defined for the image plane
as shown in Figure 3.2. An image center (ul, v1) is quantitized with a
resolution of (Au, Av). If i and j are the discrete indices of the u and
v axes respectively, the discrete axes can be defined as u = iAu and
v = jAv. Therefore, an image center (ul, 1) in the frame would have
integer indices i = l, j = 1
S' j = (roundoff to integer is implied). The
Au Av
digitized location, then, has an error of +Au and +Av for the u and v
axes respectively. The discrete axes have the same origin as their con
tinuous counterparts and have ranges corresponding to the size of a
picture frame. Without loss of generality, the coordinates were chosen
such that the u axis is in the same direction as the pipe's x axis and
the v axis is in the z direction if the camera and prism are aligned
properly (see Appendix A). Then, images typically move in the increasing
u direction (due to mean flow rate) and the top section of the prism
generates images in the top half of the front image plane; the bottom
section generates images in the bottom half. This forces top images to
have positive v (and j) values and bottom images to have negative values
while the u (and i) value will always be positive. The frame boundary
is defined by limits on (i, j):
u range: 0 i i m i 0 u r i Au
max max
v range: j mi j < j : j i Av v j max Av
min max mm max
where j max is positive and j min negative. Then the optical axis inter
max mnu
i Au
sects the front image place at (u v ) = ( max 0).
2
As mentioned earlier, the sample period is T seconds so frame k
contains images at time t = kT. Time is assumed to begin at zero for
frame one (t = 0) and increase to a final time t = k T where k
o f max max
is one minus the number of frames of data available. Therefore, the
range of k is: 0 k s k The data from an experiment and hence the
max.
input of the PFM may then be represented by the binary array
1 if image present
S =
ik 0 otherwise.
where (i, j, k) have the limits previously discussed. Array S can be
considered as the output of a fictitious sensor that performs all the
preprocessing; each (i, j) plane representing one film frame at time k.
A particle in the field of view (therefore having at least one
image) undergoes continuous motion that is sampled by the camera. Thus
it is recorded as a sequence of image locations which are defined as an
image path. In general, an image path P for particle p that is n samples
long may be written as the sequence:
Image Path P = {S.. S. ..."'' Sijn (1)
p 1j1 132 ijn
where any Sik = 1, k = 0, 1, 2, ..., n. Path P is then a subset of all
ijk p
image locations.
When a particle generates two image paths (stereo views), ideally the
two paths can be used to determine the particle's location in the pipe.
The two paths are defined as invertible image paths. One path will neces
sarily be in the top region of the front image plane (j > 0); the other
in the lower region (j < 0). When a particle generates only one image
path, its position cannot be determined and the image path is termed non
invertible. Several problems occur when the errors in the data acquisition
system are accounted for. Image locations have errors that make invertible
image paths as defined above not necessarily invertible. The inversion
requires tracing rays from the image locations into the pipe. Ideally
the two rays from a stereo image will intersect at the particle's posi
tion but the measurement error prevents this. In addition, images may
be inadvertently generated or lost. (See Chapter IV for examples.) This
can occur for a number of reasons, but it means that some images may not
represent particles and some particle's image may be missing thereby
adding spurious images that need to be eliminated and leaving gaps in
paths that need to be filled. Finally, there is the possibility of overlapping
images. This means that an image may belong to more than one path
prohibiting any possible onetoone relationship between images and
particles.
The PFM does not determine invertible image paths. It finds all
image paths in the data leaving the determination of conjugate paths to
a later procedure. Therefore, (1) represents the output of the PFM.
The PFM's task has now been set. Its input and output have been
stated precisely in this section. The following sections build the
internal procedures that the PFM uses to accomplish its task.
Image Path Attributes
Since the input of the PFM is strictly array S, it is important to
know the characteristics of image paths that it is expected to contain.
One must be careful not to give image paths characteristics found in the
threedimensional particle paths. The PFM can only operate on what it
can "see," i.e. the twodimensional projections. It is interesting to
note that very little is known about the input data. This is discussed
further in Chapter IV when the initialization problem is considered. What
is of interest here are the characteristics of an image path: What inform
ation is available; what is important; what can be determined:
From the definition of an image path given in the last section,
fixed attributes can be identified as follows.
1. An image path consists of a subset of S.
2. An image path will be in either the top or bottom part
of the film plane.
3. The time index of an image path increases monotonically.
Variable attributes include:
1. Image dynamic state (motion, velocity, acceleration) changing
with time.
2. Image location perturbed by measurement noise.
3. Confusion (local image density) varying along a path.
4. A conjugate image path may or may not exist.
5. Sample rate uncertainty.
Fixed attributes are hard and fast rules that define what constitutes an
image path. In any array S, there are many possible image paths, each
with its own qualities expressed by the variable attributes. To determine
the true image paths, these quantities must be measured and compared to
the internal path concept. The PFM, then, must filter the true paths
from the set of all possible paths (candidates). To develop a procedure
to accomplish this, the quality of the variable attributes must be
considered.
The dynamics of the image are, of course, closely tied to the
dynamics of the particle (see Appendix A). If a perfect dynamic model
of the particle existed, there would be no need for the experiment; the
structure of turbulence could be easily simulated. Since this is not
the case, it is necessary to model the particle (and image) motion. The
approach taken here uses a model that is known to be simpler than the
actual process. This is done to allow the PFM to be flexible and follow
images for different Reynolds number flows without having to modify the
model. The way that the PFM uses the model error to aid in tracking the
images will be discussed shortly. First, the image path model is iden
tified. Using the results of qualitative observations made of the
particle motion presented in Chapter II, an appropriate model for a
particle path is a constant acceleration process. This results in paths
with smooth trajectories and if the acceleration is small, they will
have small curvature. Image paths are assumed to have similar char
acteristics but only be twodimensional. If r(t) is the position of an
image, it may be written as
r(t) = u(t)i + v(t)j
where i and j are unit vectors in the +u and +v directions. The con
tinous path is assumed to be twice differentiable so velocity and
acceleration of the image may be written as
i(t) = v(t) = u(t) i + #(t)j
_(t) = a(t) = u(t) i + v(t)j
The components of these vectors describe the dynamic state of an image
and are used to make the image state vector,
u(t)
v(t)
t u(t)
x(t) = u(t)
u(t)
v(t)
At some initial time, t = 0, the initial state vector is
u(O)
v(O)
x(0)
(O(O))
Assuming a(t) is constant, we have
i(t) = v(t)
_(t) = a(t)
i(t) = 0
i+Fx
where
The state transition matrix 4(t), for this stationary system, is then
found (see e.g. Gelb, 1974). By definition,
P(t) = eF(t)
which is expanded to yield
4(t) = I + Ft +
22 3 3
+2 2.3
2 2.3
t2/2 0
0 t2/2
0 0
+ 0
0 0
0 0
0 0
fl
0
0
0
0
0
0
0
0
0
0
0
1 0 t 0 t2/2 0 1
0 1 0 t 0 t2/2
0 0 1 0 t 0
0 0 0 1 0 t
00 0 0 1 0
000 0 0 1
For the sampled system with sample period T, we can then write
{xk+ = (T) {x}
where
1 0 T 0 T2/2 0
0 1 0 T 0 T2/2
(0 0 1 0 T 0
0 0 0 1 0 T
0 000 1 0
0 0 0 0 1
This is a recursive state transition relation for an image based on an
assumption of constant acceleration. There is no process noise so this
is a deterministic process. Actually, the acceleration will be a random
process due to the measurement noise so the state becomes stochastic.
The expected range of the states are estimated from the qualitative
observations in Chapter II. The image position, of course, must remain
within the frame boundaries. The image's velocity, u, will be limited
by the fluid's maximum axial velocity as discussed in Appendix A. The
velocity, 7, will be smaller than u. Recall from Chapter II, that the
image paths were found to form at most a 10150 angle with the x axis of
the pipe. Hence v can be expected not to be more than 25% of u. The
accelerations, u and v are expected to be small and slowly changing. Note
that a zero acceleration reduces the image motion to straight lines.
Hence this model is expected to be adequate for short sections of an
image path. With a dynamic model of an image path and expected state
values, the PFM can begin to sort the proper image paths from all the
candidates.
In addition to the dynamic model, something is known of the measure
ment noise. It is assumed that the measurement can be represented as a
linear combination of image state and additive white gaussian noise, i.e.
z = H 4+ w
zk :I
where zk is the measurement vector, H is the output matrix and wk is a
white Gaussian process with mean zero and covariance R.* Since the
measurement consists of the (u, v) location of an image,
1 0 0 0 000
0 1 0 0 0 0
The combined errors due to the optics and the digitizer (manual or
automatic) are assumed to be gaussian and uncorrelated. The covariance
matrix is then written as
uR 0
R = R u
R= R
This is recognized to be only approximate since there are other identifi
able sources of error that could be individually modeled. This is not
done because, for the present, the measurements have large gaussian noise
*The use of w as the measurement noise should not be confused with some
conventions where w is process noise.
32
components (see Appendix C) and, in addition, a simple model is adequate
for the tracking problem where suboptimal estimation is accepted. It
should be noted that one source of error that was ignored and is not
particularly easy to model, was the frame alignment in the digitization
procedure. This error does reduce the performance demonstrated in
Chapter V, but the alignment problem will be significantly reduced in
future work (Lindgren, 1977) making it allowable to ignore this error for
the present.
The last three variable image attributes are more difficult to model.
Confusion comes from high local image densities causing less confidence
in the resulting decision. It is uncertain whether or not special pro
cedures could be applied to regions with high confusion. One possible
approach would be to use the conjugate image paths and perform three
dimensional tracking. In the pipe, the particles might be more separated
making the correct decision more obvious. This approach would require
knowledge of the conjugate paths and a trialanderror search of all
combinations of the candidate images would be necessary to find them.
(Recall that conjugate images are not yet known to the PFM.) This was not
incorporated in the current work and would require more consideration
before being implemented. Any benefit could easily be outweighed by the
added cost since the paths traced cannot be checked for exact accuracy,
only compared to what a human would choose. Furthermore, finding
conjugate paths is not a simple matter. Due to the measurement errors
and the characteristics of the optics, conjugate images do not have the
i location. Finding paths that have a maximum similarity of horizontal
33
(i axis) positions works for random measurement noise but does not con
sider the optical shifting and scaling that may take place (see Appendix
A). Hence for high densities where confusion is large, conjugate paths
would also be difficult to find, leaving the PFM with more work and
possibly no more information.
The path attributes, fixed and variable, allow a candidate path to
be rated as to its possibility of being a true image path. The next
section shows how the attribute metrics are applied.
Following Image Paths: The Decision Process
To begin the process, it is necessary that a partial image path
(a segment) be identified. This is referred to as the initialization
problem and is treated in Chapter IV. Then there are three possible sub
problems that can be defined: the forward, backward, and central problem.
The first two are labels given respectively to the extention of a segment
forward and backward in time. The backward problem is of interest
because it might be used to resolve some cases where abnormally large
measurement noise occurred or high confusion exists. The central problem
is a combination of the forward and backward cases. This involves con
necting two segments to make one long segment. Only the forward problem
is considered in the present work.
To consider the forward problem in detail, it is assumed that a
segment of path P for particle p has been followed up to time kl. This
is denoted as P k. Using this segment, a prediction is made for the
location of the next image in the path:
location of the next image in the path:
k= H() = Hk1
where xk (+) is the estimate of the state of the last image of segment
ki
P k and xk() is the predicted state of the next image in the segment.
This relation is obtained by taking the expectation of the measurement
equation and using the state transition equation. It is assumed that
the desired images) that extend segment P [kl will be in a window, W,
constructed around the estimate zk. The images in the window for the PFM
k
are defined as candidate images,
{g : (S = 1) and (i, j) is in W}
ijk
q = 1, 2, ..., qk
where q is the (i, j) location of image q in W at time k. Connecting
k
segment Ppk1 to a candidate image forms a forward link (simply called
a link) from time k1 to time k. Since there is uncertainty as to which
link or links are true links (those that accurately represent the real
particle paths), a link probability is defined as:
Sk = Pr{P connects to candidate q}
pq p k1 
Figure 33 summarizes these definitions. At the center of the window is
the prediction zk. Neighboring images (the i's) are shown to link to the
k
prediction by the candidate link vectors (v's) with probabilities
pq
The result of the decision process is described by the decision state
vector apk = {a(1), a(2), ..., a(q )} where a(q) is the decision of true
or false for the link to image q. The i subscript represents different
possible decision states. The number of different decision states depends
LINK PROBABILITY k
pl
1
7 CANDIDATE
Akl CANDIDATE ERROR
 PREDICTION
2kqk
II
I
k p
pq
j
il 
Sk
4^
IMAGE CANDIDATES AND ERRORS
FIGURE 33
on the number of candidate images, qk.* For qk = 1, (only one candidate
link), there are only two possible decision states: the link is true
(ai = {1}) or not (a2 = {0}). For q = 2, there are four possible deci
sion states: either both are false (a = {0,0}), one link or the other
is true (a = {1,0}, a = {0,1}), both links are true (a = {1,1}).
2 3 4
In general, the decision states have three categories (called Hypotheses):
HO: No true links  path stops
H1: One true link  path continues normally
H2: More than one true link  path branches
Figure 34 summarizes the link classes and decision state hypotheses,
where co is defined as the false link class, and cl, the true link class.
Also given is the number of states in each hypothesis. For example, let
qk = 3. Then there are three candidate links and eight different decision
states. Of these states, one is H0, three are H1, and four are H2. If
pk
all decision states, a were of equal probability, H2 would have the most
i
chance of occurring. Of course this is not the case as will be indicated
shortly. In general, for qk candidate images, there are 2 k possible
pk Aqk
decision states, pk i = 1, 2, ..., 2 of which one type is HO and q
are H1. The remaining are of H2. Each state is mutually exclusive making
the three hypotheses mutually exclusive. Therefore
3 2 k
SPr{H.} = 1 and Pr{a.} = 1
1 1
i=l i=l
As qk increases, the PFM must discriminate between an increasing number
of decision states making the work increasingly difficult.
*Decision states are the different vectors, a1, whose components are
link states a(q), q = 1, ..., qk"
K'
Link Classes
Decision State
Hypotheses
qk H H1 2
1 1 1 0
2 1 2 1
3 1 3 4
4 1 4 11
5 1 5 26
Number of
in the
Example for qk = 3
Link
Class
C0:
c :
H :
H1:
H2:
False link
True link
no links true
one link true
more than one link true
Decision States Resulting
Hypotheses H0, H1, H2
Decision
State
'1':
'O':
true link: c1
false link: c0
LINK CLASSES,
DECISION STATES
FIGURE 34
AND HYPOTHESES
a1 S2 a3 4 26 6 a7 48
0 1 0 0 1 0 1 1
0 0 1 0 1 1 0 1
0 0 0 1 0 1 1 1
To formalize the decision process, we consider a simpler example
as a guide. Assume that an object has feature vector C. Two possible
classifications exist for the object: co and cl. If there is no cost
for a correct classification and unit cost for an error, it can be
shown that the minimum error rate Bayes decision rule is:
Decide cl if Pr{cI } > Pr{col[}.
This says that cl is chosen if the probability of cl given the feature
vector C is greater than the probability of cO given C. (For a derivation
of this result, see Duda and Hart, 1973.) For the PFM, there are several
possible decision states to choose from. The feature vector becomes a
matrix of feature vectors and the problem is termed a compound decision
problem.
The possible decision states for the forward problem may be written
as (the subscripts p and k are dropped for simplification):
a., i = 1, 2, ..., n
i
where n = 2 Representing the feature matrix as E, the minimum error
rate Bayes decision rule becomes:
Decide a. if Pr{a. I} > Pr{a. i} V. (2)
i J J # i
where ac. represents one state (a vector of qk link decisions). Bayes rule
is used to calculate the required probabilities. This rule relates the
probabilities required in the decision rule (a posteriori probabilities)
to the a priori probabilities. This is written as
Pr{ia.} Pr{a,}
Pri 15E} = 
Pr{S}
where ai is one particular decision state. Formulating the decision
process in such a manner is of little direct use. Some assumptions must
be made to reduce the calculation of these probabilities to a more
tenable form. The term Pr{5ia.} is the probability of the features, E,
given a decision state, a.. Here, it is reasonable to assume that the
features are independent of one another. We can write
qk
Pr{ ia.} = Pr{(_ a(q)}
q=1
where the symbol n is used to express the product, and Pr{q la(q)} is
q
the probability of the feature vector q occurring for a classification
a(q) of the link to image q. Note that a(q) represents a choice of
either class co or cl. The term Pr{E} is the probability of the features
for all classifications. This term acts as a normalization constant and
can be ignored for the present situation. Finally, Pr{a.} is the
1
a priori probability of the decision state a.. As discussed previously,
1
the decision states have categories H0, HI, and H2. It will be assumed
that all a 's classified as H1 have equal probability. The same assumption
i
*
is made for H2. Of course, only one ai can be classified as HO. It
is therefore possible to compute the probabilities required for the deci
sion rule. In order to calculate these probabilities, the feature
vector must be specified and modeled.
This is the weakest assumption. The hypothesis H2 contains cases of
of single, double and higher branches which could be given different prob
abilities. However, a properly chosen window size will keep qk small
and the possibilities of higher order branches low.
Note on alternate decision rules. The Bayes decision rule is very
general and complicated. It is reasonable to ask; is this necessary?
Considering other alternatives, this decision rule does appear very
reasonable. An alternate rule might be to arbitrarily choose one
image in the window and assume the error generated would not be
significant to the results. A second possibility would be to choose the
image in the window that is closest to the predicted location.
The first rule might be acceptable if the window were very small.
This would require small measurement error and a better dynamic image
path model. The same would also be necessary for the second rule to
be viable. Improving the resolution and reducing the measurement error
of the system are certainly desirable. However, improving the dynamic
model can only be done after more is known of the structure of the
turbulence. This is what is being sought in flow visualization research.
In addition, the closest image to the prediction, because of confusion,
may not be the correct image. No matter how good the estimator, there
will be some variance. This manifests itself as a region of uncertainty
around the prediction and results in a finite probability of an image
intervening between the correct image and the prediction. Therefore,
the minimum error rate rule (2) based on the estimator error statistics
is used and its performance is expected to be statistically better
than these alternate procedures.
*
This rule would result from assuming the feature vector to be gaussian
and have zero mean. This would not necessarily be valid.
The Image Feature Vector
The feature vector, required by the decision process, should
be as simple as possible, but contain sufficient information for
decisions to be made at an acceptable error rate. With this as the
goal, the feature vector chosen is derived from the relative location
vector (candidate error vector), kq Recall that this was defined as
S = qq Hx ()
kq k k
(see Figure 3.3). This looks very similar to the residual in the
Kalman filter and the connection will be made shortly. First, the
Kalman filter used for estimation and prediction, will be presented.
Then the feature vector will be modeled and used in the final forms
of the decision rule.
Estimation and Prediction. The Kalman filter has been extensively
covered in the literature (see Kalman, 1960, Bryson and Ho, 1969,
Jazwinski, 1970, and Gelb, 1974) so only the results will be presented.
Refer to Figure 3.6 at the end of this chapter for a block diagram of
the decision and estimation system.
As described previously, the state x has a recursive transition
equation
k kI
where D is the stationary transition matrix. The measurement equation
was given as
z = H w+ wk
where wk is a zero mean, uncorrelated gaussian noise process with
K
covariance R. If x is the state estimate, its error is written as
k = k xk*
It can be shown that minimization of the cost function
T
J = E{ x Qx }
where Q is any positive semidefinite matrix and E{ } is the expectation
operator, for a linear, unbiased, estimator yields the following Kalman
filter relations (as in Gelb, 1974):
Prediction
Equations Pk = P (+) Tk
(+) = k( + Kk k xk
Update Pk(+) = {I KH P()
Equations
Kk= Pk() H HPk()H + R}1
Here, the symbols () and (+) are used to indicate the estimate before
and after a measurement is taken. The matrix P is the covariance of
the error, x, and K is the Kalman gain. The initial values for the system
are x () and P (). The first measurement updates these estimates to
o o
x (+) and P (+). A prediction can then be made and the process repeats.
o o
This is an optimum filter in that it gives an estimate with minimum
error covariance.
For the ideal case where the process is perfectly modeled by
the state transition equation, the covariance of the error becomes very
small making the gains, K, very small which essentially makes the output
independent of the measurements. The residual
u
= fk] =
k
:k v H
Vk
which is the difference between the measurement and the predicted
state estimate becomes a white noise (uncorrelated, gaussian) process
and hence contains no information. For the PFM, it is known that the
dynamic model, as for most real systems, is not precisely accurate.
This leads to divergence of the estimate from the actual state.
The problem of divergence is treated by Price (1968). To eliminate
divergence, solutions typically involve improving the system models
(using more states) or using a filter with only a finite memory. The
latter technique forces the system to ignore input in the distant past
and keeps the covariance matrix from becoming very small. This causes
the filter to "track" the input. Another alternative is to add arbi
trary process noise which also keeps the covariance from going to zero.
All of these techniques have good and bad points and their usefulness
is dependent on the application. The estimator that is used in the
present work is the finite fading memory filter. (See Tarn and
Zaborszky, 1970, Sachs and Sorenson, 1971, and Miller, 1971.)
The idea behind the finite fading memory filter is that new
measurements are weighted more heavily to make the filter "forget"
earlier measurements. This is appropriate for the PFM since it is
the tracking function that is desired as opposed to a truly optimal
estimator. The error covariance matrix for measurements that occurred
at time j is increased by
kj
R. = s R kzj.
J
The normal covariance, R, is a constant. The present time is k and
s 2 1. The resulting filter equations are the same except for the
covariance which is
P'() = s 0 P' T
k ki
where the prime indicates that this is a different covariance matrix
than normal. The value of s is empirically determined. The residual
vk has different characteristics as s is changed. This is demonstrated
in Chapter V. Essentially, s selects how much memory the filter has.
With s = 1.0, the filter reduces to the normal Kalman filter and all
past points affect the prediction and state estimate. This is not
desirable for the PFM because the residual may diverge. As s increases,
the past has less affect on the estimate causing the residual to always
remain finite. By doing this, the residual becomes more consistent
and therefore modelable as is shown in the next section. Intuitively
it seems apparent that if the residual remains finite, never crossing
zero (i.e. has a nonzero mean), a current value would be correlated to
the distant past. By using values of s other than one, the autocorre
lation of the residual becomes small between a present value and past
value, eventually becoming dependent only on its next to last value.
This is the assumption used in order to model the residual.
Modeling the Residual
The residual, presented in the last section, is a discrete sto
chastic process that "contains information" as to the difference between
the image dynamic model and the real process. Use of the residual to
modify parameters is a type of adaptive filtering which is a form of
machine learning. An example of a system that uses adaptive techniques
to modify the process noise is found in Jazwinski (1970) and to change
the system dynamic model in McAulay and Denlinger (1973). Other examples
are also found in Gelb (1974). However, the PFM does not use these
techniques. It operates at a more basic level because it must first
identify its next measurement. That is, it must choose which image in
the next frame i:ost likely belongs to the particle being followed.
To model the residual, several assumptions are necessary. First,
it is assumed that the residual can be characterized as a Markov chain.
Therefore, we can write
Pr{l k1, k2..'" } = Pr{Ikk1
That is, the joint conditional probability of k given all past residuals
is dependent only on k_1, not the entire history of the chain. In
addition, it is assumed that the statistics for each path':s residual
are stationary and applicable to all transitions on all paths (i.e.
u v
ergodic). Finally, it is assumed that the elements vk and vk of the
residual vector v are independent.
Using the assumptions made for the residual, the feature vector
for candidate q is specified as
= }T
kq 1kq' k
In words, the feature vector of a candidate image is a vector composed
*
of candidate residual, vkq, and the last residual yk. The matrix of
The last residual, v ,was the candidate residual selected at time
k1 and hence does not Aave a q subscript.
features, E, is then
= ^l k2' "'' kq k
This represents all the features used to make the forward following
decision for path P kl'
To actually make a decision, the decision rule is written out
explicitly. Recall that the probability of a decision state given the
feature is
Pr{'*a.} Pr{a.}
Pr{a E} = 
Pr{E}
Considering Pr{2a} = Hk Pr{e.la(j)}
q=l
where Pr{ qla(j)} = Pr{vkj, uk la(j)}Pr{rv v, lla(
because the u and v residuals are assumed independent. The probabilities
in this last relation are joint conditional probabilities of the last two
u and v residuals given the candidate's link state a(j).
At this point, it is useful to present an explicit example for the
situation qk = 2. There are four (2 ) link states given as:
Link Class
a(l) a(2)
a 0 0
1
2 1 0
Decision '1' = true
a 0 1
State '0' = false
a 1 1
4 _____
47
We can then write (while reverting to vector notation for the residuals)
Pr{~ljE} Pr{EII} Pr{~l} = Pr{kl,.l C o} Pr{vk2' 1'1 } Pr
Pr{a2E} = Pr{Ea2} Pr{a2} = Pr{vk, Ikl c} Pr{k2'kl c} Pr{2}
(3)
Pr{aB E} = Pr{(a3} Pr{ 3} = Pr{vkl'kI co} Pr{4k2,'kJ1co} Pr{(a
Pr{aIE} = Pr{Ea4} Pr{(4} = Pr{ykl,4 k cl} Pr{vk2' klci) Pr{a4}
To interpret this, the quantities on the left are referred to as the a
posteriori probabilities of the decision state, i.e. the probability of
the decision state given the features. These are only proportional to the
right hand side because Pr (E) has been eliminated from the denominator.
(Recall that this is just a scale factor.) The Pr(a ) term is the a
priori decision state probability and is conditioned by the likelihood
of the feature given the state, Pr(E Ia). Recall that the decision state,
aR, is the vector composed of all link class decisions. Then, for example,
i
with a, = {co, co}, we have
Pr{Bla Pr{v ,v klco} Pr{v v C}
Pr{ = kl'kl Pk2'k1l}
since features are assumed independent. The quantity Pr{.kl,_k_lco} is
the joint probability of candidate residual Yk, and the last residual k1
given a class co (false) link. As another example, consider
Pr{Ea2} = Pr{yl,kllcI} Pr{k2',k_ co}.
Here, Pr{vkl,klcl} is the joint probability of the candidate residual
1
and the last residual given a true link to candidate image k. After
making a few observations, a useable relation will result. Terms like
Pr(v kqv lc} = Prt{ Ico} refer to the probability of the feature
vector E = {(vkq ,vk} given a false link condition. It can be reasoned
q
that an image with a confusing false link can occur anywhere in the
window with an equally likely chance. Therefore, any one particular
feature vector has a probability of occurrence that is equal to one
divided by the total number of possible feature vectors. (See Chapter V
for a specific example.) Terms such as Pr{. Jc1} are the probability of
J
a true link. to the image q with feature vector E Figure 35 shows
K q
1 2
an example for the case qk = 2. Two candidate images, 1k and k are
found in the window around the prediction zk. Candidate residuals are
v and v. A past residual of is assumed to have occurred at
kl k2 k1
time kl. The sketch shows contours of a possible residual joint prob
ability Pr{vq ,U k cI} = constants L1, L2, or L3. If, for example,
k1 = d, then the joint transition probabilities Al andA2 are found, i.e.
A1 = Pr{v, dlcl}
A2 = Pr{vk2 dc}
as shown. The joint probability shown is unimodal but, in general, no
constraints are placed on its form so it could be multimodal. Since the
residual statistics are assumed stationary, it can be expected that this
quantity could be learned from examples of correct image path sequences.
Before discussing the learning aspect further, the a priori probabilities,
Pr{a) need to be examined in more detail.
The decision state, a., was previously classified according to three
hypotheses: the path stops, continues normally, or branches. It can be
S: prediction 1
1 2 2k
k, : candidate images
~V' k2: candidate residuals
u
u k1j
u
k2 2
Ik2
2
IMAGE PLANE CONFIGURATION
IMAGE PLANE CONFIGURATION
v
u_
v
"k2
U
v Candidate
Vq Residual
A
A2
I
L3 k Pr{vq ukl11 c1
= constant
u
vk1
Last
Residual
I k2
EXAMPLE OF JOINT TRANSITION PROBABILITY
FIGURE 35
expected that the chances are much greater that the path will continue
normally rather than stop or branch if there are no missing images and
the image density is sufficiently small to limit:the number of overlaps.
As shown in Appendix A, the overlap expectation is a function of particle
density and system resolution. The chance that a path stops is negligible
(if not near the edge of a frame) except in the case of data that were
input using the graphics tablet. For tablet data entry images are
easily omitted and will cause the performance to suffer. In essence,
the a priori probabilities define expectations of what the decision
results should be. They weight the calculations of the a posterior prob
abilities and, hence, directly affect the selected decision state.
The a priori probabilities vary as the number of candidates changes.
When there are no candidate images, qk = 0, it is assumed that the path
stops. Therefore, Pr{H0} = 1.0. For the case of only one candidate,
it is assumed that the path will most likely continue, linking to the
candidate (Pr{H1} = .99, Pr{Ho) = .01). For qk = 0 or 1, branching can
not occur so Pr{H2) = 0. When qk = 2, branching can occur so Pr{H2}
is just the probability of image overlap calculated in Appendix A. These
probabilities were found to be negligible for the case of two or more
images. Therefore, when qk = 2, decision state 4 = {cl,c1} is highly
unlikely. (Thus branching is assumed negligible.) If qk = 3, one state
is {clC]cl} which has a very small likelihood. Other states such as
{clcico} or {clcocl} would have the same probability as aq = {cici}. It
is clear from the analysis presented in Appendix A that the current data
available is a rather sparce situation since the chance of overlap and
confusion are so small. This is good for the PFM because the expected
number of decisions to be made is small and reduce to picking the one
most likely link. It is important to remember that as the density or
lighting situation changes, the probabilities of the hypotheses, HO, HI,
and H2, will change. Returning to the example for q = 2. Define
PO = Pr{ lIco}, q = 1, 2, so (3) can be written
Pr{allE} (P0)2Pr{H0}
Pr{c2IE} AIP0 Pr{H1I
2
Pr{ 3} POA2 Pr{H}1
2
Pr{aI} AIA2 Pr{H2}
where Al and A2 are the probabilities Pr{l ci} and Pr{ 2lc2}. Now,
it is assumed that Pr{H0} and Pr{H2} are very small compared to Pr{H1}.
Therefore, the decision amounts to a choice between states a2 and a .
Bayes minimum error rate decision rule (2) reduces to:
Choose a2 if A > A2 and a3 otherwise.
When A and A2 are very small, the probabilities for a a 3 and a are
small allowing the possibility of a Therefore, the decision rule becomes:
Choose al if POPr{Ho} > APr{Hi} and POPr{Ho} > A2Pr{Hi}
2 2
This just says that when neither link has a probability greater than the
chance of a confusing image being present, the path should stop decisionn
state a ). Finally, when A1 and A2 are both large, state a becomes
theoretically possible. The decision rule would be
Pr{H1}
choose a. if A Pr{H2} > POr{H
2
and
Pr{H1)
A2Pr{H2} > POr{H
2
This completes the example. A decision can be made as to which is the
most probable state. It is evident that the a priori assumptions are
very important in the stoppath and branch decisions but do not enter
in the normal condition where one link of the candidate links is chosen.
Similar results can be obtained for other values of qk.
Learning. An important aspect of the PFM is the determination of the
residual transition probabilities. It has been shown that this reduces
to determining the joint probability of the feature vector given a true
link, Pr{ Jcl}. Used here is a nonparametric procedure of identifying
the feature probabilities. (See Duda and Hart, 1973.) To determine these
feature likelihood, the PFM is given path examples that have been pre
viously identified by manual effort, or alternatively, by previous output
(which has been verified by an operator) of the PFM. An identified image
path is a sequence of images that have true links. Therefore, no deci
sion has to be made, the filter can be given the sequence as input. As
the filter follows this path a residual will be generated
V1' I 2' 3""'Vkl' yk'
This sequence contains examples of residual transitions for true links.
By using several paths (the training sample), many example transitions
are generated. The residuals are quantified and the transitions v to
k are counted to form a histogram of transition probability. There
are two histograms that result:
Pr{v v cl} and Pr{q v, v Cl
kq' k kq' ki
These are used by the PFM to make decisions on new data as discussed pre
viously in this section. The decisions made by the PFM are then dependent
on the training sample. As the characteristics of the training sample
change, the PFM's decisions will change and thereby adapt to the char
acteristics of the flow being studied.
Control of the PFM
With the capability of making a decision as to the link classes, the
PFM can follow image paths. At each step, the decision state probabilities
can be calculated and the decision state with the largest probability is
chosen. The three possible results are simply the state hypotheses as
presented earlier. It is expected that the normal response will most
often be taken (i.e. only one link found; Hi). When H0 is decided, the
path can simply be stopped. When H2 is decided, branches have been identi
fied. They can be saved for following after the primary path (the path
with largest link probability) has been found. Dealing with branches is
very difficult; many possible combinations of links occur. For example,
if a path branches and then comes back together or, in general, any case
where a path is split, there is a question as to whether this represents
the true case or is a result of confusion. The only real answer to such
a question is to do everything possible to avoid branches. It is worth
while to introduce the term isolated path. This is a path whose link
state decisions were all for the case of qk = 1. That is, no branch
decisions were necessary; the window around the estimate contained only
one image. This is what would be considered a "nice" path. If the
measurement noise is small enough and the estimator fairly accurate,
these paths are a reality. Data consisting of isolated paths would have
a high confidence level. Any other situations of confusion and overlap
would have a lower confidence level. As particle density is increased,
the possible confidence will decrease and more operator interventions
may be necessary to resolve confusing cases. The limit to the particle
density, then, is where the cases become so obscure that the operator
cannot be confident in his decision. The limit to the abilities of the
PFM to follow paths is clear. When the link state probabilities are
very close to each other (the decision boundaries), the decision is not
obvious; the PFM is indecisive. However, for reasonable measurement
noise and density, the PFM control decisions are expected to be correct
more often than not. By monitoring the probabilities likelihoodd) of
the decisions made, an operator can be warned of possibly bad decisions
and poor performance.
Measuring Performance
Perhaps the most frustrating aspect of the particle following
problem is the inability to know precisely what images belong to the
same particle. This means that any measure of error or performance of
the PFM is subject to interpretation. Obviously isolated paths can be
quickly identified by a human observer as well as the PFM. But when
confusion occurs, the true image path is only conjecture. Both the
55
human and the PFM must make decisions, and both will make mistakes. The
determination of performance of the PFM becomes a relative comparison,
i.e. does the path found by the PFM agree or disagree with the human's
belief? The human, of course, has the ability of taking more factors
into account than the PFM. However, the PFM is expected to be more con
sistent in applying its concept of what an image path is. If the PFM is
made more sophisticated, it would be expected that its decisions would
agree to more of an extent with the human. However, there is a tradeoff.
If a semiautomatic procedure is used that has a good performance and
only a few links are wrong and need human intervention, it might be best
to avoid making a more costly and complicated PFM. This can only be
judged with a great deal of experience using the program. For the
present, the performance of the PFM is judged against results given by
Jackman. No procedures were developed to handle special cases so the
error rate is not artificially decreased. Nevertheless, the performance
is fairly good as will be shown in Chapter V. When automatic data
acquisition techniques are used and the variance of an image's location
is reduced (probably by an order of magnitude), it can be expected that
the present PFM will provide satisfactory results. It is conceivable
that when errors are detected, the path in question can be discarded.
Alternatively, it can be kept if the choice is between two images that
would not alter the description of the flow field significantly. The
choice of how to deal with these "differences of opinion" allows the
experimenter to select his experimental accuracy and confidence. It
may be discomforting to know that the performance cannot be stated
quantitatively, but at least an acceptable level of confidence is
realizable.
Summary
The PFM has the task of identifying image paths (patterns) in the
data array Sijk. The approach taken here assumes that an initial path
segment has been identified which the PFM is to extend forward in time.
Referring to Figure 35, candidate images, J, q = 1, 2, ..., qk, are
found in a region (the window) around a prediction, zk, of the location
of the next image of the particle being followed. A decision must be
made as to which candidate images) is (are) most likely the correct
imagess. To do this, the probability of the transition vk1 to kq
is assumed known. The images with the most likely transition prob
abilities are selected as true links. If more than one image is selected,
branching occurs. The PFM continues following only one branch, saving
the remaining branches to be followed at a later time. If no images
pass the requirements, or no images are in the window, the path is halted.
The normal situation is for only one image to be selected with index q .
This defines the measurement, zkq*. A Kalman filter is used to generate
the prediction, and its residuals are used as features in the decision
process. The sequence z k = k k + 1, ..., k + n 1 is the
k o o o p
identified path, where k is the initial frame of the path and n is the
number of frames spanned by the path.
LU
I.
V)
V)
C)
0
Z r4
Z *
i z
0
LL
I
LL<
I
LL
I
J
Q
CHAPTER IV
PFM: INITIALIZATION AND IMPLEMENTATION
Presented in this chapter is a discussion of the procedures used
to implement the Particle Following Machine (PFM). The material pre
sented is not vital to the understanding of results given in Chapter V,
but will aid in the comprehension of where the results came from.
Three logically separate routines are outlined: the initialization
program (INIT), the Particle Following Program (PFP), and the learn
ing program (LEARN). Listings of the source codes are given in
Appendix D. Before presenting the details of the programs, it is
necessary to present the background reasoning to the initialization
problem. The PFM, presented in Chapter III, assumes that an initial
segment has been identified and precedes to follow the image path
using its path concept. Finding initial segments must be done with
much less information and is therefore rather difficult. Once the
initialization problem has been presented, the discussion of the
programs is preceded with a short description of the data structures
used by the programs.
Initialization
The Initialization Problem. If the initialization of the particle
paths were simple, the whole problem could be reduced to a series of
initialization steps. However, the problem of finding initial path
sequences is not simple. Required here is a scheme that can take an
arbitrary image location and find its closest successors or predecessors
58
representing the path. This must be done without any knowledge of the
path's attributes. Especially important and lacking is any information
concerning the particle's velocity (speed and direction). The data
contains particles traveling at different speeds and many different
directions. The fact that there is a strong general velocity in the
x direction is useful, but as the turbulence increases, it becomes
less dominant. Difficult cases to consider here are the possible
backward motion of a slow particle due to measurement error and the
possible close proximity of a slow and fast particle. In the latter
case, the slow particle has a much higher spatial sampling frequency
that may possibly obscure the faster path whose spatial frequency is
much lower.
Considering this problem in terms of spatial sampling gives some
insight to the difficulties involved. The particles have different
velocities so when sampled at a constant time rate will travel different
directions and distances. Since, in the initialization problem, there is
no information available concerning the particle velocity, the problem
becomes one of identifying the spatial sampling rate of paths of arbi
trary shape. To make matters much worse, large measurement noise
exists in data currently available making the spatial sampling rate
highly variable. If this were not the case, Fourier transform tech
niques could be utilized to advantage. Here, spatial transforms of
the data would indicate strong sampling frequencies and directions which
could be identified through filtering. The fact that paths are not
straight causes some problems, but as discussed earlier, the paths
are roughly straight over five time samples allowing one to look for
straight line segments. With such a small number of points in the
assumed line, however, measurement error obscures any trend in
the spatial frequency.
The beauty of using spatial transforms is that they utilize all
possible information available in the initialization step. This in
formation consists of knowing that in any region of the data array S,
there are straight lines that have different directions and different
spatial sampling rates. Therefore, the result of the initialization
step should be groupings of the images by path segments. The number of
groups are unknown, but each group contains up to a certain number of
images (the number of time samples considered in the region) whose
time indices are monotonically increasing. The reason for the unknown
number of groups is that the spatial dimension of the region under
consideration is arbitrary and hence may sever image paths leaving
a segment not starting at the time boundary of the region. This view
of the problem is that of a constrained clustering problem which could
also apply to particle following. The difference is that the particle
follower makes use of path history. Therefore, initialization is a
special case of the particle follower. That is, the probability of
an image in the next time sample connecting to one in the current
time is uniform; all images have equal chance. This is in essence
Jackman's initialization assumption. He then used a constant velocity
criterion to locate the next image. If only one image was located,
he assumed he had found the start of a path, otherwise he arbitrarily
chose one. This approach also returns to the problem of considering
all image path possibilities. As the image density increases, the
number of possible paths increases and the possibility of confusion
increases, thus making this approach far less desirable.
An Initialization Procedure. The final form of the initialization
program is an implementation of heuristic search procedures. A
particle in the first frame to be analyzed is chosen as the base or
reference. A window is constructed around this location with a size
sufficient to contain the next four images of the particle represented
by the reference image. All possible threeimage paths are then found
and a cost calculated for each. The path with the smallest cost is
possibly a valid path. If other paths have costs within 10% of the
minimum cost, they are retained for further processing. Threeimage
paths are not typically adequate for initial path segment identification;
there is too high a chance that three arbitrary images will have an
accidentally low cost. Therefore, a fourth and fifth frame of data
is considered. A new candidate path is formed by connecting each of
the lowcost paths identified in frames one through three to each image
in the window at frame four. The last three images of each candidate
are used to calculate a threeimage path cost as previously described.
The cost from the first three and the last three images are summed
to make a cost for the fourimage paths. Those paths with costs within
5% of the minimum are considered likely initial segments. Now, frame
five images are linked to the candidates. Again, the last three
images are used to calculate a cost which is added to the total cost
of the path calculated up to frame four. Finally, those paths having
costs within 3.3% of the minimum are selected as initial path segments
to be used as input to the particle following program. The criteria
for closeness changes from 10% to 3.3% because the costs are being
added. This causes the total cost to increase, so, by decreasing the
percentage, the criteria remains roughly the same. Note that this
technique does not guarantee that the paths found are the overall min
imum cost paths because only three images are used at a time. This
was done in lieu of using all possible fiveimage paths which would
be very time consuming and expensive.
To be more specific, we define the reference image position
(in frame k ) as r = (i j ). A region surrounding this point in
space and time (called a window) is defined as
i i < i
w w
min max
Window R = j m w
w w w
min max
k k k + 4
o o
The size of the window is selected to include five samples of the
fastest particle. The qualitative observations discussed in Chapter II
indicate that images typically travel in the +u direction with small
inclinations to the v = 0 axis. Allowing for the possibility of a
noisy, slow particle backing up, i is selected slightly less than
w.
min
i Then i is set at roughly six times the maximum expected
max
velocity. Finally, j and j are set to allow the fastest image
max min
to be inclined approximately twentyfive degrees up or down. With
the window boundaries fixed, the image data set can be searched
for all images within the window. The resulting list of images is
defined as the candidate image list. Assuming that there is a true
path segment among all possible paths that can be constructed from
the candidates, one approach would be to calculate the cost of each
path as a function of its five images. Then the minimum cost path
would be the path most likely to be the true path. In order to save
computation time, a slightly modified procedure is used.
A cost function is defined for any threeimage path. This cost
function is a heuristic whose specification was guided by qualitative
observations as well as trial and error. Basically, the cost is designed
to be minimal for a straight, evenly sampled path with a spatial sampling
rate (absolute speed) that is reasonable.
Consider a window containing n1 images in frame k + 1, n
images in frame k + 2, etc. Then, nI links connect the reference
o 1
image at O to r (images p in frame k + l,p = 1,2,...,nl).
q
Each of these links could connect to any of the n2 images E2
in the third frame (time k + 2) making the total number of possible
threeimage paths, n1 x n2. For each twoimage link, a first dif
ference is defined as
k to k + 1 : v r r V
o o 0 1 O p
k + 1 to k + 2: v r r V
o o 1 2 1 pq
The first difference is a vector quantity that is a first approximation
to velocity. The cost function, J is then defined for each twolink
pq
(threeimage) path segment using the reference image, image p and image q.
The form for the cost is
V v
J 1 
J = +
pq C C
The first term represents a change in magnitude of velocity, i.e. a
tangential acceleration. The second term finds the change in direction
(cosine of the angle) between the first and second link with the
result shifted to be zero when the directions are identical. This
term is an approximation to the direction of normal acceleration.
The reasoning behind the choice of this cost function comes from the
qualitative results that indicated that a short path segment was a
straight line, i.e. zero normal and tangential acceleration. Mini
mization of the cost yeilds a path with minimum tangential and
normal acceleration. By selecting a value of the scale factors C1
and C and a maximum allowed cost, J ,a the performance of the
2' max
initialization procedure can be adjusted to be specific for most of
the correct initial segments.
Threeimage segments are the shortest segments for which the path
accelerations can be calculated. It was found that there is a good
chance that incorrect candidate paths have a smaller cost than the
true path, especially where image density is high. Therefore, a limit,
C3, is placed on the maximum speed, I I, l'vql etc. Limiting the
allowed angle is not desirable since slow particles can have large
fluctuations in direction due to noisy measurements; the observation
of small inclinations being valid only for particles with nominal
speed.
65
Extending the search to five images reduces the number of errors
even further. To accomplish this, the search procedure looks for low
cost paths in the first three frames. Then, using frame four (k0 + 3),
each of the path segments in the reduced candidate list, is extended
r
to all images r~, r = 1,2,...,n and the first difference is found
from
qr = r r V
2 3 2 qr
Then, the cost of the last threeimage segments is calculated by
pq qr
1 q
+ r + 2
pqr
1 2
where p and q are the image indices for the candidate paths and r is
the index for images in frame four. The cost, J is added to the
pqr
candidate costs, J and the low cost paths are determined by finding
those paths whose costs are within 5% of the minimum. This procedure
s
is then repeated using the images r4, s = 1,2,...,n in frame five.
Here, the first difference is
rs s r
v = r r V
.3 4 3 r,s
and the cost becomes
qr rs
v v
2 3
1 
rs, qq qr I rs
3 2 2 3
J =+
qrs C
2 3
where q and r are indices of images in the candidate paths and s is
the index of images in frame five. This cost is added to the candidate
costs to give a total fiveimage path cost for a given candidate path:
J =J +J + J
tot pq pqr qrs
Those paths whose overall cost is within 3.3% of the minimun cost are
the final initial path segments and are used as input to the particle
following program. Chapter V shows the results of using this initial
ization procedure and discusses some of the special problems encountered.
The Data
Since this work focused on the theory of particle following,
new experimental data was not taken from the flow apparatus. Instead,
the data analyzed by Jackman (1976) was utilized. This consisted of
55 frames of raw image locations as generated by manual entry via
the data tablet as well as the manually identified invertible paths.
This data was collected using the small prism for a flow with Re = 3500
at 200C. A film rate of 25.7 frames per second was used to expose
TriX film through a 100mm lens fitted with a closeup attachment.
Images by Frame. Data, as entered via the tablet, consists of the
image locations in the film plane (i, j, k). At these locations,
array S, (see Chapter II), has the value '1'. The image locations
were entered systematically, but not necessarily in any sorted order.
Therefore, this data was sorted by i value which allows somewhat faster
windowing. Each frame of data contains different numbers of images
which, for the present purposes, are best stored as a vector of all
images and addressed via an index to the first image of each frame.
Three vectors are then used to represent the framebyframe data:
NPART, KINDEX and DATA. (See Figure 41.) NPART and KINDEX must have
dimension at least as large as the number of frames needed. (Not all
DATA
'jI1 J2 _2
INDEX
Frame 1
Pointer
N DATA
Frame 2 # of images # of images # of images
Pointer in Frame 1 in Frame 2 in Frame 3
Frame 3
Pointer
FRAMEBYFRAME DATA
FIGURE 41
ID START FRAME LENGTH
il J '2 J2 '3 J3
IMAGE LIST
ID I START FRAME LENGTH
I l i 2 J2 i3 J3
IMAGE LIST
INVERTIBLE PATH DATA
FIGURE 42
frames are necessarily used.) DATA has dimension at least as large
as two times the number of images since each image location has an i
and j value. The numbers obtained from the tablet were integers but
when used by these programs, they were stored as real numbers. This
allows addressing of the data as a complex array equivalenced to the
array DATA. Then, when convenient, both i and j values of an image
may be retrieved by only one index. The framebyframe data were
rotated and translated to common axes and, without loss of generality,
the origin was set at the lower right hand corner of the tablet as
opposed to the image plane coordinates used in Chapter III. This
merely shifts the j values making them always positive. The portion
of this data set used by any program is currently stored in core which
simplifies the programming and speeds up execution.
Invertible Path Data. The second major data set consists of the in
vertible paths identified by Jackman. These are the path sequences
P Each sequence has a length, a start frame number and the image
locations. This data set is accessed sequentially, path by path and
is not saved in core. Its structure is shown in Figure 42. The
paths are not in any particular order, but each path has an identifier,
ID, that was assigned by Jackman. Conjugate paths were the only
paths put in the data set, but they were not necessarily located
next to each other. When comparing the output of the PFM to the
paths in this data set, it was necessary to build an index that located
a path by its first point. However, because this data was manually
identified and punched, many errors occurred (they were typically small
since gross errors were corrected) making it impossible to compare
the results automatically. In addition, the paths in this data set
include generated points that were used to fill in missing data points.
The PFP was not built to handle this case since new automatic data
entry procedures will not have this problem (Lindgren, 1977).
The Initialization Program
The FORTRAN program written to perform the initialization procedure
is outlined in Figure 43. The main program, INIT, primarily handles
the selection of a reference image, various counters, and outputting
of the results. Subroutine REDDAT establishes the data arrays contain
ing the image locations, frame start pointers and frame lengths as
discussed previously. Options allow the selection of printing inter
mediate steps, plotting the results and storing or recalling the data
in the directly usable form generated by REDDAT. Input data to INIT
also includes selection of the time and space window boundaries. Sub
routine WINDOW searches the data for images in a specified window.
Data is required to have been sorted by i values which speeds up the
search significantly. Subroutine TRAK determines the minimum cost
paths. The FORTRAN source code is listed in Appendix D and contains
additional documentation of these routines.
The Particle Following Program
The implementation of the particle following machine is the
Particle Following Program (PFP). This routine takes the initial path
segments identified by routine INIT and attempts to follow them for
a specified number of frames. A simplified flow chart for the PFP
is given in Figure 44.
READ OPTIONS, INPUT DATA
READ IMAGE DATA
FOR FRAMES ONE TO NTS
PLOT IMAGE DATA FOR
FRAMES NTSO NTS
SELECT A REFERENCE IMAGE
I S ,
FIND ALL IMAGES IN WINDOW
BASED ON REFERENCE
FIND MINIMUM PATH(S)
FOR REFERENCE
NO
DONE?
YES
WRITE RESULTS
I
PLOT PATH SEGMENTS
END
SUBROUTINES
(REDDAT)
(WINDOW)
(TRAK)
SIMPLIFIED FLOWCHART OF "INIT"
FIGURE 43
SUBROUTINES
(PFPRD)
(LRNMAT)
(FILTER)
(UPDTE)
(WEST)
(DECIDE)
UPDATE )
END
SIMPLIFIED FLOWCHART OF"PFP"
FIGURE 44
This program begins by reading the options which consist of
writing or recalling the image data in a directly useable form (as
in INIT), and writing intermediate results. Subroutine PFPRD is
functionally identical to Subroutine REDDAT used by INIT. However,
larger data arrays are specified since PFP needs at least as many
frames of data as images in the desired path. LRNMAT reads the joint
probability matrices used in the decision process. These matrices
are the output of the program LEARN which is discussed in the next
section.
After the preliminaries, an initial segment is read. This con
sists of five sequential image locations assumed to start in frame
one. (Only paths starting in frame one were considered for the demon
stration of the particle follower.) The path state is then initialized.
The location of the first image is assumed to be the path's estimated
location. The initial velocity is calculated from the first difference
of the first two images in the initial segment, and the acceleration is
assumed to be zero. These values specify the initial state x ().
The initial covariance matrix is specified as
10
10
5
P () =5
1
1
which reflects fairly large uncertainty in the initial position, less
for the initial velocity and even less for the initial acceleration.
Initialization of the filter in this manner forces the first two images
to be used and their resulting error to be zero. Other possibilities
exist such as forcing all five images in the initial segment to be
used and calculating an approximate initial state from them. By
utilizing only the first two images, some errors in the initial paths
can be eliminated, and some incorrect initial paths can be eliminated
entirely. (Recall that the initial paths are not guaranteed to be
correct.)
Once the initial estimates are fixed, the Kalman filter sub
routine FILTER is initialized. This is a multiple entry routine which,
on initial call, sets the $, H, and R matrices used in the Kalman
filter equations. Entry UPDTE, calculates a new estimate using the
latest measurement, and entry EST, calculates the predicted state
(the estimated state prior to a measurement). Once the updated initial
state has been determined, PFP enters a loop that attempts to follow
the path as far as possible.
The following loop consists of making a prediction by calling EST,
windowing the estimate with WEST, choosing the next image with DECIDE,
and finally, updating the state estimate using the new image (measurement)
with UPDTE. Windowing the PFP is twodimensional since only the images
surrounding the prediction, at the same time, are required. The actual
search for candidate images is the same as for subroutine WINDOW used by
INIT. Subroutine DECIDE, calculates the candidate residuals, and using
the joint probability matrices, selects the candidate image that would
produce the most probable error. A check is made to determine whether
the decision is isolated and a message is printed if it is not. In
addition, a simple check is performed to determine whether the
choice made was different from the minimum error choice (choosing the
image closest to the predicted location). If it was different, another
message is printed. Finally, if the probability of the chosen image
is less than the minimum allowed, DECIDE sets a flag to stop the PFP
from following the path further. After a candidate image has been
chosen for use as the next measurement, the UPDTE entry to the filter
routine is called to update the state estimate. The followingloop
continues until the path is stopped by subroutine DECIDE, or the speci
fied number of images have been identified for the path. The results
are then printed and the next initial segment is used. After all
initial segments have been followed, PFP is finished. The PFP source
is given in Appendix D with additional documentation.
Program LEARN: Determination of the Joint Probabilities
Routine LEARN is very similar to PFP. it takes a given path
(identified by earlier manual effort), and applies the Kalman filter.
The transitions of the residuals are quantified and counted, thereby
developing histograms of the joint probability to be used in PFP. A
simplified flow chart is given in Figure 45.
A path is read by the program and initial state estimates are
made in precisely the same manner as in PFP. The same routine FILTER
is used to propagate the state estimate over time. After the path has
been "followed" (the measurements in this case are known beforehand),
the residuals are added into the probability histograms by routine SAVE.
After all the paths to be analyzed have been studied, the resulting
joint probability matrices are output. The source and additional
documentation for LEARN is given in Appendix D.
SUBROUTINE
(FILTER)
UPDATE )
(EST)
(UPDTE)
(SAVE)
SIMPLIFIED FLOWCHART OF "LEARN"
FIGURE 45
CHAPTER V
RESULTS OF PERFORMANCE TESTS
Initialization
Program INIT was used to find initial path segments composed of
five images with the first image in frame one. The data used was col
lected by Jackman (1976). Identification of initial path segments
is dependent on the values of parameters C1, C2 and C3 used to calculate
the cost of threeimage segments. To observe how these parameters af
fect the cost, consider a threeimage segment (il, j1), (i2, j2), and
(i3, j3). Define vector A as a vector from image one to two, and
vector B as a vector from image two to three. Then the cost function
may be written as
Total = T + N
where JT = IIA
CI
A B
1 1
IA
JN
C 2
The cost is then a sum of two costs, JT and JN which can be interpreted
as instantaneous tangential and normal acceleration costs respectively.
Figure 51 shows the costs JT and J as functions of IBI IAI and
76
JT: COST OF TANGENTIAL ACCELERATION
JN: COST OF NORMAL ACCELERATION
TANGENTIAL AND NORMAL ACCELERATION COSTS
FIGURE 51
A B
II I, (= cosine of the angle between A and B) for various values of
C1 and C2. As the parameters C1 and C2 increase, the cost of a specific
difference in speed or angle decreases.
This, in effect, decreases the cost of the vector velocity varia
tion between A and B. For the results given here, the parameter values
found to give best performance were as follows:
C = 40.
C = 1.0
C = 30.0
Recall that the last parameter, C3, is a hard upper limit on A and B.
The upper limit on cost, J was arbitrarily set to one. Therefore,
max
many combinations of J and J exist that produce costs less than Jmax
T N max
Figure 52 shows the contours of Jtotal. The contour, Jtotal = 1.0,
is the limit, so any point within this contour would represent an
allowed threeimage segment.
Figures 53 and 54 show the initial path segments as found by
INIT. The locations of the symbols '1' through '5' represent locations
of images in frame one to five. (This figure is just the superposition
of the first five frames of data.) A proper path segment consists of
five sequentially numbered images. Figure 53 shows data from the
"top" section of the frames, i.e. one stereo view. Table 51 describes
labeled segments in these two figures. Important to note here, is that
the output of INIT contains several different kinds of errors. Most
errors are a result of data entry from the tablet. When using the
tablet for data entry, numerous images were accidentally omitted or
79
Lj
C) LL
Fj
Li Z
'LI) (O ~ I
LL
I
Li
Li
tL
U
Lr
I
( 1
^ "'irKj
I in bi
vl
C0i
T
j a
_I
'ii
0 3 .4
"\ "
II.
'^ : r '
01 1
IT' r of
^ 01) ^^m
r., n1 Id
U5,
* ~ j
I
^ ri
42 7
"ji Ti "
.1" ?~3
'_1
j;
I
ZUJ
UJ CL
LO
C0
UJ
LLL
F C)
 t
z
I
LF'
Un
4 id
:3. 1
* 1 "' ,1
:~i z i
ryi j 'i1
a
** t: fl J
Li _I N1
L1 I .,
riv3
:1 K
_4t
_I'1 1 i
~"3
(V rn Lj
(11 m
t' l mj.
b/ ,1 \ N
&
I J^^ J 1
'.
7. li
LA
1 1
IJs
Ri '
(Cf
r.J
II2
IiI
1
j S
S J 1_
.I.
Zr Lb I, Y
1 (Xi
U. I
Z r I1.
Uir itI.
U, 11I
LT J
~L~rr
LUC
CD0 C
U
z.
Lw
C:)
LL

J0
o
2II
I
TABLE 51
KEY TO FIGURES 53 AND 54
NOTES (SEE FIGURE 5.3)
NOTES (SEE FIGURE 5.4)
LABEL
LABEL
A NORMAL FIVEIMAGE SEGMENTS
B TABLET MISPLACED IMAGE ONE
C MISSING IMAGE IN FRAME TWO
D PATH SEGMENT BEGAN OUT OF FIELDOFVIEW
E IMAGE INCORRECTLY STARTED, "JUMPS" TO
ALTERNATE PATH
F IMAGE MOVES OUT OF FIELDOFVIEW
G TABLET MISPLACED IMAGE FOUR
I,K INITIAL PATH COST TOO HIGH
J ISOLATED IMAGE ONE
L DARK PATH TOTALLY WRONG "CUTSACROSS"
PATHS
M FIRST THREE IMAGES CORRECT
N DOUBLE ENTRY OF IMAGE ONE
P UPPER IMAGE ONE LINK INCORRECT
Q IMAGE THREE OF LOWER PATH INCORRECT BUT
CORRECTABLE
entered twice. In addition, a failure of the tablet to correctly
locate an image could occur because of operator error. This caused
some image locations to be displaced to the left as seen in Figures
53 and 54. These errors are not considered in the overall perform
ance because they are expected to have been eliminated by technique
modifications as discussed previously. The significant errors in
these results occur when the image density is high and incorrect initial
paths are selected. Some of the errors can be corrected by the particle
following program (see last section, Performance of the PFP). Another
error is not starting a path for an image in frame one. This usually
results from the cost of the path being too high and can be corrected
by increasing J However, this must be done carefully because many
max
erroneous paths may also be allowed. The parameter values chosen re
present a tradeoff in this regard. Very few "good" images in frame
one were lost.
Performance of INIT. Table 52 contains the results of initialization
of frameone images. The different types of errors that occur are
categorized by whether or not an image in frame one was "started,"
i.e. an initial segment was found starting with a frameone image.
An initial segment can have one or more incorrect members, be in
error because of a double entry of the image's location, or be the
result of an inadvertently added frameone image. If an initial
path is not found for an image, there are a number of possible reasons.
The cost of a correct path may be too large, images may be missing or
inadvertently added, the tablet could have caused a displacement of an
image, or the correct path does not have five images in the field of view.
Except for the first of these errors, too high a cost, these errors
are a result of using the data tablet and are unavoidable. Of the
254 images in frame one, 47 were classified as images with unavoidable
data entry errors and were disregarded. Of the remaining 207 images,
148 were started correctly, 31 were started with at least the first
two images correct, 11 were started incorrectly (totally wrong) and
15 were not started because their cost was too high. These results
are expressed as percentages in Table 52. Note that initial segments
with the first two images correct are sufficient to start the PFP and
may possibly be correctable. Initial segments that are totally wrong
may be followed by the PFP but would be expected to be aborted as
short (up to fiveimage) paths since the chance of these inadvertent
paths being longer is very small.
TABLE 52
INITIALIZATION PERFORMANCE ON 207 FRAMEONE IMAGES
IMAGES STARTED
CORRECT FIVEIMAGE 72%
SEQUENCES CORRECT OR
87% CORRECTABLE
INCORRECT BUT WITH 15%
THE FIRST TWO IMAGES
CORRECT
TOTALLY WRONG 5%
IMAGES NOT STARTED
COST TOO HIGH 8%
Particle Following
Learning Results. Program LEARN generated the joint probabilities
as described in Chapter III from a training sample of fifty arbi
trarily selected image paths that had been identified by Jackman. The
resulting probability distributions are dependent on the finite fading
memory parameters. Figures 55 through 58 show the distributions for
s = 1.0, 1.2, 1.5, and 2.0. Using other arbitrary groups of fifty
paths gives reasonably similar results. These graphs are interpreted
in the same manner as the example discussed in Chapter III. In general,
these results are fairly similar for the different values of s. As s
is increased, there is a slight change in Pr {v vql _} to change from
positively correlated to negatively correlated. The Pr{v qIv 1} also
becomes negatively correlated. Some distributions are seen to be
multimodal, but only slightly, All of the distributions have non
zero means and show some dependence (correlation). Checks were made
of the crosscorrelations between u and v residuals and no significant
dependence was found. Furthermore, correlations were not found to be
significant between kq and v v etc. Therefore, the assumptions
kq k2 k3
made in Chapter III about the residual appear reasonable.
The probabilities for the u and v transitions are seen to be dif
ferent in character. One possible explanation for this comes from
considering the odd image shapes. Analysis of the tablet error in
Appendix C assumed images to be circular. As shown in Appendix A,
particle movement causes the images to be lengthened in the u direction.
Therefore, the variance of the u location of an image should be larger
than for the v location. This effect shows itself in the u residual
i .:~. ...
J l 1_ 1.i IJU'1
'4
I H'fi I II ]
1HFFI] "I *
00'[ 1
'ii
.III
I
C)
0
I<
I
L J
'1
I,
0
LA
.r a
Q:
[l
IJ
I
II
__1
S
2
a
= :i
2
r
0
IC
z
01
I i
I I
0I 13
U U I 
O
1
I 

CT;>
0
I~II
 I
00I L  ii1 II
1 j rl r I
>
S
a
! .\ L 
4
0
o
Lii
)
2i
I1
0
(I
1 IlL
o LI 3
]tritim jy I ;G N t
II
.
LUi
ZD
(
z
I
C3
, '
I I
II
,
V)
LJ
Q
F(
U
h
LT
oirl ii
^ ^
r,F
I
I 
vV\\\
S ,, i
/II ,
." ,,,' i
,~ .
II I I
I;II, i,
3 ^
s.
CL
"" .
in
:
?i '
II I I i LI 1
2
> 
2
S.
_
II
z.
LLU
I4
0
CO
SD
Oj:
U I i i '
i o
Ln
H I _I
II
iC)
I I

I
J
I
l
I)
u ;
CI
I
2'

4'
o
II
0I
z
0
:I
i
iO U 1 1
liii11
i
1
/ ~~
 I ~I.
r
I
>
C
Q.
LF  
I "S 011 '
:*j A L" .N
4
S.
0
0
I lli 11 H
lJPl If I IX i3
r4
o
II

Cr:
0
I
z
0
<_)
00'1:
...[.
u
LII
i <,

IC
a
I ,
II
C/)
<:
*J
n11 
U I
Li i t' 1I
u I
1 I
P
