[JS]
[DE]
[USE] plot
[PLUGIN] ml.plugin
[PLUGIN] math.plugin
[PLUGIN] image2.plugin
[PLUGIN] popup.plugin
[PLUGIN] R.plugin
[PLUGIN] file.plugin
[PLUGIN] help.plugin
[PLUGIN] button.plugin
[TITLE] R - Einführung in die Bildverarbeitung und Analyse
[AUTHOR] Stefan Bosse
[VERSION] 11.2023

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

[FUN]
```js
Plot.create=function (data,options) { console.log(options); drawPlot(0,Plot(data,options))}
new HelpUI({
  welcome:'Enter search pattern (with optional wildcards *)',
  label:'HelpBot',
  callback:function (text) {
    var result = R.help(text)
    return result?result:'Sorry, found nothing helpful.'
  }
})
```

[TOC]

## Vorwort

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

[R] R Set-up { lines:2; height:5 }
```R
use math,plot,imager,geometry
dev.new(width=500)
print("Libraries opened.")
print(R.version)
```

>! Es wird in diesem Kurs der R-Dialekt r+ verwendet. R+ ist eine Reimplementierung von R in JavaScript und läuft direkt in Web Browsern und node webkit Applikationsprogrammen. Es ist keinee Softwareinstallation erforderlich.

- 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](http://edu-9.de/uploads/Lehre/asd3k/tutorials/9781785880629.pdf)

## 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?)`

[R] Randomisiertes Bild { lines:20; height:5 }
```R
im1 = matrix(runif(10000)*255,100,100,mode='uint8')
plot(im1)
```

[R] Geometrisches Bild 100 × 100 Pixel als Matrix (Hintergrundwert ist hier 10) { lines:20; height:20 }
```R
im1 = matrix(10,100,100,mode='uint8')
  p = runif(2)*100
  a = runif()*360
  w = 20+runif()*10
  h = 10+runif()*10
  I = 80+runif()*100
  draw.circle(im1,I,center=p,angle=a,width=w,height=h,filled=TRUE)
  p = runif(2)*100
  a = runif()*360
  w = 20+runif()*10
  h = 10+runif()*10
  I = 80+runif()*100
  draw.rectangle(im1,I,center=p,angle=a,width=w,height=h,filled=TRUE)
plot(im1)
```

[TODO]
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.

[R] Synthetisches Porenbild { lines:20; height:20 }
```R
im.pores = matrix(10,100,100,mode='uint8')
function pores.generate (im,N) {
  for (i in 1:N) {
    # TODO
  }
}
pores.generate(im.pores,20)
plot(im.pores)
```

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)

[R] Auscchnitt für den Algortihmus zur Vermeidung von überlappenden Ellipsen durch Suche von bbox Überlappungen (Intersections) { lines:20; height:20 }
```R
  bbox.list = list()
  # Random
  p = runif(2)*100
  a = runif()*360
  w = 20+runif()*10
  h = 10+runif()*10
  I = 80+runif()*100
  s = max(w,h)
  b = {left=p[1]-s/2,top=p[2]-s/2,width=s,height=s,center=p,angle=a,major=w,minor=h}
  # ...
  nooverlap = all(map(bbox.list,
          function (b1) { isFALSE(bbox.intersect(b,b1)) }))
  print(nooverlap)
  if (nooverlap) {
    draw.circle(im1,I,center=p,angle=a,width=w,height=h,filled=TRUE)
    push(bbox.list,b)
  }
  plot(im1)
```

>! Da randomisiert neue Kandidaten erzeugt werden kann man obiges Beispiel in eine Zählschleife  einbinden, in der die Anzahl der Versuche eine nicht überlappende Ellipse zu erzeugen begrenzt wird.

[TODO]
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.

[R] Synthetisches Porenbild { lines:20; height:20 }
```R
im.pores = matrix(10,100,100,mode='uint8')
function pores.generate.overlapfree (im,N) {
  for (i in 1:N) {
    # TODO
  }
}
pores.generate(im.pores.synth,20)
```

## 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*~1~,*x*~2~⟩ in *k* Intervalle der Breite Δ=(*x*~max~-*x*~min~)/*k* unterteilt, d.h., {⟨*x*~1~,*x*~1~+Δ⟩, ⟨*x*~1~+Δ,*x*~1~+2Δ⟩, .. , ⟨*x*~1~+(*k*-1)Δ,*x*~2~⟩}.

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.

[R] Histogram Plot mit *breaks*=*k* Intervalle oder alternativ Rückgabe einer Liste { lines:20; height:20 }
```R
use math,plot
im1 = matrix(runif(10000)*255,100,100,mode='uint8')
hist(im1,breaks=20)
hist(im1,plot.fill=TRUE,xlim=[20,200])
im1.h = hist(im1,breaks=20,plot=FALSE)
hist(im.pores,plot.fill=TRUE,xlim=[10,200],breaks=30)
```

[TODO]
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?

[INPUT]

## 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.

[R] Additives Rauschen zum Bild hinzufügen (hier nicht realistisch gleichverteiltes Rauschen) { lines:20; height:10 }
```R
gamma=0.1
im.noise = matrix(runif(10000)*255*gamma,100,100)
im.pores.noisy = im.pores + im.noise
plot(im.pores.synth)
```

[R] Multiplikatives Rauschen zum Bild hinzufügen (hier nicht realistisch gleichverteiltes Rauschen) { lines:20; height:10 }
```R
gamma = 0.1
im.noise = matrix(1+(runif(10000)-0.5)*gamma,100,100)
im.pores.noisy = im.pores * im.noise
plot(im.pores.synth)
```

[TODO]
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 Rauschen wird über *gamma* festgelegt. Erzeuge ein verrauschtes Porenbild mit γ=0.1.

[R] Synthetisches Rauschen { lines:20; height:20 }
```R
function image.noise (im,gamma,noisegen,additive) {
  if {additive) {
  
  } else {
  
  }
}
plot(image.noise(im.pores,0.1,runif,FALSE))
```

## 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.
- Ein relational 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.

[R] Canny Edge Detektor (*t*~1~ und *t*~2~ sind wählbare und problemabhängige Thresholdparameter) { lines:20; height:10 }
```R
use math,imager
im.pores.edges = cannyEdges(im.pores,t1=100,t2=150)
im.pores.edges.marked = as.pixset(im.pores.edges>1,as.matrix=TRUE)
print(nrow(im.pores.edges.marked))
plot(im.pores.edges)
```
[TODO]
Bewerte das Ergebnis des Canny Edge Detektors ohne und mit verschiedenen großen Rauschen. Wann werden Artefakte markiert (also Rauschen), wie muss man die Parameter *t*~1~ und *t*~2~ verändern um die Artefakte zu reduzieren. Wo ist die Grenze erreicht?.

[INPUT]

## Clustering und Einhüllende

- Die Koordinatenliste aller markierten Pixel ist als eine Punktwolke 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)

[R] Gruppierung von Punkten nach Nachbarschaftsbeziehungen  (*t*~1~ und *t*~2~ sind wählbare und problemabhängige Thresholdparameter) { lines:20; height:20 }
```R
use math,imager,dbscan
# Clutser canny edge created marked pixel corrdinates
im.pores.edges.clusters = dbscan(im.pores.edges.marked,
                                 as.coordinates=TRUE,
                                 eps=2,minPoints=2)
im.pores.edges.hulls = im.pores.edges.clusters
# Check hull points
im.temp = matrix(0,nrow(im.pores),ncol(im.pores))
for (pxl in im.pores.edges.hulls) {
  draw.pixset(im.temp,pxl,255,as.polyline=FALSE)
}
plot(im.temp)
```

>! Das Problem: Die Gruppen enthalten eventuelle nicht zur Einhüllenden einer Pore gehörende Punkte oder die Zahl der Punkte ist unnötig hoch. Für einen Ellipsenfit benötigt man die relevanten einhüllenden Punkte (i.A. sortiert)

- Wir ignorieren zunächst diese 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

>! Die Ellipsenanpassung liefert die Parameter der allgemeinen Ellipsengleichung. Mit `normalize=TRUE` werden die geometrischen Ellipsenparameter als Liste geliefert. Der Winkel ist aber um 270° verschoben gegenüber der `draw.circle` Funktion.

- Alle Ellipsenparameter werden in einer Datentabelle eingetragen

[R] Ellipsenfit für alle Gruppen { lines:20; height:20 }
```R
use math,imager,regression
im.pores.ellis=data.frame()
for (pxl in im.pores.edges.hulls) {
  e = ellipse(pxl,normalize=TRUE)
  e$angle = e$angle-270+360
  push(im.pores.ellis,e)
  draw.circle(im.temp,100,center=[e$x,e$y],angle=e$angle,
              width=e$w,height=e$h,filled=FALSE)
}
plot(im.temp)
print(im.pores.ellis)
```

>! Durch Rauschen verursachte Artefakte können als Ergebnis NaN (Not a Number) Werte liefern!

[TODO]
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?

[INPUT]

[R] Hüllenberechnung nach dem Graham Algorithmus { lines:20; height:10 }
```R
im.pores.edges.hulls = map(im.pores.edges.clusters, function (pixset) {
  chull(pixset)
})
# Check hull points
im.temp = matrix(0,nrow(im.pores),ncol(im.pores))
for (pxl in im.pores.edges.hulls) {
  draw.pixset(im.temp,pxl,255,as.polyline=FALSE)
}
plot(im.temp)
```

## 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)

[R] Statistische Analyse { lines:20; height:10 }
```R
use math,imager
im.pores.ellis=im.pores.ellis[!is.nan(im.pores.ellis$x),]
im.pores.stats = {
  xm = mean(im.pores.ellis$x),
  ym = mean(im.pores.ellis$y),
  wm = mean(im.pores.ellis$w),
  hm = mean(im.pores.ellis$h),
  am = mean(im.pores.ellis$angle)
}
print(im.pores.stats)
```

>! Wichtig: Vorher müssen  alle Datenzeilen mit NaN Werten entfernt werden (erste Anweisung)

[TODO]
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).

[INPUT]

---

[BUTTON] Hilfe { action:post; label:Absenden; style:"color:red" }
```
{
  url:'edu-9.de:28888',
  // url:'localhost:28888',
  form:['Name','Email','Pin','Frage'],
  email:{from:'$Email', to:'sbosse@uni-bremen.de', name:'$Name'},
  subject:'Hilfe Kurs ML $TITLE',
  message:'$Name: $Frage',
  attachments:[{filename:'$FILE.json',content:'$CODE'}],
  pin:1827,
}
```

[BUTTON] Einreichung (Assignment #01-46158 ) { action:post; label:Absenden; style:"color:green" }
```
{
  url:'edu-9.de:28888',
  // url:'localhost:28888',
  form:['Name','Email','Pin','Kommentar'],
  submit: { from:'$Email', to:'sbosse@uni-bremen.de' },
  assignment:'01-46158',
  name : '$Name',
  comment : '$Kommentar',
  attachments:[{filename:'$FILE.json',content:'$CODE'}],
  pin:1827,
}
```

[BUTTON] Prüfen { action:post; label:Laden; style:"color:browm" }
```
{
  url:'edu-9.de:28888',
  // url:'localhost:28888',
  form:['ID','Pin'],
  load: { id:'$ID' },
  pin:[1827,9223],
}
```

[BUTTON] Bewerten (Lehrer) { action:post; label:Absenden; style:"color:blue" }
```
{
  url:'edu-9.de:28888',
  // url:'localhost:28888',
  form:['ID','Marking','Pin','Remarks'],
  submit: { id:'$ID', from:'sbosse@uni-bremen.de', name:'$Name' },
  marks : '$Marking',
  remarks : '$Remarks',
  attachments:[{filename:'$FILE.json',content:'$CODE'}],
  pin:9223,
}
```

---

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

Inhalt.

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

Vorwort

Literatur

Synthetische Bilder

Das Histogramm

Rauschen

Kantendetektion der Poren

Clustering und Einhüllende

Ellipsenanpassung als Minimierungsproblem

Statistische Analyse

Vorwort

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

R Set-up

use math,plot,imager,geometry
dev.new(width=500)
print("Libraries opened.")
print(R.version)

▸

[]

✗

≡

Es wird in diesem Kurs der R-Dialekt r+ verwendet. R+ ist eine Reimplementierung von R in JavaScript und läuft direkt in Web Browsern und node webkit Applikationsprogrammen. Es ist keinee Softwareinstallation erforderlich.

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

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?)

Randomisiertes Bild

im1 = matrix(runif(10000)*255,100,100,mode='uint8')
plot(im1)

▸

[]

✗

≡

Geometrisches Bild 100 × 100 Pixel als Matrix (Hintergrundwert ist hier 10)

im1 = matrix(10,100,100,mode='uint8')
  p = runif(2)*100
  a = runif()*360
  w = 20+runif()*10
  h = 10+runif()*10
  I = 80+runif()*100
  draw.circle(im1,I,center=p,angle=a,width=w,height=h,filled=TRUE)
  p = runif(2)*100
  a = runif()*360
  w = 20+runif()*10
  h = 10+runif()*10
  I = 80+runif()*100
  draw.rectangle(im1,I,center=p,angle=a,width=w,height=h,filled=TRUE)
plot(im1)

▸

[]

✗

≡

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.

Synthetisches Porenbild

im.pores = matrix(10,100,100,mode='uint8')
function pores.generate (im,N) {
  for (i in 1:N) {
    # TODO
  }
}
pores.generate(im.pores,20)
plot(im.pores)

▸

[]

✗

≡

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)

Auscchnitt für den Algortihmus zur Vermeidung von überlappenden Ellipsen durch Suche von bbox Überlappungen (Intersections)

  bbox.list = list()
  # Random
  p = runif(2)*100
  a = runif()*360
  w = 20+runif()*10
  h = 10+runif()*10
  I = 80+runif()*100
  s = max(w,h)
  b = {left=p[1]-s/2,top=p[2]-s/2,width=s,height=s,center=p,angle=a,major=w,minor=h}
  # ...
  nooverlap = all(map(bbox.list,
          function (b1) { isFALSE(bbox.intersect(b,b1)) }))
  print(nooverlap)
  if (nooverlap) {
    draw.circle(im1,I,center=p,angle=a,width=w,height=h,filled=TRUE)
    push(bbox.list,b)
  }
  plot(im1)

▸

[]

✗

≡

Da randomisiert neue Kandidaten erzeugt werden kann man obiges Beispiel in eine Zählschleife einbinden, in der die Anzahl der Versuche eine nicht überlappende Ellipse zu erzeugen begrenzt wird.

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.

Synthetisches Porenbild

im.pores = matrix(10,100,100,mode='uint8')
function pores.generate.overlapfree (im,N) {
  for (i in 1:N) {
    # TODO
  }
}
pores.generate(im.pores.synth,20)

▸

[]

✗

≡

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.

Histogram Plot mit breaks=k Intervalle oder alternativ Rückgabe einer Liste

use math,plot
im1 = matrix(runif(10000)*255,100,100,mode='uint8')
hist(im1,breaks=20)
hist(im1,plot.fill=TRUE,xlim=[20,200])
im1.h = hist(im1,breaks=20,plot=FALSE)
hist(im.pores,plot.fill=TRUE,xlim=[10,200],breaks=30)

▸

[]

✗

≡

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.

Additives Rauschen zum Bild hinzufügen (hier nicht realistisch gleichverteiltes Rauschen)

gamma=0.1
im.noise = matrix(runif(10000)*255*gamma,100,100)
im.pores.noisy = im.pores + im.noise
plot(im.pores.synth)

▸

[]

✗

≡

Multiplikatives Rauschen zum Bild hinzufügen (hier nicht realistisch gleichverteiltes Rauschen)

gamma = 0.1
im.noise = matrix(1+(runif(10000)-0.5)*gamma,100,100)
im.pores.noisy = im.pores * im.noise
plot(im.pores.synth)

▸

[]

✗

≡

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 Rauschen wird über gamma festgelegt. Erzeuge ein verrauschtes Porenbild mit γ=0.1.

Synthetisches Rauschen

function image.noise (im,gamma,noisegen,additive) {
  if {additive) {

  } else {

  }
}
plot(image.noise(im.pores,0.1,runif,FALSE))

▸

[]

✗

≡

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.
Ein relational 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.

Canny Edge Detektor (t₁ und t₂ sind wählbare und problemabhängige Thresholdparameter)

use math,imager
im.pores.edges = cannyEdges(im.pores,t1=100,t2=150)
im.pores.edges.marked = as.pixset(im.pores.edges>1,as.matrix=TRUE)
print(nrow(im.pores.edges.marked))
plot(im.pores.edges)

▸

[]

✗

≡

Aufgabe. Bewerte das Ergebnis des Canny Edge Detektors ohne und mit verschiedenen großen Rauschen. 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 Punktwolke 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)

Gruppierung von Punkten nach Nachbarschaftsbeziehungen (t₁ und t₂ sind wählbare und problemabhängige Thresholdparameter)

use math,imager,dbscan
# Clutser canny edge created marked pixel corrdinates
im.pores.edges.clusters = dbscan(im.pores.edges.marked,
                                 as.coordinates=TRUE,
                                 eps=2,minPoints=2)
im.pores.edges.hulls = im.pores.edges.clusters
# Check hull points
im.temp = matrix(0,nrow(im.pores),ncol(im.pores))
for (pxl in im.pores.edges.hulls) {
  draw.pixset(im.temp,pxl,255,as.polyline=FALSE)
}
plot(im.temp)

▸

[]

✗

≡

Das Problem: Die Gruppen enthalten eventuelle nicht zur Einhüllenden einer Pore gehörende Punkte oder die Zahl der Punkte ist unnötig hoch. Für einen Ellipsenfit benötigt man die relevanten einhüllenden Punkte (i.A. sortiert)

Wir ignorieren zunächst diese 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

Die Ellipsenanpassung liefert die Parameter der allgemeinen Ellipsengleichung. Mit normalize=TRUE werden die geometrischen Ellipsenparameter als Liste geliefert. Der Winkel ist aber um 270° verschoben gegenüber der draw.circle Funktion.

Alle Ellipsenparameter werden in einer Datentabelle eingetragen

Ellipsenfit für alle Gruppen

use math,imager,regression
im.pores.ellis=data.frame()
for (pxl in im.pores.edges.hulls) {
  e = ellipse(pxl,normalize=TRUE)
  e$angle = e$angle-270+360
  push(im.pores.ellis,e)
  draw.circle(im.temp,100,center=[e$x,e$y],angle=e$angle,
              width=e$w,height=e$h,filled=FALSE)
}
plot(im.temp)
print(im.pores.ellis)

▸

[]

✗

≡

Durch Rauschen verursachte Artefakte können als Ergebnis NaN (Not a Number) Werte liefern!

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?

Hüllenberechnung nach dem Graham Algorithmus

im.pores.edges.hulls = map(im.pores.edges.clusters, function (pixset) {
  chull(pixset)
})
# Check hull points
im.temp = matrix(0,nrow(im.pores),ncol(im.pores))
for (pxl in im.pores.edges.hulls) {
  draw.pixset(im.temp,pxl,255,as.polyline=FALSE)
}
plot(im.temp)

▸

[]

✗

≡

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)

Statistische Analyse

use math,imager
im.pores.ellis=im.pores.ellis[!is.nan(im.pores.ellis$x),]
im.pores.stats = {
  xm = mean(im.pores.ellis$x),
  ym = mean(im.pores.ellis$y),
  wm = mean(im.pores.ellis$w),
  hm = mean(im.pores.ellis$h),
  am = mean(im.pores.ellis$angle)
}
print(im.pores.stats)

▸

[]

✗

≡

Wichtig: Vorher müssen alle Datenzeilen mit NaN Werten entfernt werden (erste Anweisung)

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).

Hilfe

Einreichung (Assignment #01-46158 )

Prüfen

Bewerten (Lehrer)

Created by the NoteBook Compiler Ver. 1.25.0 (c) Dr. Stefan Bosse (Fri Dec 08 2023 17:05:31 GMT+0100 (Central European Standard Time))