Abderrahman Ouayed
UNIVERSITEIT VAN AMSTERDAM Bachelor Informatica 10 augustus 2008
MATHEMATISCHE MORFOLOGIE: SNELLE ALGORITMEN VOOR EROSIES EN DILATIES
Begeleiders: Rein v. d. Boomgaard Dick van Albada
1
UNIVERSITEIT VAN AMSTERDAM BACHELOR PROJECT 2007-2008 NAAM: ABDERRAHMAN OUAYED
BEGELEIDERS: REIN V. D. BOOMGAARD DICK VAN ALBADA
PROJECT TITEL:
MATHEMATISCHE MORFOLOGIE: SNELLE ALGORITMEN VOOR EROSIES EN DILATIES
2
SNELLE ALGORITMEN VOOR EROSIES EN DILATIES
INHOUD Doel van het project ...................................................................................................... 5 1.1Geschiedenis en toepassingsgebieden:............................................................ 6 A-Binaire Morfologie ............................................................................................... 7 B-Elementaire Morfologische transformaties....................................................... 7 C-Grijswaarden morfologie ................................................................................. 10 2.1- MIN-MAX filters................................................................................................ 12 3.1- Naïef algoritme .................................................................................................. 14 3.2- Maxline (Pitas/Mike Brookes) ........................................................................... 14 3.3- HGW algoritme .................................................................................................. 16 3.4- Gil en Kimmel algoritme ................................................................................... 18 3.5- Lemire algoritme ................................................................................................ 18
IV- Optimaliseren Maxline, HGW en Lemire algoritmen.......................................... 21 4.1- Platform en algemene zaken voor optimalisatie ........................................ 21 4.1.1- Platform: ...................................................................................................... 21 4.1.2- Compiler flags: ............................................................................................ 21 4.1.3- Datastuctures: Array.................................................................................. 21 4.1.4- Datastructuur: double-ended queue (deque) ...................................... 22 4.1.5- Macros: inline functies ................................................................................ 23 4.1.6- If -Else vermijden ......................................................................................... 23 4.2- Optimalisatie HGW ............................................................................................ 24 4.3- Optimalisatie Lemire ......................................................................................... 25 V-TEsten en meetresultaten ........................................................................................ 27 5.1- Meetmethode: ............................................................................................... 27 5.2- Testprotocol: ................................................................................................... 27 5.3- Test resultaten: vergelijking van verschillende optimalisatie technieken .. 27 5.3.1- Geoptimaliseerde HGW versies ................................................................ 28 3
5.3.2- Geoptimaliseerde Lemire implementaties.............................................. 29 5.4- Algoritmen prestatietesten voor verschillende soort data .......................... 29 Één dimensioneel representatie van een 2D beeld ........................................ 32 Conclusie ................................................................................................................... 34 VII- Verder Werk ............................................................................................................ 36 VIII- Appendix ................................................................................................................ 37 IX- Reverences .............................................................................................................. 41
4
DOEL VAN HET PROJECT
Erosie en dilatie zijn de basisoperatoren van de mathematische morfologie in beeldverwerking [8][14]. In de loop der tijd zijn veel algoritmen bedacht om prestaties van deze operaties te verbeteren en ze zijn meestal lastig te begrijpen en niet eenvoudig om te implementeren. Bovendien zijn ze niet even geschikt voor elk type beeld als het gaat om efficiëntie. In 2006 kwam D. Lemire [12] met een zeer simpel algoritme voor Max/Min filters dat alle voorafgaande sequentiële algoritmen qua snelheid overtreft. Het presteert naar behoren voor alle soort beelden of tijdreeksen. Naar aanleiding zijn “Streaming” Min/Max filter algoritme zal ik in dit project een paar algoritmen met elkaar vergelijken en nagaan wanneer welk algoritme gebruikt kan worden. Ik zal mij beperken tot het behandelen van de Maxline (Pitas)[1][4][13], Van Herk [5][6][17] en Lemire algoritmen[12]. Deze worden hier verderop geoptimaliseerd en getest op verschillende beelden en data. Aan het einde van dit project wordt het Lemire “streaming” algoritme toegepast op twee-dimensionale grijsbeelden. In dit document zal ik een korte introductie van Mathematische Morfologie van zowel binaire als van (gray-scale) grijswaardes beelden geven.
SLEUTELWOORDEN: Mathematische morfologie, Binaire en (Gray-scale) grijsbeelden, Max/Min filter, erosie, dilatie, opening, (closing) sluiting, algoritme, C++.
5
I-
MATHEMATISCHE MORFOLOGIE IN BEELDVERWERKING
1.1 GESCHIEDENIS EN TOEPASSINGSGEBIEDEN:
Morfologie stamt van het Griekse woord morph en logos en betekent ‘de studie van de vorm’. Binnen het vakgebied van de beeldverwerking wordt mathematische morfologie (MM) gebruikt om de geometrische structuur van het beeld te analyseren en te bestuderen en zo op lokaal niveau te veranderen. MM is in de jaren 60 uitgevonden door Jean Serra en Georges Matheron. Zij hielden zich bezig met het bestuderen van automatische systemen die beelden uit de mineralogie en petrografie analyseerden. Pas na het verschijnen van artikelen van Haralick/Sternberg/Zhuang, PAMI 1987[18] werd MM wereldwijd veel toegepast.
1.2 Toepassingen van MM MM behandelt de geometrische structuur van een beeld op pixel niveau door middel van een zogenaamde structurerend element (SE) [14]. De keuze van een vorm en een maat van het SE bepaalt hoe de omgeving rondom ieder pixel wordt behandeld. De MM heeft de weg gevonden naar talloze gebieden die zich bezig houden met beeldverwerking en wordt op allerlei terreinen gebruikt zoals: • Medisch gebied: ondanks steeds betere middelen/technieken zoals MRI en CT-scanners is er plaats voor geautomatiseerde visuele inspectie van beelden. • Militair gebied: wapensystemen worden vaak uitgerust met computers die zorgen voor een intelligent waarnemingniveau. • Industrie: inspectie van producten met schade of verdachte vormen in productieketens. • Astronomie: het optimaal filteren van ruis die wordt veroorzaakt door kosmische stralen uit zonbeelden verkregen met satellieten [20].
1.3 Inleiding MM MM is gebaseerd op verzamelingsleer en maakt gebruik van elementaire operaties zoals vereniging, intersectie, complement. Verder wordt er ook gebruik gemaakt van een aantal geometrische transformaties zoals translatie, homothetie en rotatie[14][8]. Oorspronkelijk werd MM toegepast op binaire beelden en dit werd Binaire Morfologie genoemd. Een binair beeld wordt gezien als een afbeelding van een verzameling coördinaten (pixels) naar een verzameling beeldwaardes. Binaire beelden bevatten alleen de waardes 0 of 1 (zwart/voorgrond of wit/achtergrond).
6
A- BINAIRE MORFOLOGIE Een kort overzicht van begrippen uit verzameling theorie die worden gebruikt in binaire morfologie. • Twee verzamelingen zijn gelijk aan elkaar: X = Y ⇔ ( x ∈ X ⇒ x ∈ Y en x ∈ Y ⇒ x ∈ X ).
•
Inclusie: als X een deelverzameling van Y is of X omvat wordt door Y:
X ⊆ Y ⇔ ( x ∈ X ⇒ x ∈ Y ). •
Intersectie: de intersectie tussen twee verzamelingen wordt aangeduid als
X ∩ Y = {x | x ∈ X en x ∈ Y }. •
Vereniging (Union) X ∪ Y = {x |x ∈ X of x ∈ Y }.
•
Het verschil: het verschil van de X op Y, genoteerd X −Y is de verzameling van elementen van X die niet tot Y behoren X − Y = {x |x ∈ X en x ∉ Y }.
•
Complement: als X een verzameling is in een domein E, dan is het complement van X in E een deelverzameling van E genoteerd , zodat = {x |x ∈ E en x ∉ X}
• •
Reflectie (Symmetrie): de reflectie ˇX van de verzameling X is ˇX = {−x |x ∈ X}. Translatie X De verzameling na translatie van X met een vector b is gelijk aan de verzameling X = {z ∈ E | z = x + b, x ∈ X}.
B- ELEMENTAIRE M ORFOLOGISCHE TRANSFORMATIES
b.1 Dilatie: Dilatie van een beeld X met een structurerend element A wordt gedefinieerd als een Minkowski som: δ# (X) = X ⊕ A = {x + a | x ∈ X en a ∈ A} De dilatie is te interpreteren als de vereniging van getranslateerde A in X i.e. verzameling punten waar A, X aanraakt.
X⊕A A
X
Figuur 1:De dilatie van een vierkant X door een circulair SE A
7
Dilatie heeft de volgende eigenschappen: • Commutativiteit. Dit betekent dat het dilateren van X met A is het zelfde als A dilateren met X. • Associativiteit: met deze eigenschap kan een structurerend element gebouwd worden door een samenstelling van een reeks dilaties. Het is handig bijvoorbeeld om een groter structurerend element te splitsen in kleinere vormen. Zo kan een 5x5 SE gebouwd worden door een dilatie van twee 3x3 SE.
⊕
=
⊕
=
Figuur 2: Decompositie grotere SE naar kleinere structurerende elementen door dilaties •
Distributiviteit over ∪: is zeer handig bijvoorbeeld bij constructie van een structurerend element in de vorm van een schijf door een vereniging te nemen van een aantal rechthoeken
U
U
=
Figuur 3: verenigen rechthoeken om een cirkel te vormen
b.2 Erosie: Erosie van een beeld X door een structurerend element A is gedefinieerd als een Minkowski substractie: ε# (X) = X ⊝ A = {x|A ( ⊆ X} = ∩)∈# X *) en het is te interpreteren als verzameling elementen van X als de origine van A langs alle punten uit X doorloopt en blijft geheel in X passen. (waar past A in X ?)
A
X⊝A
Figuur 4: Erosie van vierkant X door een S.e. A De associativiteit van de dilatie kan ook gebruikt worden bij de erosie om een groter structurerend element uit te splitsen in kleinere delen . Dit is bekend als de “chain rul” : Als A = A1⊕ A2⊕ A3 dan is X ⊝ A=X ⊝ (A1⊕ A2⊕ A3) = ((X ⊝A1) ⊝A2) ⊝A3
8
b.3 Opening: De “opening” van een beeld wordt verkregen door het dilateren van een geërodeerd beeld met het zelfde structurerende element
X ∘A = ( 0 ⊝ A ) ⊕ S Openen lijkt enigszins op erosie omdat het kleine heldere gebieden kan verwijderen, maar in het algemeen is de opening minder drastisch dan erosie. Het effect van de operator moet voorgrondgebieden bewaren die een gelijksoortige vorm heeft als het structurerend element, of die het structureerde element volledig kan bevatten, terwijl alle andere gebieden van voorgrondpixels elimineert. Bijvoorbeeld het openen met een SE in de vorm van een schijf maakt contouren glad, breekt smalle landengtes, elimineert kleine eilanden en scherpe pieken.
A
X ∘A
Figuur 5: Openen van een vierkant X met een circulaire S.e. A Het resultaat is een vierkant met gladde hoeken
Figuur 6: Filteren van vormen die lijken op het gekozen SE. door te openen.
b.4 Closing (Sluiting): De sluiting van een beeld wordt verkregen door het eroderen van een gedilateerd beeld met hetzelfde structurerende element
X • A = ( X ⊕ A ) ⊝S Sluiten lijkt op de dilatie in zoverre dat het neigt om de grenzen van (heldere) voorgrondgebieden in een beeld te vergroten. De sluiting bewaart de achtergrondgebieden die kunnen het structureerde element volledig bevatten, terwijl het alle andere gebieden van de achtergrond elimineert. Het sluiten door een SE in de vorm van een schijf maakt contouren glad, smelt smalle onderbrekingen, lange en dunne golven, elimineert prikken, vult hiaten binnen de contouren.
A X
X•A
Figuur 7: sluiten van twee aan elkaar geplakte vierkantjes.
9
C-
GRIJSWAARDEN MORFOLOGIE
Grijswaarden (Gray-scale) morfologie is een uitbreiding van de binaire morfologie. De beeldwaarden zijn niet alleen 1 en 0 maar kunnen alle grijswaarden hebben tussen zwart en wit (i.h.a waarden van 0 tot 255). Zo worden de binaire beelden een speciaal geval van grijsbeelden. De operaties in binaire morfologie kunnen in grijswaarde morfologie als volgt vervangen worden: •
Inclusie ⊆ vervangen door kleiner of gelijk ≤
• • •
Vereniging ∪ door maximum ⋁ Intersectie ⋂ door minimum ⋀ De translatie van F met a is dan F) (x)=F(x-a)
C.1 Dilatie Uitgaande van de definitie van de dilatie van binaire beelden wordt de functionele dilatie voor grijsbeelden F⊕ A = {F) | a A} of (F⊕ A) (x) = Max{ F(x-a) | a A }. In woorden uitgedrukt is de dilatie van de functie F door een structurerend element A de grootste waarde van F binnen A met A gepositioneerd in punt x. 2D grijswaarden beeld 1
3
4
5
6
5
3
2
1
1
3
6
4
8
9
2
4
7
6
7
7
7
8
9
4
5
6
6
7
7
3x3 SE
1 3 4 5 6 5 9 3 4 5 6 5 9 9 4 5 6 5 9 9 9 5 6 5 9 9 9 7 3 2 1 1 3 6 3 2 1 1 3 6 3 2 1 1 3 6 3 2 1 1 3 6 9 9 9 9 4 8 9 2 4 7 4 8 9 2 4 7 4 8 9 2 4 7 4 8 9 2 4 7 9 9 9 9 6 7 7 7 8 9 6 7 7 7 8 9 6 7 7 7 8 9 6 7 7 7 8 9 6 7 7 7 4 5 6 6 7 7 4 5 6 6 7 7 4 5 6 6 7 7 4 5 6 6 7 7 4 5 6 6 Figuur 8: Dilatie van 2D beeld met 3x3 SE. Het maximum wordt bepaald binnen ieder 3x3 raster en het resultaat wordt geplaatst in de oorsprong van het SE. De twee rijen en colommen (randpixels) zijn intakt gebleven. De randpixels eisen een speciale behandeling, zie verder in hoofdstuk VI.
6 3 4 8 7
5 6 7 9 7
10
C.2 Erosie Analoog aan de bovenste definitie is de erosie van F met A gegeven door het minimum te nemen van waardes van F binnen Amet A is gepositioneerd in x:
F ⊝ A = ⋀ {F) | a ∈A} of (F ⊝ A) (x) = min{ F(x+a) | a ∈ A }. 2D grijswaarden beeld 1
3
4
5
6
5
3
2
1
1
3
6
4
8
9
2
4
7
6
7
7
7
8
9
4
5
6
6
7
7
3x3 SE
1 3 4 5 6 5 1 3 4 5 6 5 1 1 4 5 6 5 1 1 1 5 6 5 1 1 1 1 6 3 2 1 1 3 6 3 2 1 1 3 6 3 2 1 1 3 6 3 2 1 1 3 6 1 1 1 2 3 4 8 9 2 4 7 4 8 9 2 4 7 4 8 9 2 4 7 4 8 9 2 4 7 4 2 2 2 4 6 7 7 7 8 9 6 7 7 7 8 9 6 7 7 7 8 9 6 7 7 7 8 9 6 7 7 7 8 4 5 6 6 7 7 4 5 6 6 7 7 4 5 6 6 7 7 4 5 6 6 7 7 4 5 6 6 7 Figuur 9: Erosie van 2D beeld met 3x3 S.e. Het minimum wordt bepaald binnen ieder 3x3 raster en het resultaat wordt geplaatst in de origine van SE. Voor randpixels zie toelichting in hoofdstuk VI.
5 6 7 9 7
Analoog aan de bovenste definities zijn de opening en de sluiting opnieuw te definiëren als
C.3 Opening ( F ∘ A ) (x) = ((F ⊝ A ) ⊕ A)(x) 1 1 1 1 1 1 1 2 4 2 2 2 6 7 7 7 4 5 6 6 Figuur 10: randpixels
6 5 4 1 1 1 6 5 4 2 1 1 6 5 4 2 3 6 1 1 1 2 3 6 1 1 1 1 3 6 1 1 4 7 4 2 2 2 4 7 4 2 2 2 4 7 4 2 8 9 6 7 7 7 8 9 6 7 7 7 8 9 6 7 7 7 4 5 6 6 7 7 4 5 6 6 7 7 4 5 De opening: De dilatie van het gerodeerde 2Dbeeld zie toelichting in hoofdstuk VI.
2 1 6 5 1 1 3 6 2 2 4 7 7 7 8 9 6 6 7 7 met een 3x3
4 2 2 5 6 6 5 6 6 6 7 7 4 5 6 SE.Voor
4 7 7 7 6
6 3 4 8 7
5 6 7 9 7
5 9 7 7 6 7 7 7 7 6 6 7 9 6 7 7 7 4 5 6 3x3 SE.Voor
7 9 7 7 6
6 3 4 8 7
5 6 7 9 7
C.4 Sluiting (F • A )(x)= (( F ⊕ A ) ⊝ A)(x) 9 9 9 7 9 9 9 9 9 9 9 9 6 7 7 7 4 5 6 6 Figuur 11: randpixels
6 5 9 9 9 7 6 5 9 7 9 7 6 5 9 7 3 6 9 9 9 9 3 6 9 9 9 9 3 6 9 9 4 7 9 9 9 9 4 7 9 9 9 9 4 7 9 9 8 9 6 7 7 7 8 9 6 7 7 7 8 9 6 7 7 7 4 5 6 6 7 7 4 5 6 6 7 7 4 5 De sluiting: de erosie van het gedilateerde 2Dbeeld zie toelichting in hoofdstuk VI.
7 7 6 9 9 3 9 9 4 7 7 8 6 6 7 met een
11
I-
MM ALGORITMEN
De op morfologieoperaties gebaseerde systemen van de beeldverwerking worden typisch gebruikt voor toezichttaken in industriële systemen[5][1][4][9], medische beeldverwerking, het optisch aflezen van grafische tekens, textuuranalyse, enz. Vele algoritmen voor diverse morfologische operatoren zijn gepubliceerd met als doel het optimaliseren van de verwerkingtijd daarvan. Maar er is slechts een klein aantal algoritmen eenvoudig te vertalen naar code en dit is niet altijd toepasbaar op elke soort beeld. Het recente algoritme van D. Lemire 2007[12] daarentegen is wel eenvoudig en is efficiënt voor alle soort beelden. In het artikel “Streaming maximum-minimum filter using no more than three comparisons per element” heeft hij bewezen dat het mogelijk is om de basisoperatoren erosie en dilatie te berekenen in niet meer dan 3 vergelijkingen per element. Verderop in dit project worden de naïeve, Pitas, Herk en Lemire algoritmen voor Max-Filter geïmplementeerd en tot slot worden de laatste twee algoritmen geoptimaliseerd. Deze algoritmen zijn in principe allemaal eendimensionale algoritmen, maar kunnen vanwege de associativiteit van de operaties worden gebruikt als bouwstenen voor ingewikkelder structurererende elementen.
2.1- MIN-MAX FILTERS Een filtervenster is een structurerend element dat wordt gebruikt om MM operaties uit te voeren. Als een filterwaarde is bepaald, wordt deze in het nieuwe beeld op de plaats gezet van de oorsprong van het structurerend element. Erosie en de dilatie in grijsbeelden morfologie impliceren het gebruik van Min- en Max filters. De eendimensionale Max/Min filters wordt als volgt gedefinieerd: Als D = {<= ,i=0,…,n} een data sample is en X een filtervenster van lengte k dan is de verzameling {>= } de Min/Max filters van D met het venster X als >= de maximum/minimum is van { <= ,…,<=?@*A } voor ieder i in {0…(n – k)}. Max/Min filter worden veel gebruikt in niet lineair signaal / beeldverwerking en vooral in “realtime” taken, vandaar de noodzaak van een snel algoritme dat de MM-operaties sneller kan uitvoeren. De reken complexiteit van Max/Min filters of het aantal vergelijkingen per data element hangt van de lengte af van het filtervenster. Een naïeve berekening van deze filter voor een data sample van lengte n met een filtervenster van lengte k vereist ( k - 1 ) . ( n - k ) = n . k – n + k vergelijkingen en het is dus van de orde O( n k ). Voor kleine structurerende elementen is het niet zo erg maar in real-time toepassingsgebieden van MM worden vele grotere reeksen beelden geanalyseerd met grotere SE’s en dat is alleen mogelijk als er een goed en efficiënt algoritme is.
2.2- De Maxline van Pitas verbeterd door Brookes Pitas(1989) [1],[2],[4] en[9]heeft met zijn maxline algoritme voor Max/Min filter de complexiteit log2(p) bereikt maar niet voor elk filtervenster. De maxline methode bereikt de beste performance, als de SE lengte van macht 2 is (BC ) en ietsje meer dan 3 vergelijkingen per element als de data sample willekeurig en onafhankelijk is gedistribueerd. Douglas(1996)[4] heeft een betere performance gehaald door de maxline te verbeteren en heeft met zijn maxlist een beter uitgebalanceerd algoritme voor verschillende soort data ontwikkeld. In maxline2 van Mike Brookes (2000)[1], werd de “worst-case performance” verbeterd door als eerste een cirkelbuffer te gebruiken om tijdelijke data op te slaan, daardoor is de databeweging verminderd en heeft als resultaat een betere performance. Ten tweede door de index van het huidige maximum op te slaan wordt in maxline2 niet altijd de hele buffer doorgezocht naar het nieuwe maximum. Dit heeft als resultaat dat MAXLINE2 efficiënter in zijn search operaties. In het volgende hoofdstuk is er een pseudocode van Maxline2.
12
2.3- Recursief van Herk algoritme (HGW) Het van Herk (HGW)algoritme(1992)[17] is het eerste algoritme dat heeft Max/Min filter kunnen
E
berekenen met een complexiteit die onafhankelijk is van de lengte k van het filtervenster (3- ).
@
Het is een driestaps algoritme dat gebruik maakt van drie buffers. Zie verder voor uitleg en implementatiecode. Het HGW algoritme is door Gevorkian(1996)[19] verbeterd door een aantal onnodige vergelijkingen te vermijden en heeft daardoor een gemiddelde complexiteit van 2.5-
G.H @
+ 1/J C
voor willekeurige gedistribueerde data. De beste verbetering van HGW tot nu toe is door Gil en Kimmel (2002) [6] uitgebracht maar in dit geval voor welke soort data dan ook en niet alleen voor i.i.d data. De aanpak lijkt op wat Gevorkian heeft gedaan maar wel uitgebreider. Gil en Kimmel hebben de onnodige vergelijkingen totaal uit de weg geruimd in elk van de drie stappen. De complexiteit is daardoor verlaagd naar een gemiddelde van 1.5+KLM(J)/J+N(1/J).
2.4- Streaming Max/Min filter : monotone wedge De “streaming “ methode van Lemire (2006)[12] is het beste algoritme tot nu toe die de Max en Min filter bij elkaar kan uitvoeren in niet meer dan 3 vergelijkingen per element en in de “worstcase” scenario wordt de grens van 3 vergelijkingen per element niet overschreden. In het algoritme is een “double ended queue” gebruikt om de indexen van de globale Max/Min per venster te bewaren en niet de Min/Max waardes zelf. Deze queue wordt monotone “wedge” genoemd want hij bevat stijgende indexwaardes (pointers naar monotoon dalende/stijgende Max/Min waardes). De globale Min/Max zit altijd aan het begin van de queue gevolgd door volgende Min/Max van de search stappen. Tabel 1: Overzicht van 1D-filter rekencomplexiteit Max en Min filter bij elkaar Algoritme Naïef Maxline HGW GAA GK Monotone wedge
Aantal operaties per element O(2k) O (2log k) (6 – 8 / k) (5-(7 /k)+2/J C ) 3 + 2 log k/k+ O( 1/k)
3
[Pitas, 1989] [van Herk, 1992, Gil and Werman, 1993] [Gevorkian ,1997] [Gil and Kimmel, 2002] [Lemire, 2006b]
13
II-
ALGORITMEN PSEUDOCODE EN IMPLEMENTATIE
3.1- NAÏEF ALGORITME
Om de algoritmen te kunnen testen op de juistheid heb ik ook het naïeve algoritme ook geïmplementeerd Dit geeft ook meteen de basiswaarde voor snelheidmetingen
HET ALGORITME: Stap 1: Plaats het venster aan het begin van het beeld Stap 2: Bepaal de Max in dat venster Stap 3: Sla de Max op in maxResults Stap 4: Zet het venster een stap verder en herhaal vanaf Stap 2 totdat de eindindex van het venster overeenkomt met de eind index van het beeld.
Figuur 12 : naïeve Max 1D filter (Min filter is analoog) functie Het algoritme maakt per venster (width-1) vergelijkingen en (imSize-width+1) keer. Per element is het aantal vergelijkingen van de orde O(width). width: de breedte van het SE imSize: lengte van het 1Dbeeld maxValues: gedilateerde beeld minValues:geërodeerde beeld
14
3.2- Maxline (Pitas/Mike Brookes) HET ALGORITME MAXLINE Image: maxResults: circulairBuffer: SeSize: p:
input data (het beeld). output beeld van het filter. circulaire buffer ter lengte van het SE. lengte van het SE. index huidige max
maxResults: max-filter resultaat en heeft lengte maxResultsSize. maxResultsSize = image.size()-SeSize+1 maxResults heeft als start index -1 en de waarde op dat index is -∞. Stap 0: set
n=0, p=0, maxResults[-1]=-∞ circulairBuffer[j] =-∞ for j = 0 up to (SeSize-1)
Stap 1: set j= n modulo ( SeltSize), circulairBuffer[j]= image[n] Stap 2: if(image[n]>=maxResults[n-1]) then p=j goto Stap 4 Else if (p == j) goto Stap 3 Else goto Stap 4 Stap 3: Foreach m = j-1 downto 0 { if (circulairBuffer[m]>circulairBuffer[p]) set p =m } Foreach m =SeltSize-1 downto j+1 { if (circulairBuffer[m]>circulairBuffer[p]) set p =m } Stap 4: maxResults[n]=circulairBuffer[p],n++ goto Stap 1
In het geval dat de inputdata onafhankelijk en identiek is verdeed (i.i.d) ,is het verwachte aantal gegevensvergelijkingen per inputwaarde onafhankelijk van filter SeltSize. Voor groot SeSize echter, wordt de rekentijd van (Pitas) MAXLINE[1][2][9] algoritme overheerst door gegevensbewegingen binnen de opslagbuffer en is asymptotisch evenredig aan SeSize. Brookes heeft deze gegevensbewegingen kunnen elimineren door een circulairBuffer te gebruiken om de inputwaarden op te slaan. Dit heeft geresulteerd in een beter prestatie van het algoritme. Het herziende MAXLINE algoritme door Brookes slaat de input sample in een circulairBuffer van dezelfde lengte als de structurerend element. In stap 1:Een nieuwe input wordt in de circularBuffer opgeslagen, de index j neemt de waardes (n % SeSize) om te zorgen dat j periodieke waardes neemt tussen 0 en SeSize-1. Stap 2 en 3 zorgen er voor dat de index p altijd wijst naar het maximum van de gegevens in de circularBuffer. De stap 3: wordt uitgevoerd om de hele buffer door te zoeken naar globale Max. De zoektocht wordt uitgevoerd in de omgekeerde chronologisch volgorde beginnend bij de data sample circularBuffer[j]. Alleen in het geval dat maxResults(n-1) correspondeert met de oudste data sample in de buffer image(n-SeSize) wordt stap 3 uitgevoerd. Door p te laten verwijzen naar het meeste recente maximum in stap 2 wordt het aantal keer dat stap 3 wordt uigevoerd geminimaliseerd. Stap 4: sla het huidige maximum circularBuffer[p] op in maxResults en waardeer n op voor de volgende ronde.
15
Figuur 13: De Maxline code
3.3- HGW ALGORITME In het HGW [17] algoritme wordt de data input gesplist in overlappende segmenten van lengte 2p-1 gecentraliseerd op
Recursief geformuleerd: Stap 1: foreach k = 1…,p-1 then set R[j,k]=max(R[j,k-1],x[j-k]) Stap 2: foreach k=1,…,p-1 then set S[j,k]=max(S[j,k-1],x[j+k])) Stap3: if ( k == 0 ) then set M[j,0]=R[j,p-1] Foreach k=1,…,p-2 set M[j,k]=max(R[j,k],S[j,p-k-1]) If ( k==p-1)then set M[j,p-1]=S[j,p-1]
16
HGW IMPLEMENTATIE: Het resultaat M[k,j] met k =0…p-1 bevat het maximum binnen ieder venster van lengte p bewegend van 0 tot p. De Max in stap-1 is (p-1) keer uigevoerd (met k=1 …p-1 en hetzelfde geldt voor stap-2, zo zijn R en S bij elkaar 2* ( p-1 ) Max operaties uitgevoerd. De Max in stap-3 is (p-2) keer uitgevoerd. Want R[j, p - 1] en S[j, p -1] zijn al bekend. Zie recursieve formule van het algoritme. In totaal zijn er dus 2*( p-1 ) + ( p-2 ) Max operaties gedaan, gemiddeld is het dus 3-4/p vergelijkingen per element. Als p groter wordt dan is het aantal vergelijkingen gelijk aan 3. Het input beeld B B1
B2
B3
B4
B5
B6
De linker Buffer R (stap1) R4=max R3,B1
R3=max R2,B2
R2=max R1,B3
B7
B8
B9
De rechter Buffer S(stap2) R1=max R0,B4
R0= S0 =B5
S1=max S0,B6
S2=max S1,B7
S3=max S2,B8
S4=max S3,B9
Het merge buffer M met het Maxfilter resultaat (stap3) M0=max S0,R4=R4
M1=max S1,R3
M2=max S2,R2
M3=max S3,R1
M4=max S4,R0=S4
….
….
….
….
Figuur 14: Voorbeeld: Maxfilter van een beeld B met een SE ter lengte 5
Figuur 15: van Herk implementatie zonder optimalisatie.
17
3.4- GIL EN KIMMEL ALGORITME
Het algoritme van Gil en Kimmel [5] is een verbetering van de HGW door een aantal onnodig vergelijkingen te vermijden in alle drie fases. Gevorkian[19] merkte al op dat het suffix S van het ene segment overlapt met het ene prefix R van het andere en dit heeft als gevolg dat een aantal prefixen en suffixen niet bepaald hoeft te worden. Wel wordt nu ieder segment door twee gedeeld en ieder deelsegment heeft lengte q=p/2 + (p modulo 2)/2. Stap 1 foreach k=0,…,q-1 S[j,k]=max(x[j],…,x[j+k]) Stap 2 foreach k=q,…,p R[j,k]=max(x[j+k],…,x[j+p]) Stap 3 If(S[j,q-1]<=R[j,q]) then R[j,q-1]=R[j,q-2]=…=R[j,1]=R[j,q] goto Stap:4 Else then S[j,q]=S[j,q+1]=…=S[j,p-1]=S[j,q-1] goto Stap:5 Stap 4 foreach k=q…p-1 S[j,k]=max(x[j],…,x[j+k]) Stap 5 foreach k=1…q-1 R[j,k]=max(x[j+k],…,x[j+p])
Het is duidelijk dat er een heel aantal vergelijkingen overbodig is geworden door dit algoritme. In stap 1 en 2 bij elkaar zijn (p-1) vergelijkingen gemaakt, in stap 3 (p-q) en in stap 5 (q-1), maar stappen 3 en 5 worden niet allebei uitgevoerd, slechts de ene of de andere. We gaan dus uit van het maximum van {q - 1,p - q }; er worden of ( q – 1 ) of ( p – q ) Max operaties gedaan. In totaal worden dus (p-1)+Max( q - 1, p - q ) vergelijkingen uitgevoerd. Uiteindelijk, per element wordt ( (p-1)+max(q-1,p-q))/p vergelijkingen gemaakt en het is gelijk aan (1.5 +c/p) waar c een constante is. Verder hebben Gil en Kimmel ook aangetoond dat de merge procedure ook efficiënter gemaakt kan worden maar ik ga hier niet verder op in. Aangezien de complexiteit van het Gil en Kimmel algoritme en het klein verschil in de efficiëntie met het van Herk algoritme is het niet rendabel om het verder te implementeren
3.5- LEMIRE ALGORITME Er is geen enkel algoritme bekend die het Max/Min filter kan berekenen in minder dan 1.5 vergelijkingen per element maar Lemire heeft met zijn algoritme “Streaming Maximum-Minimum Filter” [12] deze limiet wel kunnen halen en voor alle soorten beelden datasignalen
HET ALGORITME: Het algoritme berekent tegelijker tijd het Max en Minfilter maar ze kunnen ook los van elkaar berekend worden. Max filter:
INPUT: an array A indexed from 1 to n INPUT: window width w >=2 U ← empty double-ended queues, we append to “back” append 0 to U foreach i =1 upto n do { if i ≥ w + 1 then { OUTPUT: front(U) as maximum of range [i − w, i) } (1) if A[i] > A[i−1] then { pop U from back (2) while A[i] > A[back(U)] do { pop U from back } } else { append i to U } if i = w + front(U) then { pop U from front } }
18
Er wordt een double-ended queue gebruikt om de index van een globaal maximum in het venster te traceren. Het eerste element van queue (front) wijst altijd naar de index van de grootste waarde binnen het venster; naar het globale maximum dus. Uiteindelijk wijzen de waardes in de queue naar (oplopend) gesorteerd vensterdeel vanaf het laatste globale maximum binnen het venster. Ieder keer wordt de waarde met de index in de back van queue vergeleken met volgende nieuwe waarde in de data array. Als het backelement van de queue wijst naar een groter element dan wordt de index van de kleinere waarde van het dataarray toegevoegd aan de queue (push_back). Anders worden alle waardes waarvan de indexen in de queue zijn opgeslagen vergeleken met de nieuwste datainput. Deze worden weggehaald als deze kleiner zijn dan de nieuwste waarde.
Voorbeeld: Filtervenster van lengte 7. Indexen : 0123456 Data input I : 5 2 0 1 4 7 6 Begin : U(0)=0 , verwijst naar het eerste element van I ,U.back = 0 2de elt 2 < I(U.back)= 5 index van 2 wordt toegevoegd aan U U = {0 ,1}, U.back = 1 3de elt 0 < I(U.back)= 2 index van 0 wordt toegevoegd aan U U={0, 1, 2}, U.back = 2 4de elt 1 > I(U.back)= 0 U.pop_back U = {0, 1} , U.back = 1 1< I(U.back)= 2 index van 1 wordt toegevoegd aan U U = {0,1,3} ,U.back = 3 5de elt 4> I(U.b ack)= 1 3 en 1 worden weg gehaald en U ={0, 4} want “4 >1,2” en “4 <5”, U.back=4. 6de elt 7> I(U.back)= 7 alle elementen in de queue zijn kleiner dan 7, ze worden weggehaald U={5}. 7de elt 6
U.front verwijst naar het grootste element van de data input binnen het venster. Aan de eind procedure van het venster verwijst de counter nu naar het laatste element in het venster. I(U.front) is nu maximum binnen het venster; de output. i == w + front(U) controleert of front van U naar het eerste element van het venster verwijst en dit is ook het globale maximum, als dit het geval is dan wordt deze verwijderd door pop_front om het venster verder te slepen naar de volgende stap. Wat de prestatie betreft is het niet echt duidelijk dat in de “worst-case” drie vergelijkingen per element worden gemaakt. • •
•
Als de data dalend is dan wordt ieder element maar 1 keer vergeleken en dat gebeurt in (1), Als de data stijgend is dan wordt ieder element ook 1 keer vergeleken in (1) nu verder naar (2) die hoeft niet uitgevoerd te worden want U.pop_back is al uitgevoerd en U is leeg. Het slechtste geval: als U.front verwijst naar het grootste element en data input verder de volgende vorm heeft: I[0] is het maximum I[2n] < I[4n] voor ieder n=1…w/2 (w vensterlengte) I[2n] > I[2n+1]
Voorbeeld:
Datainput Vergelijkingen
I= 9121314151617 012 1313131313
Gemiddeld per element is het dus 23/13=1.77 vergelijkingen voor de Max filter in dit geval, gecombineerd met Min filter is het in het ergste geval 3 vergelijkingen per element, en als de data monotoon dan worden er ongeveer 2 vergelijkingen per element gedaan.
19
Figuur 16: Eerste Lemire implementatie
20
IV- OPTIMALISEREN MAXLINE, HGW EN LEMIRE ALGORITMEN
4.1- PLATFORM EN ALGEMENE ZAKEN VOOR OPTIMALISATIE 4.1.1- PLATFORM:
De programma’s zijn in C++ geprogrammeerd en gecompileerd met GCC-G++ compiler versie 3.4.4 in CYGWIN omgeving en onder Windows Vista 32bit. Het werkstation heeft processor AMD Turion(tm)64 X2 Mobile technologie TL-58 1,90GHz met 3G RAM. Er is een aantal belangrijke hoofdpunten die van belang zijn om deze programma’s minder traag te krijgen. Kennis van het platform en de datastuctuur keuzes spelen hierin een belangrijke rol.
4.1.2- COMPILER FLAGS:
Uit een aantal testen is gebleken dat de optimalisatie vlag -O3 de beste resultaten oplevert voor alle programma’s. De reden is dat -O3 optimaliseert niet alleen voor de snelheid zoals -O2 maar ook voor intensieve loops en inline functies. Dat is maar goed ook. Want alle programma’s zitten vol loops en een aantal inline functies. Verder heeft GNU voor de g++ compiler voor AMD processoren specifieke compiler opties ontwikkeld om betere optimalisaties te verkrijgen voor een specifieke machine. In de gebruikte machine, Turion(tm)64, geven de vlaggen -march=athlon-mp of beter nog -mtune/march=athlon64 respectievelijk samen met –O3 de beste resultaten [15]. Testen Flags –O3 en –O2: meetresultaten zijn in msec voor 1D-image input met 1.000.000 “floating points“ en de Max-minfilter van lengte 4,5,6
lengte 4 5 6
$ g++ -O3 -mtune=athlon64 mmFilter.cpp -o mmFilter vanHerk1 vanHerkF vanHerk3 lemire 0.078 0.094 0.047 0.047 0.062 0.094 0.047 0.047 0.063 0.108 0.048 0.031
lemire_nopop 0.047 0.047 0.047
lengte 4 5 6
$ g++ -O2 -mtune=athlon64 mmFilter.cpp -o mmFilter vanHerk1 vanHerkF vanHerk3 lemire 0.265 0.172 0.156 0.093 0.234 0.202 0.141 0.093 0.25 0.203 0.14 0.094
lemire_nopop 0.093 0.094 0.093
4.1.3- DATASTUCTURES: ARRAY
Image Input/output buffers zijn tot nu toe 1D arrays. De snelste manier om data te bewaren is een standaard array. Ook al is deze ietsje lastiger te beheren dan een vector datastructuur in C++. Hoewel de vectors gebruiksvriendelijk zijn, zijn de arrays sneller en daar gaat mijn voorkeur voor want met behulp van Template functie in C++ is het mogelijk om een array zo te definiëren dat het eenvoudig wordt om te manipuleren. Daarmee wordt een array gedeclareerd voor elke datatype.
21
Array Template: een simpel manier om een array te manipuleren template
class ArrayVector{ T* F; int SIZE; public: ArrayVector(int n= 10) { SIZE = n ; F =new T[n]; } Inline void clear() {delete F;} Inline int size(void) {return SIZE;} Inine T& operator [] ( int index) { return F[index]; } ~ArrayVector(){delete T;} };
4.1.4- DATASTRUCTUUR: DOUBLE-ENDED QUEUE (DEQUE)
In het Lemire algoritme is er gebruik gemaakt van een STL container double-ended queue om indexen op te slaan. In deze databuffer vindt veel dataverkeer plaats. Het is dus van belang dat deze databuffer ook sneller gemanipuleerd kan worden. Mijn voorstel is een array versie van een double-ended queue. Het is een array implementatie van een circulairequeue cirQueue die gebruik maakt van “bitmasking” zodat de data tussen front en back de enige data die zichtbaar blijft (met dank aan Dick van Albada). Alleen functionaliteiten die nodig zijn worden geïmplementeerd: push_back(int), pop_back(), pop_front(), front(), back() en size() en resize(int). Push_back(i) plaats i achter het array en verhoog back_index met 1 Pop_front() / Pop_back() schuif front_index()/back_index een stap naar voren / achter Bij herhaaldelijk pop_back() is beter een resize() functie te gebruiken. Als het aantal pop_back operaties bekend is, zetten wij de index back_index meteen op zijn plaats door één instructie namelijk back_index =(front_index + aantal pop_back). Merk op dat de lengte van de cirQueue altijd een macht van twee is (zie de Appendix). Elementen kunnen gewoon als bij array benaderd worden. In deze simpele test is te zien dat de cirQueue wel ruim 45% efficiënter is dan de STL container deque voor het meeste operaties. $ ./queuTest 10000000 deque push_back 0.156 secs deque pop_back 0.078 secs deque pop_front 0.063 secs cirQueue push_back 0.077 secs cirQueue pop_back 0.047 secs cirQueue pop_front 0.047 secs Voorbeeld: Queu m(32); // declareer een Queu van lengte 32. m.push_back(5); m[0]=5=m.back()=m.front().
22
4.1.5- MACROS: INLINE FUNCTIES
Alle kleine functies zoals in het cirQueue doen niet zoveel dus deze kunnen als inline gedeclareerd worden. Deze zijn enigszins vergelijkbaar met macro definities in C. Op deze manier kan de compiler beter optimaliseren waar het nodig is [15]. Als de compiler heeft gekozen om een functie “in te lijnen” dan wordt er een kopie van gemaakt in plaats van de functie aan te roepen. Het stack frame en de “return” functie kunnen op deze manier worden vermeden. In de van-Herk implementatie wordt veel gebruik gemaakt van max en min deze zijn ook weer met inline en static gedeclareerd.
4.1.6- IF -ELSE VERMIJDEN
Sprongopdrachten, zoals door de compiler gebruikt om een if – else statement te vertalen zijn relatief duur. Soms is het mogelijk deze opdrachten te vermijden. Als voorbeeld kan de maxfunctie, normaal geschreven als (x>y)?x:y; ook worden geïmplementeerd zonder conditionele “branches” en dat kan gunstig zijn als er veel Max aangeroepen moeten worden gedaan zoals in het geval van HGW algoritme. Als het zo is dat is het beter om Max te implementeren als volgt: Inline static uint Max( const int a,const int b) { return (uint) (0.5 * (abs ((a-b)) + (a-b) ) + b);} Of in het geval van integers type, binaire shift te gebruiken: static inline int Maxb ( int a, int b) {
a -= (a-b) & ((a-b)>>31); return a; } als a-b>=0 dan is b&((a-b)>>31)=0 dus a-=0 (is a) anders is het gelijk aan b dan is a-=a-b ie a=a-(a-b)=b bit shifting x>>31 : als x een unsigned integer is (positief of null) dan alle “geshifte” bits naar rechts zijn gedumpt en de linkere bits zijn naar 0 gezet. Dus een 32 bit unsigned integer wordt 0 na 31 shift naar rechts vandaar y & (x>>31)= y AND 0000000…000 = 0 In het geval van een signed integer wordt (meestal) het teken-bit naar rechts ingeschoven. Bij een negatieve signed integer worden dus alle bits aan de linker kant 1 en zo wordt y&(x>>31) = y & 1111…1111=y . Dit gedrag is niet portable.
23
4.2- OPTIMALISATIE HGW Zie code in de volgende pagina. Regels (67 en 112) in Figuur 17 In het geval van vergelijkingen tussen integers is MINB de beste manier. Zonder if-else expressie en met bit masking lijkt het beter te werken. static inline int MinB( int a, int b) { b = a-b; a -= b & (~(b>>31)); return a; }
1) 2) 3)
4)
regel (57 en 103 Figuur 17 ) seSize<<1 = 2*seSize , een bit shift operatie is sneller dan vermenigvuldiging. Regel (60 en 105 Figuur 17) integer (j+seSize-1) wordt vaker gebruikt dus eenmaal berekenen en opslaan. De voorafgaande verbeteringen hebben wat minuscule tijdwinst opgeleverd maar het beste resultaat is natuurlijk het uitrollen van de hoofd for-loop. De eerdere implementatie was natuurlijk een letterlijke uitschrijving van het algoritme, maar met stapjes van ieder keer seSize heeft aanzienlijk tijdwinst opgeleverd in vergelijking met de eerste code in Figuur 15. Voor een structurerend elementen kleiner dan 5 is het beter om de for-loops volledig uit te rollen. In regel 62 en 107: het maakt wel wat uit of de buffer R met een oplopend index wordt geschreven in het geval van grotere structurerende elementen. Maar als gaat om kleine buffers gaat dan is het niet zo erg. Voor het benaderen van arrays is het beter het in oplopende index te scannen want in het algemeen het is de manier hoe de data caches de arrays optimaliseren.
Figuur 17 : herkMaxValues2 eerste optimalisatie, herkMaxValues Finale optimalisatie .
24
4.3- Optimalisatie Lemire POGING 1: U.size() in regels (234,253) figuur 16, kan wel vermeden worden door een counter te nemen die ieder keer wordt opgehoogd bij een push en verminderd bij een pop_back en pop_front. Dus in plaats van (U.size()> 0); i.e. ((mask&backIndex)-(mask&frontIndex)>0); een simpele (counter>0) operatie.
POGING 2: De IF statement in ( Figuur 16 regel 230,249) kan weg gehaald worden door de eerste globaal Max/Min bij de start uit te rollen.
Figuur 18: uitrol eerste venster
POGING 3: Het is mogelijk om de pop_back niet iedere keer uit te voeren maar pas na de while loop, zo kan het aantal te verwijderen waardes uit de queue in één keer door de functie resize(n) weggehaald worden; resize knipt de queue bij index n en het kost evenveel tijd als een pop_back (zie void resize(int n) in de Appendix ).
25
Deze aanpak heeft als nadeel dat er een aantal operaties bij komen dat de geboekte winst negatief wordt beïnvloed. De reden is dat in plaats van if((image[i] >= image[U.back()])) uit te voeren wordt if((image[i] >= image[U[size-1]])) uigevoerd (zie eerder implementatie). Dit heeft als gevolg dat ieder keer U[size-1] wordt aangeroepen (Size-1 wordt berekend en samen met Qarray[(frI+(size-1))&mask]) en dit is dit is duurder dan gewoon U.back() = Qarray[mask & (bkI-1)]. Ondanks de neveneffecten van resize() is het mijn beste Lemire implementaties tot nu toe.
Figuur 19: Lemire Streaming implementatie zonder pop_back en uitrol eerste venster.
26
V-
TESTEN EN MEETRESULTATEN
Om de programma’s te testen heb ik het volgende testprotocol gebruikt:
5.1- MEETMETHODE: De metingen zijn in seconden uitgedrukt en het zijn gemaakt door middel van clock_t : met clock() methode voor en na het uitvoeren van de functie aanroep (te meten functie) en het verschil daartussen (verstreken CPU tikken) door CLOCKS_PER_SEC te delen wordt de tijd in seconden teruggegeven. De gemeten waardes zijn de totale duur van Max en Min filter. Klok resolutie Verder om te zorgen dat meetresultaten meer nauwkeurig zijn, is het noodzakelijk dat de metingen een aantal keer worden herhaald en vervolgens het gemiddelde daarvan te nemen. Door een blok van N metingen achterelkaar te doen wordt de nauwkeurigheid met factor N verbeterd, in plaats van met wortel √Z in het geval zou zijn bij het middelen van N aparte metingen. Hieronder staat uitgeschreven hoe de metingen worden uitgevoerd.
Figuur 20: Meetmethode
5.2- TESTPROTOCOL: 1234-
Vergelijken van iedere algoritmeversie met elkaar voor een gegeven datainput (100.000 data monsters ) voor verschillende structurerend element lengtes. Vergelijken van het beste algoritme op verschillende data soorten. Vergelijken van verschillende algoritmen met elkaar voor een gegeven datasoort. Data soorten zijn dalend, stijgend, randomwalk, ruis en periodiek. • •
• •
Ruisdata: de data is willekeurig gedistribueerd en is totaal onafhankelijk van elkaar. Bijvoorbeeld data[k] = random(0…1)*0.5 voor iedere k = 0…dataSize Randomwalk: Ieder data sample is gelijk aan de vorige data sample plus of minus een willekeurige waarde. Bijvoorbeeld : data[0]= 0 en data[k]= random(0…1)*(0.5) + data[k-1] voor k= 1… dataSize. Periodieke data: de data is afhankelijk van elkaar in een periodiek vorm. Bijvoorbeeld: data[k] = sinus(2*Pi *j/ periode) voor k = 0 … dataSize. Dalend (stijgenddata): bijvoorbeeld data[k] < data[k-1](data[k]>data[k-1]).
27
5.3- Test resultaten: vergelijking van verschillende optimalisatie technieken
5.3.1- GEOPTIMALISEERDE HGW VERSIES
1-
De executietijd voor eerste implementatie (herk1 zie figuur 15) stijgt lineair met de lengte van het SE. De executietijd van HerkFinal neemt juist af met toenemende lengte van het structurend element. HerFinal wordt in gemiddeld 0.0046 sec uitgevoerd. Herk1 is zeer traag want het zijn drie loops ingebed in een hoofd for loop waarin voor een beeld met lengte N en een structurerend element met lengte P worden er (N-P+1)*(3*(P-1)) instructies uitgevoerd de complexiteit is dus van de O(NP) (niet getoond)
2-
Door de hoofd-loop in stapjes van P te nemen is er een aanzienlijke tijdwinst geboekt. Voor 100.000 datapunten en voor filter lengte van 1 t/m 100 de tweede implementatie heeft gemiddeld 0.0053 seconden gescoord en de finale implementatie 0.0044 seconden waar R in oplopend orde wordt gescand en dat heeft als gezorgd een stabieler verwerkingstijd
8 vanHerk optimalisaties 7
msec
6
5
4
herkMax/MinValues2 hekMax/MinFinal
3
2
1 10
20
30
40
50 Se lengte
60
70
80
90
Figuur 21: Herk2 v.s. HerkFinal(finaal opti.)
3-
4-
Voor kleinere SE (<5) zijn de algoritmen trager want het aantal vergelijkingen wordt groter (3-4/p)en aan de andere kant binnen de hoofd-for-loop vindt er veel branches plaats. Gemiddelde HerkFilter operatie duurt 0.0046 seconden. Het was al verwacht want de complexiteit van het algoritme is niet afhankelijk van de Se lengte.
28
5.3.2- GEOPTIMALISEERDE LEMIRE IMPLEMENTATIES
1-
In de grafiek is duidelijk te zien dat de verbeteringen aanzienlijk tijdwinst hebben opgeleverd. Zo is de tweede implementatie met een arrayQueue 17% sneller geworden dan de eerste implementatie waar deqeue wordt gebruikt. En als pop_back() helemaal vermeden is in het geval van de finale implementatie waar resize() wordt gebruikt is de prestatie met 41% verbeterd. Zoals eerder is geconstateerd. Is het algoritme zeer goed in dalende en stijgende data gevolgd door periodiek data. De ruis en random walk data zijn het ergste geval voor het algoritme verder blijft de verwerkingstijd als verwachtconstant als SE groter wordt.
2-
.
x 10 4.6
4.4
4.2
4
Lemire optimalisaties
ms e c
3.8
3.6
3.4
3.2
3
Lemire: Deque Lemire: ArrayQueu Lemire: ArrayQueu met resize()
2.8
2.6 10
20
30
40
50 SE lengte
60
70
80
90
Figuur 22: geoptimaliseerde Lemire implementatie. De meet resultaten voor verschillede soort data
5.4- ALGORITMEN PRESTATIETESTEN VOOR VERSCHILLENDE SOORT DATA
Hier onder een aantal grafieken met de prestaties van de geoptimaliseerde versies van Herk en Lemire naast het naïeve en maxline algoritmen. Per data worden er metingen getoond van hoe lang ieder programma er overdoet om Maxminfilter te berekenen. De data input is een array van 100.000 waarden en de lengte van SE neemt waarden tussen 1 en 100. 1-
De geïmplementeerde algoritmen zijn pas interessant als SE groter is dan 5. Voor kleinere SE’s is de Naïeve implementatie (Figuur 12) sneller want het is het hoofd voorloop die voornamelijk wordt uitgevoerd en dat geldt ook voor de andere implementaties maar wel met extra boolean instructies en meer dan 2 loops er bij.
2-
Zo als bedoeld is het Maxline algoritme alleen geschikt voor i.i.d data maar lang niet optimaal als Lemire of Herk algoritmen. Voor andere soort data(niet willekeurige data) stijgt de gebruikte tijd lineair met de lengte van SE.
29
0.06
DalendData
Lemire Herk Maxline Naïef
0.05
seconden
0.04
0.03
0.02
0.01
0
0
10
20
30
40
50 Se lengte
60
70
80
90
100
Figuur 23: Prestaties verschillende algoritmen voor dalend data
3-
4-
Zoals verwacht is Lemire implementatie beter dan alle andere algoritmen maar wel met een klein gemiddelde 0.017 sec voorsprong op Herk. Theoretisch gezien is het Lemire algoritme beter maar de implementatie daarvan is niet efficiënt genoeg geworden door bijzondere push en pop operaties en dat is nadelig voor de prestatie. In de implementatie van het Herk algoritme is eenvoudige datastructuur gebruikt.
0.06
Lemire Herk Maxline Naïef
0.05
PeriodiekData
seconden
0.04
0.03
0.02
0.01
0
0
10
20
30
40
50 Se lengte
60
70
80
90
100
Figuur 24: Prestaties verschillende algoritmen voor periodiek data Voor periodieke data is Lemire implementatie ook weer de beste maar net als in stijgend/dalend data is het verschil met Herk implementatie heel klein.
30
0.08
0.07
Lemire Herk Maxline Naïef
0.06
RandomWalkData
Seconden
0.05
0.04
0.03
0.02
0.01
0
0
10
20
30
40
50 Se lengte
60
70
80
90
100
Figuur 25: Prestaties verschillende algoritmen voor randomwalk data 0.06
RuisData Lemire Herk Maxline Naïef
0.05
seconden
0.04
0.03
0.02
0.01
0
0
10
20
30
40
50 Se lengte
60
70
80
90
Figuur 26: Prestaties verschillende algoritmen voor Ruis data Voor ruis-data is LemireFinal enigszins minder gaan presteren dan HerkFinal. Want dat is ook weer voor het Lemire algoritme de worst-case en het is het gemiddelde van 3 vergelijkingen per element (0.071sec) en wat betreft Herk algoritme ietsje minder 3-4/p (0.064secs). Naïef
31
VI-
Twee-Dimensionale implementatie van het Lemire Max/Min filter
Om 1D Max-Min filter uit te breiden naar 2D Max-Min filter maken wij gebruik van de associatieve eigenschap van de dilatie en de “chain rule” van de erosie. Een SE in de vorm van een rechthoek (H,B) kan een dilatie zijn van een horizontale SE van breedte B door een verticale SE van hoogte H d.w.z. Als SE =(S⊕1 ⊕ S2) dan is A⊕SE= (A⊕S1)⊕S2, hetzelfde geldt voor de erosie want weer door de “chain rule” eigenschap hebben wij AΘ(S1⊕S2)= (AΘS1)ΘS2. In de praktijk, stel dat wij een beeld willen dilateren met rechthoekig SE dan dilateren wij ten eerste het beeld lijn per lijn met 1D SE S1 en daarna het resultaat transponeren, weer dilateren met de verticale SE S2 en het uiteindelijke resultaat terug zetten door nog een transpositie.
Voorbeeld: Gegeven een 3x3 beeld en een 2x2 s.e
1) Eerst het s.e. splitsen in twee delen rij en colom
2)
Het beeld dilateren rij per rij met
3) Het resultaat uit 2) transponeren voor de dilatatie met de s.e. colom
4)
Het resultaat uit 3) weer transponeren en dat is uiteindelijk de max filter
5)
Eind resultaat
Figuur 27: voorbeeld implementatie van 2D dilatie van een grijswaarde beeld
32
ÉÉN-DIMENSIONALE REPRESENTATIE VAN EEN 2D BEELD Het eerste probleem bij het toepassen van 1D filter op een 2D beeld is efficiëntie, er moet twee keer getransponeerd worden en dat is niet handig. Wij kunnen dit omzeilen door een 2D beeld te representeren als 1D-array. Hier worden rijen en kolommen benaderd door twee te definiëren functies R en C. Stel dat A een 2D beeld met N rijen ([= )i=1…N en C kolommen Als [A = {\]A ,…,\^A } , [C = {\]C ,…,\^C } …,[_ = {\]` ,…,\^` }, dan is 1D representatie van A de 1D array B = [[A , [C , … , [` ]. De twee genoemde functies zijn als volgt gedefinieerd:
R(r,i, C) = r * C + i : geeft element i uit Rij r terug. C(c,i, C) = i * C + c : geeft element I uit colom c terug [16]. Randcondities: padding methoden Tot nu toe doen alle behandelde algoritmen niets aan de problemen bij de rand. Als het filtervenster bij de rand van het beeld aankomt, blijven er ongedefinieerde pixels over die een speciale behandeling nodig hebben. Er zijn veel mogelijkheden om de rand op te vullen zonder dat het resulterende 1D beeld (seSize-1) kleiner is dan de oorspronkelijk data input. Het kan bijvoorbeeld door randpixels te herhalen of reflectie van het beeld aan de rand. Ik heb zelf voor de onderstaande methode gekozen.
Toepassen op Grijsbeelden: Ik heb MATLAB gebruikt om een beeldbestand te lezen en te converteren naar grijsbeeld ( UINT8,16 of 32). Vervolgens wordt het grijsbeeld opgeslagen in een C-header file in de vorm van een 1D-array.De 1D-array wordt geladen en verwerkt. Het resultaat is weer een 1D array dat MATLAB kan lezen en tonen. Het is een hele omweg. In MATLAB kan een C++ programma zodanig worden geïmplementeerd dat het door MATLAB Mex-compiler wordt vertaald en meteen inzetbaar is in MATLAB. Zie de code in de Appendix.
Voorbeelden toepassingen op 2d grijsbeelden met 10x10 SE.
50
100
150
200
(a)
250
(b)
(c)
(d) (e) Figuur 28: (a)input beeld, (b) dilatie, (c) erosie, (d)opening en (e) sluiting door gebruik van eigen 2d implementatie.
33
CONCLUSIE Algoritmen voor beeld verwerking zijn vaak moeilijk te vertalen naar efficiënte code. En als dat eenmaal is gelukt, blijft de vraag of het geprogrammeerde algoritme inderdaad aan de verwachte efficiëntie en snelheid voldoet zoals beschreven wordt in de theorie. De Lemire Implementatie overtreft in het algemeen de ander implementaties qua snelheid als de data monotoon of periodiek is omdat het vensterfilter wordt niet in zijn geheel doorgelopen zoekend naar het nieuwe globale maximum of minimum waarde. Wat het niet het geval is als de data i.i.d is. Zo is bijvoorbeeld, in het geval van ruisdata presteert Lemire implementatie minder goed want het is ook weer de worst-case scenario van het algoritme waar 3 vergelijkingen per element moeten uitgevoerd worden tegen 3-4/p vergelijkingen per element in het geval van Herk algoritme. In de grafiek is duidelijk te zien dat Herk algoritme nog steeds beter kan gebruikt worden als het gaat om het filteren van beelden van ruis want het algoritme is minder gevoelg voor i.i.d data. En wat betreft de implementaties van de algoritmen is van het vanHerk algoritme ook sneller geworden dankzij de eenvoud van de gebruikte datastructuren en omzeilde “branche” operaties terwijl het Lemire algoritme een complexere datastructuur en operaties gebruikt zoals push, pop back en front.
x 10 6.5
6
Lemire: monotoon en periodiek Lemire: ruis data vanHerk: monotoon en periodiek vanHerk: ruis data
msec
5.5
5
4.5
4
3.5 0
10
20
30
40
50 SE lengte
60
70
80
90
Figuur 29: Prestaties van het Lemire en van Herk algoritmen. Uiteidelijk we kunnen onze algoritme samenstellen uit drie verschillende algoritmen 12-
Als SE is kleiner dan 5 dan gebruik het naïeve algoritme En als SE groter dan 5 dan a. Als de data niet i.i.d is dan is lemire algoritme de beste b. Anders wordt van Herk algoritme gebruikt
34
x 10 7
6
msec
5
4
3 Lemire vanHerk Naief 2
1
2
3
4
5
6
7
8
9
10
SE lengte
Figuur 30: Prestaties van het Lemire, van Herk en naïeve algoritmen voor ruis data met S.e. kleiner is dan 10. x 10 7 Naive vanHerk Lemire 6
msec
5
4
3
2
1
2
3
4
5
6
7
8
9
10
SE lengte
Figuur 31: Prestaties van het Lemire, van Herk en naïeve algoritmen voor monotone data met s.e. kleiner is dan 10.
Wij kunnen nu als keuze maken dat bij kleine SE de naïeve mmfilter te gebruiken en voor grotere SE als de datainput i.i.d dan is van Herk
35
VII- VERDER WERK
1-
Het Lemire algoritme kan ook gebruikt worden voor kleinere SE’s Maar dit keer met het uitrollen van het algoritme voor verschillenden kleinere structurerend elementen en voor grotere SE’s is de vraag of het mogelijk is om de loops in het algoritme met grotere stappen te implementeren
1-
Lemire 2D Min/Max filter
a- Met 2D Max/Min filter worden er twee keer zoveel vergelijkingen per element gemaakt als in 1D versie. Uit een paar testen is gebleken dat voor kleinere 2D-SE (5x5) duurt 2D Max/Min filter tenminste 5 keer langer. Kan het beter?
b- In het Herk algoritme is er aangenomen dat de lengte van het beeld een multiple moeten zijn van de lengte van SE, verder moet SE oneven lengte hebben. Deze beperkingen heeft het Lemire algoritme echter niet. Wel is er een algemeen probleem bij beide algoritmen dat ze de filteroperaties met een willekeurige hoek niet ondersteunen in het geval van 2Dfilter. P.Soille [16] heeft al aangegeven hoe dat opgelost kan worden voor het Herk algoritme en ik stel voor hetzelfde te doen met het Lemire algoritme. Om een 2D beeld te scannen door (lijn) een SE in een willekeurige hoek kunnen een aantal pixels overgeslagen worden of meerder keer meegenomen worden in de berekening. Dat komt door de disreet representatie van een lijn in een discreet vlak (beeld). Om dat te kunnen vermijden moet een discreet lijn aan een aantal voorwaarden voldoen P.Soille [16]. • Voor bepaalde hoeken mag een discreet lijn niet meer dan één pixel hebben per rij of kolom. • Om het hele beeld te scannen door een disreet lijn in een bepaalde richting, moet de start positie goed gekozen worden.
36
VIII- APPENDIX
1- ARRAYQUEU class Queu { private: int *ar,frI,bkI,mask; public : Queu(int b){ // b mustbe 2 4 16 32 ... frI=0; bkI=0; ar =new int[b]; int n = 0; while (b) { n++; b >>= 1; } b = 1 << n; mask = b-1; } inline int &operator []( int); //~Queu() {bkI=0;frI=0;delete ar;} inline int size() { return (mask&bkI)-(mask&frI); } inline int front(){ return ar[mask & frI]; } inline int back(){ return ar[mask & (bkI-1)]; } inline void push_back(int i){ ar[mask&bkI]=i; bkI += (1& mask) ; } inline void resize(int i){ bkI =(frI+i)&mask; } inline void pop_front(void){ frI+=(1&mask); } Inline void pop_back(void){ bkI-=(1&mask); } }; inline int & Queu::operator []( int i){ return ar[(frI+i)&mask]; }
37
2-
ARRAY IMPLEMENTATIE VAN EEN VECTOR
template class FloatArrayVector{ T* F; int SIZE; public: FloatArrayVector(int n= 10) { SIZE = n ; F =new T[n]; } inline int size(void){return SIZE;} T& operator []( int index) { return F[index]; } }; 3-
MAX en MIN : verschillende implementaties inline static uint max( const int a,const int b) {return (uint) (0.5 * (abs ((a-b)) + (a-b) ) + b);} inline static uint min(const int a ,const int b) {return (uint)(-0.5 * (abs((a-b)) + (a-b) ) + a); } static inline int maxi( int a, int b) { b = a-b; a -= b & (b>>31); return a; } static inline bool isBig(int a,int b){ b=a-b; return true&(~((b)>>31)); } static inline int mini( int a, int b) { b = a-b; a -= b & (~(b>>31)); return a; } inline double max2(double a,double b) { if(a > b) return a; return b; } inline double min2(double a,double b) { if(a < b) return a; return b; }
38
4-
Object gerionteerd implementatie van min/max filter class minmaxfilter { public: virtual DataArrayType& maxvalues()=0; virtual DataArrayType& minvalues()=0; virtual ~minmaxfilter(){}; }; void display(minmaxfilter& image1) { display(image1.maxvalues()); display(image1.minvalues()); } bool isEqual(minmaxfilter& image1, minmaxfilter& image2 ) { return isEqual(image1.maxvalues(), image2.maxvalues()) & \ isEqual(image1.minvalues(), image2.minvalues()); }
5-
Lemire 2D Max Filter void lemiremax(DataArrayType &image,DataArrayType &outPutU1,int a){ DataArrayType maxValues(a- width+1); int Rows= a/Cols; //====== Dilation of the rows===== Queu U(width); for(int row =0;row image[rowelt(row,(i-1),Cols) ]) { U.pop_back();--size; while(U.size() > 0) { if(image[rowelt(row,i,Cols)] <= image[rowelt(row,U.back(),Cols)] ) break; U.pop_back();--size; }} U.push_back(i);++size; } int tmp= Cols-width;//image.size()-width; for( i = width; i image[rowelt(row,i-1,Cols)]) { U.pop_back();--size; while(U.size()>0) { if(image[rowelt(row,i,Cols)] <= image[rowelt(row,U.back(),Cols)]) break; U.pop_back();--size; } } U.push_back(i);++size; if(i== width+U.front()) { U.pop_front();--size; }
39
} maxValues[rowelt(row,tmp,Cols)] = image[rowelt(row,U.front(),Cols)]; U.reset(width); } //==== Dilation of the columns==== Queu U2(high); DataArrayType outPutU(maxValues.size()-high+1); for(int row =0;row maxValues[colelt(row,(i-1),Cols) ]) { U2.pop_back();--size; while(size > 0) { if(maxValues[colelt(row,i,Cols)] <= maxValues[colelt(row,U2.back(),Cols)] ) break; U2.pop_back();--size; } } U2.push_back(i);++size; } int tmp= Rows-high;//maxValues.size()-high; for( i = high; i maxValues[colelt(row,i-1,Cols)]) { U2.pop_back();--size; while(U2.size()>0) { if(maxValues[colelt(row,i,Cols)] <= maxValues[colelt(row,U2.back(),Cols)]) break; U2.pop_back();--size; } } U2.push_back(i);++size; if(i== high+U2.front()) { U2.pop_front();--size; } } outPutU[colelt(row,tmp,Cols)] = maxValues[colelt(row, U2.front(),Cols)]; U2.reset(high); } outPutU1= outPutU; }
6-
P en C functies int rowelt(int r,int i,int C) { return (r) * C + (i) ;}// element i uit rij r // C aantal colo // R aantal rijen // c=[1:C] en i=[1:R] int colelt(int c,int i,int C) {return (i) * C + c ;}// element i uit colom c
40
IX- REFERENCES
[1]. Brookes, D.M. “Algorithms for Max and Min Filters with Improved Worst-Cases Performance”, IEEE Trans on Circuits & Systems -II, Vol 47(9), pp930-935, Sept 2000. [2]. D. Coltuc and I. Pitas, “On fast running max-min filtering,” IEEE Trans.Circuits Syst. II, vol. 44, pp. 660–664, Aug. 1997. [3]. D. Lemire, 2006. Streaming maximum-minimum filter using no more than three comparisons per element Source ISSN:1236-6064 . [4]. Douglas, S. C., 1996. An efficient algorithm for running max/min calculation.ISCAS’96 2. [5]. Gil, J. Y., Kimmel, R., 2002a. Data filtering apparatus and method. US Patent Number 6,952,502. [6]. Gil, J. Y., Kimmel, R., 2002b. Efficient dilation, erosion, opening, and closing algorithms. IEEE Trans. Pattern Anal. Mach. Intell. 24 (12), 1606–1617. [7]. Gil, J.,Werman, M., 1993. Computing 2-d min, median, and max filters. IEEETrans. Pattern Anal. Mach. Intell. 15 (5), 504–507. [8]. Henk J.A.M Heimans,2000 . Mathematical Morphology,ISBN: 1586030566. [9]. .I. Pitas, “Fast algorithms for running ordering and max/min calculation,”IEEE Trans. Circuits Syst., vol. 36, pp. 795–804, 1989. [10]. J. Bresenham Algirithm for computer control of digital plotter. IBM System journal,4:2530,1965. [11]. Keogh, E., Ratanamahatana, C., 2005. Exact indexing of dynamic time warping. Knowledge and Information Systems 7 (3), 358–386. [12]. Lemire, D., Brooks, M., Yan, Y., 2005. An optimal linear time algorithm for quasi-monotonic segmentation. In: ICDM’05. pp. 709–712. [13]. Pitas, I., 1989. Fast algorithms for running ordering and max/min calculation. IEEE Transactions on Circuits and Systems 36 (6), 795–804. [14]. R. v.d. Boomgaard (1992) Mathematical Morphology: Extensions towards Computer Vision. PhD thesis, University of Amsterdam. [15]. Richard M. Stallman ,2003,Using the GNU compiler collection voor GCC4.3.0 ,GNU press,Free Software Foundation [16]. Soille, P., Breen, E.J. and Jones, R. (1996) Recursive implementation of erosions and dilations along discrete lines at arbitrary angles . IEEETrans. Pattern Anal Mach Intell, 18(5), 562—567. [17]. Van Herk, M., 1992. A fast algorithm for local minimum and maximum filters on rectangular and octagonal kernels. Pattern Recogn. Lett. 13 (7), 517–521. Hugues Talbot ,2006. Image Analysis and Processing Mathematical Morphology [18]. {h.talbot}@esiee.fr
[19]. D.Z. Gevorkian,J.T. Astola, and S.M. Atourian, 1997, “Improvement Gi l Weerman Algorithm for running Min Max Filters”,IEEE T.P.A.M.I vol.19 99.526-529 [20]. http://www.aanda.org/index.php?option=article&access=standard&Itemid=129&url=/art icles/aa/ref/2006/38/aa4852-06/aa4852-06.html
41
