[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 mit ML (Teil 2)
[AUTHOR] Stefan Bosse
[VERSION] 01.2024

# Merkmalsextraktion in Radiographie und CT Bildern mit ML (Teil 2)

[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.'
  }
})
```

Ziel ist die Markierung von versteckten Poren in Druckgussbauteilen (Platten) mittels X-ray Radiographie Bildern. Da es in gemessenen Bildern von realen Proben keine zuverlässige Kenntnis und Charakterisierung der Poren a-priori gibt werden synthetische Bilder mit Ground Truth Markierungen (von Poren und deren ROI Boundingboxen) verwendet um ein datengetriebene Prädiktormodell zu trainieren.

Das Prädiktormodell ist ein einfacher Pixelklassifikator. Ein Pixelklassifikator soll entscheiden ob ein Pixel eines Bildes zu einer bestimmten Klasse gehört, hier Pore oder nicht. Dazu wird ein bereich um das zentrale Pixel betrachtet (Maskierung, z.B. 20 × 20 Pixel Maske). Es wird ein einfaches Convolutional Neural Network als Musterdetektor verwendet.

Eingabe des CNN ist eine Ausschnitt des gesamten Bildes, Ausgabe ist hier ein skalarer Wert 0..1 der eine Vorhersage der Zugehörigkeit des zentralen Pixels des Ausschnitts zu einer Pore gibt.

## Daten

In dieser Übung werden Radiographie Bilddaten von einer von Aluminium Gussbauteilen mit prozessbedingten Poreneinschlüssen verwendet. Die Bilder repräsentieren die Dichteverteilung orthogonal zur Oberfläche 8also in die Tiefe gehend, Zentralprojektion). Die Proben stammen von Fraunhofer IFAM, Bremen (Dirk Lehmhus), die Messdaten vom Fraunhofer IFAM, Bremen (Dirk Lehmhus) mit einer industriellen Röntgenanlage (Pixelgröße des Detektors ist 200μm, typische Vergrößerung M=2).

1. Synthetische X-ray Radiographie Daten (Simulation) mit Ground Truth Labeling (ROIs); Bilder 1 und 2;
2. Reale X-ray Radiographie Messdaten, hier Bild 3.

[DATA] PlateWRanElliPores1:data/xray-PlateWRanElliPores1.pgm
[DATA] PlateWRanElliPores2:data/xray-PlateWRanElliPores2.pgm
[DATA] BT226004:data/BT-226-004.pgm
[DATA] ROIS1:data/PlateWRanElliPores1.json
[DATA] ROIS2:data/PlateWRanElliPores2.json

[FUN]
```js
console.log('Adding imge data ..')
PlateWRanElliPores1_Image = ImageUtils.pgm.decode(Buffer.from(PlateWRanElliPores1,'binary'))
R.addValue('image.PlateWRanElliPores1',Object.assign(PlateWRanElliPores1_Image,{
  type   : 'cimg',
  value  :  PlateWRanElliPores1_Image.data,
  format : 'GREY',
  bits   : 16
}))
PlateWRanElliPores2_Image = ImageUtils.pgm.decode(Buffer.from(PlateWRanElliPores2,'binary'))
R.addValue('image.PlateWRanElliPores2',Object.assign(PlateWRanElliPores2_Image,{
  type   : 'cimg',
  value  :  PlateWRanElliPores2_Image.data,
  format : 'GREY',
  bits   : 16
}))
BT226004_Image = ImageUtils.pgm.decode(Buffer.from(BT226004,'binary'))
R.addValue('image.BT226004',Object.assign(BT226004_Image,{
  type   : 'cimg',
  value  :  BT226004_Image.data,
  format : 'GREY',
  bits   : 8
}))
R.addValue('rois1',R.fromJS(ROIS1))
R.addValue('rois2',R.fromJS(ROIS2))
```

[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.PlateWRanElliPores1,
  image.PlateWRanElliPores2,
  image.BT226004
}
rois = {
  rois1,
  rois2
}
```

[R] Daten { lines:10; height:10 }
```R
logg(typeof(images))
logg(typeof(rois))
logg(summary(images[[1]]))
logg(summary(images[[2]]))
logg(summary(images[[3]]))
```

## Aufbereitung der Daten

>! Zunächst visuelle Inspektion der Bilder:

[R] Plot eines Bildes { lines:20; height:5 }
```R
index = 1
image = images[[index]]
# image = images.cropped[[index]]
plot(image,auto.scale=TRUE,keep.aspect=TRUE)
```

>! Boundingboxen als Umrandung des ROI (Platte im X-ray Bild) können manuell ermittelt werden (Verwendung von `plot(...,crosshair=TRUE)`)

[R] Beschneidung der Bilddaten (Cropping) { lines:10; height:20 }
```R
# Image cropping 
# Crop bboxes for images #1/2 and #3
images.cropped=list()
parameter.cropboxes={
  bbox(left=345,top=100,right=680-1,bottom=930-1),
  bbox(left=345,top=100,right=680-1,bottom=930-1),
  bbox(left=340,top=100,right=620-1,bottom=920-1)
}
for(index in 1:3) {
  images.cropped[[index]]=crop.bbox(images[[index]],
                                    bbox=parameter.cropboxes[[index]])
  logg(summary(images.cropped[[index]]))
}
print("=> images.cropped")
```

Zum trainieren des CNN werden annotierte (gelabelte) Beispieldaten benötigt. Annotierte Bilder werden aus den ROI Listen der synthetischen (simulierten) Bilder gewonnen.

Die ROI Listen stammen von einem CAD Modell einer synthetischen Platte mit randomisiert verteilten Poren. Mittels X-ray Bildsimulation wurden die radiographischen X-ray Bilder erzeugt. Die Bilder besitzen als Koordinatensystem Pixeleinheiten, die ursprünglichen ROIs aus dem CAD Modell sind in physikalischen Einheiten (mm) angegeben und müssen in Pixelkoordinaten  umgerechnet werden. Dies geschieht nachfolgend.

[R] Transformation der ROI Label Listen { lines:20; height:30 }
```R
rois.t = list()
# roi coordiante systeme origin: (0,0) [mm]
# width: 40
# scale to pixel coordinates (parameter.cropbox)
parameter.width=40
parameter.height=99
function transformPoints(roi) {
  roi=copy(roi[[1]])
  roi$x=floor((roi$x+parameter.width/2)*parameter.scale)
  roi$y=floor((roi$y+parameter.height/2)*parameter.scale)
  print(paste(roi$x,roi$y))
  roi$path=map(roi$path,function (p) {
    { 
      x=floor((p$x+parameter.width/2)*parameter.scale),
      y=floor((p$y+parameter.height/2)*parameter.scale)
    }
  })
  roi
}
for(i in 1:2) {
  print("===== Image #"+i+" ======")
  parameter.cropbox = parameter.cropboxes[[i]]
  parameter.scale=parameter.cropbox$width/(parameter.width)
  roi = map(rois[[i]],transformPoints)
  rois.t[[i]]=roi
}
print("=> rois.t")
```

Die Umrechnung ist etwas "tricky", daher besser die Übereinstimmung der ROI Boundingsboxen mit den Poren im Radiographiebild nachfolgend prüfen.

[R] Ausgabe ROI { lines:10; height:10 }
```R
index.image = 1
index.roi = 1
logg(rois.t[[index.image]][[index.roi]])
```

>! Die ROIS (x,y,path) stammen aus der CAD Modellierung der synthetischen Platten. Es fehlen noch rechteckige Bereiche ohne Poren, d.h., Hintergrund. Diese können manuell ausgewählt werden (s.o.)

[R] Background ROIs { lines:10; height:10 }
```R
rois0={
  {
    bbox(20,20,50,150),
    bbox(270,200,50,150)
  },
  {
    bbox(100,520,50,150),
    bbox(270,200,50,150)
  }
}
print("=> rois0")
```

[R] Plot eines annotierten synthetischen Bildes (1/2) { lines:30; height:20 }
```R
use geometry
index=1
image = as.cimg(images.cropped[[index]],mode='uint16')
roil   = rois.t[[index]]
roil0  = rois0[[index]]
for(roi in roil) {
  draw.pixset(image,roi$path,value=65000)
}
for(roi in roil0) {
  draw.rectangle(image,bbox=roi,value=65000)
}
plot(image,auto.scale=TRUE,zoom=1)
```

Jetzt müssen die Trainings- und Testdaten aus den vorherigen ROI Listen extrahiert werden.

- Es sollte eine Balance zwischen Poren- und Hintergrundbeispielen erzielt werden
- Das Ergebnis der Prädiktion könnte eine erhöhte Falsch-Positiv oder Falsche-Negative Rate entstehen. Die Anteile von positiven und negativen Beispielen kann variiert werden, je nachdem welche Rate man reduzieren will (bzw. welchen Fall man optimieren möchte). Ist es wichtiger möglichst alle Pixel von Poren zu markieren und ein erhöhtes Rauschen in Kauf zu nehmen, oder ist es besser das Rauschen zu reduzieren?

- Die Bildsegmente werden aus den ROIs gewonnen. Vereinfacht können sie durch rechteckige Boundingsboxen approximiert werden. Nimmt man zentrale Pixel aus dem Zentrum wird es i.A. immer eine Überdeckung mit den tatsächlichen Polygonen 8der Porenumrandung) geben. Ansonsten muss man Punkte aus der Boundingbox auf Einschluss in die Porenhülle (also das Polygon) explizit testen (und Punkte außerhalb der Porenhülle verwerfen).
- Hier wird der einfache Ansatz gewählt.

[R] Erzeugen der Bildsegmente für die Merkmalsklassen P und N { lines:10; height:100 }
```R
parameter.segment.size = 20
h2 = parameter.segment.size/2
w2 = parameter.segment.size/2
segments.p = list() # Subimages containing feature
segments.n = list() # background
async for(index in 1:2) {
  image  = images.cropped[[index]]
  roil   = rois.t[[index]]
  roil0  = rois0[[index]]
  for (roi in roil) {
    print(paste(index,roi$x,roi$y))
    # get segments around around roi center position
    for (dx in -2:2) {
      for(dy in -2:2) {
        # TODO: use hull search
        # convert uin16 to float32 and normalize to [0,1]
        segment = as.matrix(image[(roi$y-h2+dy):(roi$y+h2-1+dy),
                                  (roi$x-w2+dx):(roi$x+w2-1+dx)],
                                  mode='float32')/50000
        push(segments.p,{label='P',image=segment})
      }
    }
  }
  # await(.Last.value)
  for (roi in roil0) {
    print(paste(0,index,roi$left,roi$top))
    # get segments around around roi center position
    # TODO: use hull search
    # convert uin16 to float32 and normalize to [0,1]
    #segment = as.matrix(image[(roi$y-h2):(roi$y+h2-1),(roi$x-w2):(roi$x+w2-1)],mode='float32')/50000
    #push(segments,{label='P',image=segment})
    xp = floor(runif(30,roi$left,roi$left+roi$width))
    yp = floor(runif(30,roi$left,roi$left+roi$width))
    for(x in xp) {
      for(y in yp) {
        segment0 = as.matrix(image[(y-h2):(y+h2-1),
                                   (x-w2):(x+w2-1)],
                                   mode='float32')/50000
        push(segments.n,{label='B',image=segment0})
      }
    }
  }
}
await(.Last.value,function () {
  logg(summary(segment))
  logg(length(segments.p))
  logg(typeof(segments.p))
  logg(length(segments.n))
  logg(typeof(segments.n))
})
```

[R] Plot eines Bildsegments { lines:30; height:10 }
```R
index = 600
# Pore
segment = segments.p[[index]]
plot(segment$image,auto.scale=TRUE,zoom=0.5)
# Hintergrund
segment = segments.n[[index]]
plot(segment$image,auto.scale=TRUE,zoom=0.5)
```

[QUESTION]
Würde man mittels visueller Inspektion der Segmentbilder Poren von Hintergrund unterscheiden können?

[INPUT]

## CNN Pixel Klassifizierer

Ein semantischer Pixelklassifizierer ist lokaler Strukturdetektor mit dem Ziel ein Eingabebild auf ein Merkmalsbild abzubilden. Dazu wird der Pixeklklassifizierer auf alle Pixel des Eingabebildes angewendet. Um lokale Strukturen erkennen zu können wird ein Bildausschnitt rund um ein zentrales Pixel benötigt, z.B. ein Segment mitd er Größe von 20 × 20 Pixel.

![#pixelcls width=60%](images/PXCL-CNN.png)

[R] Erzeugen des CNN Pixelklassifizier Modells { lines:10; height:30 }
```R
use cnn
parameter.segmentSize=20
parameter.features=["P","B"]
parameter.filterA=8
parameter.filterB=8
m=cnn(
  input.size=[parameter.segmentSize,parameter.segmentSize,1],
  layers={
    {type='conv',kernel.size=[5,5],filter=parameter.filterA,stride=1,padding=2,activation='relu'},
    {type='pool',kernel.size=[2,2],stride=2},
    # {type='conv',kernel.size=[5,5],filter=parameter.filterB,stride=1,padding=2,activation='relu'},
    # {type='pool',kernel.size=[2,2],stride=2},
    {type='softmax',num.classes=length(parameter.features)}
  },
  labels=parameter.features #optional
  # trainer={method='sgd',learning.rate=0.01,momentum=0.2,batch.size=10}
  # trainer={method='adadelta',l1.decay=0,l2.decay=0.001,batch.size=1}
  # trainer={method='adadelta',l1.decay=0,l2.decay=0.001,batch.size=10}
)
logg(summary(m))
```

Der Pixelklassifizierer besteht aus:

1. Einer ersten Filterschicht (Convolution) mit 4-8 Filtern (Kernels) der Größe 5 × 5 pixel. Jedes Filter ist mit 25 Parametern unabhängig konfigurierbar und trainierbar. Diese Schicht wird direkt auf die kleinen maskierten Eingabilder als Teilausschnitt des großen Bildes angewendet.
2. Einer Poolingschicht. Im Grunde wieder eine Filterschicht, aber nicht dynamisch parametrisierbar sonder mit einer statischen Funktion, hier die Maximumsfunktion. Diese wählt aus der Filtermaske den größten Wert als Wert für die Zielmatrix aus.
3. Optional weitere Convolution-Pooling Schichtenpaare.
4. Eine neuronale vollständig verbundene Schicht, d.h. jedes Neuron dieser Schicht (mindestens eines oder eines pro Zielklasse!) wird mit **allen** Elementen der Matrizen der letzten Ausgabeschicht verbunden. Die Verbindungen sind gewichtet und parametrisierbar.

## Trainings- und Testdaten

- Zunächst müssen die annotierten Bilddaten aufbereitet werden.
- Es wird eine Datentabelle mit einer x- und einer y-Spalte erzeugt.
- Die x-Spalte enthält die Bildsegmente, die y-Spalte die Labels

[R] Erzeugen der gesamten Datentabelle { lines:10; height:30 }
```R
segments = c(segments.p,segments.n)
data = data.frame(
  x=pluck(segments,'image'),
  y=pluck(segments,'label')
)
logg(typeof(data))
```

- Die Klassen *B* (Background) und *P* (Pore) sollten gleichermaßen in den Trainings- und Testdaten enthalten sein. Die Zusammensetzung z.B. einer Spalte einer Datentabelle kann mit der `table` Funktion einfach bestimmt werden (s.u.).

>! Die Trainingsdaten sollten in der Reihenfolge randomisiert werden (Shuffling). Dies kann einfach mit der `sample` Funktion erreicht werden.

[R] Erzeugen der Trainings- und Testdatentabelle (Sampling) { lines:20; height:5 }
```R
samples.training = 1000
samples.test     = 1000
training = sample(data,samples.training)
test     = sample(data,samples.test)
logg(typeof(training))
logg(table(training$y))
logg(table(training$y))
```

## Training

[R] Training des Modells { lines:10; height:30 }
```R
parameter.epochs=1
xy=training
print(table(as.vector(xy$y)))
train(m,
  xy,
  formula=y~x,
  epochs=parameter.epochs,
  callback=function (result) {
    cprint(result)
})
await(.Last.value,function (result) {
  print('==== Finished ========')
  cprint(result)
})
```

[TODO]
Führe das Training durch. Wähle zunächst 10 Epochen. Dann notiere den letzten Klassifikationsfehler (*error*) und den "Verlust" (*loss*). Führe nach 10 Epochen nachfolgenden Test mit den Trainings- und Testdaten durch. Gibt es Unterschiede, und warum könnte der fehler des nachfolgenden Tests (mit den Trainingsdaten) vom letzten ausgegeben Fehler des Trainings abweichen? Führe weitere 10 Trainingsepochen durch (usw.).

[INPUT]

## Test

- Nach einigen Trainingsdurchläufen sollte das Modell mit den Trainings- und Testdaten exemplarisch getestet werden.
- Alternativ einen weiteren randomisierten Ausschnitt aus der großen Datenmenge (Trainings- und Testdaten sind hier nur Untermengen)

[R] Test des Modells { lines:10; height:30 }
```R
xy0=sample(data,100)
# oder
# xy0=training
# oder
# xy0=test
print(table(as.vector(xy$y)))
result = predict(m,x[100])
print(result)
result = predict(m,xy0$x)
result = cbind(result,y0=xy0$y)
print(paste('Error',(1-sum(result$value==result$y0)/100)*100,'%'))
```

[QUESTION]
Was gibt die Prädiktion (`predict`) zurück? Welche Bedeutung hat der Score Wert?

[INPUT]

## Inferenz

- Um das Modell auf ein gesamtes Bild anzuwenden muss jeweils ein maskierendes Fenster (Segment) über das Bild geschoben werden und für jedes Segment das Modell berechnet werden. 
- Das Ergebnis  der Berechnung ist dann ein Pixelwert im Zielbild (Featurebild)
- Skalierung (Normierung) der Intensitätswerte nicht vergessen! Gleiche Skalierung wählen wie bei Trainingsdaten.

[R] Anwendung des Modells auf ein gesamtes Bild { lines:10; height:30 }
```R
index = 1
image = images.cropped[[index]]
w = width(image)
h = height(image)
# margins
margin = parameter.segment.size
function predictSegment(image,x,y) {
  sx = parameter.segment.size
  sy = parameter.segment.size
  segment = image[(y-sy/2):(y+sy/2-1),
                  (x-sx/2):(x+sx/2-1)]
  x = as.matrix(segment,mode='float32')/50000
  result = predict(m,x)
  if (result$value=='P') 
    result$score
  else
    0
}
logg(predictSegment(image,50,50))
image.feature = matrix(0,h-2*margin,w-2*margin,mode='float32')
async for (y in margin:(h-margin)) {
  print(y)
  for(x in margin:(w-margin)) {
    i=x-margin+1
    j=y-margin+1
    image.feature[j,i]=predictSegment(image,x,y)
  }
}
await(.Last.value,function (result) {
  plot(image.feature,ylim=[0,1],keep.aspect=TRUE)
})
```

[TODO]
Wende die Inferenz auf ein synthetisches und das eine gemessene Bild an. Wie sind die Ergebnisse?

[INPUT]

[R] Raum für weitere Experimente { lines:10; height:30 }
```R
# TODO
```

Merkmalsextraktion in Radiographie und CT Bildern mit ML (Teil 2)

Eingabe des CNN ist eine Ausschnitt des gesamten Bildes, Ausgabe ist hier ein skalarer Wert 0..1 der eine Vorhersage der Zugehörigkeit des zentralen Pixels des Ausschnitts zu einer Pore gibt.

Daten

Synthetische X-ray Radiographie Daten (Simulation) mit Ground Truth Labeling (ROIs); Bilder 1 und 2;
Reale X-ray Radiographie Messdaten, hier Bild 3.

DATA: Variable PlateWRanElliPores1 Type: Binary DATA: Variable PlateWRanElliPores2 Type: Binary DATA: Variable BT226004 Type: Binary DATA: Variable ROIS1 Type: { x : number, y : number, path : { x : number, y : number } [17], name : string } [1][60] DATA: Variable ROIS2 Type: { x : number, y : number, path : { x : number, y : number } [17], name : string } [1][60]

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.PlateWRanElliPores1,
  image.PlateWRanElliPores2,
  image.BT226004
}
rois = {
  rois1,
  rois2
}

▸

[]

✗

≡

Daten

logg(typeof(images))
logg(typeof(rois))
logg(summary(images[[1]]))
logg(summary(images[[2]]))
logg(summary(images[[3]]))

▸

[]

✗

≡

Aufbereitung der Daten

Zunächst visuelle Inspektion der Bilder:

Plot eines Bildes

index = 1
image = images[[index]]
# image = images.cropped[[index]]
plot(image,auto.scale=TRUE,keep.aspect=TRUE)

▸

[]

✗

≡

Boundingboxen als Umrandung des ROI (Platte im X-ray Bild) können manuell ermittelt werden (Verwendung von plot(...,crosshair=TRUE))

Beschneidung der Bilddaten (Cropping)

# Image cropping 
# Crop bboxes for images #1/2 and #3
images.cropped=list()
parameter.cropboxes={
  bbox(left=345,top=100,right=680-1,bottom=930-1),
  bbox(left=345,top=100,right=680-1,bottom=930-1),
  bbox(left=340,top=100,right=620-1,bottom=920-1)
}
for(index in 1:3) {
  images.cropped[[index]]=crop.bbox(images[[index]],
                                    bbox=parameter.cropboxes[[index]])
  logg(summary(images.cropped[[index]]))
}
print("=> images.cropped")

▸

[]

✗

≡

Zum trainieren des CNN werden annotierte (gelabelte) Beispieldaten benötigt. Annotierte Bilder werden aus den ROI Listen der synthetischen (simulierten) Bilder gewonnen.

Transformation der ROI Label Listen

rois.t = list()
# roi coordiante systeme origin: (0,0) [mm]
# width: 40
# scale to pixel coordinates (parameter.cropbox)
parameter.width=40
parameter.height=99
function transformPoints(roi) {
  roi=copy(roi[[1]])
  roi$x=floor((roi$x+parameter.width/2)*parameter.scale)
  roi$y=floor((roi$y+parameter.height/2)*parameter.scale)
  print(paste(roi$x,roi$y))
  roi$path=map(roi$path,function (p) {
    { 
      x=floor((p$x+parameter.width/2)*parameter.scale),
      y=floor((p$y+parameter.height/2)*parameter.scale)
    }
  })
  roi
}
for(i in 1:2) {
  print("===== Image #"+i+" ======")
  parameter.cropbox = parameter.cropboxes[[i]]
  parameter.scale=parameter.cropbox$width/(parameter.width)
  roi = map(rois[[i]],transformPoints)
  rois.t[[i]]=roi
}
print("=> rois.t")

▸

[]

✗

≡

Die Umrechnung ist etwas "tricky", daher besser die Übereinstimmung der ROI Boundingsboxen mit den Poren im Radiographiebild nachfolgend prüfen.

Ausgabe ROI

index.image = 1
index.roi = 1
logg(rois.t[[index.image]][[index.roi]])

▸

[]

✗

≡

Die ROIS (x,y,path) stammen aus der CAD Modellierung der synthetischen Platten. Es fehlen noch rechteckige Bereiche ohne Poren, d.h., Hintergrund. Diese können manuell ausgewählt werden (s.o.)

Background ROIs

rois0={
  {
    bbox(20,20,50,150),
    bbox(270,200,50,150)
  },
  {
    bbox(100,520,50,150),
    bbox(270,200,50,150)
  }
}
print("=> rois0")

▸

[]

✗

≡

Plot eines annotierten synthetischen Bildes (1/2)

use geometry
index=1
image = as.cimg(images.cropped[[index]],mode='uint16')
roil   = rois.t[[index]]
roil0  = rois0[[index]]
for(roi in roil) {
  draw.pixset(image,roi$path,value=65000)
}
for(roi in roil0) {
  draw.rectangle(image,bbox=roi,value=65000)
}
plot(image,auto.scale=TRUE,zoom=1)

▸

[]

✗

≡

Jetzt müssen die Trainings- und Testdaten aus den vorherigen ROI Listen extrahiert werden.

Es sollte eine Balance zwischen Poren- und Hintergrundbeispielen erzielt werden
Das Ergebnis der Prädiktion könnte eine erhöhte Falsch-Positiv oder Falsche-Negative Rate entstehen. Die Anteile von positiven und negativen Beispielen kann variiert werden, je nachdem welche Rate man reduzieren will (bzw. welchen Fall man optimieren möchte). Ist es wichtiger möglichst alle Pixel von Poren zu markieren und ein erhöhtes Rauschen in Kauf zu nehmen, oder ist es besser das Rauschen zu reduzieren?
Die Bildsegmente werden aus den ROIs gewonnen. Vereinfacht können sie durch rechteckige Boundingsboxen approximiert werden. Nimmt man zentrale Pixel aus dem Zentrum wird es i.A. immer eine Überdeckung mit den tatsächlichen Polygonen 8der Porenumrandung) geben. Ansonsten muss man Punkte aus der Boundingbox auf Einschluss in die Porenhülle (also das Polygon) explizit testen (und Punkte außerhalb der Porenhülle verwerfen).
Hier wird der einfache Ansatz gewählt.

Erzeugen der Bildsegmente für die Merkmalsklassen P und N

parameter.segment.size = 20
h2 = parameter.segment.size/2
w2 = parameter.segment.size/2
segments.p = list() # Subimages containing feature
segments.n = list() # background
async for(index in 1:2) {
  image  = images.cropped[[index]]
  roil   = rois.t[[index]]
  roil0  = rois0[[index]]
  for (roi in roil) {
    print(paste(index,roi$x,roi$y))
    # get segments around around roi center position
    for (dx in -2:2) {
      for(dy in -2:2) {
        # TODO: use hull search
        # convert uin16 to float32 and normalize to [0,1]
        segment = as.matrix(image[(roi$y-h2+dy):(roi$y+h2-1+dy),
                                  (roi$x-w2+dx):(roi$x+w2-1+dx)],
                                  mode='float32')/50000
        push(segments.p,{label='P',image=segment})
      }
    }
  }
  # await(.Last.value)
  for (roi in roil0) {
    print(paste(0,index,roi$left,roi$top))
    # get segments around around roi center position
    # TODO: use hull search
    # convert uin16 to float32 and normalize to [0,1]
    #segment = as.matrix(image[(roi$y-h2):(roi$y+h2-1),(roi$x-w2):(roi$x+w2-1)],mode='float32')/50000
    #push(segments,{label='P',image=segment})
    xp = floor(runif(30,roi$left,roi$left+roi$width))
    yp = floor(runif(30,roi$left,roi$left+roi$width))
    for(x in xp) {
      for(y in yp) {
        segment0 = as.matrix(image[(y-h2):(y+h2-1),
                                   (x-w2):(x+w2-1)],
                                   mode='float32')/50000
        push(segments.n,{label='B',image=segment0})
      }
    }
  }
}
await(.Last.value,function () {
  logg(summary(segment))
  logg(length(segments.p))
  logg(typeof(segments.p))
  logg(length(segments.n))
  logg(typeof(segments.n))
})

▸

[]

✗

≡

Plot eines Bildsegments

index = 600
# Pore
segment = segments.p[[index]]
plot(segment$image,auto.scale=TRUE,zoom=0.5)
# Hintergrund
segment = segments.n[[index]]
plot(segment$image,auto.scale=TRUE,zoom=0.5)

▸

[]

✗

≡

Frage. Würde man mittels visueller Inspektion der Segmentbilder Poren von Hintergrund unterscheiden können?

CNN Pixel Klassifizierer

Erzeugen des CNN Pixelklassifizier Modells

use cnn
parameter.segmentSize=20
parameter.features=["P","B"]
parameter.filterA=8
parameter.filterB=8
m=cnn(
  input.size=[parameter.segmentSize,parameter.segmentSize,1],
  layers={
    {type='conv',kernel.size=[5,5],filter=parameter.filterA,stride=1,padding=2,activation='relu'},
    {type='pool',kernel.size=[2,2],stride=2},
    # {type='conv',kernel.size=[5,5],filter=parameter.filterB,stride=1,padding=2,activation='relu'},
    # {type='pool',kernel.size=[2,2],stride=2},
    {type='softmax',num.classes=length(parameter.features)}
  },
  labels=parameter.features #optional
  # trainer={method='sgd',learning.rate=0.01,momentum=0.2,batch.size=10}
  # trainer={method='adadelta',l1.decay=0,l2.decay=0.001,batch.size=1}
  # trainer={method='adadelta',l1.decay=0,l2.decay=0.001,batch.size=10}
)
logg(summary(m))

▸

[]

✗

≡

Der Pixelklassifizierer besteht aus:

Einer ersten Filterschicht (Convolution) mit 4-8 Filtern (Kernels) der Größe 5 × 5 pixel. Jedes Filter ist mit 25 Parametern unabhängig konfigurierbar und trainierbar. Diese Schicht wird direkt auf die kleinen maskierten Eingabilder als Teilausschnitt des großen Bildes angewendet.
Einer Poolingschicht. Im Grunde wieder eine Filterschicht, aber nicht dynamisch parametrisierbar sonder mit einer statischen Funktion, hier die Maximumsfunktion. Diese wählt aus der Filtermaske den größten Wert als Wert für die Zielmatrix aus.
Optional weitere Convolution-Pooling Schichtenpaare.
Eine neuronale vollständig verbundene Schicht, d.h. jedes Neuron dieser Schicht (mindestens eines oder eines pro Zielklasse!) wird mit allen Elementen der Matrizen der letzten Ausgabeschicht verbunden. Die Verbindungen sind gewichtet und parametrisierbar.

Trainings- und Testdaten

Zunächst müssen die annotierten Bilddaten aufbereitet werden.
Es wird eine Datentabelle mit einer x- und einer y-Spalte erzeugt.
Die x-Spalte enthält die Bildsegmente, die y-Spalte die Labels

Erzeugen der gesamten Datentabelle

segments = c(segments.p,segments.n)
data = data.frame(
  x=pluck(segments,'image'),
  y=pluck(segments,'label')
)
logg(typeof(data))

▸

[]

✗

≡

Die Klassen B (Background) und P (Pore) sollten gleichermaßen in den Trainings- und Testdaten enthalten sein. Die Zusammensetzung z.B. einer Spalte einer Datentabelle kann mit der table Funktion einfach bestimmt werden (s.u.).

Die Trainingsdaten sollten in der Reihenfolge randomisiert werden (Shuffling). Dies kann einfach mit der sample Funktion erreicht werden.

Erzeugen der Trainings- und Testdatentabelle (Sampling)

samples.training = 1000
samples.test     = 1000
training = sample(data,samples.training)
test     = sample(data,samples.test)
logg(typeof(training))
logg(table(training$y))
logg(table(training$y))

▸

[]

✗

≡

Training

Training des Modells

parameter.epochs=1
xy=training
print(table(as.vector(xy$y)))
train(m,
  xy,
  formula=y~x,
  epochs=parameter.epochs,
  callback=function (result) {
    cprint(result)
})
await(.Last.value,function (result) {
  print('==== Finished ========')
  cprint(result)
})

▸

[]

✗

≡

Aufgabe. Führe das Training durch. Wähle zunächst 10 Epochen. Dann notiere den letzten Klassifikationsfehler (error) und den "Verlust" (loss). Führe nach 10 Epochen nachfolgenden Test mit den Trainings- und Testdaten durch. Gibt es Unterschiede, und warum könnte der fehler des nachfolgenden Tests (mit den Trainingsdaten) vom letzten ausgegeben Fehler des Trainings abweichen? Führe weitere 10 Trainingsepochen durch (usw.).

Test

Nach einigen Trainingsdurchläufen sollte das Modell mit den Trainings- und Testdaten exemplarisch getestet werden.
Alternativ einen weiteren randomisierten Ausschnitt aus der großen Datenmenge (Trainings- und Testdaten sind hier nur Untermengen)

Test des Modells

xy0=sample(data,100)
# oder
# xy0=training
# oder
# xy0=test
print(table(as.vector(xy$y)))
result = predict(m,x[100])
print(result)
result = predict(m,xy0$x)
result = cbind(result,y0=xy0$y)
print(paste('Error',(1-sum(result$value==result$y0)/100)*100,'%'))

▸

[]

✗

≡

Frage. Was gibt die Prädiktion (predict) zurück? Welche Bedeutung hat der Score Wert?

Inferenz

Um das Modell auf ein gesamtes Bild anzuwenden muss jeweils ein maskierendes Fenster (Segment) über das Bild geschoben werden und für jedes Segment das Modell berechnet werden.
Das Ergebnis der Berechnung ist dann ein Pixelwert im Zielbild (Featurebild)
Skalierung (Normierung) der Intensitätswerte nicht vergessen! Gleiche Skalierung wählen wie bei Trainingsdaten.

Anwendung des Modells auf ein gesamtes Bild

index = 1
image = images.cropped[[index]]
w = width(image)
h = height(image)
# margins
margin = parameter.segment.size
function predictSegment(image,x,y) {
  sx = parameter.segment.size
  sy = parameter.segment.size
  segment = image[(y-sy/2):(y+sy/2-1),
                  (x-sx/2):(x+sx/2-1)]
  x = as.matrix(segment,mode='float32')/50000
  result = predict(m,x)
  if (result$value=='P') 
    result$score
  else
    0
}
logg(predictSegment(image,50,50))
image.feature = matrix(0,h-2*margin,w-2*margin,mode='float32')
async for (y in margin:(h-margin)) {
  print(y)
  for(x in margin:(w-margin)) {
    i=x-margin+1
    j=y-margin+1
    image.feature[j,i]=predictSegment(image,x,y)
  }
}
await(.Last.value,function (result) {
  plot(image.feature,ylim=[0,1],keep.aspect=TRUE)
})

▸

[]

✗

≡

Aufgabe. Wende die Inferenz auf ein synthetisches und das eine gemessene Bild an. Wie sind die Ergebnisse?

Raum für weitere Experimente

# TODO

▸

[]

✗

≡

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