www.martinwerner.de
Kapitel 4
Raster-Daten
381
www.martinwerner.de
Motivation
Die Übungsaufgabe, in der eine Menge aus Polygonen A mit Labels versehen werden soll, indem eine andere Menge gelabelter Polygone B mit diesen Polygonen geschnitten wird und als Label des Polygons das Label genommen wird, mit dem die Polygone die meiste Fläche haben, ist in Geometrie sehr aufwändig zu lösen, zum Beispiel
Für jedes Polygon in A
Finde zu einem Polygon aus A alle scheidenden Polygone aus B
Berechne den Schnitt und dann die Fläche für jeden Schnitt
Summiere die Flächen der gleichen Labels
Weise das Label zu
www.martinwerner.de
A litte inaccuracy saves a lot of calculation…
Eine einfachere Lösung des Problems wird geliefert, indem man den komplexen Schnitt zweier Polygone durch eine geeignete Menge an Punkt-in-Polygon-Tests ersetzt.
Dazu benötigen wir eine Menge an Punkten, die Flächeninhalt mit einer gewählten Genauigkeit repräsentieren kann.
Eine Möglichkeit besteht in einem regulären Gitter: Setze Punkte mit
383
www.martinwerner.de
Gitter
Ein n-dimensionales reguläres Gitter ist eine Menge an Punkten,
deren Abstände in jeder Koordinate fest sind. Diese Abstände heißen Gitterabstände.
Ein georeferenziertes Gitter ist ein Gitter mit Koordinaten in einem Referenzsystem (Koordinatensystem)
www.martinwerner.de
Gitter Labeln
Werte, die einem regulären Gitter zugeordnet werden, sind in natürlicher Weise als Matrix aufzufassen:
j
i ⋯
⋮ ⋱ ⋮
385
www.martinwerner.de
Gitter Labeln
Werden bei festem Gitter jedem Gitterpunkt mehrere Werte zugeordnet, spricht man von Bändern.
www.martinwerner.de
Beispiel: Wellenlängen in optischen Bildern
+ +
387
www.martinwerner.de
Zur Übungsaufgabe
1. Wähle Gitterabstand.
2. Berechne Gitter / Raster
3. Fülle Raster-Layer mit Klassen aus Datensatz B (z.B. int8) 4. Fülle Raster-Layer mit ID aus Datensatz B (z.B. int64)
5. Reduziere auf ID und zähle
Bem.: Das ist durch ein Ablaufen des Rasters einfach möglich.
www.martinwerner.de
Raster- und Vektorkarten
Vektorkarten
Liste von geometrischen Primitiven (Line, …)
+ beliebig genau
+ verlustfrei skalierbar
+ Datengröße abhängig von Komplexität - Komplexe Flächen- und
Rasterkarten
• Information pro Flächenelement / Gitterpunkt
+ konstante Genauigkeit
+ Flächen- und Kollisionsanfragen einfach
+ Systeme können für Bildbearbeitung optimiert sein (MMX, SSE, CUDA etc.)
~ Datengröße abhängig von Fläche
389
www.martinwerner.de
Struktur einer Vektorkarte
www.martinwerner.de
Struktur einer Rasterkarte
391
www.martinwerner.de
Geometrie in Rasterkarte (horizontal, vertikal)
www.martinwerner.de
Geometrie in Rasterkarte (Diagonal)
393
www.martinwerner.de
Schneiden sich diese zwei Linien?
www.martinwerner.de
Pixel werden
entsprechend ihres
Flächenanteils gefüllt oder gefärbt.
Für Farben einfach (weil vom menschlichen Auge korrekt wahrgenommen).
Für z.B. Raster mit
diskreten Werten (unsere Lösung: Sub-Pixeling / Anti-Aliasing
395
www.martinwerner.de
Raster-Darstellung in der echten Welt
gdalinfo germany-cloudfree-simple.tiff Driver: GTiff/GeoTIFF
Files: germany-cloudfree-simple.tiff Size is 67235, 94080
Coordinate System is:
LOCAL_CS["WGS 84 / Pseudo-Mercator", GEOGCS["WGS 84",
DATUM["unknown",
SPHEROID["unretrievable - using WGS84",6378137,298.257223563]], PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]], AUTHORITY["EPSG","3857"],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]]]
Angaben zu - Format
- Anzahl Pixel
- Koordinatensystem (WKT)
www.martinwerner.de
Raster-Darstellung in der echten Welt
Origin = (641594.910993173718452,7370990.488857678137720) Pixel Size = (15.545613327413848,-15.545613327413848)
Metadata:
AREA_OR_POINT=Area Image Structure Metadata:
INTERLEAVE=PIXEL Corner Coordinates:
Upper Left ( 641594.911, 7370990.489) Lower Left ( 641594.911, 5908459.187) Upper Right ( 1686804.223, 7370990.489)
- Ursprung in Geokoordinaten - Größe der Pixel - Fläche oder Punkt - Strukturinformation - Geokoordinaten
(berechnet) - Bänder
397
www.martinwerner.de
Bemerkungen
Rastergeometrie im Detail:
Ein GeoTIFF-Raster besteht aus einer rechteckigen Matrix mit Werten pro Band. Eine Geotransformation besagt, wie sich mit einem Schritt nach rechts (bzw. nach unten) in dieser Matrix die Weltkoordinaten ändern. Diese Veränderung ist konstant.
Elt. Bedeutung
GT[0] X-Koordinate in CRS oben links (für Pixel 0,0)
GT[1] X-Offset in X-Richtung für jeden Pixel-Schritt nach Rechts GT[2] X-Offset in Y-Richtung für jeden Pixel-Schritt nach rechts
(Drehung, oft 0)
GT[3] Y-Koordinate in CRS oben links (für Pixel (0,0))
GT[4] Y-Offset in X-Richtung für jeden Pixel-Schritt nach unten (Drehung, oft 0)
GT[5] Y-Offset in Y-Richtung für jeden Pixel-Schritt nach unten
www.martinwerner.de
Berechne Weltkoordinaten für beliebigen Pixel
Mit der Geotransformation bestimmen die folgenden Gleichungen die Weltkoordinaten eines Pixels mit (ganzen) Koordinaten i,j:
= + ⋅ + ⋅
= + ⋅ + ⋅
399
www.martinwerner.de
Einfache Rasteroperationen
gdal_info:
Zeige Raster-Informationen an gdal_warp:
Zusammensetzen von Teilrastern
Ändern des Gitters / der Projektion (Warping) gdal_translate
Ändere Format
Subsetting / Resampling / Interpolation gdal_rasterize
Übersetze (OGR-)Geometrie in Raster gdal_grid
Erzeuge Raster aus Werten an irregulären Punkten gdal_calc
Einfaches Rechnen mit Rasterlayern (numpy)
www.martinwerner.de
Kapitel 4.1
Rasterisierung
401
www.martinwerner.de
Rasterisierung
Übersetzen von
Geometrieobjekten in Raster
Beispiel: Punkte (mit Radius)
Berechne Pixelkoordinaten in Fließkomma
Berechne anteilige Fläche pro Pixel
Modifiziere Pixel (Blending
Computergrafik / Visual Computing)
Ist der Radius klein genug, kann die Berechnung der Fläche entfallen
www.martinwerner.de
Rasterisierung: Linien
Einfachste Möglichkeit: Interpoliere und Raster Punkte
= + ⋅
Problem: Welche Schrittweite? Tradeoff zwischen Aufwand und Wahrscheinlichkeit, alle Pixel auch zu treffen.
403
www.martinwerner.de
Bresenham-Algorithmus
Beobachtung: Nach Setzen einen Pixels kann der nächste Pixel nur ein Nachbarpixel sein
Einschränkung auf Oktant: Linien, die nach unten-Rechts gehen
Folge:
Der Y-Schritt pro X-Schritt ist kleiner 1 1 Pixel pro X-Schritt zu setzen
Nur Linien in diesem Oktant
www.martinwerner.de
Welcher Y-Schritt
Allgemeine Liniengleichung
−
− = −
−
/
/ /
405
www.martinwerner.de
Optimierung
Der Bresenham-Algorithmus berechnet jetzt nicht die echte Y-
Koordinate, sondern verfolgt die ganzzahlige Y-Koordinate und den Fehlerterm (also die Differenz des echten Y-Wertes und des aktuellen ganzzahligen Y-Wertes. Sobald dieser Fehler größer wird als 1, wird er um 1 verringert und die Y-Koordinate angepasst.
Explizit also Pseudo-Code also:
error=0, int y = y0, float deltaerr = abs( (y1-y0)/(x1-x0)) for x from x0 to x1
putpixel(x,y)
error := error + deltaerr if error ≥ 0.5 then
y := y + sign(deltay) error := error - 1.0
www.martinwerner.de
Beispiel: Bresenham-Algorithmus für Linien (2)
Eine Fallunterscheidung unter Ausnutzung der Symmetrie wird jetzt verwendet, um den Basisalgorithmus auf alle 8 Oktanten abzubilden.
Merke: Der Bresenham-Algorithmus hat immer eine schnelle Richtung und eine langsame Richtung. In die schnelle Richtung (X oder Y)
macht man immer einen Schritt, in die langsame Richtung nur ab und zu.
407
www.martinwerner.de
Bresenham-Algorithmus in Python
− 1
Fehler übersteigt 0.5 Durch Y-Schritt reduziert sich auf − 1 und kann dann negativ sein!
www.martinwerner.de
Beispiel: Bresenham-Algorithmus (Anti-Aliasing)
Das Fehlerglied kann für Anti- Aliasing verwendet werden
es kommen in die langsame Richtung nur der richtige Pixel nach Bresenham und die beiden vertikalen Nachbarn in Frage
Ist der Fehler groß, beeinflusst es den Pixel darüber
409
www.martinwerner.de
Einfaches Beispiel
Mit diversen Vereinfachungen (z.B.
ist die Liniendicke von der Steigung abhängig), können wir eine einfache Anti-Aliasing-Regel umsetzen.
Wäre der Fehler 1 oder -1, so wäre der zentrale Pixel nicht getroffen und müsste Gewicht 0 bekommen
www.martinwerner.de
Linie mit und ohne Anti-Aliasing (einfaches)
411
www.martinwerner.de
Bresenham-Algorithmus für Kreise
Bei einem Kreis ist im ersten Oktant die schnelle Richtung die Y- Richtung. Alle anderen Oktanten entstehen durch symmetrische Operationen (verschieben, spiegeln) der Koordinaten
Abbildung: Wikipedia, abgerufen 2019
Wir betrachten den Kreis um Ursprung
+ =
Für einen (negativen) Schritt in X-Richtung ändert sich die rechte Seite wie folgt:
− 1 + − + = −2 + 1
Analog für Y: + + 1 − − = 2 + 1
Damit kann der Fehler fortgeschrieben werden und ein Schritt ausgeführt werden ,sobald der Fehler zu groß wird.
www.martinwerner.de
Rasterisieren von Polygonen
Zum Rasterisieren von Polygonen gibt es zwei gängige Methoden mit ihren Vor- und Nachteilen
Angenommen, ein Punkt-im-Polygon-Test ist verfügbar, dann kann man für jeden Pixel einen solchen Test mit dem gegebenen Polygon durchführen. Ein Pixel wird gesetzt, wenn er im Innern des Polygons liegt.
Vorteil: Genauigkeit und Kompatibilität mit Projektionen
Nachteil: Langsam, Rundungsartefakte möglich
413
www.martinwerner.de
Kapitel 4.1.2
Punkt-Im-Polygon-Test
www.martinwerner.de
Beispiel: Punkt in Polygon
Gegeben sei ein Polygon durch seine Eckpunkte in Umlaufrichtung und ein Punkt. Liegt dieser Punkt im Innern des Polygons?
415
www.martinwerner.de
Beispiel: Punkt in Polygon
Gegeben sei ein Polygon durch seine Eckpunkte in Umlaufrichtung und ein Punkt. Liegt dieser Punkt im Innern des Polygons?
Jordanscher Kurvensatz:
Jede stetige, schnittpunktfreie Kurve (Jordankurve) teilt den Raum in zwei disjunkte Mengen, deren Rand die Kurve selbst ist. [Beweis
ziemlich aufwändig…]
www.martinwerner.de
Beispiel: Punkt in Polygon
Gegeben sei ein Polygon durch seine Eckpunkte in Umlaufrichtung und ein Punkt. Liegt dieser Punkt im Innern des Polygons?
Jordanscher Kurvensatz:
Jede stetige, schnittpunktfreie Kurve (Jordankurve) teilt den Raum in zwei disjunkte Mengen, deren Rand die Kurve selbst ist. [Beweis
ziemlich aufwändig…]
Algorithmus: Jordan-Test
Werfe ausgehend vom Punkt einen Strahl in irgendeine Richtung und zähle die Schnittpunkte. Fälle:
417
www.martinwerner.de
Außerhalb
www.martinwerner.de
Innerhalb
419
www.martinwerner.de
Nicht entscheidbar
www.martinwerner.de
Rasterisieren mit Punkt-im-Polygon-Test
Der Algorithmus zum Rasterisieren ist nun ziemlich einfach:
Für jeden Pixel, berechne die Koordinaten im Koordinatensystem des Polygons und führe den Jordan-Test durch.
Optimierungen:
Es reicht, horizontale Strahlen zu schicken.
Man kann das Kreuzprodukt verwenden, um zu testen, ob ein Schnitt vorliegt.
421
www.martinwerner.de
Kapitel 4.1.3
Scanline-Algorithmus für Polygone
www.martinwerner.de
Rasterisieren von Polygonen – Scanlines
Grundidee:
Wir durchlaufen das Raster in eine Richtung, typischerweise Zeilenweise (deswegen Scanline-Algorithmus).
Für jede Zeile berechnen wir
die Schnittpunkte mit dem Polygon
Ordnen diese nach der X-Koordinate
Zeichnen von links nach rechts „Runs“ ein
423
www.martinwerner.de
Implementierung
Da die Kanten, die eine Zeile schneiden, meist auch die Kanten der darunter liegenden Zeile schneiden, werden diese vernünftig
organisiert und nicht jedesmal neu berechnet.
Sortiere alle Kanten des Polygons nach minimaler Y-Koordinate
Berechne aus dieser Kantenliste eine lister „aktiver“ Kanten („Active Edge Table“, AET)
Speichere Steigungsinformationen um den Wechsel der Scanline ohne Berechnung durchführen zu können
Ignoriere horizontale Linien
Behandele Punkte geeignet (doppelt zählen…)
www.martinwerner.de
Polygon-Füllen
1. Edge-Table (sortiere nach Y- Anfang der Kanten, danach nach X)
Ignoriere horizontale Kanten
1 2
3
4
425
www.martinwerner.de
Polygon-Füllen
2. Initialisiere Active Edge Table Nutze Anfang der Liste, um die erste Scanline zu suchen.
Fülle Active Edge Table, Bestimme Schnittpunkte im Raster und sortiere nach X
Jeder Schnittpunkt schaltet um von Füllen zu Nicht-Füllen und umgekehrt.
Zeichne Scanline
www.martinwerner.de
Polygon-Füllen
3. Zeilenschritt
- Aktualisiere Edge Table
- Entferne Kanten, die enden
- Füge nächste Kante ein, wenn wir die Y- Koordinate erreicht haben
- Verarbeite Schnittpunkte wie
427
www.martinwerner.de
Polygon-Füllen
3. Zeilenschritt
- Aktualisiere Edge Table
- Entferne Kanten, die enden
- Füge nächste Kante ein, wenn wir die Y- Koordinate erreicht haben
- Verarbeite Schnittpunkte wie gehabt
www.martinwerner.de
Polygon-Füllen
3. Zeilenschritt
- Aktualisiere Edge Table
- Entferne Kanten, die enden
- Füge nächste Kante ein, wenn wir die Y- Koordinate erreicht haben
- Verarbeite Schnittpunkte wie
429
www.martinwerner.de
Polygon-Füllen
www.martinwerner.de
Polygon Füllen
Vorteile vom Scanline-Algoirthmus
Speicher-Lokal im Raster (zeilenweise von links nach rechts
Nachteile
Sortieren ist komplex
Auflösungsdifferenz wird nicht ausgenutzt Weitere Alternativen:
z.B. Rasterisiere Randlinie (dicht!) und Fülle im Inneren
431
www.martinwerner.de
Bemerkung
Scanline-Algorithmen sind eine gängige Klasse in der Computational Geometry
Formuliere Algorithmen so, dass sie in einer Richtung monoton (hier mit steigendem Y) verarbeitet werden können
www.martinwerner.de
Kapitel 4.2
Vektorisieren
433
www.martinwerner.de
Vektorisierung
Vektorisierung ist der Versuch, eine im Raster gefasste Geometrie wieder als Vektor-Geometrie darzustellen
Vektorisierung ist komplex, da z.B. nicht klar ist, welche Objekte zusammen gehören.
Wir betrachten nur einen vereinfachten Fall: Ein einzelnes Polygon ist mit Bresenham-Linienalgorithmus gemalt worden, hat daher keine Löcher. Wie rekonstruieren wir die Geometrie?
www.martinwerner.de
Möglichkeit: Direkt als Umlauf
Starte mit blauem Pixel
Suche unter den 9 Nachbarn gegen den Urzeigersinn von der Umkehrung des eigenen Schrittes ausgehend den ersten Pixel
Konstruiere Liniensegment
435
www.martinwerner.de
Strukturverbesserung
Das sind nun zu viele Ecken, einer pro Pixel
Ecken in dieselbe Richtung (180°) können entfernt werden
Polygon war mit großer Wahrscheinlichkeit nicht im Gitter, in dem nur 45° dargestellt werden können
Dieses Problem löst man mit Linienvereinfachung:
Gegeben eine Linie (einen Polygonring), finde einen Ring mit weniger Punkten, der geometrisch und/oder visuell so wenig wie möglich abweicht.
www.martinwerner.de
Kapitel 4.3
Line Simplification
437
www.martinwerner.de
Kartographische Vereinfachung
Beispiel: Douglas-Peucker-Algorithmus: Gegeben ein Fehlermaß pro Punkt.
So lange Summe aller Fehler zu groß, unterteilte Objekt am Punkt mit größtem Fehler und fahre mit beiden Hälften fort.
www.martinwerner.de
Kartographische Vereinfachung
439
www.martinwerner.de
Kartographische Vereinfachung
Beispiel: Douglas-Peucker-Algorithmus: Gegeben ein Fehlermaß pro Punkt.
So lange Summe aller Fehler zu groß, unterteilte Objekt am Punkt mit größtem Fehler und fahre mit beiden Hälften fort.
www.martinwerner.de
Im Detail / in Python
Rekursion in der Programmierung
Die Fähigkeit eines
Unterprogramms, sich selbst aufzurufen.
Vorteil:
Einfache Formulierung vieler Algorithmen
Nachteil:
441
www.martinwerner.de
Beispiel: Line Simplification mit Douglas Peucker
Wir erweitern den Parameter um den Bereich der Verarbeitung
Array mit den Punkten
Fehlerschranke
Anfangs und End-Index
Rückgabe der Vereinfachung: Indizes der noch enthaltenen Punkte
Algorithmus
Berechne den Fehler für jeden Punkt im Bereich der Verarbeitung, der entsteht, wenn der gesamte Bereich durch ein Liniensegment vom
Anfangs zu Endpunkt ersetzt wird. Bestimme Index des maximalen Fehlers. Ist der Fehler größer als die gegebene Fehlerschranke, rufe rekursiv auf mit den Bereichen Anfangspunkt bis Punkt maximalen Fehlers und Punkt maximalen Fehlers bis Endpunkt.
www.martinwerner.de
Kartographische Vereinfachung
Der Algorithmus terminiert?!
Falls der Fehler die Eigenschaft hat, dass der Fehler der Eingabe verglichen mit sich selbst Null wird und der Fehler nicht-negativ ist, terminiert diese Vorschrift – im Worst Case mit der Eingabe
Ist der Algorithmus korrekt?
Der Algorithmus ist korrekt, weil zum Zeitpunkt der Terminierung kein Segment entsteht, in dem das Fehlermaß überschritten wird.
443
www.martinwerner.de
Welchen Fehler?
Wir verwenden im klassischen Douglas-Peucker-Algorithmus den orthogonalen Fehler, also den Abstand eines Punktes aus der Ursprungstrajektorie zum Kandidaten-Liniensegment.
Dieser ist lokal und generalisiert (wenn an keiner Stelle im Innern eines Segmentes der Fehler der Vereinfachung gegenüber dem Original größer ist als die Fehlerschranke, dann gilt dies auch für die zusammengesetzte vereinfachte Trajektorie, da der Fehler an den Stützpunkten Null ist (sie liegen ja auf der
Ursprungstrajektorie!)
www.martinwerner.de
Implementierungsbeispiel
void simplify(TrajectoryType &u, size_t istart, size_t iend, DistanceType max_error) {
// simplify segment and recurse
#ifdef DEBUG_douglas_peucker_simplify
std::cout << istart << "," << iend << std::endl;
#endif
if ((istart -iend) < 2) return; // there is nothing I can add.
size_t the_index=0;
DistanceType the_contribution = largest_contribution(u,istart,iend,the_index);
if (the_contribution > max_error) {
booleanMap[the_index] = 1;
Diese Implementierung verwendet eine (Klassen-)globale
445
www.martinwerner.de
Implementierungsbeispiel
DistanceType largest_contribution(TrajectoryType &u, size_t istart, size_t iend,size_t &the_index) {
// find the point which contributes most DistanceType the_contribution = -1;
typename TrajectoryType::value_type s = u[istart];
typename TrajectoryType::value_type e = u[iend];
for (size_t i=istart+1; i < iend; i++) {
DistanceType contribution = (*d)(s,e,u[i]);
if (contribution > the_contribution) {
the_contribution = contribution;
the_index = i;
} }
return the_contribution;
}
Aus libtrajcomp – https://www.github.com/mwernerds/trajcomp
Hier wird eine Point-Line-Distanzfunktion aufgerufen, die verändert werden kann. Ist in der Regel die orthogonale Distanz.
www.martinwerner.de
Implementierungsbeispiel
std::vector<size_t> indizes(TrajectoryType u, DistanceType max_error,
element_segment_distance<typename TrajectoryType::value_type,DistanceType> &d) {
this->d = &d;
booleanMap.resize(u.size());
for (size_t i=0; i < u.size(); i++) booleanMap[i] = 0;
booleanMap[0] = booleanMap[u.size()-1] = 1;
numericMap.resize(u.size());
simplify(u,0,u.size()-1,max_error);
// now build the simplified trajectory std::vector<size_t> ret;
for (size_t i=0; i < booleanMap.size(); i++)
Initialisiere die Property Map, berechne die
447
www.martinwerner.de
Implementierungsbeispiel
TrajectoryType operator()(TrajectoryType u, DistanceType max_error,
element_segment_distance<typename TrajectoryType::value_type,DistanceType> &d) {
auto which = indizes(u,max_error,d);
TrajectoryType ret;
for (auto w:which) ret.push_back(u[w]);
return ret;
}
Aus den Indizes kann man auch wieder eine Trajektorie bauen, indem man die Punkte zusammensetzt.
www.martinwerner.de
Beispiel: Douglas-Peucker (vereinfacht, in Python)
449
www.martinwerner.de
Beispiel: Douglas-Peucker (vereinfacht, in Python)
www.martinwerner.de
Douglas-Peucker Hauptprogramm
451
www.martinwerner.de
Kapitel 4.4
Rasteroperationen / „Map Algebra“
www.martinwerner.de
Typische Rasteroperationen
Bei Operationen zwischen zwei Rastern unterscheidet man zunächst
Local
Der Wert des Ausgaberasters ist berechnet aus den Werten an denselben Orten des Eingaberasters
Focal
Der Wert des Ausgaberasters ist berechnet aus den Werten in einer (rechteckigen) Umgebung des Ortes im Eingaberaster
Zonal
453
www.martinwerner.de
Eingabedaten für die Raster-Operationen
Local Focal Zonal Global
www.martinwerner.de
Beispiel: Lokale Operationen
Lokale Operationen sind gebräuchlich für zwei oder mehr (Stacks) Raster mit gleicher Größe und Projektion.
Zum Beipsiel:
Arithmetic operations (addition, subtraction, multiplication, division)
Statistical operations (minimum, maximum, average, median)
Relational operations (greater than, smaller than, equal to)
455
www.martinwerner.de
Beispiel: Fokale Operation
Gauss‘sche Unschärfe
Berechne Pixel als Mischung der 9 umgebenden Pixeln, z.B. mit Gewichtsmatrix
1 1 1 1 4 1 1 1 1
An jeder Stelle wird der Pixel erzeugt, indem die Pixel aus dem Quellbild gelesen werden, mit der Matrix gewichtet werden und die Summe gebildet wird, die dann den Wert des Pixels ergibt.
Ein ICAML-Tutorial zu diesem Thema:
https://icaml.org/canon/data/images-videos/Python_Convolution/Python_Convolution.html
www.martinwerner.de
Zonale Operation
Zonale Statistik wird oft verwendet, um Zusammenhänge über geographische Gebiete zu modellieren.
z.B. Durchschnittstemperatur in einem Bundesstaat Zonen: Bundesstaaten
Raster: Temperatur
Ergebis: Raster, in dem die Durchschnittstemperatur der Zonen hinterlegt wird
z.B. Durchschnittstemperatur über Wasser
457
www.martinwerner.de
Globale Operationen
Eine globale Operation ist beispielsweise das Euklidische Distanztool.
Für jede Zelle wird die Distanz zur nächsten Quelle berechnet.
Abb. aus http://desktop.arcgis.com/de/arcmap/10.3/tools/spatial-analyst- toolbox/understanding-euclidean-distance-analysis.htm
Quelle: Raster mit 1 Pixel pro Stadt Ergebnis: Distanzkarte
www.martinwerner.de
Beispiel: Morphologie-Operatoren
Morphologie basiert auf Strukturelementen (SE), die zur lokalen Analyse eines Bildes verwendet werden
Der Einfachheit halber betrachten wir zunächst nur Morphologie auf binären Bildern (schwarz und nicht-schwarz)
Ein Strukturelement ist dabei ein zusammehängender, in der Regel kleiner Graph / kleines Bild-Element, z.B.
459
www.martinwerner.de
Grundoperatoren der Morphologie
Erosion: Ein Pixel des Resultats wird gesetzt, wenn das SE
vollständig in das Bild passt (also die gesetzten Pixel Teilmenge der im Bild gesetzten Pixel ist)
= { ∶ ⊆ }
A
oder
A) Jeder Pixel, für den alle Nachbarpixel gesetzt sind
B) Jeder Pixel, für den die Pixel links, rechts, drüber, und drunter gesetzt sind C) Jeder Pixel, der Teil einer diagonalen Linie ist
B C
www.martinwerner.de
Grundoperatoren der Morphologie
Dilatation: Ein Pixel des Resultats wird gesetzt, wenn das SE einen nichtleern Schnitt mit den gesetzen Pixeln des Bildes hat
= { ∶ ∩ ≠ ∅}
oder
461
www.martinwerner.de
Eigenschaften von Erosion und Dilataion
Verschiebungsinvariant
Die Erosion (Dilatation) einer Verschiebung ist die Verschiebung der Erosion (Dilatation)
Nicht Idempotent E*E != E / D*D != D
Monoton
Gilt I<J für die Bilder, so gilt E(I) < E(J) und D(I) < D(J)
Dual
Die Negierung der Erosion der Negierung ist die Dilatation und umgekehrt.
www.martinwerner.de
Kapitel 5
Deep Learning in der Geoinformatik
463
www.martinwerner.de
Outline
A Short History of Deep Learning
Deep Learning Elements
Neurons
Neural Networks
Back Propagation and Gradient Descent
Some Basic Deep Learning Architectures
Dealing with Point Clouds
And now? How would I?
www.martinwerner.de
Kapitel 5.1
A Short History of Deep Learning
465
www.martinwerner.de
History of Deep Learning
https://pbs.twimg.com/media/DuE4LnRWs AEpHjQ.jpg:large
www.martinwerner.de
History of Deep Learning
Nothing of this is really new. It is an old and established discipline.
The current hype comes from several factors
Advances in computational performances (GPUs, TPUs)
Creation of Huge Datasets
(Smaller) Advances in Stochastic Gradient Decent
Novel Ideas about Regularization
Novel Ideas for Capacity (Weight) Reduction
Convolutional Neural Networks
But, Deep Learning is not very powerful per se:
467
www.martinwerner.de
Kapitel 5.2
Deep Learning Elements
www.martinwerner.de
Kapitel 5.2.1
Neurons and Neural Networks
469
www.martinwerner.de
Biology-Inspired but simplified
S
Input 1 Input N
w1 wN
Activation
www.martinwerner.de
Linear Neuron
The simplest Neuron is a linear one.
This means
Activation Function is linear
A bias term is added
Then, we can write the output as +
..
For simplicity, the bias is often made an
S
Activation
471
www.martinwerner.de
Let us learn something
Lets assume two inputs to the neuron and f(x)=y the activation function.
Question: What can we represent in this way:
Answer: Lets calculate a bit (with explicit bias)
+ = + +
Now, for binary classification, we need a simple decision rule. What about (output > 0)
Then, we can learn sets that have the structure
+ = + + ≥ 0
This is easily seen to be a split along a line in space. Lets do this
www.martinwerner.de
This is a typical linear neuron decision
473
www.martinwerner.de
However, XOR is impossible to represent with a single neuron
0 1
0 1
There is no line that separates the two colors!
www.martinwerner.de
Solution: Add multiple layers (MLP)
x y
a b
o
This architecture has a bias term for all hidden nodes (a) and the output node which is hidden.
That is, there are nine weights!
Each of the early neurons decides a) Above the line
b) Below the line
1 a
B
475
www.martinwerner.de
The first scientific! AI winter
(the term AI winter is used for periods of cut funding as well)
Now, for a long time, no real progress was made. People got frustrated, left the field. The frustration points were:
Finding optimal weights is NP-complete – exponential runtime
While solving XOR is possible with a MLP, it is impossible to train, because the expected output of the inner connections is unknown.
Many people turned away from this part of machine learning
Dates are difficult to assign as related machine learning techinques are still evolving:
Starts about the time that the implications of the unsolvability of XOR for general intelligence become clear
Challenge Problem has been identified: train MLP
Ends about the time where multilayer perceptrons are successfully trained
Challenge Problem has been fully solved without avoiding it.
www.martinwerner.de
The solution to the MLP problem
477
www.martinwerner.de
Back Propagation
Where do the weights come from?
Finding the optimal weights is NP-complete (that is, as hard as the TSP;
Blum and Rivest, 1992)
Fortunately, we can find a sufficient set of weights through back propagation (e.g., Rumelhart et al. (1985))
First, we compare the output of a forward pass with the expected value.
Then, we slightly adjust each of the weights backwards in the network by a very small amount.
We do this over and over again (training)
We do so, because the error function we chose is differentiable and
sufficiently smooth such that the local direction of error reduction is sensible globally (which need not be the case)
www.martinwerner.de
Backpropagation Details
Forward Pass
All units within a layer have their values set in parallel
Next layer only after first layer has completely been computed
Layer Function needs to
Have bounded derivative only
However, linear aggregation of the input before applying one non-linear function simplifies learning procedure
Total Error Function
479
www.martinwerner.de
Backpropagation Rules (*)
Let us fix a single case c. Then
= −
Now, let denote the activity of a unit in the forward pass. Then use the chain rule
= ⋅
Now, with an activity function of = we can calculate and substitude the second factor:
== ⋅ 1 −
This means, we know how the total input of node changes the total error for this case. But as the total input is a linear sum of the inputs, we can compute
www.martinwerner.de
Backpropagation Rules (*)
= ⋅ = ⋅
And analogously, we can calculate derivatives for y
= ⋅ = ⋅ w
Now, we have seen how to calculate for any unit in the penultimate layer when given information from the last layer
This can be iterated backwards such that the derivatives become known
481
www.martinwerner.de
Fully Explained
It is a very good idea to spell out this for the XOR problem. You can follow the following article (using different names than here)
https://medium.com/@14prakash/back-propagation-is-very- simple-who-made-it-complicated-97b794c97e5c
One way of thinking about back-propagation is that it is a major factorization of the derivative into things that we can calculate as numbers!
= ⋅ ⋅
https://medium.com/@14prakash/back-propagation-is-very-simple-who-made-it- complicated-97b794c97e5c
www.martinwerner.de
Now allows many architectures
Classical Networks
Input, a few hidden layers, an output
Difficulty: expressivity (number of layers) vs. trainability (number of parameters)
Convolutional Neural Networks and Pooling
Input an image, Layers are now calculating some local convolution of the image and dimensionality is reduced by pooling, that is taking only a subset of the data points.
Less Weights (only once for the convolution kernel which is swiped over the image, not for every pixel)
Recurrent Networks
483
www.martinwerner.de
Second (scientific) AI Winter
Now, Backpropagation can train deep networks and, therefore, XOR,but
Not enough processing power (no GPUs, for example)
Lack of Datasets (big and annotated datasets, because in real-world scenarios you would need those)
Overfitting (mainly, because you need to choose a sufficiently expressive architecture but don‘t have enough data to train)
Vanishing Gradient Problem
During learning, you multiply a lot of very small numbers which eventually get too small for sensible learning on finite accuracy machines
People turned away, because practical examples of deep networks were not brought to significant success, especially as other techniques became very powerful including support vector machines
www.martinwerner.de
Breakthroughs
Training tricks
ImageNet Dataset (2009, 16 million annotated images)
Visibility through ILSVRC (1 million images, 1,000 classes) 2013: AlexNet trained on ImageNet using two GPUs
Dropout
Rectified Linear Units (ReLU) instead of sigmoid or tanh activations
485
www.martinwerner.de
In computer vision
CNNs:
Use a focal operation (Filter each pixel with a set of shared matrices) and sample taking
maximum or average value in filter block
Results:
Errors drop significantly year by year
Architectures get deeper and deeper
Trainable with tricks
Some results from the golden years of CNNs follow
www.martinwerner.de
CNN: Übersicht…
487
www.martinwerner.de
ILSVC over the early years
Slide from a taken from the Internet…
www.martinwerner.de
2015
In 2015, Microsoft Research Asia won with a 150 layer network
Almost superhuman performance (3.5 % error, later even improved)
GoogLeNet 2014 had 22 layers
Is the next AI winter just around the corner?
We have been successful in image regognition, speech, and translation.
But we rely on excessive datasets that we cannot generate in
489
www.martinwerner.de
Kapitel 5.3
Some Basic Deep Learning Architectures
www.martinwerner.de
Architectures
Perceptron (P)
491
www.martinwerner.de
Architectures
Feed Forward (FF)
www.martinwerner.de
Architectures
Deep Feed Forward (DFF)
Architectures
Recurrent Neural Network (RNN)
Long / Short Term Memory (LSTM)
Gated Recurrent Unit (GRU)
www.martinwerner.de
Architectures
Auto Encoder (AE)
495
www.martinwerner.de
Architectures
Deep Convolutional Network (DCN)
www.martinwerner.de
Architectures
Deconvolutional Network (DN)
497
www.martinwerner.de
Architectures
Deep Residual Network (DRN)
www.martinwerner.de
Exkurs 5.4
Deep Learning with Point Clouds
499
www.martinwerner.de
A first Success Story: PointNet
www.martinwerner.de
Why?
501
www.martinwerner.de
Classical Point Cloud Treatment
Extract hand-crafted features (e.g., structure tensor + friends)
Should be invariant for certain transformations
Can be global or local
Usually need a context definition (for pure 3D points)
Including Deep Feed-Forward Architectures!
Volumetric CNNs
Step towards a voxelgrid and use (learned) 3D convolutions
Multiview CNNs
Render several perspective views of the point clouds and feed them to a CNN
Limited to aspects represented by 2D aspects (e.g., classification, but not completion)
www.martinwerner.de
Central challenge
Point Clouds are Unordered Collections of Points
and there is no sensible ordering function Model Functionalities Needed
Classification outputs a score for each candidate class
For Scene Understanding / Segmentation, the model outputs scores for each point and each candidate class
503
www.martinwerner.de
A first glance at the architecture
www.martinwerner.de
PointNet Architecture
Based on three main properties, assertions and their consequences
The order of the points shall not matter
Nearby things shall be able to interact with each other
The system should become invariant under rigid transformation including rotation, translation, and flip
505
www.martinwerner.de
Treating the Order of Points
To make a model invariant under the order of input points can be done basically in three ways:
Sort input into a canonical order,
However, no order exists that preserves data locality completely
Treat the input as a sequence and train with all permutations of the input
However, it has been shown that order matters still.
Excessive training times (There are n! permutations)
Use a simple, symmetric function to aggregate information from each point
Okay, lets go for it…
www.martinwerner.de
The symmetric function
It would be easy to use addition or multiplication as they are perfectly commutative. But more flexibility is needed and a trainable (learnable) function is preferred.
Therefore, f is a function mapping the point cloud to a single real number
507
www.martinwerner.de
Local and Global Information Treatment
For now, we just transformed the whole point cloud into a single feature vector
…
We can now just train any machine learning system like a SVM or a MLP on this very result
However, this can only rely on global information
But, we will need a combination of local and global information
This is done in the Segmentation Network
www.martinwerner.de
The Segmentation Network
It concatenates 64 per point features with 1024 global features for a matrix of nx1088 of features
509
www.martinwerner.de
Invariance w.r.t. rigid transformations
The remaining piece is how to achieve invariance under rotation, translation etc.
Idea: Predict an affine transformation matrix (T-Net) and apply this transformation to the input points
These mini-networks have the same structure as the global
network: point independent feature extraction, max pooling, and fully connected layers
This can as well be applied again to the feature space.
But beware, it is a large matrix and difficult to optimize
Therefore, a constraint makes it almost orthogonal by adding to the loss
www.martinwerner.de
Why PointNet? Because it looks nice and works in practice
511
www.martinwerner.de
Yes, it works…
www.martinwerner.de
… but it is also theoretically sound!
513
www.martinwerner.de
Exkurs 5.4.1
Extensions for PointNet++
www.martinwerner.de
Extension towards real-world problems
PointNet uses a single Max-Pooling layer, which means that all features are single-scale
Point Clouds have varying sampling density, especially with fixed sensors
PointNet++ is based on a hierarchical grouping analyzing larger and larger extracts of the point cloud
Implemented as Compression: At each and every step, a point set is abstracted to a point set with fewer points
515
www.martinwerner.de
Only the ideas
Sampling Layer
Iterative Farthest Point Sampling (FPS)
Iteratively add the farthest point from the input to the current set
Grouping Layer
Assign some neighboring points using a
ball query
Pro: same scale, Con: different number of elements
kNN
Pro: same number of elements, Con: different scale
Ball query preferred as PointNet can deal with varying inputs
Many additional tricks
See https://arxiv.org/pdf/1706.02413.pdf
www.martinwerner.de
Kapitel 5.4.2
Wie macht man das wirklich („hands-on“)
517
www.martinwerner.de
Lets go for it
Run a computer / container with tensorflow
I am runnning NVIDIA‘s optimized tensorflow container (need an account at NVIDIA container registry)
Optimized by NVIDIA for DGX-1 Familiy
On 8 interconnected V100 GPUs (256 GB total memory)
Trains about 2 hours to 88.8 % accuracy on point cloud classification for ModelNet 40 dataset
Inside the container (or in the Dockerfile)
apt-get install libhdf5-dev (for HDF5 file support)
pip install h5py
git clone https://github.com/charlesq34/pointnet
python train.py
Automatically downloads dataset
Runs a few epoochs and outputs results
www.martinwerner.de
Hands-On
root@ede2a32eccac:/workspace/pointnet# python train.py --2019-01-17 06:41:18--
https://shapenet.cs.stanford.edu/media/modelnet40_ply_hdf5_
2048.zip
Resolving shapenet.cs.stanford.edu (shapenet.cs.stanford.edu)...
171.67.77.19
Connecting to shapenet.cs.stanford.edu
(shapenet.cs.stanford.edu)|171.67.77.19|:443... connected.
HTTP request sent, awaiting response... 200 OK
519
www.martinwerner.de
Results look like
[…]
eval mean loss: 0.549058 eval accuracy: 0.886769
eval avg class acc: 0.860618
**** EPOCH 249 ****
[…]
----1---
eval mean loss: 0.546670 eval accuracy: 0.888393
eval avg class acc: 0.858817
(after less than two hours including data download on a NVIDIA V100 GPU, e.g., DGX-1)
www.martinwerner.de
Per-class performance and visualization of errors
$> pip install scipy
$> pip install image # for PIL
$> pip install matplotlib # for visualizations
$> python evaluate.py –visu
This now creates output of erroneous classifications in the dump folder and gives per class performance results. Looks like
airplane: 1.000 bathtub: 0.860 bed: 0.980 bench: 0.700 bookshelf: 0.900
521
www.martinwerner.de
Some Results
535
www.martinwerner.de
