NASA Technical Memorandum 101306

Application of Optical Correlation Techniques to Particle Imaging Velocimetry

Mark P. Wernet Lewis Research Center Cleveland, Ohio and Robert V. Edwards Case Western Reserve University Cleveland, Ohio

Prepared for the Conference on Sensors and Measurement Techniques for Aeronautical Applications cosponsored by the AIAA, NASA, and AFWAL Atlanta, Georgia, September 7-9, 1988
I '
L

NASA

MAY

2 5 i369

UNGLEY RESEAKCH CE14TE2

Application of Optical Correlation Techniques to Particle Imaging Velocimetry

Mark P. Wernet NATIONAL AERONAUTICS AND SPACE ADMINISTRATION Lewis Research Center Cleveland, Ohio 44135

Robert V. Edwards Case Western Reserve University Cleveland, Ohio 44106

Abstract Pulsed laser sheet velocimetry yields non-intrusive measurements of velocity vectors across an extended 2-dimensional region of the flow field. The application of optical correlation techniques to the analysis of multiple ex osure laser light sheet photographs can reduce andfor simplify the data reduction time and hardware. Here, Matched Spatial Filters (MSF) are used in a pattern recognition system. Usually MSFs arc used to recognize characters or assembly line parts. In our application, the MSFs are used to identify the iso-velocity vector contours in the flow. The patterns to be recognized are the recorded particle images in a pulsed laser light sheet photograph. Measurement of the direction of the particle image displacements between exposures yields the velocity vector. In this work, the particle image exposure sequence is designed such that the velocity vector direction is determined unambiguously. A global analysis tcchnique is used in comparison t o the more common particle tracking algorithms and Young's fringe analysis tcchniqnes.

*
CU

spacing is inversely proportional t o the velocity magnitude. Measurement of the fringe spacing and orientation gives the local velocity vector. By repeating tlie process over a grid of points on the transparency, a 2-D velocity vector map is constructed. For low number density double pulsed photo raphs, different processing schemes are applied. These ow number density processing techniques are referred to as Particle Imaging Velocimetry (PIV). Due to tlie low particle density, individual particle displacement vectors call be determined, instead of the average displacement of a group of particles as in laser speckle. Particle tracking algorithms employing image processing techniques can be used to reduce the recorded PIV transparencies. Alternatively, the point-wise laser beam readout technique can be used to determine the local velocity vectors. Hence, different processing schemes may be applied depending on the recorded particle number density.

'i

'1'11~t.cc:hnicl~lc rccortling r~~llll.il)lc of speckle pal.l.ern translations that occur durin surface motion has beconle a useful tool to measure in-p ane displacements of solid srlrfaccs.l-2 Current, interests are focused on a.pplying this t.ccl~niq~~e rneaslrrcmcnl, of particle l o tllc tlisl)lil.cc~r~a~t.s Illlid flow fiel(ls.~-RI)ol~l)le i l l 2-11 1)11Isir1g a laser light sheet perpendicular to the direction of observation yields a record of the particle positions a t of 11t11nl)cr two ir~stances 1,irrlc over I,hc entire plane. '1'11~ density and size of particles in the fluid determines the data analysis approach. For a large number density of recorded particle images, the laser speckle data analysis procedure is used. This technique is called Laser Speckle Velocimetry (LSV). In LSV, the recorded particle photographs are usually analyzed by illuminating small re ions by a low power laser beam. The far field diffraction pattern of the illumination beam contains speckly Young's fringes. The fringes are oriented perpendicular to the local velocity vector, and the fringe

The point-wise laser beam readout of a recorded particle image transparency is a tedious and time consuming process. The point-wise analysis of the transparency is similar to laser anemometry techniques. In their present state, both PIV and LSV techniques are similar in'nature to the point-wise measurements obtained using laser anemometry, which arc time consuming and have positioning constraints. For a typical 100cm2 region of a transparency analyzed on a l m m grid, there are on the order of 104 poir~tsto be analyzed. At each point, the far field fringe lrallcrn spacing ant1 orientation must be dctermincd. 'I'his cxl~a~lsting arlalysis technique is most tlaba efficiently handled by an automated procedure. Ilowcvcr, even an automated procedure employing spccializcd array processors can still rccluire I~untlredsof hours of processing time. For all of the effort that has been expended to analyze the particle image transparencies and determine the fringe spacing and orientation at w l O 4 points, there still exists a 180' directional ambiguity. The fringe patterns generated by the readout beam are symmetric, hence, the velocity vector is known only t o *180. More information must be included a t the recording stage of the process to resolve the direction unambiguously.

In this paper, we describe a 2-D optical parallel processing technique for analyzing recorded PIV transparencies with unambiguous velocity vector determination. The techniques of Matched Spatial Filtering used in character recognition systems are employed to identify individual particle image records. The directional ambiguity dilemma is eliminated by using a n asymmetric laser light sheet exposure sequence. The optimal exposure pattern must yet be determined. Three unevenly spaced pulses yield enough information for determining the flow direction, but also can result in many false correlations. Other less ambiguous patterns are being investigated. The analysis procedure can be done optically using a magnetooptical spatial light modulator. T h e light modulator is used to display Hartley Binary Phase Only Filters (HBPOF) used in an optical correlator. The BPOFs are written on the light valve through a computer interface. Light valves with 128x128 elements are currently available. Alternatively, the photographs to be analyzed can be digitized and the MSF process performed digitally using Fast Fourier Transforms. This approach requires supercomputing facilities, or specialized array processors, but does not suffer from the limited spatial resolution of the commercially available light valves. PIV and LSV Double pulsed velocimetry enables the instantaneous recording of the in-plane velocity vectors across an entire 2-D flow field, see figure 1. The fluid under study is seeded with small particles. Illumination is supplied by a thin sheet of double pulsed laser light. The particle positions a t these two instances of time are recorded on a photographic plate. A switched or pulsed laser source supplies the light pulses. The pulse lengths, dt, are usually on the order of 20-30 nanoseconds.8 The short pulses are needed to record the images of the micron sized particles before they move enough to distort the s image. The pulsing interval At, which can range from j to seconds, is chosen based on two criteria: 1) to limit thc maximum particle in-plane displacements to a clesircd value, and 2) to limit the out-of-plane displacement to less than the light sheet thickness Az. The first criterion ensures that the recorded imagcs arc not so far apart on the imagc planc that a practical data analysis proccdurc cannot bc usctl. 'I'hc sccontl rcst.riction mcrcly i n s ~ ~ r c s that most of thc particles are within thc illr~mination volume during both p ~ ~ l s c s .
y.

obtained, but discrete particle ima es may be recorded instead. In either case, one must c oose the best technique for reducin the data. Although there are some techniques whic attempt particle tracking, most approaches rely on the beam readout technique. The beam readout technique can be used for both the PIV and LSV recording regimes. The photo raphic transparency is divided up into a grid ok points A typical photograph may contain 104 grid points. At each grid point a low power laser beam is used to illuminate a small re ion, (xl rnrn diameter). The particle images within tRe illuminated region act as paired sources producing Young's fringes in the far field. The fringe pattern spacing is inversely proportional to the average particle displacement within the illumination region. The orientation of the fringes is perpendicular to displacement direction. The Young's fringe patterns are digitized by a CCD camera and processed in a computer. The processing techniques vary depending on the available hardware, but the primary operation is to Fourier transform the fringe pattern. The Fouricr transformed fringe pattern is used to compute the power spectral density of the fringe pattern. The power spectral density is a symmetric function containing three peaks, see figure 2. The displacement of the side peaks from the origin peak is proportional to the vclocity. The beam readout technique does not directly provide any information on the sign of the velocity vector direction. The estimated vector a t each interrogation region is known only to *180.

A similar approach is to compute a select number of 1-D autocorrelations instead of the 2-D autocorrelation. This significantly reduces the number of computations. The 1-D autocorrelations are obtained by integrating the 2-D interrogation regions along the x- and y-axes. The integration can be done either numerically in the computer after digitization by a 2-D CCD array, or optically using cylindrical lenses before digitization by two 1-D CCD arrays.10 The particular approach will depend on available hardware, although the optical integration will be faster. After tlie 1-D projections of the 2-D interrogation region are obtained, 1-D autocorrelations are computed to determine the s- and y-components of the local vclocity vector. Thcsc 1-D analysis techniques work well when only a single particlc image pair is present in the interrogation region. When lnorc than two imagcs arc present, sornc spurior~speaks can arise in the autocorrclations, which can lead to i~~tlcl.crrninatc vclocity vcct.ors.
A basic problcm with all of these PIV ant1 LSV data reduction tcchniqucs is tlic directional ambiguity. Thc particlc tracking algorithms can dcterminc thc velocity magnitude, but a &1SO0angle ambiguity still rcniains. 'I'wo modified PIV tcchniqucs which address thc tlircctional ambiguity problem arc the image shifting t,cchniqr1cll~l2, holographic image scparation.13 and
111 1,110 irn;lgc shifl,ing I,ccllniquc, thc irnagc ficltl is shift.al t)ct,wccn thc two illumination pulscs, such t.hat. the second image of cvcry particle image pair is rcct,ilinearly displaced from the first image. A mirror ~nountcd a galvanornctcr scanner in the recording on system optical path enables the image shift between pulses. Thc magnitude and direction of the image shift is set such that even the most negative velocities are still positively displaced. This image shifting aids in the particle tracking algorithms since it sets a bound on tlie most probable position for the second particle image. For example, if the shift is a positive upward shift to the right, the tracking program knows to look for t,hc sccond

Thcre are major tliffcrcnccs bctwccn thc application of laser spccklc to surfacc displaccmcnts and particle tlisplaccmcnts. 7'11c main tliffcrcncc is that t,hc light scattering propcrtics of small particlcs suspcndctl in a f l ~ ~ iarc tlisl.inct,ly different than t,hc l i Ilt scat,t.cring tl cl~aractcrist.ics a solitl sorfacc. 'I'hc f rlitls arc of i l l ~ ~ ~ ~ ~ i ~ ~ a llighl,lshool, of li11i1.c I)y ;I. , ( ! ( 1.hick11oss Ax. Ilcncc! sc;~t,l,orir~g oc:c~~~.s ;L vol~lrnctlisl,ribr~Lior~ scat.l.cring frorn of sitcs instead of a surfacc distribution. The scattering pa.rticlcs arc ~lsoally small (0.1pm t o 10pm). The nr~mbcr density and particle sizc vary dcpcnding on thc fluid and the application.9

For speckle patterns to occur, the number of scattering sites per unit volume must be high enough so that the diffraction limited images overlap in the ima e plane.5,8 If the number density of particles is low, or i the scattering properties of the particles are insufficient for exposing thc photograph, laser speckle may not be

particle image up and to the ri ht of the first image. The sign of the velocity is assigned ased on the magnitude of the measured displacement. Displacements larger than the rectilinear ima e shift are positive velocities, while displacements smafler than the shifted distance are negative.

The holographic ima e separation technique permits the unambiguous identification of the first and second particle images. The double pulsed particle images are recorded holographically on a photographic plate. Two separate reference beam arrangements are used for each exposure. In this way, the first and second positions of all the particles can be viewed separately by subsequent reconstruction using the two reference beams. The separation of the particle images permits the computation of the cross correlation function which gives the displacement and unambiguous direction information. All of these techniques discussed are either computationally intensive or have limited applicability, and yet very few offer unambiguous direction determination. There is a strong need for a fast data reduction technique which also provides unambiguous direction information. Matched Spatial Filtering Pattern recognition is a widely investigated optical processing technique. The purpose of a pattern recognition system is to determine the presence and location of a reference pattern in an input scene. The primary operation in a pattern recognition system is the correlation of an input scene with a reference pattern. A matched spatial filter is the optimum filter for pattern recognition since it maximizes the ratio of the peak signal energy to the mean square noise energy in the output of the processing system. The signal is any object of interest in the scene, and the noise is the rest of the image. A coherent optical processing system, shown in figure 3, is capable of performing optical pattern recognition. The correlation is performed by multiplying thc Fourier transforms of the input scene and a refcrencc pattern, followed by an inverse Fourier transform. Lcnscs arc used to perform thc optical Fourier transforms. In a coherent optical processor, the complex conjugate of the Fourier transformed reference pattern is stored on a transparency. Sincc thc transparency containing the refcrericc pattern is placed in the Fourier planc, it is a spatial filtcr. Thc storcd rcfercncc pattcrn is a cornplcx clcct,ric ficltl ;~ntlis sirnilar t,o i\ I~olograrn. 111[,his t,t!xb, ;\ spatial filtcr will bc any transparency or wavcforrn modulating device placed in the filter plane of an optical correlat,or, or sirn~~lat,cd opt.ical corrclator systtrn.

fact that the spectrum of a translated signal is only modified by a linear phase factor. The phase factor contains the necessary information t o image the signal probability of occurrence a t the proper position in the output plane. However, the system is sensitive to the relative orientation of the input signal with respect to the filter. The input scene can be scanned for all possible orientations of the signal by rotating the filter to all possible orientations. More formally, the input scene transparency placed at plane PI is Fourier transformed by lens L2. Let the input scene transparency contain the signal function s(x,y). The Fourier transform of the input scene transparency S(p,q) is defined as:

The electric field, S(p,q), is imaged onto a transparency containin a matched spatial filter at plane P2, which is also calle% the filter or Fourier plane. The MSF contains F*(p,q) and (unavoidably) some other terms. Propagating an electric field through a transparency is equivalent to multiplying the complex transmittance function contained in the transparency, F*(p,q), by the electric field passing through the transparency. The multiplication operation performs the spatial filtering. The filtered signal S p,q).F*(p,q) is then inverse Fourier transformed by lens 3. The resultant electric field is observed a t the output plane P3, and is given by:

When S(p,q is located in the input scene and the hlSF contains F* p,q) = S*(p,q), we see that the electric field at plane P2 is pure real, since the conjugate phases cancel. Thc cancellation of the phase gives a uniform illumination a t plane P2 (a planar wavefront), which when Fourier transformed by Lg,gives a dclta function rcsponsc. This brigllt spot occurs on the o ~ ~ t p plane ut corresponding to the location of s(x,y) in the input transparency. Thc MSF that we have assumed here is a sirnplificd vcrsion of the actual casc. In reality, tile complcx waveform F*(p,q) can only bc recorded on a transparency I)y 11si11g 11ifl1 frc(ll~e~~cy r carrier. l'hc 11igl1l r e q ~ l c ~ ~ c y cilrricr is easi y obtaincd usin a planar off-axis reference beam.14 Using a re erence beam to record the complex waveform incrcascs the spaccbandwidth rcq~iircmcnts the I)l~otographic-recording on 1.ransparcncy. Thc interfcromctrically recorded MSF actually contains:

When thc t.ranSl)i\rcn(:y conl,aining the input sccnc is placcd a t thc input ~ ) l i ~ olcthc cohcrcnt opttical n correlator, thc output, plane contains the li lit intensity distribution from thc convolution of the re crcncc pattcrn with the input sccnc. If the input scene transparency contains the reference pattern, the output of the system essentially rcprescnts the probability that thc desired signal has occurred in the input scene. A bright spot in the output planc indicatcs a high probability of occurrence a t that position, a dim spot indicates a low probability. This process simultancously detects signals with s i m i l a ~ ~ e n t a t i o regardless of their position in ns the input scene. This space invariance results from the

where R p, ) represents the an off-axis reference electric field, anb ~ l p , q represents the filter pattern spectrum. ) The equation above contains three separable terms. The first two terms together make up a dc spot. The third term contains the filter spectrum multiplied by a linear

phase factor. The fourth term, which is the term of interest, represents the complex conjugate of the filter spectrum multiplied by the another linear phase factor. The output electric field obtained on plane P3 is:

which yields a t the output plane of the optical correlator:

+ roe [s*f*d(x-b,

y-c)]

A spatial filter obtained by combining a MSF and an amplitude transmittance filter in this manner will aid in revealing the properties of the POF. The Fourier spectra of finite, real valued objects fall off with increasing frequency. The amplitude only filter preceding the MSF will have a reciprocal response behaving similar t o a high pass filter. High pass filtering is a standard image processing technique t o emphasize edge information. The correlation peaks obtained using POFs in optical correlators are very sharp, since the high spatial frequency information decorrelates very quickly. This high frequency emphasis also makes t h e W F more sensitive to scale and rotation changes than a conventional MSF. Although the P O F can produce sharper correlation peaks with lower sidelobes, the MSF is still the optimum filter in the presence of noise.

where b and c are constants determining the reference beam incidence angle. In equation 5, the first term, which contains two parts, represents a dc spot that occurs on-axis (x=O,y=O) in the output plane. The Fourier transform of the collimated reference beam R(p,q) is a point r,. 'l'he second term represents the convolution of the filter function with the input scene (s*f). This term appears off-axis at position (x=b,y=c). The third term, which is the one of interest, is exactly the expected output from the optimum filter. This term represents the cross correlation of the input scene with the filter function (s*f) and appears on the output plane at (x=-b,y=-c). When s=f a maxima occurs on the output plane. Phase Onlv And Binary Phase Onlv Filters Since Vander Lugt's original introduction of MSFs there has been considerable effort devoted to the study of optical correlators. Analysis of the Fourier transform of continuous tone pictures has revealed that the phase information is considerably more important than the amplitude information.15 Classical Vander Lugt MSFs use both amplitude and phase information. The amplitude portion of the MSF alters the amplitude transmission of the spatial filter. This amplitude modulation decreases the li ht transmitting efficiency of the MSF. A purely Phase 8 n l y Filter (POF) has 100% transmission efficiency. The incident light is not attenuated, but the phase is modulated across the plane. Ilencc, POFs may offer an alternative to classical MSFs for opt.ica.I corrola1,or systcn~s rcquirir~g hifill triu~slnissio~~ offici(:~~cy. Srlc:ccsslr~lrc:sr~ll.sr~sir~g 1:s i l l ol)l,ical c:orrd;~l,or 1'0 experiments have been reported.lG~17,18The performance of a P O F ca.n be dcscrihed mat,liematically by using a con~bir~al,ion a classical MSI: ant1 an aml)lit~~tle ol l , r ; \ , r ~ s ~ r ~ ifilt,(>r. l,(!t, ~,~I C MSl? r(!s~)ollsofr111(:t,io11 t,t,i~~~ I ~ I)(! written in polar lorrn:

Phase only filters can be generated by bleaching a previously developed photographic plate, or by using a plate with a dichromated gelatin emulsion. However, restricting the continuous phase filter function to binary values offers several advantages. A binary process is easy to control and the two phase levels can be set very accurately. Another advantage of using a binary phase only filter is the elimination of the dc spot in the output of an optical correlator system. These binary filters can be digitally generated on a computer, which eliminates the complex optical preprocessing required in hlSF synthesis, or in the synthesis of continuous POFs. Binary Phase Only Filters (BPOF) are generated in a two step process. First the discrete Fourier transform F(p,q) of the function to be searched for, f(x,y), is calculated. Then the BPOF, which will be referred to as FBpoF(p,q)is made by setting the phase at each point (p,q) t o one of two levels, depending on the sign of the cosine of the phase an le in F(p,q). A technique proposed by Psaltisls for generating BPOFs is:
+l i f Re [ F(p,q)] > 0

FBPOF(P,~) =

-1

otherwise

The Psaltis technique is equivalent to taking the cosine transform FC(p,q) of f(x,y) and then determining FBWF(P,~) by: +1 i f Fc(p,q) > 0 FBWF(P,~) = -1 otherwise

'I'~I(! [)(!rfor~~~:i,~~(:(! 1~1'01~a11 opl,i~i~.l of 1,11(! ~ I I (,ort-(!li~l,~r (:;I.III)e tlclnonstra~l,etIOy col~sitlerillgt l ~ c cot~t.in~~orls cosinc r.riu~sforlr~ f(x,y). 'I'he cosinc transforrn can I)e written: or

where for real f(x,y):

where I F(p,q!l is the magnitude of F(p,q), and Hp,q) is the phase of E p,q). A P O F can be obtained by inserting an amplitude transmittance filter containing the inverse amplitude of the MSF (11 I F(p,q) I ) in front of the MSF. When these two filters are placed in series, the two amplitudes cancel, producing unity amplitude transmittance. The phase only response is then:

where S { )denotes the Fourier transform operation and * denotes complex conjugation. Equation S shows that whenever Fdp,q) is used in an optical correlator, the input scene will be correlated with both f(x,y) and f(-x,-y).

The optically recorded MSF had to be recorded using a spatial carrier in order to get the terms of interest off-axis, away from the dc spot. Since no spatial carrier is required for separating the outputs, the space-bandwidth product of the recording medium is substantially reduced. So we see that the BPOF yields two terms off axis, without the use of a spatial carrier, and without any optical processing. The BPOFs discussed above do not always perform well for cases where the input scene contains large amounts of clutter or objects. The BPOFs can cause the objects to interfere with each other, which causes a variation in the correlation peak amplitudes on the output plane. An alternative to the BPOF is the Hartley-BPOF, or HBPOF. The HBPOF is constructed in a similar two step process as the BPOF, except that the I-Iartlcy transform is used instead of the Fourier transform. The Hartley transform can be defined for real valued input functions a s 9

incident over the two-dimensional area of the SIGHT-MOD can be locally modulated. Dcpending on the orientation of the output polarizer, both binary amplitude modulation, or binary phase modulation can be obtained. The binary phase modulation operational mode permits the construction of BPOFs on the SIGHT-MOD. The SIGHT-MOD is electrically addressable through a computer interface. The device is capable of displaying several frames per second. The largest matrix size available is a 128x128 elcment device. As pointed out.before, BPOFs do not require high space-bandwidth products in order to obtain good performance. Through the development of these high speed, electrically addressable, binary phase modulators, coherent optical correlator systems can be constructed which cmploy BPOFs at the filter plane.19,20,28-40 Such a coherent optical correlator system is capable of generating the interrogation patterns inside the computer, and then displaying the BPOFs on the spatial light modulator. The input scene is interrogated by displaying successive filter response functions on the light modulator. Hence, the true potential of pattern recognition systems sensitive to object shape and rotation can be obtained using 2-D spatial light modulators employing BPOFs for signal detection. Ontical Correlations For Particle Ima.~inrr Velocimetrv

where F(p,q) is the Fourier transform of the real valued filter function f(x,y). The HBPOF is generated by:
$1 i f H(p,q) 2 0

FIBPOF(P,~) =

-1

otherwise

The IIBPOFs seems to perform very much the same as the BPOFs, with the exception of having much lower variation in correlation peak amplitude.20 Snatial Liaht Modulators

A coherent optical correlator requires a spatial filter for proper operation. A matched spatial filter is the most common type of filter used in an optical correlator. Making the MSFs requires much care, preparation, and trial and error. Yet each MSF is capable of only recognizing a single object pattern. Processing the input scene using optically recorded MSFs requires many transparencies. The BPOFs can be synthesized entirely on a digital computer. The exact object function is not required to be Fourier transformable since the discrete Fourier transform contains sufficient information for generating the BPOFs. The cosine of the phase of the original Fouricr transformcd function is analyzcd a t cach point, and assigned one of two discrete phasc values. The FlPOFs arc more scnsitivc to scale and rotation changcs than conventional continuous, complcx valucd MSFs.16 'L'hcsc propcrl,ics alorlc (lo not rnakc the usc of 1lI'OT;'s in rohcrent optical processors attractive. lIowcvcr, with the advcnt of high spcctl two-dimensional spatial light motlulators, RPOFs can bc readily synt,hcsizctl in the c.or~ll)r~l,cr cncotlctl o111,o ligltl, n~otl~rl;tl,or I I S i~l l a~ltl I,hc for coherent opt,ical corrclator syslcms.
'I'llcrc arc s~vcralI,y~)es spatial light 111otluli~1,ors of av,zilal)lc. Sornc o l tl~csc optically atltlrcssal)lc ant1 arc use liquid crystal light valves.21-23 Other elcctronically addressable liquid crystal displays arc used in small portable tclevisions.24-26 Semetex corporation manufactures a SIGHT-MOD magnetooptic spatial light modulator.27 The SIGHT-MOD uses the Faraday effect to selectively rotate the polarization of thc transmitted linearly polarizcd coherent light. The SIGIIT-MOD consists of a matrix of small elements, or pixels. Each pixel can operate independently, so that the light

Optical correlator systems have been employed to identify a wide class of objects including: alpha numeric characters, assembly line parts, human faces, military aircraft, ships, tanks, and multiple exposure particle image records.l6,18-20,41142 In our application, an optical correlator will be used to identify recorded particle images from a multiple-pulse laser light sheet technique. The pulsed laser sheet technique, also called Particle Irnaging Velocimetry (PIV), records the in-plane particle positions a t the instants the light sheet is pulsed. Typically, two pulsed light sheets are used, yielding two records of each particle in the light sheet. The individual particle displacements between exposures is proportional to the in-plane velocity component of the particle. This two pulse recording scheme results in a 160' directional ambiguity when traditional point-wise Young's fringe analysis or autocorrelation techniques are used to reduce the data. When the photographic transparency containing the recorded particle images is placed a t the input plane of a coherent optical correlator system, the cntire pllotograph can be proccsscd in parallel. In addition t o speeding up the data reduction process, the opl,ical corrclator can bc used to determine the vclocity vcctor directions unambiguously. The directional ambiguity is removed by using the appropriate pulsing scheme.

A pulsing scqncnce yieldin a non-symmetric pattern ci1.11 trsrvl 1,o( : l i ~ t ~ i ~ ~,II(! firodoll;Ll il~t~l)ig~ril.y. I)(! ~i~k: '1'11~ ~)roblcrnis t,o dctcrminc thc optimum pulsing scqucnce which givcs good signal discrimination on the output plane, unambiguous velocity direction information, and a low number of false correlations. A triple pulse light sheet scheme can be used with non-symmetric time periods between the pulses. The asymmetry enables the particle direction to be unambiguously determined. For example, let Atl2 = 2.At23 where Atlz is the time interval is the time interval between pulses 1 and 2, and bctwecn pulses 2 and 3. The corresponding recorded

particle image sequence for a left to right flow direction is:

This exposure sequence yields unambiguous direction information. However, this is not the traditional exposure technique used in PIV and LSV approaches. Typically, the time intervals are equally spaced. The asymmetric spacing in the recorded particle image triplet above produces a fringe pattern with 3 frequency components. The phase information from the 3-pulse pattcrn contains the direction information. The traditional PIV and LSV data reduction techniqncs rlsc the squared modulus of the fringe pattern to obtain the velocity information. The squared modulus of the fringe pattern retains no phase information. Hence, the conventional PIV and LSV techniques are unable to make use of the 3-pulse encoded phase information for directional ambiguity removal. A coherent optical processor is able to use the phase information from the 3-pulse pattern to determine the velocity direction unambiguously. Another recording scheme could consist of a 5-pulse pattern. The ratio of the spacings is 2:2:1:1. The recorded pattern is shown below for a left to right flow:

computer, although the simulation is computationally intensive. Two simulation programs were written. One computer simulation was written on a high performance personal computer (pc) to analyze the output plane light intensity distribution of the optical correlator. The second simulation program was written to simulate the sequential scanning of the input scene for successive velocity vectors. This second simulation program was executed on a CRAY-XMP supercomputer. The 2-D velocity vector maps were obtained from the CRAY simulations. . The programs can be easily modified to sclect diffcrcnt types of spatial filtcrs (MSF, POF, BPOF, and IIBPOF) and/or different types of interrogation patterns (3-pulse . The MSFs used in the computer simulations Or are define as the complex conju ate of the filter pattern spectra. The ease in modifying t e program permits the analysis of different filtering techniques without the trial and error process of optically recording the filters. Thc optical Fourier transforms performed by lenses in the optical correlator are obtained in the simulation by using Discrete Fourier Transforms (DFT). The DFT closely approximates the actual Fourier transform. The Fast Fourier Transform (FFT) is a computationally efficient DFT and was used in the computer simulations. The pc simulation program will be discussed first to gain an understanding of the expected output light intensity distribution from the optical correlator. Then the CRAY simulations will be discussed along with the program output 2-D velocity field map.

""""2

fl

The 5-pulse pattern can be generated by shuttering a continuous pulsed laser source. The 3- and 5-pulse exposure sequence particle image photographs can be reduced using an optical correlator. Individual particle displacement records for each particular velocity magnitude and flow direction are used to make the spatial filters in the PIV optical correlator. Spatial filters are made to cover the expected range of velocity magnitude and flow direction. Using conventional MSF techniques would require making individual MSFs for each expected velocity vector on the input scene. This could easily escalate into an unmanageable number of MSFs to store and manipulate. Sequentially positioning large numbers of MSFs into the frequency plane of the optical correlator would be cxtrcmcly tedious and time consuming. Alternatively, an electronically addressable spatial light modulator can be used to display MSFs in the frcquency plane of an optical correlator. Tlle relatively low spatial frequency content of the Fourier spcctra of the 3- and 5-pulse exposure sequences arc wcll suited for the BPOF tcchniqlrc. The BPOFs can be r~scdinstead of MSFs since wc do not, have a detailed signal ncstctl in a cornplcx backgrountl. Our noise is the random location of intlcpcntlcr~t, particle imagc rocortls on t11c i n l ) ~ ~sccnc t, trn~~sl).zrcat:y. IlI'OFs call 1)csy~lt,l~c~si;zc~tl1 1 ~ 'I'lrc ill 1. col~~putcr witlrout illly optical proccssin autl ikrc ~ci\dily encded onto a magnetooptic spatial l i g k mmhlator thror~ h a cornputcr interface. A spatial ligl~t, rnotl~~lal,or placei in the frequency plane of the PIV optical correlator can sequentially display the BPOFs which have been previously generated and stored in the computer. A CCD array camera placed at the output plane in the optical correlator digitizes the light intensity distribution for analysis in the computer. Computer Simulation of O~tical Correlators The optical correlation system discussed for analyzing PIV photographs can easily be simulated on a digital

The pc simulation of the optical correlator system uses a 64x64 element input scene, zero padded43 to 128x12s elements. The zero padding eliminates the wraparound error, but also increases the CPU time. A flow sheet for the simulation is given in figure 4a. The program queries the user for the number of particle image recordings in the input scene. The data are entered in the same format regardless of whether the recording scheme is the 3- or 5-pulse. The data for each particle image recording on the input transparency are entered as: the (x,y) coordinates of the mean position of the recorded image on the input scene; S1-the spacing between images 1 and 2 (or the spacing between exposures 1, 2, and 3 for the 5-pulse case); S2-the spacing between images 2 and 3(or the spacing between exposures 3,4, and 5 for the 5-pulse case); and the angular orientation of the particle image record, 8. For the triple pulse case, the program determines the coordinates of each particle image and places a value of unity at each location in the input scene array, see figure 5. The total power from each particle image triplet is 3.0. For the 5-pulse case, the program again places a value of unity at each particle image position, scc figure 6. The energy content for the 5-pulse pattcrn is 5.0. After all of the input scene data havc been cntcrcd, the 2-D FFT of the input scene is cornputctl. Ncxt, the spatial filter pattcrn is cornputctl. 'I'llc pararnc!t.crs tlcscribing thc spatial liltcr arc st.orctl as constants within the program. The coordinates of the particle image positions in the filter pattern array are cornpr~tctl, value of unity is placed at each position a within the array, and the filter pattern is Fourier transformed. The complex conjugate of the Fourier transformed filter pattern is the MSF, and can be stored in a complex array in the simulation. Recall that the optical recording of an MSF required an off-axis reference beam to record the complex waveform. The MSFs are generated and stored very easily in the computer simulation. At this point in the program, the spatial filter can be modified to some type other than the classical MSF. When optically recording a MSF for use in an actual

optical correlator, one tries to maximize the transmittance of the MSF. The higher the transmittance, the more optically efficient the MSF. The MSF optical transmittance calculated in the computer simulation is maximized by normalizing the computed MSF such that the maximum transmittance is unity. The MSF is scanned and the maximum transmittance value, F,,,, is determined. All points in the MSF are divided by the maximum amplitude transmittance of the = filter, Fhs~(p,q) Ft(p,q)/FmjlX. The resulting continuous valued MSF conta~ns transmittance values from 0 to 1. MSFs are never 100% optically efficient, but maximizing the transmittance results in more power a t the output plane. The POFs are produced by computing tllc MST: antl then normalizing each point in the filter by the absolute magnitude at that point, FpOF(plq)= F*(p,q)/ I F(p,q) I. The magnitudc normalization produces a filtcr with unity magnitude a t all points, which gives 100% optical efficiency. The Binary Phase Only Filters (BPOF) are generated by examining the sign of the cosine of the phase of the MSF a t each point and assigning a value of +1 or -1. The resulting filter pattern is both real and symmetric, and has no attenuation since it has unity magnitude at all points. The real and symmetric properties of the BPOF are similar t o the properties of a cosine transform of the filter function.19 The BPOF contains both the filter function spectrum and the complex conjugate of the filter function spectrum. After the appropriate filter spectrum has been generated, a point-by-point complex multiplication is then performed on the F F T of the input scene and the complex conjugate of the filter spectrum. This simulates the spatial filtering operation which occurs a t the frequency plane in an optical correlator. The filtered signal is then inverse FFT'd back to the space domain. The resulting image is the output electric field from the simulated correlator. The electric field is a complex quantity which cannot be directly measured, but the intensity of the image can be measured. The output image intensity is calculated by multiplying the output electric field by its complex conjugate. Only a portion of the output image, which contains 128x128 elements, contains useful information. The input image and spat,ial filter pattern were zero patltlcd to avoid wraparound error. The output plane only information. Tllc contains 64.64 clcrncnts of u s c f ~ ~ l appropriate choice of t,hc inpr~t scene ant1 filtcr pattern within the 128x128 arrays yields the output image within t,hc first 64x64 clcrnont,s of t,he or~t,pr~t, array, see figl~rc 7. The MSF and POI: filters arc complex functions which arc symmetric across the 128x128 array. The protluct of 1,hcsc fl~nct,ior~s I,II(! i11p11t wil.11 S(:CII(! rollowed I)y 1110 inverse F F T operation produces only a 64x61 subarea of 11scf111 it~rornlI,1011. ;I ' 'I'hc situation is sligl~tlytliffcrcnt when 11I'OI:s arc ~ ~ s a tdthe filt,cr plant. 'L'lic 13POF is a reill a11t1 c symmetric function which will give a symmetric result after being inverse FFT'd. The output plane of an optical correlator using a BPOF a t the frequency plane will possess some symmetry. The symmetry manifests itself as an invertcd copy of the filter pattern at the conjugate position of the filter pattern. The inverted copy produces the convolution of the input scene with an invertcd copy of the filter pattern. This information is located in the region 65-128, 65-128) on the output plane, see figure 7. T is subarea of the output plane is usually ignored in tllc computer simulation.

The output plane intensity distribution from the pc computer simulation is displayed in an isometric 3-D bar graph. The 3-D bar graph accurately represents the discrete nature of the data. The plotting routine is written in Fortran 77 using the Halo44 graphics primitives library. These plots enable easy visual evaluation of the output intensity distribution. A sample output for a MSF in the optical correlator is given in figure 8. The graphs typically cover the range x:1-64, and y:1-64, which contains the desired output information. . The total energy output for the MSF is approximately half of the input plane energy. The optical efficiency of the MSF is not 100%, so some attenuation is expected. Thc POFs give 100% optical efficiency, therefore, the output plane energy is identical to the input plane energy. The BPOFs are also 100% efficient, however, the output plane information is spread over the two subareas (1-64, 1-64) and (65-128, 65-128 . The total energy is conserved over the entire output p ane, approximately equally divided between the two subareas.

The optical correlator simulation on the CRAY-XhIP computer is very similar to the pc simulation. The flow chart for the simulation is given in figure 4b. The CRAY simulation uses 256x256 element arrays, which contain 128x128 points of useful information. The input scene data for the CRAY simulation is obtained by multiple calls to the subroutine used for generating the particle image record. After all the input particle images have been written to the input scene array, the input scene is Fourier transformed. The Fourier transformed input scene is stored in a complex array. The simulation now enters a loop where the filter pattern parameters are sequentially varied. The input scene data must be searched over a range of vector magnitudes and angles. There are three loops used to generate the filter patterns. The first loop selects the spacing S1. The second loop picks the spacing S2. The value of S2 ranges from S1/2 to S1/2 - 1. This accounts for any pixel rounding off errors in the digitized PIV transparency. The third loop selects the angle range on 8. Due to the limited pixel resolution, the number of unique velocity vectors which must be searched varies with the magnitude of the vclocity vector. For small interrogation vector magnitudes, the angle increment between unique vectors is large. A unique vector is defined as a vector which occupies a t least one new pixel location comparetl to the previous vector. For a given vclocity vcctor magnitutle, the program increments the angle 8 by lo, and calculates the coordina1,cs of tllc particle images in thc filtcr pattcrn. If the ncw coortlinatcs are not different from the oltl coortlinatcs tllc angle is incremcnted again. This process is rrpcabctl ~rnl.ila 11nit111c particlc image rccortl is fonnd. 'I'hc I I I I ~ ~ I Ifilt,(!r 1);1ttcrnscovering ~ I I Cra.ngc of C expected vclocity magnitudes and angles are sequentially gcl~cral.ctl,Ipolrricr I,ransformctl, point-wise multiplied by 1,11(! I:ol~rit!r l,r;~~~sfor~r~etl ,C C I I ~ ~II~IIS I arra,y, a.11tI1,11c11 1,Ilc filtcrctl signal is inverse Fourier transformed. This procctlurc of gcncratir~gfiltcr patterns antl taking tllc product of thc input scene spectrum and the filtcr spectrum simulates the sequential display of BPOFs on a light valve in an actual coherent optical processor. For each filter pattern generated, the output plane intensity distribution is calculated and scanned for correlation peaks. Each point in the output plane is compared to a threshold level. When a peak exceeds the threshold level, threc quantities are stored in three independent arrays. The elements of the arrays correspond to each pixel on the output plane in the simulated correlator. The current filtcr a.ngle and magnitudc settings are stored in

r,

two different arrays corresponding to the indices of the detected correlation peak. In addition, the correlation peak amplitude is stored at the detected peak position in another array. Before a detected peak's parameters are stored, the current peak amplitude is compared with the value stored in the peak amplitude array. If the new value exceeds the stored value, the new data are stored. If the stored peak value is greater than the detected peak value the stored values are preserved. This comparison procedure ensures that the most probable velocity vector at a particular location is stored, where the probability is indicated by the correlation peak amplitude. After all the filter patterns have been used, the vector magnitude and angle st.ora e arrays contain the tlctected 2-D velocity vector map o the input scene. The arrays are written to a file and downloaded to a pc computer for graphing.

The diagonal terms are the matched terms, which cancel to yield the dc portion of the filtered signal. This dc component will become the correlation peak on Fourier transformation to the output plane. The off-diagonal terms are the sidelobe information. The filtered signal can be rewritten as:

Ex~ected u t ~ u Plane Information O t An analysis of the matched filtering operation for a PIV optical correlator will indicate the expected output response and signal to noise ratio of the output. The finite energy in the input signal should limit the correlation peak amplitude. Also, the correlation length of the interrogation signal should be related to the energy in the sidelobes. The calculated results will be compared to the pc simulations. The pc simulation requires 5.5 seconds to compute a 128x128 real to complex single precision FFT. The complex to complex inverse FFT takes 11 seconds to compute. The entire pc simulation program take 23 seconds to perform the optical correlation operation. 3-Pulse Exposure Seauence For the triple pulse exposure scheme, the particle image is recorded at three instances in time. The time intervals between pulses are chosen to make the spacings non-symmetric. Our choice is to make the time interval ratio 2:l. Assuming a noise free background on the input scene, the functional form of the 2-D, 3-pulse signal is:

The output electric field is obtained by inverse Fourier transforming the filtered signal:

The intensity is given by:

however, we will assume no y- dcpcntlcncc to rctlucc the bookkeeping:

Using a definition of SNR as the ratio of the correlation peak power to the sidclobe power, we find SNR = 1.5. The output plane from the computer simulation of the 3-pulse MSF case is shown in figure 8. The calculated SNR is identical to the estimate obtained from the pc simulation. The computer simulation output plane total energy is 1.667. The MSF is approximately 50% cfficicnt.

'I'llc triple pulse 1'01: cnsc sllor~ld yield n larger SNll than the MSF case. Again, let the input scene contain s(x) as defined in equation 12, and define the POF as: 3-Pulse MSF The MSF for this case is just the complex conjugate of the input signal spectrum, F*(p) = S*(p). The filtering operation for the 3-pulse MSF case is:

The filtering opcration for the "matchcd" case is: S(P)'FPOF(P) =

' F # W
=

= I S(P) I

(21)

the SNR = 6.0 will approach the computer simulation value of 5.81 as more higher order terms are included in the approximation. The higher order terms produce more off-axis sidelobes on the output plane. The input scene total power is conserved in this 3-pulse P O F simulation, since the total output energy is 3.0. 3-Pulse B P O F A direct analytical study of the BPOF or HBPOF is not feasible. A 3-pulse BPOF computer simulation result is shown in fi ure 10. The distribution of sidelobes is not in a linear fas ion as in the P O F case. The sidelobes are distributed randomly around the correlation peak. The B P O F SNR is 2.35, which is apprqximately half the value obtained with the continuous POF. This is due to the symmetry of the BPOF. The BPOF spectrum contains both the interrogation pattern spectrum and the conjugate of the interrogation pattern spectrum. Hence, for the matched case (F(p,q) = S(p,q)), the input scene is cross correlated with the filter function S(p).S*(p) and also autocorrelated with itself Sz(p). This means that half of the input scene energy is used in the autocorrelation and the other half is used in the cross correlation. This is the reason for the ~ 5 0 % reduction in the SNR for the B P O F case over the P O F case. The total output plane energy from the sinlulation in figure 10 is approximately 1.5. The total energy in the 64x64 array is only half of the input scene energy. Figure 11 shows the entire output plane of 128x128 points. The autocorrelation term a t the conjugate position is readily observed. Also, the total energy is now accounted for. The B P O F is 100% optically efficient, but distributes the sidelobe energy in a different form than the continuous POF.

= IS(P).S*(P)I~
The optical correlator output electric field is given by:

From the 3-pulse MSF case we know the form of S(p).S*(p), which when inserted into equation 22 gives: L(x) = 1

j( 3 + 2cos[p(x3-x*)] + 2cos[p(x3-x2)]
-m

+ 2 c o ~ [ ~ ( x ~ x l ) l )d*p. e ~ ~ ~

(23)

There is no known closed form solution for this inte ral. The kernel of the transform can be rewritten in the orm (l+x)+ and then approximated by

Fi

(l+x)+

1 1 + zsx - 1

1 + --ax3 - ... 16

(24)

Using the first two terms in the approximation we find:

3-Pulse HBPOF The HBPOF performance is very similar to that of the BPOF. The HBPOFs are phase only filters, so we expect 100% optical efficiency. The HBPOFs also contain the same conjugate symmetry as the BPOF. The HBPOF output will contain the desired filter spectrum plus the complex conjugate of the filter spectrum. The difference between the HBPOF and the BPOF is the sensitivity of the correlation peak amplitude to the filter pattern orientation angle 8. The BPOF and IIBPOF give identical results when 8=0, or 90, as in figure 10. IJowever, as 8 varies between O0 and 45' the B P O F correlation peak amplitude decreases, reaching a minimum a t 6'=45', see figure 12. The energy loss in the correlation peak appears as an energy ain in the sidclobes, which rctluccs thc SNR for t c BPOF. Figure 13 shows the output plane light intensity distribution for the IIBPOF under the identical conditions as figurc 12. The HBPOF correlation pcak valuc and SNR arc csscntially I,l~c same as for the 8=0 case. At 8=1S0, the BPOF corrclation pcak amplitude has decreased 6% over the valuc a t 8=0, whilc the IIBPOF correlatioll pcak amplitude increased by 0.8%. Both filters display equal performance a t 8=0, but as angle 8 increases the B P O F degrades in performance. The HBPOF maintains its pcrformance throughout the filter pattern angle orientation range. This feature is important for an optical correlator system which requires a minimum threshold level for si nal detection. We want the correlation peaks t o e as large as possible to give better noise discrimination. For the class of objects we are tlenling with the IIRPOF gives sr~periorperformance.

where denotes an approximate rcprcsentation ol the actual function. A closed form solution for the Fourier transform in equation 25 exists and the result is:

The intensity is:

The cstirnated SNR using this approximation is SNR = 6.0. The output from the computer simulation for this 3-pulse P O F case, shown in figure 9, gives a SNR = 5.81. The approximation appears t o be pretty close for this case. The information we are interested in is why the P O F gives smaller sidelobes than the MSF. The signal in the P O F case is of the form (l+x):. The square root operation rctluccs the sidelobe arnplitotles relative to the correlation peak amplitude. The estimated valuc of

5-Pnlse Patterns The 5-pulse output plane intensity distributions are very similar to the 3-pulse cases discussed above. The main differences are an increase in both the correlation peak amplitude and the number of sidelobes. An anticipated benefit of the 5-pulse sequence is a reduction in the number of false correlations due t o the random distribution of the recorded particle images. The requirement of 5 particle images being properly oriented is more stringent than for 3 particle images. The differences will be borne out in the CRAY simulation results.

The 2-D velocity vector map from the 5-pulse 2:2:1:1 case is shown in figure 17. There are 79 total vectors detected in this output map. There are 74 correctly identified vectors, 5 missed vectors, and 5 falsely identified vectors. These results are very good, 87% of the input data vectors have been correctly identified. The lower number of falsely identified vectors for the 5-pulse case is attributed to the lower probability of the random particle locations satisfying the required 5 particle image positions in the pattern. Of the 74 correctly identified velocity vectors, there were no errors in determining the velocity magnitude. The HBPOFs are very sensitive to scale and rotation. The detected angles did not exactly agree with the input data. This is not a detection error, but is a quantization error. The limited resolution of the arrays used reduces the angular resolution. The program skips some angles because they do not produce unique velocity vector patterns. The variation between the detected angle and the actual angle occurs when the actual angle lies between two unique velocity vectors as determined by the computer. This is an inherent limitation due to the discrete nature of the process. Hi her resolution arrays will decrease these quantization ehects. Conclusions The results of a computer simulation of optical correlation techniques applied to the reduction of particle imaging velocimetry transparencies has been presented. An asymmetric pulse exposure scheme allows the velocity vector direction to be unambiguously determined. The 5-pulse exposure sequence yielded a high quality 2-D velocity map with few false correlations. Hartley binary phase only filters were used to simulate the operation of a spatial light modulator placed a t the frequency plane in the optical correlator. Although the simulations were intended to model the operation of an optical correlator employing a spatial light modulator, the digital evaluation of the transparencies is also a viable technique. For the digital evaluation, the input scene transparency can be digitized and supplied t o the comp~ltcrprogram. The required processing time varies directly with the range of velocity vector magnitudes and angles which must be searched. Future work will involve optically recording a 5-pulse exposure PIV transparency and processing the transparency using both the optical correlator employin a spatial light modulator, and a digital analysis on t e computer.

CRAY Simulation Results The output data from the CRAY simulation program is a 2-D velocity vector map of the input particle image transparency. The input scene chosen for the simulation is a four quadrant swirling flow. The velocity vector map for this flow field is shown in figure 14. The four swirl pattern was enerated using a solid body rotation in each e quadrant. ~ f velocity varies with radius. The vector positions were determined using a random number generator. The input scene contains 79 velocity vectors in a 128x128 array. The angle range encompasses a full 360'. The particle image separations (Sl) range from 2 pixels to 12 pixels. The input scene contains the images of the particles in the swirling flow a t the instances in time that the light sheet is pulsed. The four swirl input scene for a 3-pulse exposure is shown in figure 15. This is the scene that the optical correlator must scan and identify the images correspondin to a unique particle. A few of the outliers are easy t o i entify. However, closer to the center of the swirls, the three images of the same particle are more difficult to identify.

The results obtained from the pc simulations were used to select the HBPOFs and to set the minimum threshold level for a tlctcction. For the 3-pulse case t,hc threshold level is 0.97. For the 5-pulse case the threshold level is 1.35. These are the minimum threshold levels. The detected correlation peaks may be larger than these minimum levels. The source of variation in the correlation peak amplitude is from the overlapping sidelobes from other partial correlations. Another source of correlation peak variation is the angular dependence of the correlation peak. This angular variation is minimizetl by using the IIBPOFs. From t,he rcq~~irctl searching ranges on S1 (10 valr~cs), S2 (2 values for each Sl), and 8 (360 angles for each S1 and S2), the numbcr of velocity vectors to be searched is 7200. The numbcr of unique velocity vcctors tlctcrmincd ant1 scarchctl by the simr~lationprogram is 1221. Thc l,Ol,al c l ' u r~ll~I,ilrl~! l)Olth 1,110:!- ikll(1 !~-])II~sc fOI' simrllations is 229 sccontls. This corrcspo~~tls to approximately 0.18 seconds per interrogation pattcrn. The results from the 3-pulse 2:l ratio spacing pattern are shown in figure 16. This 2-D velocity map contains many false and missed correlations. There are 48 correctly identified vcctors, 32 missed vectors, and 69 falsely identified vectors. These results are unacceptable. The random location of the particle images results in many false correlations for the 3-pulse pattern.

References
1)

Lccndcrtz, J . A.;"lnterferomctric displacement measurement on scattering surfaces utilizing spcckle effect", J. Phys. E 3, p214, 1970. Incichcn, B., Elgin, P., Dandliker, R.;"Hybrid optical a.nd clcctronic image processing for s1,rain measurements by speckle photography", Appl. Opt. 19, ~ 2 1 9 1 ,1980. Barker, D. B., Fourney, M. E.;I1Measuring fluid velocities with speckle patterns", Optics Lttrs 1, pp 135-137, 1977. Dudderar, T. D., Simpkins. P. G.;I1Laser s ~ ~ e c k l e in a fluid mediumll, haturc 2'70, pp 45-47, 1977.

2)

3)

4)
'

5)

Meynart, R.;I1Equal velocity frin es in a Raylei h-Bernard flow by a specile method", Appl. b p t . 19, p1385, 1980. Meynart, R.;ItSpeckle velocimetry study of vortex pairing in a low-Re unexcited jet", Phys. Fluids 26, p2074, 1983. Meynart, R.;"Instantaneous velocity field measurements in unsteady gas flow by speckle velocimetry", Appl..Opt. 22, p535, 1983. Adrian, R. J.; Yao, C. S.;"Development of pulsed laser velocimetry (PLY) for measurements of fluid flow" Eighth Biennial Symposium On Turbulence, Sept. 26-28 1983, U. Missouri Rolla, J. L. Kakin and G. K. Patterson, Eds. Adrian, R. J., Yao, C. S.;"Pulsed laser technique application t o liquid and gaseous flows and the scattering power of sccd materials", Appl. Opt. 24, p44, 1985.

Mallory, W. R.;I1Dynamic Spatial Filter for Optical Signal Processin Using a Liquid Crystal Light Valve", Intern. oc. for optical Engineering, A83-22794 08-60, Aug 25, 1981. Michaelson, J. D., Sawchuk, A. A.;I1Non-Linear Optical Processing Using Liquid Crystal Light Valves", Soc. of Photo-Optical Instrum. Engineers, A81-32464 14-35, 1980. Mok, F, Diep, J., Liu, H. K., Psaltis, D.;IoReal-Time Computer Generated Hologram by Means of Liquid Crystal Television Spatial Light Modulator", Opt. Lttrs, vol 11, Nov. 1986. Liu, H. K., Davis, J. A., Lilly, R. A.;"Optical-Data-Processing Properties of a Liquid Crystal Television Spatial Light Modulator", Opt. Lttrs, vol 10, Dec. 1985. Tai, A. M.;I1Low-Cost LCD Spatial Lilght Modulator With High Optical Quality , Appl. Opt., vol 25, 1986. Ross, W. E., Snapp, I<. M., Anderson, R. H.;"Fundamental Characteristics of the Litton Iron Garnet Magnetooptic Spatial Light Modulator", Proc. Soc. Photo-Engineers 388 55-64, 1983. Flannery, D., Keller, P., Cartwright, S., Loomis, J.;"Characterization of Improved Binary Phase-Only Filters in a Real-Time Coherent Optical Correlation", SPIE vol 753, Acousto-Optic, Electro-Optic, and hfagneto- Optic Devices and Applicat, 1987. Cottrell, D. M., Davis, J. A., Schamschula, M. P., Lilly, R. A.;"Multiplexing Capabilities of the Binary Phase-Only Filter", Appl. Opt., vol 26, pg 934, 1987. Gregory, D. A.;I1Real-Time Pattern Recognition Using a Modified Liquid Crystal Television in a Coherent Optical Correlator", Appl. Opt., vol 25, pg 467-6 Keller, P., Flannery, D., Cartwright, S., Loomis, J.;"Performance of Binary Phase-Only Correlation on Machine Vision Imagery", SPIE vol 728, Optics, Illumination, and Image Sensing lor hlachine Vision, 1986. Davis, J. A., Day, T., Lilly, R. A., Liu, Ii. I<.;"Multichannel Optical Correlator/Convolver Utilizing the Magnetooptic Spatial Light Modulator", Appl. Opt., vol 26, pg 2479-83, 1987. Davis. J. A., Ileissenberger, G. M., Lilly, R. A., Cottrell, D. M., Brownell, M. F.;"High Efficiency Optical Reconstruction off Binary Phase-Only Filtersllughes Liquid Crystal Light Valve", Appl. Opl., vol 26, pg 929-33, 1087. Flannery, D. L., Biernacki, A. M., Loomis, J. S., Cartwright, S. L.;"Real-Time Coherent Correlator Using Blnary Magnetooptic Spatial Light Modulators atand Fourier PLanes", Appl. Opt., vol 25, pg 466, 1986. Cottrell, D. MI Davis. J. A., Schamschula, M. P., Lilly, R. A.;"Multiplexing Capabilities of the Binary Phase-Only Filter, Appl. Opt., vol 26, pg 934-37

6)

7)

8)

9)

10) Yao, C. S., Adrian, R. J.;"Orthogonal compression and I-D analysis technique for measurement of 2-D displacements in pulsed laser velocimetry", Appl. Opt. 23, pl687, 1984. 11) Adrian, R. J.;I1Image shiftin technique to resolve directional ambiguity in dou~lf+pulsed velocimetry", Appl. Opt. 25, p3855, 1986. 12) Landreth, C. C., Adrian, R. J., Yao, C. S.;"Double pulsed particle ima e velocimeter with directional resolution for comp ex flows" Tenth Biennial Symposium On Turbulence, Sept. 22-24, 1986, U. Missouri Rolla.

13) Coupland, J. M., Pickering, C. J. D., Halliwell, N. A.;I1Particle image velocimetry: theory of directional ambiguity removal using holographic image separation", Appl. Opt. 26, p1576, 1987. 14) Goodman, J. W.; Introduction to Fourier Optics, McGraw-Hill Book Pub., New York, 1968. 15) Oppenheim, A. V., Lim, J. S., Proc. IEEE 69, pg 529, 1981. 16) IIorner, J. L., Gianino, P. D.;"Phase Only Matchetl Filtering", Appl. Opt, vol 23, pg 812-16, March 1984. 17) Gianino, P. D., IIorner, J. L.;"Atlditional Properties of the Phase-Only Correlation Filter", Opt. Eng., vol 23, pg695-97, December 1984. tion 18) IIorner, J. L., Leger, J. R.;I1Paltern Recogl~i With Binary P h a s c O n l y Filters", Appl. Opt., vol 24, pgG09-1 I, March 1985.
19)

l'sall.is, D., I'nck, E. C., Vcnkatcsli, S. S.;"Ol)l.icnl Irnage Corrcl;~tiol~ Will1 a IJinary S]);ll,inl 1,igllt Modulator", Opt. Eng., vol 23, pg 698, 19S4. Cottrell, D. M., Lilly, R. A., Davis, J. A., Day, T.;"Optical Correlator Performance of Binary Phase-Only Filters Using Fourier and Martley Transforms", Appl. Opt, vol 26, pg 3755, 1987.

20)

21) Jost, S. R.;I1A New Light-Modulated Liquid Crystal Light Valve"; J. of Appl. Phys., vol 49, October, 1978.

36)

Mills, S. A., Ross, W. E.;I1Dynamic Magneto-optic Correlator: Real Time Operation", SPIE vol 753, Acousto-Optic, Electro-Optic, and Magneto-Optic Devices and Applications, pg 54-63, 1987.

40)

Wernet, M. P., Edwards, R. V.;"Real Time Optical Correlator Using a Magnetooptic Device Applied to Particle Imaging Velocimetry", Appl. Opt., vol 27, pg 813-15, 1988. Vander Lugt, A.;"Signal Detection by Complex Spatial Filtering", IEEE Trans. Inform. Theo y, VOI IT-10, pp139-145, 1964. Horner, J. L., Bartlet, H. O.;"Two-Bit Correlation", Appl. Opt., vol 24, pg 2889, 1985. + Gonzalez, R. C., Wintz, P.;Digital Image Publishing Co., 1977. Processing, Addiso:-Wesley Halo Graphics I<ernel and Device Drivers Functional Description Manual, Version 2.26a, Media Cybernetics, Inc, 1986.

37) Javidi, B., I<uo, C. J.;"Joint Transform 1ma.g~ Correlation Using a Binary Spatial Light Modulator a t the Fourier PLane", Appl. Opt., vol 27, pg 663-65, 38) Duthie, J . G., Upatnieks, J.;"Compact Real-Time Coherent Optical Correlators", Opt. Eng., vol 23, pg 7-11, January/February, 1984. Ross, W. E., Psaltis, D., Anderson, R. H.;"Two-Dimensional hlagnetc-Optic Spatial Light Modulator for Signal Processing", Opt. Eng., vol22, July/Augustl983.

41)

42)
43) 44)

39)

u;F:fD
FIELD LIGHT SHEET 7,

/I

CYLINDRICAL

Figure 2:

Light intensity distribution after Fourier transformation of a speckle fringe pattern. The power spectrum is symmetric. The average particle displacement is the distance (xo,yo) along the line adjoining the central peak with the sidelobes. The velocity vector is indeterminate to *180.

COLLIMATING

INPUT SCENE

SPATIAL FILTER (REFERENCE SPECTRUM)

CCD ARRAY

I I -I POINT SOURCE OF MONOCHROMATIC PLANE LIGHT

FREQUENCY PLANE

IMAGING LENS

OUTPUT PLANE

Optical arrangement of a coherent optical correlator. The PIV transparency is placed at the input plane (PI). The spatial filter is placed at the frequency plane (Pz). The correlated output is observed on plane p3.

CRAV SlMJL(lT1OH FLOU CHART

PC SIMULBTIOH FLOW CHART

D a t a Entry: X.Y.Sl.S2.0

Spectrum w i t h I n p u t Scene

FFT and S t o r e I n p u t Scene Spectrum

I
1
I

I n p u t Scene Array Data E n t r y V i a Subroutine Calls FFT and S t o r e I n p u t on Vector b g n i t u d r s and R n g l r s b d i f y F i l t e r Pattern

+
FFT F i l t e r P a t t e r n F i l t e r e d Signal Search Output P l a n e

( ~ o d i f y i l t e r Pattern F

Spectrum With I n p u t Scene Spectrum

I n u e r s e FFT F i l t e r e d Signal

Generate I s o m e t r i c

A
End Program

Figure 4a:

PC simulation flow chart.

Figure 4b:

CRAY simulation flow chart.

TOTAL ENERGY = 3.0

16

32

48

64

Figure 5:

3-pulse pattern (2:l spacing ratio).

TOTAL ENERGY

5.0

I6

32

48

64

Figure 6:

5-pulse pattern (2:2:1:1 spacing ratio).

On11 PL(ICEIEWT V I T I I H TIE SI)IIL~TIOH 1RRAVS

12'

m
I n p u t Scene

.-F i l l e d subareas i n t h e I n p u t Sccnc and F i l t e r Pattern arrays i n d i c a t e the data l o c a t i o n w i t h i n t h e arrags. The Cross hatched r e g i o n i n t h e Output Plane a r r a y i n d i c a t e s t h e l o c a t i o n -6 the Correlated r e s u l t . 11.0 in6icat.d i s t h e conjugate i u g e l o c a t i o n f o r t h e BPOF and MBPOFI. Conjugate l u g e

HBPOF

Output Plane

Figure 7:

Relative placement of data within the 128x128 arrays used in the computer simulations.

P A AMPLITUDE = 1.000 EK SIDELOBE AMPLITUDE = 0.111 S R = 1.50 N TOTAL ENERGY = 1.667

16

32

48

64

Figure 8:

3-pulse MSF output plane light intensity distribution for 8=0.

P A AMPLITUDE = 2.560 EK SIDELOBE AMPLITUDE = 8.59E-2 S R = 5.81 N TOTAL ENERGY = 3.000

Figure 9:

3-pulse POF output plane light intensity distribution for 0=O0. Note the lower sidelobes compared to the MSF case.

PEAK AMPLITUDE = 1.037 SIDELOBE AMPLITUDE = 5.109E-2 SNR = 2.35 TOTAL ENERGY = 1.479

Figure 10:

3-pulse BPOF output plane light intensity distribution for 6'=0. sidelobes are randomly distributed about the correlation peak.

The

PEAK AMPLITUDE = 1.037 SIDELOBE AMPLITUDE = 5.109E-2 SNR = 0.53 TOTAL ENERGY = 3.000

Figure 11:

3-pulse BPOF total (128x128) output plane light intensity distribution with conjugate image readily observable. The total input scene energy is conserved.

= 0.9758 SIDELOBE AMPLITUDE = 7.44E-2 SNR = 1.86 TOTAL ENERGY = 1.50

PEAK AMPLITUDE

10

32

10

04

Figure 12:

3-pulse BPOF output plane for 6'=45'.

PEAK AMPLITUDE = 1.046 SIDELOBE AMPLITUDE = 0.1163 SNR = 2.30 TOTAL ENERGY = 1.501

Figure 13:

3-pulse IIBPOF output plane light intensity distribution for 6'=45O. Note the more rapid dccay in sidclobe energy than for the BPOF case with b'=45'. Also, the correlation peak amplitude is very close to the value for the BPOF case with 6'=0.

Figure 14:

2-D velocity vector map of input data to CRAY simulation. There are 79 total vectors.

Figure 15:

Pixel plot of the input data to the CRAY simulation for the 3-pulse exposure case. This plot represents the input scene transparency in the PIV optical correlator.

Figure 16:

Output velocity vector map for the 3-pulse HBPOF simulation. There are many false and missed correlations.

Figure 17:

Output velocity vector map for the 5-pulse HBPOF simulation. The >pulse exposure yields fewer false or missed correlations.

Nat~onal Aeronaut8cs and Space Adrn~n8slrallon

Report Documentation Page

1. Report No.

3. Recipient's Catalog No.

NASA TM-101306
4. Title and Subtitle 5. Report Date

A p p l i c a t i o n o f O p t i c a l C o r r e l a t i o n Techniques t o P d r t i c l e Imaging V e l o c i m e t r y
I

6. Performing Organization Code

I
I

7. Author(s)

Mark P. Wernet and Robert V. Edwards

8. Performing Organization Report No.

E-4297
10. Work Unit No.

9. Performing Organization Name and Address

N a t i o n a l A e r o n a u t i c s and Space A d m i n i s t r a t i o n Lewis Research Center Cleveland, Ohio 44135-3191

1
11. Contract or Grant No. 13. Type of Report and Period Covered

N a t i o n a l A e r o n a u t i c s and Space A d m i n i s t r a t i o n Washington, D.C. 20546-0001

1
14. Sponsoring Agency Code
I

Techni c a l Memorandum

Prepared f o r t h e Conference on Sensors and Measurement Techniques f o r A e r o n a u t i c a l A p p l i c a t i o n s cosponsored by t h e AIAA, NASA, and AFWAL, A t l a n t a , Georgia, September 7-9, 1988.
16. Abstract

Pulsed l a s e r sheet v e l o c i m e t r y y i e l d s n o n - i n t r u s i v e measurements o f v e l o c i t y v e c t o r s across an extended 2-dimensional r e g i o n o f t h e f l o w f i e l d . The a p p l i c a t i o n o f o p t i c a l c o r r e l a t i o n techniques t o t h e a n a l y s i s o f mu1 t i p l e exposure l a s e r l i g h t sheet photographs can reduce a n d l o r s i m p l i f y t h e d a t a r e d u c t i o n t i m e and hardware. Here, Matched S p a t i a l F i l t e r s (MSF) a r e used i n a p a t t e r n r e c o g n i t i o n system. U s u a l l y MSFs a r e used t o r e c o g n i z e c h a r a c t e r s o r assembly l i n e p a r t s . I n o u r a p p l i c a t i o n , t h e MSFs a r e used t o i d e n t i f y t h e i s o - v e l o c i t y v e c t o r contours i n t h e f l o w . The p a t t e r n s t o be r e c o g n i z e d a r e t h e r e c o r d e d p a r t i c l e images i n a p u l s e d l a s e r l i g h t sheet photograph. Measurement o f t h e d i r e c t i o n o f t h e p a r t i c l e image displacements between exposures y i e l d s t h e v e l o c i t y v e c t o r . I n t h i s work, t h e p a r t i c l e image exposure sequence i s designed such t h a t t h e v e l o c i t y v e c t o r d i r e c t i o n i s determined unambiguously. A g l o b a l a n a l y s i s technique i s used i n comparison t o t h e more common p a r t i c l e t r a c k i n g a l g o r i t h m s and Young's f r i n g e a n a l y s i s techniques.

17. Key Words (Suggested by Author(s))

18. Distribution Statement

Optical correlation P a r t i c l e imaging v e l o c i m e t r y

Unclassified - Unlimited S u b j e c t Category 35

19. Security Classif. (of this report)

20. Security Classif. (of this page)

21. No of pages

22. Price*

Unclassified
NASA FORM 1626 OCT 86

Unclassified

22

AO 2

'For sale by the National Technical Information Service, Springfield, Virginia 22161

National Aeronautics and Space Administration

FOURTH CLASS MAIL ADDRESS CORRECTION REQUESTED

Postage and Fees P a ~ d Nal~onal Aeronaut~cs and Space A d r n ~ n ~ s l r a l ~ o n NASA 451

NASA