[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 - Merkmalsextraktion in Radiographie und CT Bildern (Teil 1)
[AUTHOR] Stefan Bosse
[VERSION] 01.2024

# Merkmalsextraktion in Radiographie und CT Bildern (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.'
  }
})
```

## Daten

In dieser Übung werden CT Bilddaten von einer radialen μCT Messung von Aluminium Gussbauteilen mit prozessbedingten Poreneinschlüssen verwendet. Die Bilder reüräsenieren rekosntruierte Schichtaufnahmen entlang der Längsachse der Platten (Abmessungen ca. 150 x 50 x 5 mm). Die Proben stammen von Fraunhofer IFAM, Bremen (Dirk Lehmhus), die Messdaten von einem Zeiss Xradia μCT Gerät, Universität Siegen, Robert Brandt.

- Die Bilder sind in der Liste `images` zusammengefasst. Sie stammen von einer Probe mit den Schichten 186, 217, 244 und 300.

[DATA] AluDC_CT_Slice_186:data/AluDC_CT_Slice_186.pgm
[DATA] AluDC_CT_Slice_217:data/AluDC_CT_Slice_217.pgm
[DATA] AluDC_CT_Slice_255:data/AluDC_CT_Slice_255.pgm
[DATA] AluDC_CT_Slice_300:data/AluDC_CT_Slice_300.pgm

[FUN]
```js
console.log('Adding imge data ..')
AluDC_CT_Slice_186_Image = ImageUtils.pgm.decode(Buffer.from(AluDC_CT_Slice_186,'binary'))
R.addValue('image.AluDC_CT_Slice_186',Object.assign(AluDC_CT_Slice_186_Image,{
  type   : 'cimg',
  value  :  AluDC_CT_Slice_186_Image.data,
  format : 'Gray',
  bits   : 16
}))
AluDC_CT_Slice_217_Image = ImageUtils.pgm.decode(Buffer.from(AluDC_CT_Slice_217,'binary'))
R.addValue('image.AluDC_CT_Slice_217',Object.assign(AluDC_CT_Slice_217_Image,{
  type   : 'cimg',
  value  :  AluDC_CT_Slice_217_Image.data,
  format : 'Gray',
  bits   : 16
}))

AluDC_CT_Slice_255_Image = ImageUtils.pgm.decode(Buffer.from(AluDC_CT_Slice_255,'binary'))
R.addValue('image.AluDC_CT_Slice_255',Object.assign(AluDC_CT_Slice_255_Image,{
  type   : 'cimg',
  value  :  AluDC_CT_Slice_255_Image.data,
  format : 'Gray',
  bits   : 16
}))

AluDC_CT_Slice_300_Image = ImageUtils.pgm.decode(Buffer.from(AluDC_CT_Slice_300,'binary'))
R.addValue('image.AluDC_CT_Slice_300',Object.assign(AluDC_CT_Slice_300_Image,{
  type   : 'cimg',
  value  :  AluDC_CT_Slice_300_Image.data,
  format : 'Gray',
  bits   : 16
}))
```

[R] R Starten { lines:2; height:5 }
```R
install.packages(["math","ml","image"])
# worker()
```

[R] R Set-up { lines:2; height:20 }
```R
use math,plot,imager,geometry,cnn
dev.new(width=500)
print("Libraries opened.")
print(R.version)
options(digits=3)
images = {
  image.AluDC_CT_Slice_186,
  image.AluDC_CT_Slice_217,
  image.AluDC_CT_Slice_255,
  image.AluDC_CT_Slice_300
}
```

## Analyse der Bilder

- Im Vordergrund stehen die geometrischen Eigenschaften der Merkmale die (semi-)automatisch durch Algorithmen extrahiert werden sollen.
- Die Merkmale sind i.A. nur latente Merkmale, z.B. die Markierung von interessanten Bereichen in Bilder (ROI), die zu einer bestimmten Klasse gehören (z.B. Poren, Defekte, Risse usw.)
- Bei Bildern untersucht man die Intensitätsverteilung, insbesonde eine ggfs. vorhandene Hintergrundhelligkeit (Offset)

[TODO]
Welchen Helligkeitswert sollten Bereich ohne Materie/Material bei CT Schichtaufnahmen besitzen, was ist hier zu beobachten?

[INPUT]

[R] Plot eines Bildes { lines:20; height:20 }
```R
use math,imager,plot
parameter.index = 1
# image.AluDC_CT_Slice_186 AluDC_CT_Slice_217
# image.AluDC_CT_Slice_255 AluDC_CT_Slice_300
image = images[[parameter.index]]
logg(summary(image))
plot(image,keep.aspect=TRUE,auto.scale=TRUE,zoom=2)
```

>@ Anstelle des klassischen Image Plots kann mit dem `crosshair=TRUE` Argument ein anderer Plotter mit einem Analysekreuz und der Angabe der aktuellen Position und relativen Intensität an dieser Stelle erzeugt werden. Achtung: Der realtive Wert 1 bezieht sich bei Autoskalierung auf den dort ermittelten Maximalwert (und 0 auf Minimum)! Besser mit `ylim=[0,65535]` normiertes Werteintervall setzen.

[R] Analyse eines Bildes { lines:20; height:20 }
```R
use math,imager,plot
parameter.index = 1
image = images[[parameter.index]]
logg(minMax(image))
hist(image,breaks=100,xlim=[0,70000])
```

## Normalisierung von Intensitätsverteilungen

Intensitätsverteilung von Bildern können auf zwei Arten normalisiert werden:

1. Globale verschiebung und Skalierung nach dem globalen minimalen und maximalen Intensitätswert ⇒ Offset und Bereichskorrektur
2. Globale Histogrammnormalisierung
3. Lokale Intensitätsangleichung, z.B. bei inhomogen beleuchteten Bildern (Mikroskop, rohe Röntgenbilder mit Intensitäsgradienten durch Quelle) ⇒ "Contrast Limited Adaptive Histogram Equalization"

>! Die Normalisierung der Intensitätswerte kann auch ein Clipping beinhalten, bezogen auf Minimal- und Maximalwert.

- Z.B. wenn aktuelle der minimale und maximale Wert in einem Bild (16 Bit, full scale Bereich [0,65535])  [3000,60000] wäre und der Eingabebereich auf [10000,50000] gesetzt wird würden alle Werte kleiner 10000 auf 0 und größer 50000 auf *I*~max~=65535 gesetzt werden! 
- Dazu kann man die `normalize(x,inputRange=[a,b],outputRange=[a,b])` Funktion aus dem *imager* Paket verwenden.
- Minimum und Maximumswerte werden mit `minMax` aus dem *math* Paket bestimmt
- Eine reine Offsetkorrekur mit `cimg-offset` führt **kein** Clippling durch und führt bei Unter- oder Überlauf zu Artefakten!

[R] Globale Normalisierung der CT Bilder mit Offsetkorrektur und Skalierung { lines:20; height:20 }
```R
images.norm=list()
# mm = {min=1000, max=60000}
for (image in images) {
  mm = minMax(image)
  cprint(mm)
  # Nur Offsetkorrektur
  image.norm = image-mm$min
  # ODER Offsetkorrektur und Skalierung
  # image.norm = normalize(image,inputRange=[mm$min,mm$max])
  push(images.norm,image.norm)
}
```

>? Ist es sinnvoll jedes Bild einzeln mit individuellen Grenzen zu normieren? Was wäre die bessere Methodik? variiere die Grenzen (inputRange von normalize), wie verändert sich das normierte Bild vom Kontrast?

[INPUT]

[R] Plot eines normierten Bildes { lines:20; height:20 }
```R
use math,imager,plot
parameter.index = 1
image = images.norm[[parameter.index]]
logg(summary(image))
plot(image,keep.aspect=TRUE,auto.scale=TRUE,zoom=2)
hist(image,breaks=100,xlim=[0,70000])
```

>! Aber: Jetzt können wir den Kontrast verbessern, doch kann es immer noch einen örtlich variierenden mittleren Intensitätswert geben (Beleuchtungsgradient).

[TODO]
Wo finden sich Veränderungen der mittleren Hintergrundhelligkeit in unseren Beispielbildern. Bemerke: Die Merkmale sind die dunklere Flecken im Materialbereich (Poren)?

[INPUT]

### Profilierung und Homogenisierung

- Um mittlere Intensitäsgradienten beurteilen und bestimmen zu können kann man gemittelte Profile in horizontaler und vertikaler Achsenrichtung verwenden (`profile(im,axis="x"|"y",normalize?,min?,max?))`. 
- Ein Profil ist ein Linienintegral. Jeder Wert eines Linienintegrals *L*(*t*) entlang einer Koordinate *t* (hier x- oder y-Achse) setzt sich aus dem mittleren Intensitäswert einer Linie orthogonal zu *L* mit dem Kreuzungspunkt *t*.
- Das Setzen der *min* und/oder *max* Argumente führt zum Ignorieren von Intensitätswerten kleiner oder größer min/max bei der Profilberechnung (Integrierung). Damit können Hintergrund und kleinere Merkmale bei der Berechnung ausgeblendet werden.

[R] Profilierung eines Bildes { lines:20; height:20 }
```R
use math,imager,plot
parameter.index = 1
image = images.norm[[parameter.index]]
logg(summary(image))
image.vprof = profile(image.norm[30:70,1:width(image)],min=30000)
plot(image.vprof)
```

- Hat man das Profil, kann man das Bild homogenisieren (geht nur einfach bei achsenausgerichteten Profilen und mittleren Intensitäsgradienten)
- Dazu wird aus dem Linienvektor (Profil) ein Bild expandiert welches dann mit einfachen arithmetischen Pixeloperationen mit dem inhomogenen Bild verrechnet wird. Dazu verwendet man einfach den `matrix` Konstruktir und gibt als werte den Profilvektor an, der dann zeilenweise dupliziert wird.

[R] Profilhomogenisierung eines Bildes { lines:20; height:20 }
```R
use math,imager,plot
parameter.index = 1
image = images.norm[[parameter.index]]
log(summary(image))
image.vprof = profile(image.norm[30:70,1:width(image)],
                      min=30000)
image.prof  = matrix(image.vprof,ncol=width(image),
                                 nrow=height(image),
                                 mode='float32')
image.profS = image.prof/60000+0.5
plot(image.profS,auto.scale=TRUE)
plot(image/image.profS,auto.scale=TRUE)
```

[R] Profilhomogenisierung aller Bilder { lines:20; height:20 }
```R
use math,imager,plot
images.nprof = list()
for (image in images.norm) {
  image.vprof = profile(image.norm[30:70,1:width(image)],
                        min=30000)
  image.prof  = matrix(image.vprof,ncol=width(image),
                                   nrow=height(image),
                                   mode='float32')
  image.profS = image.prof/60000+0.5
  image.nprof = image/image.profS
  push(images.nprof,image.nprof)
}
```

### Globale Histogramnormalisierung

Die Histogrammnormalisierung (Histogrammspreizung) ist eine  Methode, die den  Kontrast  eines  Graustufenbildes automatisch zu erhöhen. Ziel der automatischen Kontrastanpassung ist es, die Pixelwerte eines Bilds so zu verändern, dass der gesamte verfügbare Wertebereich abgedeckt wird. Dazu wird das aktuell dunkelste Pixel auf den niedrigsten und das aktuell hellste Pixel auf den höchsten Intensitätswert abgebildet und alle dazwischenliegenden Pixelwerte  linear  verteilt.

- Homogenisierung impliziert die Abbildung einer Verteilung (des gegebenen Histogramms) auf eine andere Verteilung (eine breitere und gleichmäßigere Verteilung der Intensitätswerte), so dass die Intensitätswerte über den gesamten Bereich verteilt sind.
- Um den Ausgleichseffekt zu erzielen, sollte die Neuzuordnung durch die kumulative Verteilungsfunktion (cdf) gegeben sein. Für das Histogramm ***H***(*i*) ist seine kumulative Verteilung ***H ***~cd~(*i*) (N: Anzahl der Pixel):

$
%%ascii
H_(cd)(i) = (sum_(0 le j < i) H(j)) / N
$

- Schließlich verwenden wir ein einfaches Remapping-Verfahren, um die Intensitätswerte des entzerrten Bildes zu erhalten:

$
%%ascii
I_(eq)(x,y)=H_(cd)(I(x,y))
$

[R] Globale Histogramnormalisierung der CT Bilder { lines:20; height:20 }
```R
images.norm.hist=list()
for (image.norm in images.norm) {
  image.norm.hist = hist.eq(image.norm)
  push(images.norm.hist,image.norm.hist)
}
plot(images.norm.hist[[1]],auto.scale=TRUE)
```

[QUESTION]
Was passiert mbei der Histogrammnormalisierung im Vergleich zu den vorherigen Normalisierungen? verwende sowohl *images* als auch *images.norm*. Wo gibt es Unterschiede? Wie sieht das Histogramm vorher und nachher aus?

[INPUT]

### Lokaladaptive Histogramnormalisierung

- Es wird unterschieden:

1. Intensitätsänderung aufgrund zu detektierender Merkmale (z.B. Aufhellung durch Materialporen)
2. Intensitätsänderung aufgrund ungleichmäßiger "Beleuchtung"

- Die globale Histogramnormalisierung kann keine lokalen Intensitätsgradienten und Unterschiede ausgleichen. 
- Bei Mikroskopaufnahmen hat man häufig eine räumlich leicht bis stark variierende mittlere Intensität entlang von Hyperboloiden durch ungleichmäßige Beleuchtung von Oberflächen. 
- Bei Röntgenradiografie Bildern har man das typische Beleuchtungsprofile der Röntgenröhre, i.A. eine 2-dim Gaußfunktion;
- Bei CT Bildern kommen variierende Helligkeiten vor allem vom Rekonstruktionsalgortihmus und polychromatischer Röntgenstrahlung.

CLAHE (Kontrastbegrenzter adaptiver Histogrammausgleich) ist ein Algorithmus zur Verbesserung des lokalen Kontrasts in Bildern und wird häufig in Anwendungsbereichen wie Unterwasserfotografie, Verkehrskontrolle, Astronomie und medizinischer Bildgebung eingesetzt.

- Die adaptive Histogrammentzerrung (AHE) unterscheidet sich von der gewöhnlichen Histogrammentzerrung insofern, als das adaptive Verfahren mehrere Histogramme berechnet, die jeweils einem bestimmten Bildausschnitt entsprechen, und diese zur Neuverteilung der Helligkeitswerte des Bildes verwendet.

- Es ist daher geeignet, den lokalen Kontrast zu verbessern und die Definition von Kanten in jedem Bereich eines Bildes zu verbessern.

- AHE neigt dazu, Rauschen in relativ homogenen Bereichen eines Bildes zu überverstärken.
 
- Eine Variante der adaptiven Histogrammentzerrung, die als kontrastbegrenzte adaptive Histogrammentzerrung (CLAHE) bezeichnet wird, verhindert dies, indem sie die Verstärkung begrenzt.

[R] Lokale Histogramnormalisierung der CT Bilder mit dem CLAHE Algorithmus { lines:20; height:20 }
```R
images.norm.clahe=list()
for (image.norm in images.norm) {
  image.norm.clahe = clahe(image.norm,invert=TRUE)
  push(images.norm.clahe,image.norm.clahe)
}
plot(images.norm.clahe[[1]],auto.scale=TRUE)
```

[QUESTION]
Was passiert bei der lokal adaptiven Normalisierung im Vergleich zu den vorherigen Normalisierungen? verwende sowohl *images* als auch *images.norm*. Wo gibt es Unterschiede? Wie sieht das Histogramm vorher und nachher aus?

[INPUT]

## Merkmalsextraktion

- Die Poren sind die Merkmale und dunklere Stellen in dem Bild
- Die Merkmale sollen zunächst in den Bilder markiert werden und ein binäres Merkmalsbild (0/1) erzeugen
- Einfachste Methode ist eine Thresholddiskriminierung → empfindlich auf Intensitätsschwankungen in Bildern
- Besser sind lokalapdative Threhsholding Verfahren

### Statische Threshold Diskriminierung

[R] Einfache Merkmalsmarkierung mit Intervallthresholding { lines:20; height:20 }
```R
use math,imager,plot
parameter.index = 1
parameter.tmin  = 30000
parameter.tmax  = 40000
image = images.norm[[parameter.index]]
image.thr  = matrix(0,ncol=width(image),nrow=height(image))
image.thr[(image > parameter.tmin) &
          (image < parameter.tmax)]=255
plot(image.thr,auto.scale=TRUE)
```
### Bestimmung einer optimalen Schwelle

Folgende Verfahren könnten angewendet werden, i.A. immer unter Verwendung von Histogrammen:

1. Manuell / Erfahrung
2. Global 
   + Nach der Otsu Methode
   + Nach der Triangulationsmethode (Dreicksbeziehung des Histogramms)
   + ...
3. Lokal
   + Lokal Adaptiv (siehe nachfolgenden Abschnitt)

In der Bildverarbeitung wird die nach Nobuyuki Otsu benannte Otsu-Methode zur automatischen Bildschwellenwertbildung verwendet. In der einfachsten Form gibt der Algorithmus einen einzelnen Intensitätsschwellenwert zurück, der Pixel in zwei Klassen unterteilt (binarisiert): Vordergrund und Hintergrund. Dieser Schwellenwert wird durch Minimierung der klasseninternen Intensitätsvarianz oder äquivalent durch Maximierung der klasseninternen Varianz bestimmt und entspricht einem global optimalen k-Mittelwert, der im Intensitätshistogramm durchgeführt wurde.

![#otsuvisu width=60%](images/Otsus_Method_Visualization.gif)
*Visualisierung  von  Otsu's Methode (Wikipedia)*

![#triavisu width=60%](images/trianglethres.png)
*Visualisierung  der Triangulationsmethode (Zack, 1977)*

[R] Bestimmung einer optimalen Threshold mit Otus und Triangulationsmethode { lines:20; height:20 }
```R
use math,imager,plot
parameter.index = 1
image = images.norm[[parameter.index]]
log(minMax(image))
image.thr.otsu = otsu(image,range=65536)
image.thr.tri = triangle(image,range=65536)
logg(image.thr.otsu)
logg(image.thr.tri)
```

[R] Binarisierung mit Otus und Triangulationschwelle { lines:20; height:20 }
```R
use math,imager,plot
parameter.index = 1
image = images.norm[[parameter.index]]
image.bin.otsu = matrix(0,ncol=width(image),nrow=height(image),mode='uint8')
image.bin.otsu[image>image.thr.otsu]=1
image.bin.tri = matrix(0,ncol=width(image),nrow=height(image),mode='uint8')
image.bin.tri[image>image.thr.tri]=1
plot(image.bin.otsu,auto.scale=TRUE,keep.aspect=TRUE)
plot(image.bin.tri,auto.scale=TRUE,keep.aspect=TRUE)
```

[QUESTION]
Sind diese Methoden hier brauchbar? Wenn nein, wo liegt das Problem?

[INPUT]

## Cropping

>! Häufig müssen Bilder auf relevante Inhalte beschnitten werden die die zu detektierenden Merkmale enthalten (Hintergrunddiskriminierung)

Cropping kann manuell oder automatische erfolgen. Automatisches Cropping verwendet verschiedene Bilddtransformationsn um den relevanten Bereichs eines Bildes detektieren zu können. Kantendetektion wird häufig verwendet (neben Intensitätsanalyse). Es können einfache Gradientenfilter oder der Canny Detektor verwendet werden. Die nachfolgende automatische Beschneidung des Bildes nutzt einen einfachen Gradientenfilter mit expliziter Hintergrundausblendung (*background*).

- Die `autocrop(img,background,foreground?,threshold?)` Funktion kann für das automatische Cropping mittels Gradientenfilter verwendet werden. Die Angabe der Hintergrundintensität und ggfs. eine Gradientenschwelle (*threshold*) sind erforderliche Angaben. 
- Achtung: Alle Werte sind auf 0..255 normiert!!!!
- Rückgabe ist eine Liste `{center:[x0,y0],cropbbox:{x,y,w,h},input:cimg,output:cimg,cropped:cimg}`

>! Das automatische Cropping kann fehlschlagen. Dann gibt nachfolgende *autocrop* Funktion kein Bild und keine *cropbbox* zurück. Die latenten Bilder *input* und *output* können zur internen Analyse herangezogen werden. Es wird dann in *cropbox* und *cropped* (Ergebnis) ein numerischer Wert zurückgegeben!

[R] Automatisches Cropping { lines:20; height:20 }
```R
use math,imager,plot
parameter.index = 1
image = images.norm[[parameter.index]]
# Achtung: Alle Werte sind auf 0..255 normiert!!!!
crop = autocrop(image,background=140,threshold=180)
logg(typeof(crop))
cprint(crop$center)
cprint(crop$cropbox)
plot(crop$cropped,auto.scale=TRUE,keep.aspect=TRUE)
```

### Adaptive Threshold Diskriminierung

- Adaptive (gradientebasierte) Threshold Detektoren können gut bei variierenden lokalen Intensitätsverteilungen eingesetzt werden um Strukturen zu diskriminieren (zu verstärken)

- Einer der Nachteile der Verwendung einfacher Schwellwertmethoden besteht darin, dass wir unseren Schwellwert *T* manuell angeben müssen. Darüber hinaus erfordert das Finden eines guten Wertes von *T* möglicherweise viele manuelle Experimente und Parameterabstimmungen, was in den meisten Situationen einfach nicht praktikabel ist.

- Um uns bei der automatischen Bestimmung des Wertes von *T* zu unterstützen, haben gibt es Methoden wie Otsu Thresholding. Während die Methode von Otsu uns bei der Bestimmung von *T* unterstützen kann, bleibt uns nur ein einziger Wert von *T*, um das gesamte Bild zu diskriminieren.

- Bei einfachen Aufnahmen mit kontrollierten Lichtverhältnissen ist dies in der Regel kein Problem. In Situationen, in denen die Beleuchtung im gesamten Bild ungleichmäßig ist, kann ein einziger Wert von T unsere Schwellwertleistung ernsthaft beeinträchtigen.

Die allgemeine Annahme, die allen adaptiven und lokalen Schwellwertverfahren zugrunde liegt, ist, dass kleinere Bereiche eines Bildes eher eine annähernd gleichmäßige Beleuchtung aufweisen.

- Die Nachbarschaft muss groß genug sein, um genügend Hintergrund- und Vordergrundpixel abzudecken, sonst ist der Wert von *T* mehr oder weniger irrelevant.

- Es ist üblich, entweder das arithmetische Mittel oder das Gaußsche Mittel der Pixelintensitäten in jedem Bereich zu verwenden (es gibt andere Methoden, aber das arithmetische Mittel und das Gaußsche Mittel sind bei weitem am beliebtesten).

- Im Arithmetischen Mittel trägt jedes Pixel in der Nachbarschaft gleichermaßen zur Berechnung von T bei. Und im Gaußschen Mittel tragen Pixelwerte, die weiter vom (x, y)-Koordinatenmittelpunkt der Region entfernt sind, weniger zur Gesamtberechnung von *T* bei.

Die allgemeine Formel zur Berechnung von *T* lautet also:

$
%%ascii
T = "mean"(I_L) - C
$

wobei der Mittelwert entweder der arithmetische oder der Gaußsche Mittelwert ist, *I*~L~ der lokale Teilbereich des Bildes *I*(*x*,*y*) mit einer bestimmten Blockgröße ist, und *C* eine Konstante ist, mit der wir den Schwellwert *T* fein einstellen können (Kompensation).

[R] AdaThreshold { lines:20; height:20 }
```R
use math,imager,plot
parameter.index = 1
# Size of a pixel neighborhood that is used to calculate a threshold value for the pixel
parameter.adasize  = 10
# Constant subtracted from the mean or weighted mean
parameter.adacomp  = 20
image = images.norm[[parameter.index]]
image.adathr = adaptiveThreshold(image,parameter.adasize,parameter.adacomp)
plot(image.adathr,auto.scale=TRUE)
```

[QUESTION]
Wo liegen die Probleme bei der adaptiven Threshold Diskriminierung? Wie ist die Abhängigkeit von den Parametern *adasize* (search size) und *adacomp* (compensation)?

[INPUT]

## Clustering

1. Es werden Listen von Pixelkoordinaten [x,y] benötigt. Ausgangspunkt ist ein Binärbild. Mittels `as.pixset(img>0,..)` werden die Koordinaten aller gesetzten Pixel des binären Thresholdbildes in eine Liste oder matrix eingetragen (pixel set)
2. Da die Ränder des Objekts und Hintergrund markiert sind, ist ein (adaptiver) Ausschnitt sinnvoll (Rechteckiger Bereich)
3. Mittels DBSCAN Clustering werden alle Punkte der Punktliste gemäß den Parametern gruppiert: Annahme: Jede Gruppe bildet eine Pore!

[R] Pixel Sets und DBSCAN { lines:20; height:20 }
```R
use math,imager,plot,cluster
parameter.adasize  = 90
parameter.adacomp  = 200
parameter.index = 1
image = images.norm[[parameter.index]]
image.adathr = adaptiveThreshold(image,parameter.adasize,parameter.adacomp)
image.adathr.crop = image.adathr[45:100,20:820]
logg(summary(image.adathr.crop))
image.adathr.pixels = as.pixset(image.adathr.crop>0,as.matrix=TRUE)
logg(nrow(image.adathr.pixels))
image.adathr.pixels.clusters = dbscan(image.adathr.pixels,
                                      eps=5,
                                      minPoints=2,
                                      as.coordinates=TRUE,
                                      as.matrix=TRUE)
logg(length(image.adathr.pixels.clusters))
plot(image.adathr.crop,keep.aspect=TRUE,zoom=2,auto.scale=TRUE)
```

[R] Cluster Analyse { lines:20; height:20 }
```R
print("====== Center of Mass / BBOX of Clusters =======")
image.adathr.crop.roi = matrix(0,height(image.adathr.crop),
                                 width(image.adathr.crop))
map(image.adathr.pixels.clusters, function (cluster) {
  c=com(cluster)
  bb=bbox(points=cluster)
  cprint([com=c,bb])
  draw.rectangle(image.adathr.crop.roi,bbox=bb)
})
plot(image.adathr.crop.roi,keep.aspect=TRUE,zoom=2,auto.scale=TRUE)
```

[TODO]
Variiere *eps* und *minPoints* und untersuche das Ergebnis auf best mögliche Übereinstimmung der Porengruppierung durch visuelle Inspektion. Mit welchen Parametern ergeben sich gute Ergebnisse, was sind negative Effekte?

[INPUT]

Merkmalsextraktion in Radiographie und CT Bildern (Teil 1)

Daten

Die Bilder sind in der Liste images zusammengefasst. Sie stammen von einer Probe mit den Schichten 186, 217, 244 und 300.

DATA: Variable AluDC_CT_Slice_186 Type: Binary DATA: Variable AluDC_CT_Slice_217 Type: Binary DATA: Variable AluDC_CT_Slice_255 Type: Binary DATA: Variable AluDC_CT_Slice_300 Type: Binary

R Starten

install.packages(["math","ml","image"])
# worker()

▸

[]

✗

≡

R Set-up

use math,plot,imager,geometry,cnn
dev.new(width=500)
print("Libraries opened.")
print(R.version)
options(digits=3)
images = {
  image.AluDC_CT_Slice_186,
  image.AluDC_CT_Slice_217,
  image.AluDC_CT_Slice_255,
  image.AluDC_CT_Slice_300
}

▸

[]

✗

≡

Analyse der Bilder

Im Vordergrund stehen die geometrischen Eigenschaften der Merkmale die (semi-)automatisch durch Algorithmen extrahiert werden sollen.
Die Merkmale sind i.A. nur latente Merkmale, z.B. die Markierung von interessanten Bereichen in Bilder (ROI), die zu einer bestimmten Klasse gehören (z.B. Poren, Defekte, Risse usw.)
Bei Bildern untersucht man die Intensitätsverteilung, insbesonde eine ggfs. vorhandene Hintergrundhelligkeit (Offset)

Aufgabe. Welchen Helligkeitswert sollten Bereich ohne Materie/Material bei CT Schichtaufnahmen besitzen, was ist hier zu beobachten?

Plot eines Bildes

use math,imager,plot
parameter.index = 1
# image.AluDC_CT_Slice_186 AluDC_CT_Slice_217
# image.AluDC_CT_Slice_255 AluDC_CT_Slice_300
image = images[[parameter.index]]
logg(summary(image))
plot(image,keep.aspect=TRUE,auto.scale=TRUE,zoom=2)

▸

[]

✗

≡

Anstelle des klassischen Image Plots kann mit dem crosshair=TRUE Argument ein anderer Plotter mit einem Analysekreuz und der Angabe der aktuellen Position und relativen Intensität an dieser Stelle erzeugt werden. Achtung: Der realtive Wert 1 bezieht sich bei Autoskalierung auf den dort ermittelten Maximalwert (und 0 auf Minimum)! Besser mit ylim=[0,65535] normiertes Werteintervall setzen.

Analyse eines Bildes

use math,imager,plot
parameter.index = 1
image = images[[parameter.index]]
logg(minMax(image))
hist(image,breaks=100,xlim=[0,70000])

▸

[]

✗

≡

Normalisierung von Intensitätsverteilungen

Intensitätsverteilung von Bildern können auf zwei Arten normalisiert werden:

Globale verschiebung und Skalierung nach dem globalen minimalen und maximalen Intensitätswert ⇒ Offset und Bereichskorrektur
Globale Histogrammnormalisierung
Lokale Intensitätsangleichung, z.B. bei inhomogen beleuchteten Bildern (Mikroskop, rohe Röntgenbilder mit Intensitäsgradienten durch Quelle) ⇒ "Contrast Limited Adaptive Histogram Equalization"

Die Normalisierung der Intensitätswerte kann auch ein Clipping beinhalten, bezogen auf Minimal- und Maximalwert.

Z.B. wenn aktuelle der minimale und maximale Wert in einem Bild (16 Bit, full scale Bereich [0,65535]) [3000,60000] wäre und der Eingabebereich auf [10000,50000] gesetzt wird würden alle Werte kleiner 10000 auf 0 und größer 50000 auf I_max=65535 gesetzt werden!
Dazu kann man die normalize(x,inputRange=[a,b],outputRange=[a,b]) Funktion aus dem imager Paket verwenden.
Minimum und Maximumswerte werden mit minMax aus dem math Paket bestimmt
Eine reine Offsetkorrekur mit cimg-offset führt kein Clippling durch und führt bei Unter- oder Überlauf zu Artefakten!

Globale Normalisierung der CT Bilder mit Offsetkorrektur und Skalierung

images.norm=list()
# mm = {min=1000, max=60000}
for (image in images) {
  mm = minMax(image)
  cprint(mm)
  # Nur Offsetkorrektur
  image.norm = image-mm$min
  # ODER Offsetkorrektur und Skalierung
  # image.norm = normalize(image,inputRange=[mm$min,mm$max])
  push(images.norm,image.norm)
}

▸

[]

✗

≡

Ist es sinnvoll jedes Bild einzeln mit individuellen Grenzen zu normieren? Was wäre die bessere Methodik? variiere die Grenzen (inputRange von normalize), wie verändert sich das normierte Bild vom Kontrast?

Plot eines normierten Bildes

use math,imager,plot
parameter.index = 1
image = images.norm[[parameter.index]]
logg(summary(image))
plot(image,keep.aspect=TRUE,auto.scale=TRUE,zoom=2)
hist(image,breaks=100,xlim=[0,70000])

▸

[]

✗

≡

Aber: Jetzt können wir den Kontrast verbessern, doch kann es immer noch einen örtlich variierenden mittleren Intensitätswert geben (Beleuchtungsgradient).

Aufgabe. Wo finden sich Veränderungen der mittleren Hintergrundhelligkeit in unseren Beispielbildern. Bemerke: Die Merkmale sind die dunklere Flecken im Materialbereich (Poren)?

Profilierung und Homogenisierung

Um mittlere Intensitäsgradienten beurteilen und bestimmen zu können kann man gemittelte Profile in horizontaler und vertikaler Achsenrichtung verwenden (profile(im,axis="x"|"y",normalize?,min?,max?)).
Ein Profil ist ein Linienintegral. Jeder Wert eines Linienintegrals L(t) entlang einer Koordinate t (hier x- oder y-Achse) setzt sich aus dem mittleren Intensitäswert einer Linie orthogonal zu L mit dem Kreuzungspunkt t.
Das Setzen der min und/oder max Argumente führt zum Ignorieren von Intensitätswerten kleiner oder größer min/max bei der Profilberechnung (Integrierung). Damit können Hintergrund und kleinere Merkmale bei der Berechnung ausgeblendet werden.

Profilierung eines Bildes

use math,imager,plot
parameter.index = 1
image = images.norm[[parameter.index]]
logg(summary(image))
image.vprof = profile(image.norm[30:70,1:width(image)],min=30000)
plot(image.vprof)

▸

[]

✗

≡

Hat man das Profil, kann man das Bild homogenisieren (geht nur einfach bei achsenausgerichteten Profilen und mittleren Intensitäsgradienten)
Dazu wird aus dem Linienvektor (Profil) ein Bild expandiert welches dann mit einfachen arithmetischen Pixeloperationen mit dem inhomogenen Bild verrechnet wird. Dazu verwendet man einfach den matrix Konstruktir und gibt als werte den Profilvektor an, der dann zeilenweise dupliziert wird.

Profilhomogenisierung eines Bildes

use math,imager,plot
parameter.index = 1
image = images.norm[[parameter.index]]
log(summary(image))
image.vprof = profile(image.norm[30:70,1:width(image)],
                      min=30000)
image.prof  = matrix(image.vprof,ncol=width(image),
                                 nrow=height(image),
                                 mode='float32')
image.profS = image.prof/60000+0.5
plot(image.profS,auto.scale=TRUE)
plot(image/image.profS,auto.scale=TRUE)

▸

[]

✗

≡

Profilhomogenisierung aller Bilder

use math,imager,plot
images.nprof = list()
for (image in images.norm) {
  image.vprof = profile(image.norm[30:70,1:width(image)],
                        min=30000)
  image.prof  = matrix(image.vprof,ncol=width(image),
                                   nrow=height(image),
                                   mode='float32')
  image.profS = image.prof/60000+0.5
  image.nprof = image/image.profS
  push(images.nprof,image.nprof)
}

▸

[]

✗

≡

Globale Histogramnormalisierung

Die Histogrammnormalisierung (Histogrammspreizung) ist eine Methode, die den Kontrast eines Graustufenbildes automatisch zu erhöhen. Ziel der automatischen Kontrastanpassung ist es, die Pixelwerte eines Bilds so zu verändern, dass der gesamte verfügbare Wertebereich abgedeckt wird. Dazu wird das aktuell dunkelste Pixel auf den niedrigsten und das aktuell hellste Pixel auf den höchsten Intensitätswert abgebildet und alle dazwischenliegenden Pixelwerte linear verteilt.

Homogenisierung impliziert die Abbildung einer Verteilung (des gegebenen Histogramms) auf eine andere Verteilung (eine breitere und gleichmäßigere Verteilung der Intensitätswerte), so dass die Intensitätswerte über den gesamten Bereich verteilt sind.
Um den Ausgleichseffekt zu erzielen, sollte die Neuzuordnung durch die kumulative Verteilungsfunktion (cdf) gegeben sein. Für das Histogramm H(i) ist seine kumulative Verteilung *H *_cd(i) (N: Anzahl der Pixel):

\[ {H}_{{{c}{d}}}{\left({i}\right)}=\frac{{\sum_{{{0}\le{j}<{i}}}{H}{\left({j}\right)}}}{{N}} \]

Schließlich verwenden wir ein einfaches Remapping-Verfahren, um die Intensitätswerte des entzerrten Bildes zu erhalten:

\[ {I}_{{{e}{q}}}{\left({x},{y}\right)}={H}_{{{c}{d}}}{\left({I}{\left({x},{y}\right)}\right)} \]

Globale Histogramnormalisierung der CT Bilder

images.norm.hist=list()
for (image.norm in images.norm) {
  image.norm.hist = hist.eq(image.norm)
  push(images.norm.hist,image.norm.hist)
}
plot(images.norm.hist[[1]],auto.scale=TRUE)

▸

[]

✗

≡

Frage. Was passiert mbei der Histogrammnormalisierung im Vergleich zu den vorherigen Normalisierungen? verwende sowohl images als auch images.norm. Wo gibt es Unterschiede? Wie sieht das Histogramm vorher und nachher aus?

Lokaladaptive Histogramnormalisierung

Es wird unterschieden:

Intensitätsänderung aufgrund zu detektierender Merkmale (z.B. Aufhellung durch Materialporen)
Intensitätsänderung aufgrund ungleichmäßiger "Beleuchtung"

Die globale Histogramnormalisierung kann keine lokalen Intensitätsgradienten und Unterschiede ausgleichen.
Bei Mikroskopaufnahmen hat man häufig eine räumlich leicht bis stark variierende mittlere Intensität entlang von Hyperboloiden durch ungleichmäßige Beleuchtung von Oberflächen.
Bei Röntgenradiografie Bildern har man das typische Beleuchtungsprofile der Röntgenröhre, i.A. eine 2-dim Gaußfunktion;
Bei CT Bildern kommen variierende Helligkeiten vor allem vom Rekonstruktionsalgortihmus und polychromatischer Röntgenstrahlung.

Die adaptive Histogrammentzerrung (AHE) unterscheidet sich von der gewöhnlichen Histogrammentzerrung insofern, als das adaptive Verfahren mehrere Histogramme berechnet, die jeweils einem bestimmten Bildausschnitt entsprechen, und diese zur Neuverteilung der Helligkeitswerte des Bildes verwendet.
Es ist daher geeignet, den lokalen Kontrast zu verbessern und die Definition von Kanten in jedem Bereich eines Bildes zu verbessern.
AHE neigt dazu, Rauschen in relativ homogenen Bereichen eines Bildes zu überverstärken.
Eine Variante der adaptiven Histogrammentzerrung, die als kontrastbegrenzte adaptive Histogrammentzerrung (CLAHE) bezeichnet wird, verhindert dies, indem sie die Verstärkung begrenzt.

Lokale Histogramnormalisierung der CT Bilder mit dem CLAHE Algorithmus

images.norm.clahe=list()
for (image.norm in images.norm) {
  image.norm.clahe = clahe(image.norm,invert=TRUE)
  push(images.norm.clahe,image.norm.clahe)
}
plot(images.norm.clahe[[1]],auto.scale=TRUE)

▸

[]

✗

≡

Frage. Was passiert bei der lokal adaptiven Normalisierung im Vergleich zu den vorherigen Normalisierungen? verwende sowohl images als auch images.norm. Wo gibt es Unterschiede? Wie sieht das Histogramm vorher und nachher aus?

Merkmalsextraktion

Die Poren sind die Merkmale und dunklere Stellen in dem Bild
Die Merkmale sollen zunächst in den Bilder markiert werden und ein binäres Merkmalsbild (0/1) erzeugen
Einfachste Methode ist eine Thresholddiskriminierung → empfindlich auf Intensitätsschwankungen in Bildern
Besser sind lokalapdative Threhsholding Verfahren

Statische Threshold Diskriminierung

Einfache Merkmalsmarkierung mit Intervallthresholding

use math,imager,plot
parameter.index = 1
parameter.tmin  = 30000
parameter.tmax  = 40000
image = images.norm[[parameter.index]]
image.thr  = matrix(0,ncol=width(image),nrow=height(image))
image.thr[(image > parameter.tmin) &
          (image < parameter.tmax)]=255
plot(image.thr,auto.scale=TRUE)

▸

[]

✗

≡

Bestimmung einer optimalen Schwelle

Folgende Verfahren könnten angewendet werden, i.A. immer unter Verwendung von Histogrammen:

Manuell / Erfahrung
Global
- Nach der Otsu Methode
- Nach der Triangulationsmethode (Dreicksbeziehung des Histogramms)
- ...
Lokal
- Lokal Adaptiv (siehe nachfolgenden Abschnitt)

Visualisierung von Otsu's Methode (Wikipedia)

#triavisu
Visualisierung der Triangulationsmethode (Zack, 1977)

Bestimmung einer optimalen Threshold mit Otus und Triangulationsmethode

use math,imager,plot
parameter.index = 1
image = images.norm[[parameter.index]]
log(minMax(image))
image.thr.otsu = otsu(image,range=65536)
image.thr.tri = triangle(image,range=65536)
logg(image.thr.otsu)
logg(image.thr.tri)

▸

[]

✗

≡

Binarisierung mit Otus und Triangulationschwelle

use math,imager,plot
parameter.index = 1
image = images.norm[[parameter.index]]
image.bin.otsu = matrix(0,ncol=width(image),nrow=height(image),mode='uint8')
image.bin.otsu[image>image.thr.otsu]=1
image.bin.tri = matrix(0,ncol=width(image),nrow=height(image),mode='uint8')
image.bin.tri[image>image.thr.tri]=1
plot(image.bin.otsu,auto.scale=TRUE,keep.aspect=TRUE)
plot(image.bin.tri,auto.scale=TRUE,keep.aspect=TRUE)

▸

[]

✗

≡

Frage. Sind diese Methoden hier brauchbar? Wenn nein, wo liegt das Problem?

Cropping

Häufig müssen Bilder auf relevante Inhalte beschnitten werden die die zu detektierenden Merkmale enthalten (Hintergrunddiskriminierung)

Die autocrop(img,background,foreground?,threshold?) Funktion kann für das automatische Cropping mittels Gradientenfilter verwendet werden. Die Angabe der Hintergrundintensität und ggfs. eine Gradientenschwelle (threshold) sind erforderliche Angaben.
Achtung: Alle Werte sind auf 0..255 normiert!!!!
Rückgabe ist eine Liste {center:[x0,y0],cropbbox:{x,y,w,h},input:cimg,output:cimg,cropped:cimg}

Das automatische Cropping kann fehlschlagen. Dann gibt nachfolgende autocrop Funktion kein Bild und keine cropbbox zurück. Die latenten Bilder input und output können zur internen Analyse herangezogen werden. Es wird dann in cropbox und cropped (Ergebnis) ein numerischer Wert zurückgegeben!

Automatisches Cropping

use math,imager,plot
parameter.index = 1
image = images.norm[[parameter.index]]
# Achtung: Alle Werte sind auf 0..255 normiert!!!!
crop = autocrop(image,background=140,threshold=180)
logg(typeof(crop))
cprint(crop$center)
cprint(crop$cropbox)
plot(crop$cropped,auto.scale=TRUE,keep.aspect=TRUE)

▸

[]

✗

≡

Adaptive Threshold Diskriminierung

Adaptive (gradientebasierte) Threshold Detektoren können gut bei variierenden lokalen Intensitätsverteilungen eingesetzt werden um Strukturen zu diskriminieren (zu verstärken)
Einer der Nachteile der Verwendung einfacher Schwellwertmethoden besteht darin, dass wir unseren Schwellwert T manuell angeben müssen. Darüber hinaus erfordert das Finden eines guten Wertes von T möglicherweise viele manuelle Experimente und Parameterabstimmungen, was in den meisten Situationen einfach nicht praktikabel ist.
Um uns bei der automatischen Bestimmung des Wertes von T zu unterstützen, haben gibt es Methoden wie Otsu Thresholding. Während die Methode von Otsu uns bei der Bestimmung von T unterstützen kann, bleibt uns nur ein einziger Wert von T, um das gesamte Bild zu diskriminieren.
Bei einfachen Aufnahmen mit kontrollierten Lichtverhältnissen ist dies in der Regel kein Problem. In Situationen, in denen die Beleuchtung im gesamten Bild ungleichmäßig ist, kann ein einziger Wert von T unsere Schwellwertleistung ernsthaft beeinträchtigen.

Die allgemeine Annahme, die allen adaptiven und lokalen Schwellwertverfahren zugrunde liegt, ist, dass kleinere Bereiche eines Bildes eher eine annähernd gleichmäßige Beleuchtung aufweisen.

Die Nachbarschaft muss groß genug sein, um genügend Hintergrund- und Vordergrundpixel abzudecken, sonst ist der Wert von T mehr oder weniger irrelevant.
Es ist üblich, entweder das arithmetische Mittel oder das Gaußsche Mittel der Pixelintensitäten in jedem Bereich zu verwenden (es gibt andere Methoden, aber das arithmetische Mittel und das Gaußsche Mittel sind bei weitem am beliebtesten).
Im Arithmetischen Mittel trägt jedes Pixel in der Nachbarschaft gleichermaßen zur Berechnung von T bei. Und im Gaußschen Mittel tragen Pixelwerte, die weiter vom (x, y)-Koordinatenmittelpunkt der Region entfernt sind, weniger zur Gesamtberechnung von T bei.

Die allgemeine Formel zur Berechnung von T lautet also:

\[ {T}=\text{mean}{\left({I}_{{L}}\right)}-{C} \]

wobei der Mittelwert entweder der arithmetische oder der Gaußsche Mittelwert ist, I_L der lokale Teilbereich des Bildes I(x,y) mit einer bestimmten Blockgröße ist, und C eine Konstante ist, mit der wir den Schwellwert T fein einstellen können (Kompensation).

AdaThreshold

use math,imager,plot
parameter.index = 1
# Size of a pixel neighborhood that is used to calculate a threshold value for the pixel
parameter.adasize  = 10
# Constant subtracted from the mean or weighted mean
parameter.adacomp  = 20
image = images.norm[[parameter.index]]
image.adathr = adaptiveThreshold(image,parameter.adasize,parameter.adacomp)
plot(image.adathr,auto.scale=TRUE)

▸

[]

✗

≡

Frage. Wo liegen die Probleme bei der adaptiven Threshold Diskriminierung? Wie ist die Abhängigkeit von den Parametern adasize (search size) und adacomp (compensation)?

Clustering

Es werden Listen von Pixelkoordinaten [x,y] benötigt. Ausgangspunkt ist ein Binärbild. Mittels as.pixset(img>0,..) werden die Koordinaten aller gesetzten Pixel des binären Thresholdbildes in eine Liste oder matrix eingetragen (pixel set)
Da die Ränder des Objekts und Hintergrund markiert sind, ist ein (adaptiver) Ausschnitt sinnvoll (Rechteckiger Bereich)
Mittels DBSCAN Clustering werden alle Punkte der Punktliste gemäß den Parametern gruppiert: Annahme: Jede Gruppe bildet eine Pore!

Pixel Sets und DBSCAN

use math,imager,plot,cluster
parameter.adasize  = 90
parameter.adacomp  = 200
parameter.index = 1
image = images.norm[[parameter.index]]
image.adathr = adaptiveThreshold(image,parameter.adasize,parameter.adacomp)
image.adathr.crop = image.adathr[45:100,20:820]
logg(summary(image.adathr.crop))
image.adathr.pixels = as.pixset(image.adathr.crop>0,as.matrix=TRUE)
logg(nrow(image.adathr.pixels))
image.adathr.pixels.clusters = dbscan(image.adathr.pixels,
                                      eps=5,
                                      minPoints=2,
                                      as.coordinates=TRUE,
                                      as.matrix=TRUE)
logg(length(image.adathr.pixels.clusters))
plot(image.adathr.crop,keep.aspect=TRUE,zoom=2,auto.scale=TRUE)

▸

[]

✗

≡

Cluster Analyse

print("====== Center of Mass / BBOX of Clusters =======")
image.adathr.crop.roi = matrix(0,height(image.adathr.crop),
                                 width(image.adathr.crop))
map(image.adathr.pixels.clusters, function (cluster) {
  c=com(cluster)
  bb=bbox(points=cluster)
  cprint([com=c,bb])
  draw.rectangle(image.adathr.crop.roi,bbox=bb)
})
plot(image.adathr.crop.roi,keep.aspect=TRUE,zoom=2,auto.scale=TRUE)

▸

[]

✗

≡

Aufgabe. Variiere eps und minPoints und untersuche das Ergebnis auf best mögliche Übereinstimmung der Porengruppierung durch visuelle Inspektion. Mit welchen Parametern ergeben sich gute Ergebnisse, was sind negative Effekte?

Created by the NoteBook Compiler Ver. 1.27.1 (c) Dr. Stefan Bosse (Tue Feb 20 2024 23:50:53 GMT+0100 (Central European Standard Time))