Bepaling van gemiddelde en spreiding van de snelheidsverdeling met de kruiscorrelatiemethode

(1)

I

1-1

I

1-I

I

1-I

l!

f

~t

T

U

Delft

TechnischeUniversiteit Delft

Faculteit der Civiele Techniek Vakgroep Waterbouwkunde

(2)

I

I-I

I

Bepaling van gemiddelde en spreiding van de snelheidsverdeling

met de kruiscorrelatiemethode J. de Graaff en R.E. Slot rapport no. 6-95, oktober 1994

(3)

I

INHOUD blz 1. INTRODUCTIE. 2. PROBLEMEN . 3. METHODS 4. PROGRAMMATUUR • 5. TESTS . 6. RESULTATEN . 7. DISCUSSIE 8. CONCLUSIES • 9. VERDER ONDERZOEK 10. LITERATUUR . 1 2 3 4 21 22 24 26 28 30 APPENDICES

Ai. FLOW DIAGRAMS

Fig.: A 1.Velocity determining program..

Fig.: A2. Measurement of ~mean and ~cr in velocity distribution program determined by the program with regard to the reference .. Fig.:A3. Determining measurements of the partieles with use of

cross correlation method.

Fig.: A4. Linked list structure of object numbers at t=t1 and object numbers at t=t2 after first selection with use of ref. area.

Fig.: A5. Creating referenee: symmetrie velocity distribution. Fig.: A6. Creating referenee: asymmetrie velocity distribution. Fig.: A7. Generating partieles.

Fig.: A8. Histogram of velocity distribution with mean and sigma.

A2. DESCRIPTION OF THE PROGRAMS A 1. MAIN PROGRAM (CRITOp2y.tip) A 1.1. Generating referenee

A 1.2. Mean and sigma values . A 1.3. Generating partieles

A 1.4. Remove edge eonneeed objeets

A 1.5. Determining positions of partieles from images A 1.6. Cross eorrelation and referenee area

A 1.7. Criterium

A 1.8. Reference area

A 1.9. Reduee number of ref. area to 1

A 1.10. Movements, velocities and velocity distribution A 1.11. Mean and sigma of the distribution

A2. Module: genpv A3. Module: genpvaa . A4. Module: movrxy1 . A5. Module: genpartde

A1 A2 A3 A4 A5 A6 A7 A8 A9 A9 A9 A10 A10 A11 A12 A13 A13 A13 A14 A15 A17 A18 A19 A2.0

(4)

hist1.c

I

pg3

I

K22

I

K22 /\23 /\24

I

/\25

I

A2. DESCRIPTION OF THE PROGRAMS, continued A6. Module: thrns256.c

A7. Module: meansigm1.c

A8. Histogram of a velocity distribution A8.1. Module: hist1.c

A8.2. Module: plverdn1

B. PROGRAM USTINGS

Files of programs descripted in this document CRITOP2y.tip genpv.c genpvaa.c genpartdc.c Rxy256.c movrxy1.c thrns256.c meansigm1.c C. FIGURES

(5)

I

1. INLEIDING

De laatste tijd wordt er veel onderzoek gedaan aan cohesief sediment. Het speelt een belangrijke rol bij het ondieper worden van havens en waterwegen, en bij andere soortgelijke problemen. Om cohesief sedimenttransport te kunnen voorspellen zijn gegevens nodig over de verdelingen van de groottes en de bezinksnelheden. Er bestaan veel methoden om de groottes van gesuspendeeerde deeltjes te bepalen,

maar de meeste zijn niet toepasbaar bij cohesieve slibvlokken,omdat die kwetsbaar zijn. Als de vlokken niet meteen bij het bemonsteren uiteenvallen gebeurt het wel in de verdere analyse b.v. de Coulter Counter of de pipetmethode. Bij analyse met de owen-buis ontstaat er nog een extra probleem:de lange duur van de analyse leidt tot verdere vlokvorming, waardoor de gemeten verdeling nog minder realistisch wordt.

Om deze problemen op te lossen worden opnames gemaakt door onderwater-camera's, die direkte informatie geven over de onverstoorde monsters. Uit één opname kunnen de vlokgroottes worden bepaald, en uit twee achtereenvolgende opnames met bekende tussenliggende tijd kunnen de bezinksnelheden worden bepaald.

Tot nu toe gebeurde de analyse van vlokopnames meestal met de hand.

Beeldbewerking met de computer is een manier om dit automatisch te doen. Omdat dit tijd bespaart kunnen meer vlokken worden geanalyseerd,wat tot representatievere verdelingen leidt.

Het onderwerp van dit rapport is het ontwikkelen en testen van een beeldbewer-kingsprogramma om het gemiddelde en de spreiding van de snelheidsverdeling te bepalen door middel van de kruiscorrelatie. Dit gebeurt met behulp van gesimuleerde opnames van deeltjes op een daarvan te onderscheiden achtergrond.

Een belangrijk nadeel van kruiscorrelatie is de zeer lange rekentijd. In dit rapport worden methoden beschreven om de rekentijd terug te brengen.

(6)

2. PROBLEEMSTELLING

De kern van het probleem is hoe dezelfde slibvlok in de eerste opname wordt teruggevonden in de tweede opname. De opnames kunnen wel 150 slibvlokken of meer bevatten.

Gedurende de tijd tussen de twee opnames kunnen de slibvlokken in het oog van de camera:

(1) - veranderen in vorm; (2) - veranderen in grootte; (3) - verdwijnen of verschijnen. (1) is als gevolg van rotatie.

(2) is als gevolg van passage door de lichtmes. Op het moment van de eerste opname treedt een slibvlok de lichtmes binnen en in de tweede opname bevindt het zich in de lichtmes met als gevolg dat gedurende de tijd tussen twee opnames het deeltje (aanzienlijk) groter is geworden, het opblaas ofwel explosie effect. Of omgekeerd bij het verlaten van de lichtmes waarbij het deeltje in de tweede opname (aanzienlijk) kleiner is geworden. (implosie effect).

(3) is eveneens als gevolg van de passage door de lichtmes. Gedurende de tijd tussen twee opnames kan een deeltje de lichtmes binnentreden, met als gevolg dat het deeltje alleen aanwezig is in de tweede opname. Het is alsof er een extra deeltje is bijgekomen. Het omgekeerde kan ook, waarbij een deeltje in de tweede opname is verdwenen.

En bovendien kunnen de deeltjes elkaar overlappen.

Hieruit blijkt dat als de tijd tussen de twee opnames groot is, het deeltje door middel van patroonherkenning zeer moeilijk is terug te vinden. Hoe groter de tussentijd, des te groter is de kans dat er een verkeerd deeltje wordt gekozen. Is de tussentijd klein, dan is de verplaatsing zo klein dat digitalisatiefouten optreden.

Is het deeltje eenmaal gekozen, dan wordt het verschil van de posities van de zwaartepunten bepaald: de verplaatsing. En omdat de tijd tussen de twee opnames bekend is de snelheid.

2

I

(7)

I

3. METHODEN

Er zijn verscheidene methoden om het deeltje in de tweede opname terug te vinden: (1) - patroonherkenning

(2) - kruiscorrelatie, referentiegebied, patroonherkenning (3) - kruiscorrelatie, referentiegebied, hoogste kruiscorrelatie (4) - kortste afstand

I

Ad (1):

De keuze is gebaseerd op de vorm die het meest overeenkomt met het deeltje van

opname 1.

I

Ad (2):

Kruiscorreleren vanuit twee 5122 binaire plaatjes, als resultaat van segmentatie,

levert een kruiscorrelatiefunctie op, die na contrast stretching voorgesteld kan

worden als een 2562 grijswaardeplaatje. Een andere naam voor de

kruiscorrelatiefunctie is correlogram. Door middel van iteratieve threshold wordt een

referentiegebied bepaald, waaruit een of meer deeltjes van de tweede opname

worden geselecteerd. En in geval van meer deeltjes wordt aan de hand van

patroonherkenning zoals in (1) een deeltje gekozen.

Ad (3):

Op dezelfde wijze als (2), in dit geval wordt als zich meer dan een deeltje in het

referentiegebied bevindt, het deeltje gekozen aan de hand van de keuze van de

plaats die overeenkomt met de hoogste kruiscorrelatie.

Ad (4):

Het deeltje in de tweede opname wordt geselecteerd aan de hand naar de kortste

afstand, een programma dat door Rijkswaterstaat wordt ontwikkeld.

I

3

I

(8)

4. PROGRAMMATUUR

De programmatuur die is ontwikkeld, getest en in dit rapport is beschreven is gebaseerd op methode: (3), geïllustreerd aan de hand van het ondestaande figuur:

t=t 1 exposure

Voor het ontwikkelen en het testen van de programmatuur beschreven in dit rapport worden de opnames en daaropvolgende segmentatie gesimuleerd met de generator, zie onderstaande figuur:

t=t 1 generator

t=t2

Een belangrijk voordeel van het gebruik van deze generator is dat de deeltjes verplaatst worden met gegeven ofwel opgelegde snelheden, dus met een opgelegde snelheidsdistributie met bijbehorende gemiddelde en spreiding, de referentie. De referentie geeft de juiste snelheidsverdeling weer, waaruit de juistheid van de programmatuur kan worden bepaald. De programmatuur van de snelheidsanalyse begint vanaf de kruiscorrelatie.

4

I

(9)

I

4.1 Simuleren opnames en segmentatie

I

4.1.1 Systeem en cameraveld

I

Het veld van de camera beslaat slechts een gedeelte van het systeem ofwel het lichtmes waarin de slibvlokken zich bevinden. Gedurende de tijd tussen de twee opnames treden slibvlokken het cameraveld binnen of gaan er weer uit. In deze

simulatie is uitgegaan van het systeem met gekozen zijden: 1000 * 1000

pixeleenheden. Voor de zijden van het cameraveld, de window, is gekozen: 512

* 512 pixeleenheden. 1000 _w_hal_e _s_y_stem

I

o o o 512

V

/

I--- wind N ~ l[) aw (camera field)

I

In de simulatie is er rekening mee gehouden dat er deeltjes zijn die gedurende de tijd tussen de twee opnames het cameraveld binnentreden of verlaten.

4.1.2 Generator

I

De generator kan in twee delen worden gesplitst: positie-snelheid en deelt

jes-generator.

De positie-snelheidsgenerator geeft de posities van de zwaartepunten van de

deeltjes weer op het moment van de eerste opname en hun snelheden, de referentie. Deze zwaartepunten en snelheden worden als invoer gebruikt voor de deeltjesgenerator . Eveneens wordt in de deeltjesgenerator de parameter tijd, de tijd

tussen de twee opnames ingevoerd. Dan worden de liggingen van de

zwaartepunten op het moment van de tweede opname berekend. De ligging van de

zwaartepunten van de eerste zowel als van de tweede opname zijn gegeven ten

opzichte van de linkerbovenhoek van het gehele systeem.

Om tot binaire plaatjes te komen zijn eerst de liggingen van de zwaartepunten ten opzichte van de linkerbovenhoek van de window bepaald en van daaruit de deeltjes.

I

4.1.2.1 positie - snelheidsgenerator .

I

Het geven van de posities en de snelheden gaat at random. De posities worden in het systeem homogeen verdeeld. De snelheden worden verdeeld in twee soorten: (1) symmetrisch en (2) asymmetrisch.

De symmetrische verdeling heeft de vorm van een klok (gaussisch). Asymmetrisch

meer van een paissionverdeling. Voor het genereren van een symmetrische

snelheidsverdeling is hierin gebruik gemaakt van een kruis- of munt model en van een asymmetrische verdeling het biedtafelmodel.

I

5

(10)

I

Symmetrische snelheidsverdeling

I

De symmetrische snelheidsverdeling wordt gegenereerd aan de hand van het kruis

-of muntmodel. Kruis of muntmodel moet worden gezien als het werpen van een gegeven aantal munten en tellen hoeveel er kruis of munt zijn. Zoals bekend is de kans dat je van een muntstuk kruis of munt werpt gelijk aan % . De kans dat van n aantal munten k kruisen wordt geworpen is gelijk aan:

P (X=k)

=( ~)( ~

r

I

P(X

=

k) : kans op k aantal kruisen

I

Dit is een binominale kansverdeling. Als steeds met n - aantal munten wordt geworpen en geturfd wordt met het aantal kruisen dan zal het histogram een klokvormige verdeling tonen.

De random generator levert integer waarden van 0 tot en met 32767 = 215-1. De kansverdeling voor al deze waarden is homogeen. Kruis is geworpen als de generator een waarde heeft gelevert in het gebied van 1,2, ... ,16383 en munt als het in het gebied ligt van 16384,16385, ... ,32767.

I

Van daaruit wordt de snelheid gegenereerd:de grootte ofwel de absolute waarde van de snelheid en de richting, de hoek ten opzichte van het verticaal (zie onderstaande fig.):

I

centerpoint partiele

-.

I

y

t

I

6

I

(11)

I

De grootte van de snelheid die in deze tests is toegepast ligt in de range van 100 tot en met 120 (pixel eenheden) en de hoek in de range van -10 en 10°. In dit

model zijn 20 munten toegepast. Er worden 2 worpen uitgevoerd, één voor de

grootte en één voor de hoek. De grootte van de snelheid

I

v

I

en van de hoek a.

worden berekend volgens:

I

v

I

= 90

+

k,

a. = -10°

+

k20

I

k.: aantal kruisen 1e worp

k2: aantal kruisen 2e worp

Bijvoorbeeld 1eworp levert 15 kruisen en de 2e worp 5 kruisen op. Dit geeft voor

de grootte van de snelheid 100 +15

=

115 en voor de hoek -10° +5°

=

_5°.

I

Asymmetrische snelheidsverdeling

I

De asymmetrische snelheidsverdeling wordt gegenereerd aan de hand van het

biedtafel model. Behalve het poissonmodel bestaat er ook een andere methode om

een asymmetrische snelheidsverdeling te kunnen opleggen, die van het poisson

model aardig benadert: het biedtafelmodel. De naam van dit model is afkomstig van

de biedtafel een onderdeel van een zeer populaire spelletjesprogramma prijzenslag

uitgezonden door RTL 4.

Aan de biedtafel zitten 4 kandidaten om de prijs van het artikel te bieden. Voor het ene artikel bieden ze een beetje te hoog en voor het ander juist iets te laag. Dit kan worden vertaald naar de opgelegde snelheden. Het steeds te hoog of steeds te laag zitten heeft de neiging om volgens poisson te verdelen. Dit gebeurt ook met de snelheden van de slibvlokken.

De grootte van het bod hangt af van het artikel. Voor een ring ga je hoger bieden

dan voor een plastic eendje dat je goedkoop in de winkel kunt kopen.

Om deze afhankelijkheid te elimineren is als eerste stap genomen de 4 bods te

normaliseren naar de prijs van het artikel:

I

bod,

=

bod/prijs artikel i: 1,2,3,4

bod, : bod van de ie kandidaat

bod, : bod van de ie kandidaat genormaliseerd naar de prijs van het artikel

De (genormaliseerde) bods van alle 4 de kandidaten zijn afhankelijk van dat ene artikel. Daarom is van deze bods het gemiddelde bepaald:

I

4 bodg = % Lbodg, i

=

1

I

7

I

(12)

t

I

Dit kan verder worden vertaald naar de snelheid.

velocity siltparticle ...= bodg

I

in this exomple: .... bodg-l - velocity - 25 pixel unita

o

1 25 2 50 3 75 4 ~ bodg 100 ~ Vsi1t partiele (pixel units)

I

centerpoint portiele

<.

x

~

I

In dit voorbeeld komt een gemiddelde 1 van de genormaliseerde bods overeen met de snelheid van 25 pixeleenheden. Is het bijvoorbeeld 0.8 dan komt het overeen met de snelheid van 0.8 * 25 = 20 pixeleenheden.

Van daaruit wordt de snelheid gegenereerd in de x en in de y component:

I

y

t

I

_I

I

8

I

(13)

I

Vertaling naar de snelheid ging in dit geval over in:

Vx = konstante, * bodg1,2

vy = konstante, * bod93•4

I

waarin:

I

bodg'.2

=

%(bodg,

+

bod92)

I

bodg'.2 : gemiddelde vld genormaliseerde bods

vld eerste twee kandidaten

bod93.4 : gemiddelde vld genormaliseerde bods

vld derde en de vierde kandidaat

I

In deze tests zijn de waarden van de konstante, en konstante, respectievelijk 10

en 25.

Genereren gebeurt uit de hand van het lezen van de gegevens van het spelletje opgeslagen in een data file.

De referentie

I

De snelheidsgenerator geeft de opgelegde snelheden aan alle deeltjes weer die in

het hele systeem bevinden. Het programma bepaalt de snelheidsverdeling van de

slibvlokken die zich in het cameraveld bevinden. Dus gekeken is naar de opgelegde

snelheidsverdeling voor de vlokken die zich in het cameraveld bevinden. Er is

rekening mee gehouden dat gedurende de tijd tussen de twee opnames sommige deeltjes het cameraveld intreden of verlaten.

De deeltjes worden geselecteerd aan de hand van de ligging van hun

zwaartepunten die zich zowel op het moment van de eerste als van de tweede

opname in het cameraveld bevinden. Onderstaande figuur geeft een voorbeeld weer

voor 4 deeltjes waar 2 en 3 worden geselecteerd en hun opgelegde snelheden als

referentie worden gebruikt.

I

camera field

/

I

first exposure second exposure

9

(14)

Dit is nodig omdat op het moment van de eerste opname deeltje 1 zich buiten het cameraveld bevindt en evenzo op het moment van de tweede opname deeltje 4. Daarom zijn uit het oogpunt van de camera de verplaatsingen van deze deeltjes onbekend, en ze kunnen dus niet als referentie worden toegepast.

I

Deze generator levert objecten op in de binaire plaatjes behorende bij de eerste zowel als bij de tweede opname.

Als invoer worden toegepast: de liggingen van de zwaartepunten en de opgelegde snelheden afkomstig van de positie-snelheidsgenerator, de tijd tussen de twee opnames en de radius van de cirkelvormige objecten. Met deze snelheden en de tijd worden de liggingen van de zwaartepunten op het moment van de tweede opname berekend ten opzichte van de linkerbovenhoek van het cameraveld. En vervolgens worden de objecten in de binaire plaatjes vastgelegd. Voor objecten zijn in deze tests tot dusverre alleen cirkelvormige ronde stippen toegepast.

I

4.1 .2.1. Deeltjesgenerator

I

Vastleggen ronde stippen op binaire plaatjes.

r----~-x

I

De pixel krijgt de waarde 1 toegekent als de positie van zijn zwaartepunt binnen de cirkel valt. (zie onderstaande fig.)

y 6.x

I

pixel (x,y)

/

(x,y) pos

I

Het criterium van deze toekenning is dat het grootste gedeelte van de oppervlakte

I

van de pixel binnen de cirkel valt.

10

I

(15)

I

Toch schuilt er bij het nemen van dit criterium nog een klein foutje. Voor objecten met een kleine straal kan het zwaartepunt van de pixel binnen de cirkel vallen, maar meer dan 50% van zijn oppervlakte valt nog buiten de cirkel:

I

pixel object

I

De oppervlakte van a en van b tezamen is groter dan van c, zodat meer dan 50% van de oppervlakte buiten de cirkel valt. Echter in dit onderzoek is deze fout niet van belang. Cirkelvormige objecten zijn slechts een ruwe benadering van wat zich in de werkelijkheid afspeelt: slibvlokken hebben vaak een totaal andere vorm. Bovendien geeft de correctie op deze fout een aanzienlijk meer rekenwerk. Verder zijn we geinterreseerd in objecten met een oppervlakte van ongeveer 100 pixeleenheden, overeenkomend met een straal van ongeveer 5 á 6 pixeleenheden.

I

11

I

(16)

Kruiscorrelatie is verreweg het meest tijdrovende rekenproces in de snelheidsanalyse. Het is een tweedimensionale kruiscorrelatie uit twee 5122 plaatjes. In de ontwikkeling van het programma is een truc bedacht om met aanzienlijk minder rekentijd ongeveer dezelfde resultaten te verkrijgen.

In discrete vorm wordt de tweedimensionale kruiscorrelatie beschreven: Rxy[k] [1]=I:I:X[i] [j] *Y[i+k] [j+l]

i j

I

,

I

4.2 Kruiscorrelatie

I

Rxy[k][1] is een tweedimensionale kruiscorrelatiefunctie, wordt ook wel een tweedimensionale correlogram, kortweg correlogram genoemd. X[i][j] en Y[i

+

k][j

+

I] zijn de twee 5122 plaatjes.

Echter, kruiscorreleren moet geschieden binnen het beeldveld van X en van Y. Hierdoor loopt de range van k en van I van -128 tot en met 127 en de range van i en van j van 128 tot en met 384. De uitdrukking van Rxy wordt dus:

I

k,1 : -128,-127, ... ,127

I

384 384 Rxy[k] [1] =

E

L

Xli] [j] Y[i+k] [j+l] i=128 j=128

Het correlogram kan worden voorgesteld als een 2562 plaatje met de oorsprong, de (0,0) coördinaat, in het midden. Het realiseren van dit correlogram vereist op z'n minst 2564 == 4.3 miljard vermenigvuldigingen. De beelden X en Y zijn binaire

images; de waarden van de pixels bedragen slechts 0 of 1. Dus de uitkomst van elke vermenigvuldiging is eveneens 0 of 1:

(i) 1 als X[i][j]

=

1 én Y[i][j]

=

1; (ii) 0 overig.

Bij voorwaarde (i) wordt bij de waarde van Rxy[k][l] 1 opgeteld. Hierdoor ben je de vermenigvuldigingsinstructie in het geheel kwijt. Van de tijdrovendheid ben je echter nog niet los. Met de HP9000 is voor dit rekenproces 2% uur voor nodig. Terugbrengen rekentijd

Gezocht werd naar de mogelijkheden om de 2 Y2 uur durende rekentijd terug te brengen. Eerst is gekeken naar de pack en unpack procedure en daarna naar de datareductie door het kleiner maken van de binaire plaatjes.

I

12

I

(17)

---- ---

--I

I

Pack en Unpack

De afzonderlijke bits in de binaire plaatjes worden gepackt, in de register gebracht waar deze met de logische operator worden verwerkt en tenslotte weer worden ge-unpackt. HP9000 is voorzien van 32 bits registers. In een lus kan 32 bits tergelijkertijd ofwel parallel worden verwerkt. Dit kan een tijdwinst opleveren met een factor 32. Hierdoor zou de rekentijd van 2 Y2uur kunnen worden teruggebracht tot 5 minuten. Echter na overleg met A. vld Meer TNO - TPD (de belangrijkste persoon bij de ontwikkeling van TCLi softwarepakket) blijkt dat het uitvoeren van de kruiscorrelatie op deze manier leidt tot een zeer gecompliceerd c-programma, voldoende om deze winst geheel teniet te doen.

Datareductie binaire plaatjes

De rekentijd kan ook kleiner worden gemaakt door eerst het binaire plaatje kleiner te maken met een factor 2 of 4 en dan pas te gaan kruiscorreleren. De tijdwinst hierbij bedraagt respectievelijk 16 en 256 maal. Het bezwaar hierin is dat door het kleiner maken van de binaire plaatjes informatieverlies optreedt. Getest is hoe groot de invloed van dit informatieverlies is op het correlogram en het daaruitvolgende referentiegebied, waaruit blijkt of het op deze manier wel geschikt is voor het bepalen van de snelheidsverdeling.

a. Kleiner maken met factor 2

Door het kleiner maken van de binaire plaatjes van 5122 naar2562 wordt de grootte

van de ranges van i,j,k en I een factor 2 kleiner. In plaats van -128 tot en met 127 worden deze voor k en I -64 tot en met 63 en voor i en j in plaats van 128 tot en met 384, 64 tot en met 192. De uitdrukking voor het correlogram wordt:

192 192

Rxy[k] [1]

=

L

E

x[i] [jfY[i+k] [j+1]

i=64 j=64 k,1 : -64,:-63, ...,63

I

Het correlogram kan worden voorgesteld als een 1282 plaatje.

Om het weer terug te brengen naar een 2562 plaatje wordt het opgeblazen met een

factor 2 (2 times blow up) en voor de daartussenliggende punten geïnterpoleerd. Door het kleiner worden van alle 4 de ranges met een factor 2 is het totaal aantal vermenigvuldigingsinstructies of het aantal uitvoeringen (i) of (ii) (page 12) gereduceerd met een factor 24

=

16. De rekentijd wordt daardoor in plaats van 2 Y2

uur naar verwachting slechts ongeveer 10 minuten.

13

(18)

I

b. Kleiner maken met factor 4

I

Hier worden de binaire plaatjes met een factor 4 kleiner gemaakt van 5122

naartZê", De grootte van de ranges van i,j,k en I zijn eveneens een factor 4 kleiner.

I

In plaats van -128 tot en met 127 worden deze voor k en I -32 tot en met 31 en

voor i en j in plaats van 128 tot en met 384, 32 tot en met 96. De uitdrukking voor

I

het correlogram wordt:

96 96

I

Rxy[k] [1] =

L L

X[i] [j]y[i+k] [j+1]

i=32 j=32

I

k,1 : -32,-31, ... ,31

Invloed in forma tie verlies door reductie van de binaire plaatjes op de kruiscorrelatie.

I

Het correlogram kan worden voorgesteld als een 642 plaatje.

Dit wordt teruggebracht naar een 2562 plaatje door het op te blazen (blowed up) met een factor 4 en vervolgens de daartussenliggende punten te interpoleren. Door het kleiner worden van alle 4 de ranges met een factor 4 is het totaal aantal vermenigvuldigingsinstructies of het aantal uitvoeringen (i) of (ii) (page 12) gereduceerd met een factor 44= 256. Dus de rekentijd wordt in plaats van 2% uur

naar schatting slechts circa 35 seconden!

De vraag luidt hoe groot het informatieverlies op het correlogram is en uiteindelijk op het referentie gebied. Om tot antwoord te komen zijn testen uitgevoerd. Bepaald werd wat voor verschillen dit zal opleveren op de correlogrammen en uiteindelijk op het referentiegebied . Deze testen zijn steekproefsgewijs uitgevoerd. Het systeem bevat 500 deeltjes. Voor radius is gekozen 2, 5 en 10 pixeleenheden. Voor meer duidelijkheid zijn ook de doorsneden in de x en in de y richting van het correlogram bepaald. Het snijpunt van de doorsneelijnen ligt op de top van het correlogram van de kruiscorrelatie zonder datareductie.

I

14

I

(19)

I

Tests

Reductie met factor 2

I

De rekentijd bedraagt 9 min 44 seconden, afgerond 10 minuten, wat is verwacht. Op pseudo grijze plaatjes van de correlogrammen (fig. 2.01 tot en met fig. 2.18) is te zien dat door reductie nauwelijks verschillen optreden. De verschillen in de posities van de top van het correlogram weergegeven in de x en in de y coördinaat bedragen hooguit 1 pixeleenheid.

In de doorsnedeplots (fig 3.01 tot en met 3.12) zijn de verschillen duidelijker waarneembaar. Deze zijn het grootst rond de flanken van de piek in het correlogram. Daar zijn de hellingen het steilst, waardoor bij zeer geringe verplaatsing van de piek grote verschillen optreden.

Voor objecten met een radius van 2 pixeleenheden variëren de verschillen tussen de -4 en 4 met hier en daar een uitschieter naar -8 en 8. En bij de flanken:

opgelegde snelheidsverdeling x - richting linkerflank rechterflank y - richting linkerflank rechterflank

I

symmetrisch asymmetrisch - 37 - 41 60 61 - 19 - 40 65 18

I

Het schijnt dat de verschillen kleiner worden als de radius van de objecten toeneemt.

Voor objecten met een radius van 5 pixeleenheden variëren de verschillen tussen de -2 en 2 met uitschieters naar -3 en 3. En bij de flanken:

I

opgelegde snelheidsverdeling x - richting linkerflank rechterflank y - richting linkerflank rechterflank symmetrisch asymmetrisch - 9 27 - 16 19 - 9 28 - 22 14

I

tussen de -1 en 1 met uitschieters naar -3 en 3. En bij de flanken:En voor objecten met een radius van 10 pixeleenheden variëren de verschillen

I

opgelegde snelheidsverdeling x - richting linkerflank rechterflank y - richting linkerflank rechterflank symmetrisch asymmetrisch - 3 -8 13 9 - 1 -9 13 7

Bij het groter worden van de straal van de objecten neemt de breedte van de piek in het correlogram eveneens toe, hierdoor neemt de steilheid van de flanken af, zodat de verschillen in de flanken afnemen.

I

15

(20)

Reductie met een factor 4

De rekentijd nodig voor de kruiscorrelatie bedraagt naar verwachting 37 s ! Een zeer groot verschil met 2% uur!

Echter reductie met een factor 4 geeft meer informatieverlies dan met een factor 2 en naar verwachting leidt het tot grotere verschillen ten opzichte van de normale kruiscorrelatie.

Op pseudo grijze plaatjes zijn de verschillen nog nauwelijks zichtbaar, behalve bij objecten met een radius van 2 pixeleenheden. De correlatieruis ziet er scherper uit. Bij opgelegde asymmetrische snelheidsverdeling is de komeetstaart in plaats van naar beneden rechts, geheel naar beneden gericht. De verschillen in de x als in de y coödinaten, die de positie van de top van het correlogram weergeven variëren tussen 1á 3 pixeleenheden. Ook in de doorsnedeplots zijn de verschillen duidelijk groter. Voor objecten met een radius van 2 pixeleenheden toont correlatieruis grove pieken met een waarde van 32 bij opgelegde symmetrische snelheidsverdeling en 51 bij asymmetrische snelheidsverdeling. Deze leiden in de verschillen tot uitschieters van rond de 30 bij symmetrische en rond de 40 á 50 bij asymmetrische snelheidsverdeling. Daarbuiten variëren de verschillen tussen 0 en 10, waarin vaak nullen voorkomen. Bij de flanken zijn de verschillen:

opgelegde snelheidsverdeling x - richting linkerflank rechterflank y - richting linkerflank rechterflank symmetrisch asymmetrisch 128 237 112 - 82 - 117 - 153

Bij asymmetrische snelheidsverdeling (y - richting) valt de piek van de factor 4 reductie binnenin de piek van het normale correlogram. (fig 1.04) zodat de afwijking een naar beneden gerichte piek is met een dal ter waarde van -153. Ook hier schijnt dat als de diameter van de objecten groter worden, de verschillen ten gevolge van de datareductie zullen gaan afnemen.

Voor objecten met een radius van 5 pixeleenheden variëren de verschillen tussen de -10 en 12 met hier en daar een uitschieter van -49 of 49.

Bij de flanken zijn de verschillen:

opgelegde snelheidsverdeling x - richting linkerflank rechterflank y - richting linkerflank rechterflank symmetrisch asymmetrisch - 48 - 33 - 49 - 43 49 40 49 60 16

I

(21)

I

En voor objecten met een radius van 10 pixeleenheden variëren de verschillen tussen de -3 en 7, en bij de flanken:

opgelegde snelheidsverdeling

x - richting linkerflank rechterflank y - richting linkerflank rechterflank

I

symmetrisch asymmetrisch - 21 - 24 25 27 - 21 - 28 23 22

I

Ook hier neemt de afwijking af ten gevolge van het het groter worden van de piekbreedte dat weer het gevolg is van het groter worden van de objecten.

I

17

(22)

I

Invloed datareduetie op de ligging van de top van het eorrelogram

I

De eoordinaten van de ligging van de top en de invloed daarop is weergegeven in

de twee onderstaande tabellen:

I

Position tops:

I

radius normal 2 times 4 times

dots distribution cross data data

correlation reduction reduction 2 symmetrie (0,33) (1,34) (1,33) 2 asymmetrie (2,6) (3,6) (5,5) 5 symmetrie (0,33) (0,34) (1,34) 5 asymmetrie (2,6) (3,7) (5,9) 10 symmetrie (0,33) (0,33) (1,34) 10 asymmetrie (2,6) (3,7) (5,9)

radius 2 times 4 times

dots distribution data data

reduction reduction 2 symmetrie (1 ,1) (1,0) 2 asymmetrie (1,0) (3,-1 ) 5 symmetrie (0,1) (1,1) 5 asymmetrie (1 ,1) (3,3) 10 symmetrie (0,0) (1,1) 10 asymmetrie (1 ,1) (3,3)

,

I

Differenee with regard to normal eorrelogram:

De invloed van de tweevoudige datareduetie op de ligging van de top van het

eorrelogram is hooguit 1 pixel groot. Van viervoudige datareduetie tot 3 pixels.

I

N=500 dt=O.3

The coordinates of the top of the correlogram are expressed in pixel units.

18

I

(23)

I

4.3 Selecteren deeltje uit tweede opname

Voorselectie

I

De voorselectie gebeurt aandehand van hetreferentiegebied. Het referentiegebied wordt bepaald met iteratieve threshold (zie eerste rapport):

EE

thre$hold

I

corr.lo~rQm binary image

with refer.ne. area

Opname 1 en opname 2 bevatten respectievelijk n, en n2 aantal deeltjes:

I

+

x

+ +

+

_x

I

+

x

+

x

+

x

+

i

x

+

x x

1 +

.

+

x

+

x

+

'.

~

_

*~_

+

x

+

__tfeference~ X

+

X

+

-.- oero

x

+!

x

! x

+

. x

+

x

+

x x

+

x x

+

x

+

x x x

+

x x

I

De (0,0) coördinaat van het correlogram is verschoven en ligt op het zwaartepunt

van deeltje i (opname 1). Van opname 2 liggen de zwaartepunten van deeltjes k,1

en m in de referentie gebied en de overigen niet, zodat k, I en m worden

geselecteerd. Met de tweede selectie wordt een van deze drie deeltjes gekozen.

I

19

I

(24)

I

Selectie

I

Van k, I en m komt de ligging van het zwaartepunt van I overeen met de hoogste

I

kruiscorrelatie, zodat deeltje I wordt gekozen.

k

m

I

top ofthe c oer etoqrcm

I

Meer dan een referentiegebied .

Bij situaties waar sprake is van gering aantal deeltjes kan thresholding meer dan één referentiegebied opleveren als gevolg van correlatieruis. Dit kan invloed hebben op de bepaling van snelheidsverdeling zodat het aantal gereduceerd moet worden tot één, het meest waarschijnlijke.

Het selecteren van dat ene is gebaseerd op de ligging van de top van het correlogram: het referentiegebied waarin die top bevindt, zie onderstaande fig.:

~

(1,

G

threshold .ot selection

0

o 0

I

correlogrom number of one

reference oreos referenee 0rea

I

20

I

(25)

I

· I

I

5.

TESTS

I

Getest is de juistheid van dit programma voor het bepalen van de gemiddelde en de sigma van de snelheidsverdeling. Bekeken wordt de afwijkingen van de gemiddelde en de sigmawaarden ten opzichte van de referentie. Gebruik wordt gemaakt van verschillende tijden tussen de twee opnamen. Wordt deze tijd zeer klein, dan zullen naar verwachting digitalisatiefouten gaan optreden. En wordt deze zeer groot dan neemt de kans toe dat het deeltje buiten het referentiegebied valt of dat er een verkeerd deeltje wordt geselecteerd. Dit betekent dat daartussen een optimum wordt verwacht. Ook kan worden verwacht dat de afwijking toeneemt, als de diameter van de desbetreffende deeltjes groter wordt, of het aantal toeneemt. In beide gevallen is de oorzaak overlap.

In deze test zijn twee verschillende snelheidsverdelingen opgelegd: symmetrisch en asymmetrisch.

In het gehele systeem worden 200, 500 en 1000 deeltjes gegenereerd. Voor asymmetrische verdeling i.p.v 1000932, omdat ik op het moment van testen met de biedtafelgegevens niet verder kwam omdat de uitzendingen zijn gestopt. Omdat het cameraveld een kwart van het hele systeem bestrijkt zullen naar verwachting in het cameraveld respectievelijk 25, 50, 125 en 250 deeltjes worden verwacht. De objecten zijn cirkelvormig met één gegeven radius. Ze zijn gelijk van vorm en van oppervlakte. De radiussen die hier zijn toegepast zijn 2, 5 en 10pixeleenheden. De middelste radius van vijf pixeleenheden komt overeen met een oppervlakte van circa 80 pixeleenheden, ongeveer de oppervlakte van deeltjes dat in praktijk het meest gangbare is. De tijd tussen de twee opnamen: Llt loopt met een stapgrootte van 0.01 van 0.01 tot en met 0.50 dimensieloze of pixel eenheden. Omdat de gemiddelde verplaatsing van de deeltjes gedurende de Llt interessanter is, is onder aan de tijdas van de grafische interpretatie de gemiddelde verplaatsing uitgedrukt in pixeleenheden weergegeven. In deze test is gekeken naar de y-component van de snelheid. De programmatuur die hier is getest is met twee en 4-voudige datareductie, daar met de HP9000 voor iedere Llt stap van 0.01 een rekenwerk is vereist van ongeveer 2 % uur, de duur van de kruiscorrelatie.

I

21

I

(26)

6.

RESULTATEN

De resultaten zijn grafisch weergegeven in fig.: 4.1 tot en met fig.: 4.8. 6.1 Symmetrische snelheidsverdeling (fig.: 4.1 - fig.: 4.5)

Gemiddelden

Vrijwel alle punten liggen netjes op de lijn van de 0%. Slechts bij 4-voudige datareductie en een radius: r van 10 pixeleenheden liggen slechts 2 punten bij de 5%. Er is een uitzondering: bij r

=

2 en wanneer het totaal aantal deeltjes in het systeem 200 (N = 200) bedraagt vallen welgeteld 6 punten ruim buiten de 0% lijn waarvan minstens 3 op -50% of lager komen te liggen.

Sigma

Behalve bij N = 200 neemt de afwijking toe met het groter worden van de radius.

Wordt het aantal deeltjes in het systeem groter, dan zal, als N gaat van 200 naar 500, !l.creveneens toenemen. Neemt N verder toe van 500 naar 1000 dan is van verdere toename van !l.crgeen sprake meer. Bij r = 2 blijkt toename van N nauwelijks invloed uit te oefenen op !l.G, hooguit 3 punten bij !l.t = 0.1O.

Bij r = 10 valt duidelijk op, dat er punten aanwezig zijn die een vloeiende lijn volgen, vervolgens een discontinu een sprong maken en daarna weer zo'n lijn gaan volgen. Bij r = 2 is het duidelijkst te zien, dat !l.crvanaf Il.pixel = 15 zeer snel omhoog gaat schieten.

Tenslotte als deze resultaten waarbij 4-voudige datareductie is toegepast worden vergeleken met 2-voudige datareductie, dan valt nauwelijks enig verschil te bespeuren. Alleen bij N = 200 en 4-voudige datareductie valt op te merken, dat meer punten beneden de 0% lijn liggen.

6.2 Asymmetrische snelheidsverdeling (fig.: 4.6 - fig.: 4.8)

Gemiddelde waarden

Vergelijking van de resultaten van twee- en van 4-voudige datareductie levert nauwelijks enig verschil op. Bij 4-voudige datareductie vallen bij r = 2 meer punten boven de 0% lijn uit, maar blijven wel onder de 20% liggen. En bij r = 10 en N = 200 lopen de punten wat stijler omhoog, wanneer de verplaatsing kleiner of gelijk is aan 2 pixels.

Bij r=5 en 10 lopen de punten van ongeveer 15 á 20 naar -10 á -15%. Bij r=2 wanneer 2-voudige datareductie is toegepast van 0 naar -10%.

De nul-doorgang toont in !l.t of !l.pixels een shift naar rechts, wanneer de radius

groter wordt. 22

I

(27)

I

sigma

Afgezien van enkele punten bij r

=

2, waarbij 4-voudige datareductie is toegepast,

is er nauwelijks enig verschil te bespeuren in vergelijking met deze waarbij 2-voudige datareductie is toegepast.

Bij r = 2 en 5 loopt .10' af van 50 naar -50%. Ook daar valt op te merken, dat er een shift te bestaat naar rechts wat betreft de ligging van de punten op de 0% lijn. Bij r

=

2 bij een verplaatsing van 2 á 3 pixels en bij r

=

5 van 4 á 5 pixels. Bij r

=

10,

afgezien N=200, tussen 5 á 7 pixels.

I

Eveneens bij de sigma vertoont de nuldoorgang een shift naar rechts bij groter worden van de straal van de objecten.

I

23

I

(28)

24

I

'

I

7. DISCUSSIE

7.1 Mogelllke oorzaken van optredende afwijkingen in gemiddelde en sigma. 7.1.1. Digltalisatiefouten

Het optreden van digitalisatiefouten is het duidelijkst waarneembaar op de plots van

l1a versus I1t (l1pixel) van symmetrische snelheidsverdeling, waar de radius van de

deeltjes twee pixeleenheden bedraagt. Dit is nog het duidelijkst bij een experiment met 2-voudige datareductie (fig.:4.4). Bij een verplaatsing van 25 pixeleenheden of meer liggen de punten op de 0% lijn, zo rond de -5 en 5%. Beneden de 25 pixeleenheden begint de lijn op te lopen, vervolgens steil omhoog en benadert de lijn I1t of I1pixel=Oasymtotisch. Bij r=5 en 10 verlopen de punten, wanneer het aantal deeltjes in het systeem gering is (N=200) ongeveer hetzelfde als bij r=2. Alleen in plaats van ongeveer 25 pixeleenheden gaat het pas omhoog vanaf 20 pixeleeneden. Bij experimenten met asymmetrische snelheidsverdeling is dit moelijker te zien. Aan het verloop van de curves zou kunnen worden afgeleid dat het ook daar vanaf I1pixel=25 steil omhoog begint te lopen. Deze punten vallen ook lager uit, wat een ander oorzaak kan hebben. Voor meer duidelijkheid zou deze test moeten worden herhaald, maar dan moet I1tlopen tot 2.0 i.p.v.tot 0.5, zodat I1pixelsloopt tot 60 i.p.v.

tot 15.

7.1.2. Overiap

De plots van l1a versus I1twaar r=5 of 10 bij symmetrische snelheidsverdeling toont duidelijke verschillen met deze van r=2. Bij N=500 of meer vallen duidelijk meer punten ver boven de 0% lijn. Ook in het gebied van verplaatsingen van ~5 pixeleenheden. De punten komen nog hoger te liggen als de radius groter wordt van 5 naar 10 pixeleenheden. Dit veschijnsel treedt eveneens op bij asymmetrische snelheidsverdeling, maar dan minder duidelijk. Verder is het bij r=10 te zien dat er punten zijn, die een vloeiende lijn volgen, vervolgens via een discontinu sprong een volgende lijn.

Door de grotere radius neemt de kans toe dat deeltjes elkaar overlappen. De figuren 1.01 tot en met 1.03 geven dit duidelijk weer. Bij r=2 is daar nauwelijks een dat overlapt, bij r=5 8 deeltjes en bij r=10 meer dan 40. Ook neemt de kans op overlap toe als het aantal deeltjes toeneemt. Dit blijkt duidelijk in de plots terug te vinden als N=500 wordt vergeleken met N=200.Echter aanzienlijk minder duidelijk wordt het als N=500 wordt vergeleken met N=1000 (sym. verdeling) of N=932 (asym. verdeling).

Verwachting is dat l1a nog groter zou worden. Dit is wel het geval bij r=10 als de symmetrische verdeling wordt opgelegd. Bij een asymmetrische verdeling is het omgekeerd. Bij r=5 valt er geen enkel verschil uit te halen.

Dat sommige punten een vloeiende lijn volgen, gevolgd door een discontinu sprong en daarna weer verder een lijn volgen zou te wijten kunnen zijn aan het feit dat deeltjes eerst los van elkaar zijn, vervolgens overlappen om daarna weer los van elkaar. Dit heeft als gevolg dat 2 zwaartepunten ineens een zwaartepunt wordt of omgekeerd. Dit heeft natuurlijk ook consequenties op de positie van deze punten. Wat

in discontinue sprongen is terug te vinden.Het gaat hier om sommige deeltjes in een groter aantal.

(29)

I

7.1.3 Buiten het referentiegebied vallen van deeltjes

De verschijnselen ten gevolge van buiten het referentiegebied vallen van de deeltjes,

is het best waar te nemen in L\crversusL\t, waarbij asymmetrische snelheidsverdeling is toegepast. Bij een grotere L\t lopen de deeltjes, die in de staart van de asymmetrische snelheidsverdeling liggen (met grotere snelheid), de kans op t=t2 buiten het referentiegebied te raken. Hoe groter de radius, des te groter het referentiegebied. De zwaartepunten van de deeltjes die anders erbuiten vallen,vallen nu binnen. Dit is te zien bij het rechtergedeelte van deL\t-as. Veel punten beneden de 0% lijn benaderen deze lijn dichter, naarmate de radius groter wordt. Deze trend is waar te nemen bij N=500 en N=1000, maar niet bij N=200.

I

(30)

8. CONCLUSIES

Kruiscorrelatie

Kruiscorreleren van twee binaire plaatjes ter grootte van 5122 pixels vraagt met de HP9000 een rekentijd van 2 % uur en is terug te brengen door datareductie toe te passen op twee binaire plaatjes.

Twee- en viervoudige datareductie brengt respectievelijk de rekentijd terug naar 10 minuten en slechts ongeveer 37 seconden.

Informatieverlies ten gevolge van tweevoudige datareductie geeft een verschil in de positie van de top van het correlogram weer van hooguit 1 pixel en van een viervoudige datareductie hooguit 3 pixels, zodat voor het bepalen van de gemiddelde snelheid direct uit de positie van de top van het correlogram, 2-voudige datareductie acceptabel is.

Gemiddelde en spreiding snelheidsverdeling.

symmetrische snelheidsverdeling:

Voor bepaling van de gemiddelde snelheid is de geteste programmatuur goed bruikbaar. De fout ligt binnen de 5%.

Wanneer kleine deeltjes met een diameter van 4 pixels worden bepaald, is het programma met 4-voudige datareductie ten gevolge van informatieverlies niet meer betrouwbaar.

Voor de bepaling van de spreiding is dit programma tamelijk goed tot goed bruikbaar. De fout ligt binnen de 10% als de verplaatsing van de deeltjes 20 pixels of meer bedraagt en de deeltjes klein zijn met een diameter van ongeveer 4 pixeleenheden. Bij verplaatsing van 20 pixelseenheden of minder gaan digitalisatiefouten optreden. Voor grotere deeltjes met een diameter van 10 tot 20 pixels is dit programma tamelijk tot goed bruikbaar als het aantal deeltjes op het scherm ongeveer 50 bedraagt. Bij een grotere aantal gaat overlap een rol spelen.

asymmetrische snelheidsverdeling

Voor de gemiddelde snelheid is dit programma tamelijk bruikbaar, als de verplaatsing 2 pixels of meer bedraagt. De fout bevindt zich binnen de 15%. Echter voor de sigma is dit programma nauwelijks bruikbaar

26

I

(31)

I

Dit kan als volgt worden samengevat:

Dit programma is goed toepasbaar voor bepaling van de gemiddelde snelheid als de snelheidsverdeling symmetrisch is en tamelijk goed bruikbaar bij een asymmetrische verdeling. Voor de spreiding kan dit programma tamelijk goed worden gebruikt als de snelheidsverdeling symmetrisch is en waarbij de verplaatsing 20 of meer pixels bedraagt. Bij een asymmetrische snelheidsverdeling is de bruikbaarheid te beperkt. Voor kleine deeltjes is dit programma met 4-voudige datareductie niet betrouwbaar, pas wanneer de diameter 10 pixels of meer bedraagt wordt het betrouwbaar.

I

27

I

(32)

9. VERDER ONDERZOEK

In het gebied van de snelheidsverdelingsonderzoek is nog aanzienlijk veel experimenteel werk te verzetten.

Er zijn andere methoden denkbaar, met name de gemiddelde snelheid direct uit de positie van de top van het correlogram en vergelijken met deze methode. Andere opgelegde snelheden waar eddys in voorkomen.

Als volgende stap het ontwikkelen van de programmatuur, gebaseerd op een combinatie van methoden, zoals kruiscorrelatie - kortste afstand.

9.1 Gemiddelde snelheid direct uit de top van het correlogram.

Met de kruiscorrelatie bestaat er een ander mogelijkheid om de gemiddelde snelheid te bepalen, direct uit de positie van de top van het correlogram. Door de x- en de y-coordinaat van de top te gaan delen door de tijd tussen de twee opnamen volgt direct de gemiddelde snelheid.

Voor sigma wordt dit moeilijker. De piekbreedte hangt samen met sigma, maar ook met de grootte van de deeltjes. De piek wordt breder bij toenemende grootte van de deeltjes. In fig.: 3.01 tot en met 3.12 is dit duidelijk te zien.

9.2 Kortste afstand

Een geheel ander methode. De keus van het deeltje in opname t =t2 is gebaseerd op de kortste afstand van het desbetreffende deeltje in opname t =t1. Rijkswaterstaat ontwikkelt op dit moment de programmatuur dat op dit principe is gebaseerd. Een handicap bij Rijkswaterstaat is, dat hun programmatuur niet is getest met gebruik van opgelegde snelheden.

9.3 Andere opgelegde snelheden

Tot nu toe is uitgegaan van twee verschillende snelheden: symmetrische en asymmetrische. Eveneens is interresant om eddys in te voeren. Wervels komt vaak in de praktijk voor.

Wat betreft de asymmetrische snelheidsverdeling. Omdat het de eerste tests zijn werd het biedtafelmodel toegepast. Beter is een gedefinieerder asymmetrische snelheidsmodel toe te passen met een parameter voor de asymmetrie.

9.4 Patroonherkenning

Selectie op basis van patroonherkenning (methode( 1), zie pag.: 3) of combinatie kruiscorrelatie - patroonherkenning (methode(2)). Voor het genereren van deeltjes kan in de richting van fractals worden gedacht.

28

I

(33)

I

9.5 Combinaties van methoden.

Een stap verder lijkt mij het ontwikkelen van programmatuur op basis van de combinatie kruiscorrelatie-kortste afstand. Het programmatuur zal kruiscorrelatie toepassen wanneer er sprake is van veel deeltjes en kortste afstand bij weinig deeltjes. Bij weinig deeltjes treedt ten gevolge van een statistisch effect sterke correlatieruis op, wat deze methode onbetrouwbaar maakt. kortste afstand daarentegen is door minder aanwezigheid van concurrerende deeltjes de kans kleiner dat er een verkeerde deeltje wordt geselecteerd. Er zal er een criterium bevinden voor de keuze van welk van de twee methoden, een interressant iets voor dit verdere onderzoek.

I

,

J

I

,

I

29

(34)

10. LITERATUUR

I

,

I

[1] - A Kelly and I Pohl,

An introduction to programming in C

The Benjamin/Cummings Publishing Company,lnc (1984) [2] - TU Delft, UvA, TPD

Image Processing for Industrial Applications, Introductory Course: 28-29 Oktober 1991 [3] - manuals of Multihouse TSI:

TCl-lmage User's Manual, Part One TCl-lmage User's Manual, Part Two TCl-lmage Programmer's Manual [4] - J. de Graaft and R.E. Slot,

Determining particle size distributions from video images by use of image processing (in preface)

30

I

(35)

I

.

I

(36)

I

.

I

_{APPENDIX A1}

I

FLOW

I

_DIAGRAMS

I

(37)

I

r I

I

II I

I

II I

I

II I I

I

I I I

I

I I I

I

I I I

I

II I

I

velocity

determ

i

n

i

ng

program

generoting relerence relerence generoting porticles porticles ---- -determining positions correlotion cross

t

2 dimensionol correlogrom r---~---, ~---~ criterium

H

t

reoble ~ not reoble concel iterotive !reshold !rom correlogrom

t

reofereneereo{s) determining one number of I. ~ reference ~ oreo{s)

t

more thon one select one ol the ref. oreos I

t

lirst selection (.)

porticle which position

corresponds with

highest cross correlotion

(

..

) movements of the porticles I I I I I I I

T

L- --J

t

porticles in rel. oreo second selection (

..

)

-determining movements determining velocities A1

FIG.: A 1

velocity determining progrom itselfs

V

_I I I I I velocities ol the porticles - _j

(38)

FIG.:

A2

measurement of

6 mean and 6(5 in velocity distribution

determined by the program

with regard to the reference

reference determined by the program determining mean and sigma determining mean and sigma determining difference with the reference difference

of mean and sigma

with the reference

( ~mean,t:.(Ï) A2

I

(39)

I

of

FIG.: A3

I

determining movements the particles with cross correlation use of method

I

positions of

porticles ot I=tl

I

positions of portici•• ot t-t2

I

correlogrom

I

(program moverxy1) reference area determining number of pixels of reference area

(pixel. of aquol 1) number pixel. eq determining coordinates of pixels equol 1 coordinot pixel. eq tronslating ret. areo with coordinotes of centerpoints of porticles ol 1-11 translate select centerpoints of porticles (1=12) with ref. area 1 centerpoi in refer ••tect centerpoint.

whleh ploee corresponds

"ith high.st ero •• correlotion seleeted determining movement ~ movement ot o porticIe A3 of uol esof uol 1 d ref. creo first selection nts of perticles ence oreo second selection porticle

(40)

FIG.:

linked list structure of object numbers at t=t1

and object numbers at t=t2 after first selection with use of ref. area

head (in program moverxyl)

toil 1

tail 2

no corresponding elements at t=t2

Vertical list: particles at t=t 1

Horizontal list: particles at t=t2 (in ref. area)

ELEMENT OF PARTICLE AT t=t1

ELEMENT OF PARTICLE AT t=t2

EJ

object number

link: pointing to the next element of porticle ot t-t2 object number link: pointing to the first element of porticle ot t-t2 link: pointing to the next element of porticle ot t=t1 A4

A4

I

(41)

I

creating referenee

symmetrie velocity distribution

(program genpv.c) rondom generator

FIG.: A5

I I I I

i

_I

/

I I I I I I I heads or tails model

velo city inx and y components

(x and y direction)

I

rondom values in the

range: 0...2" - 1 --f--- --t: I I I I I I I I IL...--"""T""----' k heods I and n-k toils I I I I determining max. value of the rondom numbers max. value of rondom numbers random volues oor 1 normolizing rondom numbers

with max. value

toss n cains rondom numbers between 0 and 1 multiply normolized rondom numbers with 1000 from number of heads follows velocity

positions of particles in system ot t=t 1

velocity size and direction

velocity converting

into x and y

components

determining positions

of porticles in system ot t=t2

position upperleft corner

camera field in system

\ position of particles

in system ot t=t2 determining

r

ositions o porticles in 1-... ;.--camera field ot t=t 1 and t=t2 select porticles if these are in camera field ot 1=11 and l=t2

velocities of the selected porticles (cenierpoints) • the reference

A5 time bet ween two exposures random volue: 14 o to 2 -1 --> 1(head) 14 15 2 to 2 -1 --> 0 (taiI)

(42)

ereating referenee asymmetrie veloeity distribution

(program genpvaa.c)

rondom generator lees biedtafel gegevens

FIG.:

A6

biedtafel gegevens

5 kolommen:4 bods + prijs ortikel

- - --, - -splits biedtafel gegevens 1e twee kolommen normaliseren voor x-richting

I

I I I I I I I I I I I I I I

7

-determining max. value of the rondom numbers max.value of

ondom numbers 3·vooren ."s:richtingkolom

normolizing rondom numbers with max.value noor prijs artikel 1· en 2e bods genormaliseerd 3·en •• bads genormaliseerd andam numbers between0 and1 normaliseren noor prijs artikel mulliply normolized rondom numbers with 1000

vermenigvuldigen mei 10 en delen door 2

positions of porticles n system ot t=tl

---

I-

--biedtafel model

vermenigvuldigen met 25 en delen door 2

snelheden in x-richting (pixelunits)

determiningpositions of porticles in system ot t=t2 time between two exposures

,

position of porticles in syslem ot l=t2 determining fositions o porticiesin t- ...;.- -camera field ot t=t 1 andt=t2 select porticles if these arein camerafield ol I~tl ond 1=12

veloeitiesof theselecled parlieles (centerpoints)•the reference

A6

position upperleftcorner comera fieldin system snelhedenin y-richting (pixelunits)

I

(43)

I

(1) position of upper left corner af camera field

(2) radius

(3) velocities (from genpv.c or genpvaa.c)

(4) time betweentwo expasures

FIG.: A7

I

Generating

partieles

I

(genpartdc.c)

I

positions of centerpoints in system ot t=t 1

I

~ determining

positions of ~

centerpoints

in system ot t=t2

--pasitions of cent

in system ot t=t

(1)

determining _determin_i_ng

-porticies in _port_i_cles _i_n

camera field (2) camera field

-,

(2)

I

(3)

I

(4) erpoints 2

I

(1 )

I

porticles in camera field

ot t=ll

porticles in camera field

ot t=t2

I

A7

I

(44)

FIG.: AB

Histogram of velocity distribution with

mean and sigma

veloeities of the partieles

reference or

determined by the progrom

determining mean and ( sigma meon ond sig making histogram

t

see olso fig. 02) ma number of particles closs histogram with numbers of particles value of meon volue of sigma size of closs AB

I

(45)

I

_{APPENDIX A2}

I

DESCRIPTION

_{THE PROGRAMS}

OF

I

(46)

I

A1. MAIN PROGRAM (CRITOP2y.tip)

A1.1. Generating reference First is used:

copyt ut

copynNNN

t,ttt: time between two exposures (At);

nn, NNN: total number of particles (whole system).

These two commands are needed because the program that generates positions and veloeities (the reference) does not accept the parameters tand nn.

The reference is generated by use of the module genpv or genpvaa, it generates positions (centerpoints) of partieles and its veloeities. The module genpv generates symmetrie velocity distributionin x and y direction and genpvaa asymmetrie velocity distribution (also in x and y direction). This module is called by use of:

genpv NNN N x1 x2 x3x4 xS x6 pwx pwy ttt

NNN: total number of partieles (whole system);

N:number of centerpoints of particles, present in camera field;

x1: x-positions of partieles, at time of first exposure (in whole system); x2: y-positions of partieles, at time of first exposure (in whole system);

x3: veloeities of partieles, x-direction (in whole system); x4: veloeities of partieles, y-direction (in whole system); x5: veloeities of partieles, x-direction (in camera field);

x6: veloeities of partieles, y-direction (in camera field); pwx: position of upper left corner of window, x-coordinate; pwy: position of upper left corner of window, y-coordinate; ttt: time between two exposures.

The variables: x1, x2, x3,x4, x5 and x6 are 1 dim.arrays of float; pwx, pwy and Ut are float variables and N and NNN are integers.

The positions of partieles in whole system are taken with respect to the upper left corner of the system. The number of centerpoints present in window (N) means that these are present both at t=t1 and t=t2=t1+At.

A1.2. Mean and sigma values

The mean and sigma values of the reference are determined by self made module

meansigm1:

meansigm1 N xS mnfX sigmfX

meansigm1 Nx6

mnry

sigmry

N: number of centerpoints of partieles,present in camera field;

x5: velocities of partieles,x-direction; x6: veloeities of partieles, y-direction;

mnrx: mean value, x-direction;

mnry: mean value, y-direction;

sigmrx: sigma value, x-direction;

sigmry: sigma value, y-direction. Ag

(47)

A1.3. Generating particles: disks with radius: ra (dots)

The generated particles are disk shaped dots with radius RAD or ra.The module does not accept parameter RAD so that first is used:

copyRAD ra

The bitplanes belonging to t=t1 and t=t2=t1+~t and are cleaned by use of:

belr bt1; belr bt2;

bt1: bitplane belonging to t=t1 bt2: bitplane belonging to t=t2

The cleaned bitplanes consist of only zeros.

The objects generated and placed on the bitplanes: bt1 and bt2 occur by use of

genpartde:

genpartdc NNN x1 x2 x3 x4 pwx pwy tttra bt1 bt2

NNN: total number of particles (whole system);

N:number of centerpoints of particles, present in window;

xt: x-positions of particles, at time of first exposure (whole system); x2: y-positions of partieles, at time of first exposure (whole system);

x3: velocities of particles, x-direction (whole system);

x4: velocities of particles, y-direction (whole system);

pwx: position of upper left corner of window, x-coordinat;

pwy: position of upper left corner of window, y-coordinat;

ut:

time between two exposures;

ra: radius of the disk shaped objects;

bt1: bitplane consisting of objects at time t=t1; bt2: bitplane consisting of objects at time t=t2.

The x and y positions of particles (x1 and x2) are centerpoints generated by the program genpv (symmetrical velocity distribution) or genpvaa (asymmetrical velocity distribution). Also direct from genpv or genpvaa are the parameters x3 and x4. The parameters pwx and pwy are equal as used in moduls: genpv or genpvaa.

A1.4. Remove edge connected objects

Sometimes if the number of the objects is low and always if it is high, it occurs that the objects are connected to the window. This phenomenon gives extra errors and have to be removed before:

temeclcon bt1; temedcon bt2

bt1: bitplane consisting objects at time t=t1 bt2: bitplane consisting objects at time t=t2 A10

I

,

I

(48)

I

To check out the particles in camera field (t=t1 and t=t2) is used: displv=1/r=bt1/g=bt2

disp: display

v=1: viewport number 1

r=bt1: contents of bitplane bt1 are displayed in [ed color g=bt2: contents of bitplane bt2 are displayed in green color That gives in viewport 1 (v=1):

objects at t=t1 (red color) objects at t=t2 (green color)

If the objects belonging to t=t1 overlap with objects of t=t2, the overlap part is to see in yellow color.

A1.5. Detennining positions of particles trom images: t=t1 andt=t2.

At first is determined the number of objects and its numbers present in bitplanes: t=t1 and t=t2. These numbers are determined with use of:

blabel bt1 Iab1im 8 nobj1 blabel bt2 Iab2im 8 nobj2

bt1: bitplane consisting of objects (t=t1);

bt2: bitplane consisting of objects (t=t2);

lab1im: grey value image consisting of objects with seq. number (t=t1); lab2im: grey value image consisting of objects with seq. number (t=t2);

8: 8-connectivety;

nobj1: number of objects at t=t1; nobj2: number of objects at t=t2.

Next the positions are determined with use of shape:

shape Iab1im Iab1 xt1 yt1 shape Iab2im Iab2 xt2 yt2

lab1im: see above;

lab2im: see above;

lab1: label numbers of objects (t=t1);

lab2: label numbers of objects (t=t2);

xt1: positions of objects at t=t1, x-coordinates;

yt1: positions of objects at t=t1; y-coordinates;

xt2: positions of objects at t=t2, x-coordinates; yt2: positions of objects at t=t2; y-coordinates.