ImageAnalStereol2009;28:35-43 OriginalResearchPaper VELOCITY FIELDCOMPUTATIONIN VIBRATEDGRANULAR MEDIAUSINGAN OPTICALFLOWBASED MULTISCALE IMAGEANALYSISMETHOD JOHANDEBAYLE1,AHMEDRAIHANE1,ABDELKRIM BELHAOUA1,OLIVIERBONNEFOY1, G´ERARDTHOMAS1,JEAN-MARCCHAIX2ANDJEAN-CHARLESPINOLI1 1´EcoleNationaleSup´ erieuredesMinesdeSaint-Etienne. La boratoiredesProc´ ed´ esen Milieux Granulaires (LPMG), UMR CNRS5148. 158,coursFauriel,42023Saint- ´Etiennecedex2,France,2SIMAP-LTPCM, INPGrenoble-CNRS-UJF.BP 75-DomaineUniversitaire. 3840 2Saint-Martin d’H` eres,France e-mail: debayle@emse.fr (AcceptedFebruary9, 2009) ABSTRACT An image analysis method has been developed in order to compu te the velocity field of a granular medium (sand grains, mean diameter 600 µm) submitted to different k inds of mechanical stresses. The differential methodbasedonopticalflow conservationconsistsin descri bingadensemotionfield withvectorsassociated to each pixel. A multiscale, coarse-to-fine, analytical app roach through tailor sized windows yields the best compromisebetween accuracy and robustnessof the results, while enablingan acceptable computationtime. Thecorrespondingalgorithmispresentedanditsvalidatio ndiscussedthroughdifferenttests. Theresultsofthe validationtestsoftheproposedapproachshowthatthemeth odissatisfactorywhenattributingspecificvalues toparametersinassociationwiththesizeoftheimageanaly siswindow. Anapplicationinthecaseofvibrated sandhasbeenstudied. Aninstrumentedlaboratorydevicepr ovidessinusoidalvibrationsandenablesexternal opticalobservationsofsandmotionin3Dtransparentboxes . At 50Hz, byincreasingtherelativeacceleration Γ, the onset and developmentof two convectiverolls can be obs erved. An ultra fast camera recordsthe grain avalanches, and several pairs of images are analysed by the p roposed method. The vertical velocity profiles arededucedandallowtopreciselyquantifythedimensionso fthe fluidizedregionasa functionof Γ. Keywords: grain motion, granular media, horizontal vibrat ions, multiscale image analysis, optical flow, velocityfield. INTRODUCTION The aim of this work is to compute the velocity field of a granular medium submitted to different kinds of mechanical stresses. In particular, it focuses on the velocity field in a vibrated granular medium of sand grains (mean diameter 600 µm) and allows assessingthedimensionsofits fluidizedregion,which correspondsto specific physicalproperties required in processengineeringforindustrialapplications. For this purpose, an image analysis method has been developed. 2D image sequences of the granular medium are firstly acquired with an ultra- rapid CCD camera. The optical flow (dense velocity field) between two images is then computed using a differential method. This technique is based on two assumptions. The first one is the light intensity conservation, implying that the luminance of the grains does not change significantly between two successive image frames at time values tandt+dt.The second one requires that each velocity vector remains constant within a small spatial window. The size of the analyzing window is related to the accuracy of the method. In order to handle small and large displacements, the optical flow is computed with an iterative and multiscale approach (pyramidal decomposition) using a coarse-to-fine analysis. The level of the pyramid multiscale analysis is related to the robustness of the method. In order to get the best compromisebetweenaccuracyandrobustness,several validation tests are achieved.For a given numerical or physicalmotion(translation,rotation),pairsofimages beforeandafterthedisplacementarecompared.Using the optimal parameters, the image analysis method is performed in the framework of a study of an horizontally vibrated granular medium submitted to different accelerations. The vertical velocity profiles are then deduced and allow to precisely quantify the dimensions of the fluidized region as a function of the vibration acceleration. 35 DEBAYLEJET AL:Velocityfieldcomputationinvibratedgranularmedia MATERIALS,EXPERIMENTAL SETUPSANDIMAGE ACQUISITION MATERIALS The granular medium is made of sand grains (more than 99% of α-quartz) without internal close porosity (2660 kg/m3). The size distribution given on Fig. 1 shows a moderate polydisperse material with a span of (d90%−d10%)/d50%=0.75 and a volume mean diameter, denoted d4,3(the ratio of the fourth to the third moments of the grain size distribution), of 600 µm. The grains are approximately spherical (Fig. 2). Fig. 1.Volumetric grain size distribution. The mean diameter is d 4,3=600 µm with a span of (d90%− d10%)/d50%=0.75 where n %in mass of the grains havediameterlowerthand n%. Fig.2.Opticalmicroscopeimageshowingtherounded shapesofthe sandgrains. EXPERIMENTAL SETUPS Threeexperimentalsetupshavebeenused. Translation bench. The translation bench is composed of a parallelepipedic open box (container) filled with sand. This container is fixed on an horizontal linear motion guide, allowing translational motion exclusively. The camera is placed above the setup.Rotation bench. This setup is composed of a plastic cap of 3.3 cm inner diameter with the open- side up. A pushpin goes through it center and fixes it on an horizontal plank, allowing rotational motion exclusively.Thecontainerisfilledwithsandgrainsand the camera is placed abovethe setup, aligned with the pushpin’sverticalaxis(Fig. 3). Fig. 3.Rotation bench: a plastic cap filled with granular material can undergo a rotational motion aroundits centerwherea pushpinisnailed. Vibration bench. This experimental setup (Raihane et al., in press) is mainly composed of a mobile table, a transparent container and sand grains. The horizontal table is linked to a marble stand by four horizontal linear motion guides ( THK, EPF7M16+55LM ). A sinusoidal motion is transmitted to the table by an electromagnetic shaker ( TIRA, S513) driven by a signal generator ( LING DYNAMIC SYSTEMS, DSC4 ) coupled with a power amplifier (TIRA, BAA120 ). The motion is driven with an accelerometer ( BRUEL & KJAER, 4371 V ). The direction of the vibrations is referred to as x-axis. The parallelepipedic container fixed on the table is composed of 8 mm thick glued plexiglas plates (Fig. 4). Fig. 4.Vibration bench: an electromagnetic shaker delivers a sinusoidal vibration at a controled frequency f. It drives an horizontal table mounted on four bearings guided by four rails. The vibration amplitude A is controlled by a piezoelectric accelerometer. 36 ImageAnalStereol2009;28:35-43 The longitudinal, transverse and vertical dimensions Lx,Ly,Lzof the containerare respectively 40, 80 and 80 mm. The initial height of the granular packingis60mmandthecameraobservesthevertical faceorthogonalto thevibration x-axis. IMAGEACQUISITION An ultra-fast CCD camera ( JAI, CMOS CV-A33 ) with a maximal frequency of 5400 frames/s and a maximal field of 494 ×660 pixels allow us to follow the movements of the grains along the walls of the box.Theanalysisoftheimagessequencesisdescribed below. VELOCITYFIELDCOMPUTATION This section is devoted to an image analysis method allowing the computation of a dense velocity field from an image sequence of vibrated granular structures. For that, an optical flow multiscale method is proposed. OPTICALFLOW BASED IMAGE ANALYSIS METHOD Optical flow (Horn and Schunck, 1981) is an estimationoftheapparentmotion(velocity)ofobjects within an image sequence. It is closely related to motion estimation. Nevertheless,the term optical flow is specifically used to describe a dense 2D apparent motion field from the projection of a 3D scene onto theimageplane.Differentmethodshavebeenreported for determiningopticalflowsuchas: –phase-based methods (Fleet and Jepson, 1990; Fleet,1992):thevelocityisdefinedin termsofthe phase behavior of band-pass filter outputs. They provide a high accuracy but are generally less efficientforlarge displacements; –differential methods (Horn and Schunck, 1981; Lucas and Kanade, 1981; Uras et al., 1988): velocity is computed from spatio-temporal derivatives of image intensities. These techniques offer a good trade-off between robustness under noiseanddensityofthe flowfields; –region-basedmatchingmethods(Burt etal., 1983; Anandan, 1989; Little and Verri, 1989): the velocity is defined as the shift yielding the best fit between image regions, according to some similarity measure. These methods provide more robustness with respect to differentiation and are generally quicker but they are less successful for sub-pixelvelocities;–energy-basedmethods(AdelsonandBergen,1986; Heeger, 1988; Barman et al., 1991): optical flow is computed using the output from the energy of velocity tuned filters in the Fourier domain. It hasbeenshownthatcertainenergy-basedmethods are equivalent to region-based or differentiation methods. In review papers, Barron et al.(1994) and Galvin et al.(1998) have evaluated different methods and concluded that the differential method proposed by LucasandKanade(1981)yieldsthebestresults. Other numerous methods combining local (using spatial constancy assumptions) and global (using smoothness assumptions) techniques have also been developed so as to propose an accurate smooth and dense optical flow field (Brox et al., 2004; Bruhn and Weickert,2005;LeBesneraisandChampagnat,2005). For the proposed application, a dense optical flow is notabsolutelynecessaryandarawvelocityfieldisfirst required. Therefore, the computation of the velocity field in this paper is based on the differential method proposedbyLucasandKanade(1981). The Lucas and Kanade (1981) method, as many differential techniques, is based on the assumption of intensityconservation( i.e.,assumingthattheintensity of the objects within the image sequence Idoes not change significantly between two successive image framesattimes tandt+dt),thatis to say: I(x,y,t) =I(x+dx,y+dy,t+dt),(1) wheret(x,y)andt(dx,dt)are vectors of the pixel location and displacement respectively;tdenoting the matrix (vector)transposition. For small displacements, a first order Taylor expansioncanbe applied: I(x+dx,y+dy,t+dt) =I(x,y,t) + ∂I ∂x(x,y,t)dx+∂I ∂y(x,y,t)dy+∂I ∂t(x,y,t)dt.(2) Combining Eq. 1 and Eq. 2, the optical flow constraintis thengivenas: ∂I ∂x(x,y,t)dx+∂I ∂y(x,y,t)dy+∂I ∂t(x,y,t)dt=0,(3) whichresultsin: ∂I ∂x(x,y,t)vx+∂I ∂y(x,y,t)vy+∂I ∂t(x,y,t) =0,(4) wherev=t(vx,vy)isthevelocityvectoratlocationand time(x,y,t). 37 DEBAYLEJET AL:Velocityfieldcomputationinvibratedgranularmedia Thus, Eq. 4 could be written as the differential system: (Ix,Iy)t(vx,vy)+It=0, (5) whereIx,IyandItdenote the space and time partial derivativeof I(x,y,t),respectively. This is a system of one equation with two unknowns vx,vyant this it can not be solved as such. This is known as the aperture problem of the optical flow algorithms. Consequently, another set of equations is needed, given by some additional constraints(Horn andSchunck,1981). To solve the aperture problem, Lucas and Kanade (1981) assume a locally constant velocity in the neighborhood of the considered pixel. Assuming that the velocity vector v=t(vx,vy)is constant within a small window Wofsizenpixelscenteredat (x,y)and ordering the pixels within as 1 ,...,n, yields an over- determinedsystem:  Ix1Iy1 Ix2Iy1...... IxnIyn /bracketleftbigg vx vy/bracketrightbigg + It1 It2... Itn =0.(6) This system can be written in the following matricial form: Av+b=0, (7) whereA,v,bdenotethespacepartialderivativematrix, the velocity vector and the time partial derivative vector,respectively. Tosolvethissystemofequations,theleastsquares method is used in the Lucas-Kanade optical flow estimation: v= (tAA)−1[tA.(−b)], (8) thatis to say: /bracketleftbigg vx vy/bracketrightbigg = n ∑ i=1I2 xin ∑ i=1IxiIyi n ∑ i=1IxiIyin ∑ i=1I2 yi −1 −n ∑ i=1IxiIti −n ∑ i=1IyiIti . (9) This estimation is reliable if the matrixtAAis invertible, i.e.,withnozeroeigenvalues.Consequently, the Lucas and Kanade (1981) method computes the vectorvatapixel (x,y)onlyifbotheigenvalues λ1,λ2 oftAAare greater than a predefined threshold value λ. Note that if λ1,λ2are both large, the point (x,y) corresponds to a corner point (Harris and Stephens, 1988).Therefore,intheproposedmethod,thevelocity vectorsare computedatcornerpoints.MULTISCALE APPROACH Thetwokeycomponentsofanopticalflowmethod are accuracyand robustness.These componentsrelate to thesizeofthewindow W: –a smallwindowis preferablefor accuracyin order not to “smooth out” the details contained in the images, –a large window is preferable for robustness to handlelarge motions. Thereisthereforeacompromisebetweenaccuracy and robustness when choosing the window size. In this way, a multiscale and iterative implementation (Bouguet, 2000) is proposed for providing both a robustandaccuratemethod. The multiscale representation is given by a Gaussian pyramid (Burt and Adelson, 1983). In this way, the optical flow is computed in a coarse-to- fine analysis. The velocity is first computed at a low resolution ( i.e., at a high level of the pyramid), in an iterative way so as to ensure the stabilization of the solution, before advancingto the next level. Then, the field is warped and upsampled at the next level of the pyramid and the optical flow is computed again. This loop estimation is performed until the highest image resolution( i.e.,atthelowestlevelofthepyramid). This iterative multiscale algorithm, denoted Motion-2D (Belhaoua, 2007), for estimating the optical flow between two image frames is therefore decomposedin thestepspresentedin Alg.1. Note that all computations are achieved at a subpixelaccuracylevel,usingbilinearinterpolation. REGULARIZATION Despite that velocity vectors are only computed for large eigenvalues oftAA(see section about optical flow), some outliers could appear.So, in orderto have a smoothed optical flow ( i.e., to discard outliers), a median filtering could be applied on the dense motion field(usingeitherthenormortheangleofthevectors). In the proposed algorithm, the median filtering is only performed if the norm or the angle of the considered vector v=t(vx,vy)is far from the median value of those within a 3 ×3 window centered on the pixellocatedat (x,y). This regularization process has been added to Motion-2D asanoption. 38 ImageAnalStereol2009;28:35-43 input:ImagesIandJ output:Velocityfield betweenimages IandJ pyramidaldecompositionin L+1levelsofthe two imageframes IandJ:I0=I,I1,...,ILand J0=J,J1,...,JL; initialization oftheopticalflow: vL+1=0; forl←Lto0do imageI+ l:warpingoftheimage Ilwith the upsampledopticalflow v∗ l+1: I+ l=Il+v∗ l+1; loopestimationofthe opticalflowbetween imagesI+ landJ: fork←1toKdo ηk l≡iterativeLucas-Kanade. ηk lis the opticalflowatlevel lfor iteration k. end opticalflowatlevel lbetweenimages I+ l andJ:ηl; finalopticalflowatlevel lbetweenimages IlandJ:vl=v∗ l+1+ηl; end estimation ofthe finalopticalflowbetween imagesIandJ: v=v0=L ∑ l=0η∗ l η∗ lis theopticalflow ηl+1upsampledatlevel l; Algorithm 1 : Algorithm Motion-2D . The velocity field is computed using an iterative multiscale imageanalysismethod. VALIDATION In general, the performance of the optical flow methodsareevaluatedonrealsequencesandsynthetic sequences for which motion fields are known. Some common evaluation criteria are average angular error, standard deviation and density of measurements (Barronetal., 1994). GENERATION OFTEST IMAGES In order to validate the implementation of the proposed algorithm in terms of accuracy and robustness, several tests have been performed on images of granular structures with imposed motion. Fourpairs ofimageshavebeenbuilt: •A−B(numerical test / translation): A picture Aof a sand packing is shifted by exactly ppixels with Aphelion™ softwareto givea picture B.•C−D(numerical test / differential translation): A pictureCof a sand packing is split into different zonesCiwhichareshiftedby pipixelseachtogive apictureD. •E−F(physical test / translation): The sand packing is shifted by dmm with the translation bench. The pictures of the initial and final states arelabeled EandFrespectively. •G−H(physical test / rotation): The sand packing isrotatedbyanangleof αdegreeswiththerotation bench.Thepicturesoftheinitialandfinalstatesare labeledGandHrespectively. VELOCITYFIELD COMPUTATION FORTEST IMAGES A−B: Various pairs of images A−Bwere created corresponding to translations between one and ten pixels. Several analyses were performed on these images by varying the level of pyramid Lfrom 0 to 4 and the window size Wfrom 3 to 25 pixels. Results of these analyzes are presented on Fig. 5. They show that there exists conditions for a displacement of ppixels to be accurately estimated. In fact, when using small window ( W9). Results obtained for pyramid level L=2 are satisfactory for all the investigated window sizes. For the following analyzes, the pyramid level L is kept equal to two and the window width Wis chosen to be equal to three grain’s mean diameter d4,3. The latter choice will be justified in the sectionResults. 39 DEBAYLEJET AL:Velocityfieldcomputationinvibratedgranularmedia Fig.5.Numericaltest: thecomputeddisplacement (in pixels) between two images A and B (for different numerical translations) is plotted versus thepyramidlevelL fordifferentwindowsizesW. C−D: Using a pyramid level L=2 and a window widthW≈3d4,3, the displacement field between imagesCandDwas calculated. For each z, the displacement is averaged along the x−axis. The resulting profile is plotted on Fig. 6. The different regions correspondingto different values of displacement ppixels are well detected and separated. Fig. 6.Differential numerical test: the mean computed and real displacements between two images C and D are plotted versus the vertical position z. The thick (resp. thin) line represents the computed (resp. real) displacement. The inset shows the different translations performed on imageC to obtainimage D. E−F: Several translation tests, from 0.1 mm to 1.5 mm, were performed. The resulting images were analyzed with parameters L=2 andW≈ 3d4,3. Table 1 shows the experimental results and calculatedresults.Table 1. Physical translation: Comparison betweenrealand p calculateddisplacementofasand packing. The image analysis method is performed with parametersL =2andW =15. d[mm]d[pixel] pcalculated 0.1 0.82 0.83 ±0.02 0.3 2.46 2.49 ±0.02 0.4 3.28 3.41 ±0.03 0.5 4.1 4.2 ±0.09 0.8 6.56 6.63 ±0.04 1. 8.2 8.29 ±0.09 1.5 12.3 12.08 ±0.3 G−H: Fig. 7 shows the result of the rotation test analysis by Motion-2D . By going away from the center of rotation (pushpin center), the displacement increases according to the law d= α×r(r:radius). Fig. 7.Displacement field from rotation: the cup filled with granular material was rotated around the pushpin. The resulting images were analyzed with parametersL =2andW =15≡3D4,3. Following these validation experiments, the algorithm is parameterized with a pyramidal level L=2levels,a numberof100iterationsanda window sizeWofapproximatelythreemeangraindiameters. RESULTS EXPERIMENTAL OBSERVATIONS In vibration experiments, a typical run consists in increasing the value of the relative acceleration Γ from zero at a fixed frequency f=50 Hz. Different behaviors were observed as a function of Γ. To facilitate the discussion, the vertical faces orthogonal 40 ImageAnalStereol2009;28:35-43 (resp. parallel) to the vibrations direction are labeled “North/South” (resp. “East/West”). The observations are thefollowing ones(Raihane etal., in press): 0<Γ<Γ1∼0.3 : The whole granular packing behaves like a compact solid which follows the motion ofthebox(atrestin theboxframe). Γ1<Γ<Γ2∼0.9: Grains simmering at the free surface; the other regions are at rest in the box frame. Γ2<Γ<4.5: AtΓ2, a transition occurs from a solid- like behavior of the whole packing to a fluid- like one in the upper region, the lower region remaining unchanged. At the free surface, grains move towards the walls from the middle of the granular medium and then, avalanche into the gap created between the granular packing and the North/South walls. On the East/West faces, two convective counter-rotating rolls are observed (Fig. 8). Only the upper region of the packing is fluidized.Byincreasingtherelativeacceleration Γ, thethicknessofthefluidizedregion Efincreases. Fig. 8.Convective rolls: sketch of convective motion in a granular packing submitted to horizontal and sinusoidal vibrations with acceleration Γ>1. In the west/east faces, twocounterrotatingconvectiverollsareobserved. At the north/south walls, grains avalanche to a distanceE ffromthe top. QUANTITATIVEDESCRIPTION Invibrationexperiments,theaimistomeasurethe velocity field of grains at the North/South faces and the fluidized region thickness Ef. For this purpose, a CCD camera records the grain avalanchesat the north face at the rate of one frame per oscillation period (50 framespersecond).Foreachacceleration,severalpairs ofimagesareanalyzedby Motion-2D (10pairsfor Γ≤ 1.8, 30 pairs else). The resulting displacement fields are averagedandthen,foreachacceleration,a vertical velocity profile is deduced (Fig. 9). For all these analyses, the pyramid level is L=2 and the windowwidthW≈3d4,3. The choiceof suchwindow width is justified by the fact that, during their fall, grains can undergo simultaneously rotation motions. Hence, the algorithm would meetdifficulties to recognizethem if the window recovers only one grain. A window width of three grain mean diameters seems to be a good compromise to avoid chaotic behavior of granular movements and to take into account the displacement gradient(Medved,2002). Fig. 9.Velocity profiles of grains at the North wall for a granular packing of initial height Hi=60 mm horizontally vibrated at three different accelerations. Profiles show a flat zone and a large peak corresponding respectively to the solid lower zoneandthe fluidizedupperzoneofthepacking. To determine the thickness of the fluidized region for each acceleration, a fluidization velocity threshold isdefined(insetofFig.9).Here,whenthethresholdof 0.06mm/sis exceeded,thegranularmediumis saidto befluidized. Fig. 10 shows the evolution of the fluidized thickness Efversus the relative acceleration Γfor a granularpackingvibratedatfrequency f=50Hz. Fig. 10. Evolution of the fluidized thickness with relative acceleration Γof a granular packing of initial height H i=60 mm. An onset of fluidization is observed for Γ=1and a saturation for high accelerations. 41 DEBAYLEJET AL:Velocityfieldcomputationinvibratedgranularmedia The phenomena occuring on East and West faces canbedescribedinthe samewaywith Motion-2D . Fig.11showsthe computeddisplacementfieldfor animagesequenceduring20s. (a) image at t=20s (b) velocityvector field (c) velocity streamlines Fig.11.Computedvelocityfieldofagranularmedium vibrated during 20 s. The vector velocity field (b) is computed between images from t =0 s to t =20 s (a). The streamlines (c) shows the direction and the magnitudeoftheflowvelocity.The velocity field shows two convective rolls in the left and right upper corners of the image. This numericalresultisinaccordancewiththeexperimental observations(Fig. 8). Finally,thisquantitativedescriptionusing Motion- 2Dprovides the velocity field in a vibrated sand granular packing and therefore enables to estimate its fluidizedthickness. CONCLUSION An optical flow based multiscale image analysis method has been developed in order to compute the velocityfieldinvibratedgranularmediaconstitutedof sandgrains.Thismethodprovidesagoodcompromise between accuracy and robustness for acceptable computation time, using specific parameters assessed from several validation tests. The velocity vectors are computed both using tailor sized windows and a coarse-to-fine analysis on 2D image sequences. The results provide a quantification of the fluidized thickness of the granular medium versus the relative accelerationofthevibration (usingafixedfrequency). REFERENCES Adelson EH, Bergen JR (1986). The extraction of spatiotemporal energy in human and machine vision. In: Proc IEEE Worksh Visual Motion. Charleston, SC, USA;May,151-66. Anandan P (1989). A computational framework and an algorithm for the measurement of visual motion. Int JComputVision2:283–310. Barman H, Haglund L, Knutsson H, Granlund G (1991). Estimation of velocity, acceleration and disparity in time sequences. In: Proc IEEE Motion Worksh. Princeton,NJ, USA; October,44-51. Barron JL, Fleet DJ, Beauchemin SS, Burkitt TA (1994). Performance of optical flow techniques. Int J Comput Vision12:43–77. Belhaoua A (2007). Analyse et mod´ elisation dynamique d’images de structures granulaires soumises ` a des vibrations. MSc Thesis, ENSM-SE, Saint-Etienne, France. BouguetJY(2000). Pyramidalimplementationofthelucas- kanade feature tracker. Tech. rep., Intel Corporation, MicroprocessorResearchLabs. Brox T, Bruhn A, Papenberg N, Weickert J (2004). High Accuracy Optical Flow Estimation Based on a theory for warping. In: Proc 8th Eur Conf Comput Vision. Prague,CzechRepublic;May,25-36,vol.4. 42 ImageAnalStereol2009;28:35-43 Bruhn A, Weickert J (2005). Lucas/Kanade Meets Horn/Schunck: Combining Local and Global Optic Flow Methods. IntJComputVision61:211–31. Burt P, Yen C, Xu X (1983). Multiresolution flow-through motionanalysis. In:ProcIEEEIntConfComputVision PatternRecogn.Washington,DC, USA; June,246-52. Burt PJ, Adelson EH (1983). The laplacian pyramid as a compactimagecode. IEEETransCommun31:532–40. FleetD,JepsonA(1990).Computationofcomponentimage velocity from local phase information. Int J Comput Vision5:77–104. Fleet DJ (1992). Measurement of Image Velocity. Kluwer AcademicPublishers,Norwell. GalvinB,McCaneB,NovinsK,MillsS(1998).Recovering motion fields: An evaluation of eight optical flow algorithms. In: Proc Brit Machine Vision Conf. Southampton,UK;September,195-204. HarrisC,StephensMJ(1988). Acombinedcornerandedge detector. In: ProcAlveyVision Conf.Manchester,UK; September,147-52.HeegerDJ(1988). Opticalflowusingspatiotemporalfilters. IntJComputVision1:279–306. Horn BKP, Schunck BG (1981). Determinig optical flow. ArtifIntel17:185–203. Le Besnerais G, Champagnat F (2005). Dense optical flow byiterativelocalwindowregistration. In:ProcIEEEInt ConfImageProc.Genova,Italy.September,137-40. Little JJ, Verri A (1989). Analysis of differential and matching methods for optical flow. In: IEEE Worksh VisualMotion.Irvine,CA, USA; March,173-180. LucasBD,KanadeT(1981). Aniterativeimageregistration technique with an application to stereo vision. In: ProcIntJointConfArtifIntel.Vancouver,BC, Canada; August,674-9. Medved M (2002). Connections between response modes in a horizontally driven granular material. Phys Rev E 65:021305. Uras S, Girosi F, Verri A, Torre V (1988). A computational approachtomotionperception. Biol Cybern60:79–97. 43