Einführung in die Bildverarbeitung und Bildanalyse mit R (Teil 1)

Vorwort

Bitte folgenden Code ausführen um notwendige Bibliotheken zu laden. Nicht erforderlich bei nativer R Software.

• R+ ist nahezu R Syntax mit einigen Ergänzungen die die Programmierung und lesbarkeit von Programmen erleichtern.

• Nicht alle R Packages sind in R+ enthalten, und es kann Abweichungen geben. Hier werden vor allem das imager, dbscan, und fft Package verwendet.

• Dennoch lassen sich die meisten öffentlich verfügbaren algorithmischen Lösungen und Hilfestellungen mit R+ umsetzen.

Literatur

• Lese "Learning R Programming" Kun Ren,2023

Ziele

• Grundlegendes Verständnis und Anwendung eines vollständigen Datenpfades, d.h. von der (synthetischen) Erstellung von Bildern mit Merkmalen bis zur Merkmalsextraktion.
• Die detaillierten Schritte sind sekundär. Wichtig ist das Verstehen des Ablaufs, der Schritte und der Verfahren (Algorithmen).

Synthetische Bilder

• Synthetische Bilder können mit Matrizen oder Bildobjekten erstellt werden.
• Matrizen sind immer monochromatische Grauwertbilder (zweidimensional)
• Eine Matrix wird mit matrix(init,nrow,ncol) erzeugt
• Alternativ kann als Initialwert auch ein Vektor (z.B. randomisiert mit runif(nrow*ncol)*scale) übergeben werden
• Es sollten für Bilder kompaktierte und typisierte Arrays verwendet werden, d.h., mit dem Zusatz mode='uint8'. Gängige Datentypen sind uint8, uint16, und uint32. bei einige Berechnungen und Operationen muss float32 verwemdet werden.
• Im Package geometry gibt es eine Reihe von Funktionen um geometrische Objekte in Matrizen und Bildern zu erzeugen, wie z.B. (x ist eine Matrix oder ein Bild):
  • draw.line(x,value,from,to)
  • draw.rectangle(x,value,center?,left?,top?,width?,height?,angle?,filled?)
  • draw.circle(x,value,center,width?,height?,radius?,angle?,filled?)

Aufgabe. Erzeuge eine Funktion pores.generate die ein künstliches Porenbild (im.pores) mit N(=20) Poren und randomisierten Parametern (wie oben) erzeugt. Zunächst lassen wir Überlappungen zu. Breite und Höhe sollen auf den Bereich [5,15] Pixel begrenzt sein.

aW0ucG9yZXMgPSBtYXRyaXgoMTAsMTAwLDEwMCxtb2RlPSd1aW50OCcpCmZ1bmN0aW9uIHBvcmVzLmdlbmVyYXRlIChpbSxOKSB7CiAgZm9yIChpIGluIDE6TikgewogICAgcCA9IHJ1bmlmKDIpKjEwMAogICAgYSA9IHJ1bmlmKCkqMzYwCiAgICB3ID0gNStydW5pZigpKjEwCiAgICBoID0gNStydW5pZigpKjEwCiAgICBJID0gODArcnVuaWYoKSoxMDAKICAgIGRyYXcuY2lyY2xlKGltLEksY2VudGVyPXAsYW5nbGU9YSx3aWR0aD13LGhlaWdodD1oLGZpbGxlZD1UUlVFKQogIH0KICBpbQp9CnBvcmVzLmdlbmVyYXRlKGltLnBvcmVzLDIwKQpwbG90KGltLnBvcmVzKQ==

Da wir nachfolgend über eine Pixelanalyse die elliptischen Poren rekonstruieren wollen ist eine Porenüberlappung nicht sinnvoll (wenn auch möglich in z.B. Schliffbildern).

• Eine Möglichkeit besteht darin bevor man die Ellipse erzeugt deren rechteckige Boundignbox zu ermitteln und nach Überlappungen der vorherigen
• Eine Boundingbox ist eine Liste mit mindestens {left,top,width,height} Elementen. Man wählt eine quadratische bbox da die Ellipse noch rotiert werden kann, d.h. width=height=max(elli.w,elli.h).
• Mit der Funktion bbox.intersect (imager) kann man die Überlappungsbbox von zwei einzelnen bboxes berechnen. Wenn sie nicht überlappen wirf FALSE zurückgegeben.
• Das wird für alle bisherigen (gespeicherten) bboxes der Ellipsen geprüft, was einen Prüfvektor unter Verwedung der map Funktion ergibt
• Die all Funktion prüft ob alle Element TRUE Werte haben, erst dann ist dir Prüfung erfüllt (Siehe nachfolgendes Beispiel)

Aufgabe. Erzeuge eine Funktion pores.generate.overlapfree die ein künstliches Porenbild (im.pores.synth) mit N(=20) Poren und randomisierten Paremtern (wie oben) ohne Überlappungen erzeugt.

aW0ucG9yZXMgPSBtYXRyaXgoMTAsMTAwLDEwMCxtb2RlPSd1aW50OCcpCmZ1bmN0aW9uIHBvcmVzLmdlbmVyYXRlLm92ZXJsYXBmcmVlIChpbSxOKSB7CiAgYmJveC5saXN0ID0gbGlzdCgpCiAgZm9yIChpIGluIDE6TikgewogICAgcCA9IHJ1bmlmKDIpKjEwMAogICAgYSA9IHJ1bmlmKCkqMzYwCiAgICB3ID0gNStydW5pZigpKjEwCiAgICBoID0gNStydW5pZigpKjEwCiAgICBJID0gODArcnVuaWYoKSoxMDAKICAgIHMgPSBtYXgodyxoKQogICAgYiA9IHtsZWZ0PXBbMV0tcy8yLHRvcD1wWzJdLXMvMix3aWR0aD1zLGhlaWdodD1zLGNlbnRlcj1wLGFuZ2xlPWEsbWFqb3I9dyxtaW5vcj1ofQogICAgbm9vdmVybGFwID0gYWxsKG1hcChiYm94Lmxpc3QsCiAgICAgICAgICBmdW5jdGlvbiAoYjEpIHsgaXNGQUxTRShiYm94LmludGVyc2VjdChiLGIxKSkgfSkpCiAgICBpZiAobm9vdmVybGFwKSB7CiAgICAgIGRyYXcuY2lyY2xlKGltLEksY2VudGVyPXAsYW5nbGU9YSx3aWR0aD13LGhlaWdodD1oLGZpbGxlZD1UUlVFKQogICAgICBwdXNoKGJib3gubGlzdCxiKQogICAgfQogIH0KICBiYm94Lmxpc3QKfQpwb3Jlcy5nZW5lcmF0ZS5vdmVybGFwZnJlZShpbS5wb3JlcywyMCkKcGxvdChpbS5wb3Jlcyk=

Das Histogramm

Ein Histogramm gibt eine diskretisierte verteilung von Werten eines Vektors, einer Matrix, eines Bildes, oder jeglicher n-dimensionaler Datenstruktur wieder.

Dabei wird der zu erwartende Wertebereich ⟨x_min,x_max⟩=⟨x₁,x₂⟩ in k Intervalle der Breite Δ=(x_max-x_min)/k unterteilt, d.h., {⟨x₁,x₁+Δ⟩, ⟨x₁+Δ,x₁+2Δ⟩, .. , ⟨x₁+(k-1)Δ,x₂⟩}.

Die Anzahl der Intervalle bestimmt die Granularität (breaks Argument). Die programmatische Rückgabe eines Histogramms (plot=FALSE) liefert eine Liste:

list [5] (
 breaks:vector,
 counts:vector,
 mids:vector,
 density:vector,
 equidist=boolean
)

Das Element density ist ein auf die Summe 1 normierter Vektor von count. Mit dem xlim Argument kann der Start und Endwert der Werteverteilung begrenzt werden.

Aufgabe. Welche Anzahl von Intervallen ist sinnvoll bzw. wie sollte man die Intervallgröße (breaks) wählen und wovon hängt diese Wahl ab? Vergleiche das erzeugte Porenbild und das Rauschbild. Ändere die runif(N,min,max) Funktion in eine rnorm(N,min,max) Funktion ab. Wie unterscheidet sich die Verteilung?

Rauschen

Echte Bilder sind verrauscht. Detekorrauschen ist i.A. gaußverteilt um einen Mittelwert des Messwertes. Wir wollen jetzt Rauschen zu unserem synthetischen Bild hinzufügen.

Aufgabe. Erzeuge eine Funktion image.noise die zum einem Bild Rauschen hinzufügt: Additiv (wie oben), jeweils wählbar mit runif(n,min,max) oder rnorm(n,min,max) Generatoren, und multiplikativ. Die Höhe des Rauschens wird über gamma festgelegt. Erzeuge ein verrauschtes Porenbild mit γ=0.1 und 0.3. n muss gleich length(im) gewählt werden. Weiterhin soll der Dynamikbereich eingeschränkt werden (Skalierung) und ein Hintergrundoffset hinzugefügt werden (nach der Skalierung). Dadurch wird Unter/Überlauf durch Rauschhinzufügung vermieden.

ZnVuY3Rpb24gaW1hZ2Uubm9pc2UgKGltLHNjYWxlLG9mZnNldCxnYW1tYSxub2lzZWdlbixhZGRpdGl2ZSkgewogIGltPWltKnNjYWxlK29mZnNldAogIGlmIChhZGRpdGl2ZSkgewogICAgaW0ubm9pc3k9aW0rKG5vaXNlZ2VuKGxlbmd0aChpbSksLWdhbW1hLGdhbW1hKSoyNTUpCiAgfSBlbHNlIHsKICAgIGltLm5vaXN5PWltKigxK25vaXNlZ2VuKGxlbmd0aChpbSksLWdhbW1hLGdhbW1hKSkKICB9CiAgaW0ubm9pc3kKfQpwbG90KGltLnBvcmVzKQpwbG90KGltYWdlLm5vaXNlKGltLnBvcmVzLDAuNSwyMCwwLjMscnVuaWYsRkFMU0UpKQpwbG90KGltYWdlLm5vaXNlKGltLnBvcmVzLDAuNSwyMCwwLjEscnVuaWYsVFJVRSkpCnBsb3QoaW1hZ2Uubm9pc2UoaW0ucG9yZXMsMC41LDIwLDAuMyxybm9ybSxGQUxTRSkpCnBsb3QoaW1hZ2Uubm9pc2UoaW0ucG9yZXMsMC41LDIwLDAuMSxybm9ybSxUUlVFKSk=

Aufgabe. Welches Rauschen hat einen größeren Effekt? Additiv oder Multiplikativ. Weclhe Art würde physikalischen Rascuhquellen entsprechen?

Kantendetektion der Poren

Um nun die Poren in einem Bild charakterisieren zu können müssen deren Konturen berechnet werden. gradientebasierte Filter sind dafür geeignet.

• Neben dem Sobel Filter ist der iterative Canny Detektor (cannyEdges(x,t1,t2)) ein geeignetes Verfahren um die Konturen der Poren in einem neuen Bild zu markieren.
• Dann brauchen wir die markierten Pixelkoordinaten [x,y]. Die Menge aller markierten Pixel (pixset) kann man mit as.pixset(Boolean matrix) extrahieren.
• Eine relationale Operation auf eine Matrix angewendet liefert eine gleichgroße Boolesche Matrix, z.B. m>1. as.pixset liefert eine Liste mit Vektoren [x,y] für alle Pixel die den Wert 1/ TRUE besitzen. Wählt man as.matrix=TRUE werden die Koordinaten der Pixel in einer Matrix geliefert.

Aufgabe. Bewerte das Ergebnis des Canny Edge Detektors ohne und mit verschieden großen Rauschen additiv und multiplikativ). Wann werden Artefakte markiert (also Rauschen), wie muss man die Parameter t₁ und t₂ verändern um die Artefakte zu reduzieren. Wo ist die Grenze erreicht?.

Clustering und Einhüllende

• Die Koordinatenliste aller markierten Pixel ist als eine Punktemenge im 2-dim. Raum anzusehen
• Die Liste ist i.A. nicht semantisch geordnet, d.h. aufeinanderfolgende Elemente gehören nicht unbedingt zu der gleichen Pore (Kontur)
• Ein Clustering Algorithmus wird nun verwendet um die Punktemenge in Gruppen einzuteilen die eine Nachbarschaft bilden
• Dazu wird hier der DBSCAN Algorithmus verwendet (siehe Vorlesungsmaterial)
• Im Idealfall gehören alle Punkte einer Gruppe zu der EInhüllenden einer Pore (kann, muss aber nicht)
• Und es kann Punkte geben die zu keiner Gruppe gehören (Noise)

• Wir ignorieren zunächst dieses Problem und verarbeiten alle Punkte einer Gruppe (daher im.pores.edges.hulls = im.pores.edges.clusters)

Ellipsenanpassung als Minimierungsproblem

• Die Algorithmik der Ellipsenanpassung ist im Vorlesungsmaterial nach zu lesen.
• Iterativ wird für jede Punktegruppe eine Ellipse berechnet
• Das Ergebnis soll als Overlay des Punktebildes visualisiert werden

• Alle Ellipsenparameter werden in einer Datentabelle eingetragen

Aufgabe. Prüfe die Übereinstimmung der angepassten Ellipsen mit dem Punktebild in Abhängigkeit mit dem Rauschanteil (variieren). Welche Beobachtungen kann man machen? Vergleiche die Ergebnisse der vollen Punktegruppe mit nachfolgender Hüllenpunkte Berechnungen (mit dem Graham Algorithmus). Gibt es Unterschiede?

Statistische Analyse

Abschließend wird das Ergebnis der Ellipsenanpassung statistisch ausgewertet. Dabei sind verschiedene Aggregate wichtig:

• Zentrumspunkt aller Ellipsenzentren (gleichmäßige Verteilung der Poren im Bild sollte den Mittelpunkt des Bildes ergeben)
• Mittelwerte und Streuung (std) der Axenlängen
• Mittelwert und Streuung des Ellipsenwinkels (Ausrichtung)

Aufgabe. Ergänze die statistische Analyse um den Median (fivenum benutzen) und die Streuung (sd) für die Axenlängen und den Winkel. Trage Ergebnisse unten ein (Copy und Paste der obigen Ausgabe).