Post on 17-Sep-2018
transcript
FORUM 3D, TECHNISCHE UNIVERSITÄT MÜNCHEN
Automatische Koregistrierung von Laserdaten aus mehreren Schräg-
ansichten städtischer Gebiete
Marcus Hebel 1, Uwe Stilla 2
1 Forschungsinstitut für Optronik und Mustererkennung, Ettlingen 2 Fachgebiet Photogrammetrie und Fernerkundung, Technische Universität München
Kurzfassung: In Schrägsicht von einem fliegenden Träger aufgenommene Laserdaten
erlauben eine gleichzeitige Erfassung von Fassaden und Dachlandschaften. Aufgrund
der durch die Schrägsicht entstehenden Abschattungen sind jedoch mehrere von ver-
schiedenen Ansichten streifenförmig aufgenommene Laserpunktwolken zu koregistrie-
ren und zu kombinieren. Die vorliegende Arbeit beschreibt die Vorgehensweise zur au-
tomatischen Vorklassifizierung der 3D Einzelpunkte und zur anschließenden Koregist-
rierung anhand erkannter Häuserdächer. Anstelle des ICP-Verfahrens kommt dabei
eine neue Methode zum Einsatz, die die korrekten Zuordnungen zwischen den Daten-
sätzen implizit herstellt. Die vorgestellte Verfahrensweise wird anhand von vier exemp-
larischen Datensätzen getestet, die das Stammgelände der Technischen Universität
München aus verschiedenen Ansichten zeigen.
1 Einleitung
1.1 Airborne Laser Scanning
Die Basis des luftgestützten Laserscanning (ALS, engl. airborne laser scanning) bildet
ein LiDAR-Messverfahren (engl. light detection and ranging), das direkte Entfer-
nungsmessungen liefert und sich im Bereich der Photogrammetrie und Fernerkundung
zur Geländeaufnahme etabliert hat. So werden ALS-Daten heute oftmals als Basis für
eine 3D-Stadtmodellierung verwendet. Typische Anwendungen können bei der Stadt-
planung, dem Tourismus und der Telekommunikation (Planung der Standorte für Mobil-
funkanlagen) gesehen werden. Eine Übersicht und umfassende Beschreibung der ALS-
Grundlagen wird z.B. in [18] gegeben. Lasermessungen weisen einige Vorteile im Ver-
gleich zu klassischen Luftbildern auf. So liegen die Aufnahmen direkt als 3D-Daten vor,
die Messungen sind unabhängig von vorliegenden Beleuchtungsverhältnissen und es
wird eine hohe relative Genauigkeit und Punktdichte geliefert, typischerweise mehr als
100.000 Punkte pro Sekunde. Heutzutage kommerziell verfügbare Laserscanner wie
der von uns benutzte Riegl LMS-Q560 sind darüber hinaus in der Lage, den gesamten
Signalverlauf des reflektierten Laserpulses aufzuzeichnen, was neue Methoden der Da-
tenauswertung ermöglicht [8]. Zusätzlich zur Erfassung mehrerer Entfernungswerte pro
Laserpuls (z.B. verursacht durch Vegetation) können Merkmale wie Intensität oder Brei-
te der reflektierten Pulse untersucht werden. Die aufgezeichneten Entfernungsdaten
werden in der Regel georeferenziert (z.B. UTM-Koordinaten), wozu typischerweise eine
Kombination aus Differential-GPS Empfänger und kreiselbasierter Inertialsensorplatt-
form (IMU, engl. inertial measurement unit) am Sensorträger zum Einsatz kommen. Für
unsere Messungen wurde das kommerziell erhältliche Navigationssystem POS AV der
Firma Applanix eingesetzt, das auf einem Hubschrauber Bell UH-1D montiert war. Ab-
bildung 1 zeigt einen Ausschnitt der über München erfassten 3D Laserpunktwolke.
Abb.1: 3D Laserdaten aufgezeichnet in 45° Schrägsicht vom Testgebiet TUM. Im Vordergrund ist die Alte Pinakothek zu erkennen, im Hintergrund das Stammgelände der Technischen Universität München
Es handelt sich bei Abbildung 1 nicht um eine photographische Aufnahme, sondern um
eine 3D Punktwolke, deren bildhafter Eindruck durch die Einfärbung der Einzelpunkte
anhand der gemessenen Intensität der rückgestreuten Laserpulse (nahes Infrarot bei
1.5 µm Wellenlänge) entsteht. Am Boden kann man die einzelnen Scanzeilen des La-
serscanners erkennen.
1.2 Problembeschreibung
Üblicherweise werden flugzeuggetragene Lasersensoren in Nadirsicht betrieben, d.h.
mit Blick senkrecht von oben auf die Szene. Die von uns verwendete Sensorkonfigura-
tion erlaubt dagegen auch die Erfassung von Fassaden. Allerdings führt diese Sensor-
anordnung zunächst dazu, dass aus Flugrichtung gesehen hinter Gebäuden Abschat-
tungen entstehen, weil der schräg eingerichtete Lasersensor dort keinen Einblick hat.
Ausgleichen lässt sich dies durch eine Mehrfachabtastung des gleichen urbanen Ge-
biets aus verschiedenen Richtungen mit anschließender Koregistrierung der einzelnen
Datensätze.
Eine Nachverarbeitung ist leider notwendig, da Ungenauigkeiten bei den Navigationsda-
ten bestehen können [10]. Vor allem wenn wie in unserem Fall eine eigene GPS-
Bodenstation für Differential-GPS fehlt, können die gemessenen Sensorpositionen um
mehrere Meter abweichen. Dies lässt sich allerdings noch nachträglich durch Korrektur-
daten des Satellitenpositionierungsdienstes der deutschen Landesvermessung
(SAPOS) berichtigen. Zusätzlich aber führen unvermeidbare Abweichungen der Anbrin-
gung des Sensors am Sensorträger zu einem systematischen Versatz der Messwerte,
den es im Hinblick auf eine Sensorkalibrierung herauszufinden gilt. Man kann hierbei
aber auf jeden Fall davon ausgehen, dass die Einzeldatensätze sich bereits grob in Po-
sition befinden. Eine Feinregistrierung ist jedoch noch notwendig, um einen konsisten-
ten Gesamtdatensatz zu erhalten.
1.3 Vorgehensweise
Der ICP-Algorithmus (engl. iterative closest point, vgl. Abschnitt 2.2) von Besl und
McKay [2] ist das Standardverfahren zur Koregistrierung von Punktmengen. Aufgrund
der beschriebenen Abschattungen kann das ICP-Verfahren nicht ohne weiteres auf die
hier vorliegenden Punktwolken angewendet werden, da der Algorithmus anfällig für
nicht-überlappende Bereiche in den Datensätzen ist [13]. Vorliegende geometrische
Beziehungen werden vom ICP-Verfahren nicht beachtet, daher können falsche Zuord-
nungen aufgrund teilweise fehlender Punkte (Verdeckungen) zu schlechten Resultaten
führen, wenn die ICP-Iterationen in ein falsches lokales Minimum hineinlaufen [14]. Um
dies zu vermeiden zielt unsere Vorgehensweise darauf ab, zunächst alle Punkte ausfin-
dig zu machen, die am vielversprechendsten für eine korrekte gegenseitige Registrie-
rung der Datensätze sind. Da die Sichtbarkeit von Häuserfassaden im Gegensatz zu
der von Dächern stark von der Anflugrichtung abhängt, müssen die dortigen Punkte als
erstes herausgefiltert werden. Ebenso treten bei einem flugzeuggetragenen Sensor be-
vorzugt Verdeckungen durch Gebäude am Boden auf, weswegen auch die Bodenpunk-
te nicht für die Koregistrierung verwendet werden sollten. Wie in [10] beschrieben, sind
auch Objekte mit einer unregelmäßigen Form wie etwa Bäume und sonstige Vegetation
hinderlich für eine erfolgreiche Zuordnung der Datensätze, da die diskreten Messpunkte
hier sehr unregelmäßig verteilt sind. Folglich sollten auch diese Punkte nicht berück-
sichtigt werden. Vor allem die Dächer von Gebäuden weisen charakteristische Merkma-
le auf, die eine gegenseitige Feinregistrierung der ALS-Datensätze zulassen.
Allerdings, selbst wenn es gelingt ausschließlich die Messpunkte auf Dachflächen zu
identifizieren, besteht dann noch ein Problem für das Standard-ICP-Verfahren. Erstens
ist aufgrund der unterschiedlichen Blickwinkel des Laserscanners auf Dachflächen die
dortige Punktdichte sehr unterschiedlich, zweitens kann man ohnehin nicht erwarten,
bei zweimaliger Abtastung der Szene die gleichen Messpunkte zu erhalten. Das hier
vorgestellte Verfahren passt daher zunächst mit einem robusten Schätzverfahren
(RANSAC) Ebenen an die Datenpunkte an, um so die Dachflächen geometrisch zu be-
schreiben. Jeweils drei nicht-parallele Dachebenen schneiden sich in einem virtuellen
3D-Punkt. Nun ist es vergleichsweise einfach, identische Dachflächen in den grobregis-
trierten Datensätzen ausfindig zu machen. Daher korrespondieren auch die jeweiligen
Schnittpunkte miteinander, sodass anders als beim ICP-Verfahren mit dieser Methode
keine Punkt-zu-Punkt Zuordnungen mehr gesucht werden müssen. Die Bestimmung
der zur korrekten Feinregistrierung führenden Transformation (Translation, Drehung)
geschieht nun anhand der einander implizit zugeordneten virtuellen 3D-Punkte.
1.4 Vergleichbare Arbeiten
In den letzten Jahren sind sehr viele Artikel zur Koregistrierung von Punktmengen veröf-
fentlicht worden. Seit Besl und McKay im Jahr 1992 den ICP-Algorithmus vorgestellt
haben, hat sich dieses Verfahren als Standardlösung für das Registrierungsproblem
etabliert. Auch unser Ansatz baut in Teilen darauf auf, obwohl wir das Problem der
Punktkorrespondenzsuche umgehen. Erwähnenswert ist, dass es in der Vergangenheit
einige alternative Konzepte zur geometrischen Ausrichtung von Punktmengen gegeben
hat, so z.B. von Pottmann et al. im Artikel „Registration without ICP“ [12]. Andere Me-
thoden zur Minimierung der Fehlerquadrate können benutzt werden, um Zuordnungen
zwischen benachbarten ALS-Streifen herzustellen [10]. Die Tatsache, dass Daten von
Laserscannern unregelmäßig im Raum verteilt sind, wird von vielen Autoren durch Ein-
führung einer Dreiecksvermaschung (TIN, engl. triangulated irregular network) ange-
gangen [9]. Wir verfolgen für die vorliegenden Daten einen geeigneteren Ansatz, indem
wir ein regelmäßiges 2D Bodengitter zugrunde legen und pro Rasterfläche nur den
höchsten vorkommenden Messpunkt betrachten. Das erlaubt gleichzeitig Punkte von
Häuserfassaden herauszufiltern. Zudem ist diese Datenstruktur sehr effizient nutzbar
für Suchoperationen und die Untersuchung von Nachbarschaften.
Es hat sehr viele Ansätze zur Verbesserung des Standard-ICP-Verfahrens gegeben,
sogar so viele, dass Rusinkiewicz und Levoy sich in [14] die Mühe machten, existieren-
de ICP-Varianten in Kategorien einzuteilen. Die beiden Autoren, die am „Digital Michel-
angelo Project“ [16] beteiligt waren, unterscheiden dabei die einzelnen Verfahrensvari-
anten danach, welcher Einzelschritt des Original-ICP-Verfahrens verändert wurde:
(1) Auswahl von Untermengen der Punktwolken
(2) Zuordnung zwischen Einzelpunkten
(3) Gewichtung von Punkt-zu-Punkt Korrespondenzen
(4) Zurückweisung einzelner Punktpaarungen
(5) Auswahl einer geeigneten Fehlermetrik
(6) Methode zur Minimierung der Fehler
In dieser Kategorisierung würde unser Verfahren am ehesten die Punkte (1), (2) und (4)
betreffen. Mit allen bekannten Modifikationen des ICP-Verfahrens wird versucht, dessen
Robustheit, Performanz und/oder Genauigkeit zu verbessern. Ein kritisch zu betrach-
tender Nachteil des ICP-Verfahrens ist die mangelnde Robustheit gegenüber Ausrei-
ßern in den Datenpunkten [4]. Wir behandeln dieses Problem, indem wir zunächst nur
diejenigen Messwerte greifen, die am ehesten korrespondierende Punkte in den jeweils
anderen Datensätzen haben sollten. Ferner wird durch die Betrachtung von Ebenen-
schnittpunkten der Ausreißeranteil erheblich reduziert. Eine andere Möglichkeit wäre ein
Herausfiltern offensichtlich falscher Korrespondenzen, etwa durch Auswertung eines
variablen Distanz-Schwellwerts, wie es in [19] getan wird. Viele Autoren (z.B. [3], [6])
ersetzen die übliche euklidische Metrik durch andere Distanzmaße, um so zusätzliche
Merkmale (Intensitätswerte, lokale Normalenvektoren) auswerten zu können.
Ein Unterpunkt unseres Ansatzes ist die Segmentierung ebener Flächen in den Punkt-
wolken. Hierfür wurden bereits die unterschiedlichsten Techniken veröffentlicht. Eine
Übersicht ist in [7] zu finden. Einige Autoren sind daran interessiert, neben Ebenen
auch z.B. sphärische, zylindrische oder kegelförmige Objekte zu segmentieren. In [13]
werden zwei Methoden beschrieben, Modelle durch Minimierung von Fehlerquadraten
an die Daten anzupassen. Dagegen verwendet man in [17] die dreidimensionale
Hough-Transformation zur strukturellen Analyse der Punktmengen. Unter den verfügba-
ren Verfahren hat vor allem der RANSAC-Algorithmus [5] einige Vorteile bei der Extrak-
tion von festgelegten Formen aus den Datensätzen [15]. Wir verwenden das RANSAC-
Verfahren zur robusten Schätzung der Ebenenparameter, außerdem liefert der sich er-
gebende Anteil sogenannter Outlier und Inlier eine wichtige Information darüber, ob
man es an der vorliegenden Position im Datensatz eher mit planaren Strukturen (z.B.
Gebäude, Straßen) oder unregelmäßig geformten Objekten (z.B. Vegetation) zu tun hat.
2 Überblick über die verwendeten Methoden
Nach dieser Einleitung und Literaturübersicht, in der bereits vom ICP- und RANSAC-
Verfahren die Rede war, werden diese beiden Methoden im Folgenden kurz vorgestellt.
Wir beschränken uns dabei auf das Wesentliche, was zum Verständnis der nachfolgen-
den Abschnitte notwendig ist. Für einen tieferen Einblick sei auf die Originalartikel [2]
und [5] verwiesen.
2.1 Random sample consensus (RANSAC)
Beim RANSAC-Verfahren (engl. random sample consensus), wie es von Fischler und
Bolles in [5] beschrieben wird, handelt es sich um eine Technik zur numerischen Schät-
zung von Parametern eines mathematischen Modells, welches einer beobachteten
Menge von Messwerten unterliegt. Die Anwendung des RANSAC-Schätzverfahren ist
vor allem dann sinnvoll, wenn die Daten einen hohen Anteil an Ausreißern enthalten,
die im Zusammenhang mit der Modellschätzung als „Outlier“ bezeichnet werden. Dem-
entsprechend bezeichnet man alle Datenpunkte, die durch das gesuchte Modell be-
schrieben werden, als „Inlier“. Im Gegensatz zur klassischen Methode, unter Verwen-
dung aller Datenpunkte eine Minimierung der Fehlerquadrate herbeizuführen (z.B. line-
are Regression), wird beim RANSAC-Verfahren eine minimale Teilmenge („Sample“)
der Daten gewählt, mit der die Parameter des Modells geschätzt werden können.
Um das für diese Arbeit relevante Beispiel zu nennen: Soll an eine Menge von dreidi-
mensionalen Datenpunkten eine Ebene angepasst werden, so benötigt man als mini-
male Auswahl drei nicht-kollineare Punkte, um die Ebenenparameter zu schätzen (Ebe-
ne durch diese drei Punkte). Wenn man davon ausgehen kann, dass die Menge haupt-
sächlich Punkte enthält, die auf einer Ebene liegen (Inlier) und andere Punkte, die da-
von abweichen (Outlier), so sieht man direkt ein, dass die Bestimmung einer Regressi-
onsebene für alle Punkte nicht zum gewünschten Ergebnis führt. RANSAC schätzt die
Ebenenparameter nur anhand der zuvor bestimmten Inlier, unter der Annahme, dass
die Wahrscheinlichkeit, bei zufälliger Entnahme von drei Datenpunkten nur Inlier anzu-
treffen, hinreichend hoch ist.
Zur Bestimmung der Ebenenparameter wird also eine zufällige Auswahl von drei nicht
kollinearen Punkten (pi, pj, und pk) aus der Datenmenge gewählt. Der Normalenvektor
n0 der resultierenden Ebene ist dann gegeben durch m = (pi – pj)×(pi – pk), n0 = m/|m|,
und mit (x–pi)·n0 = 0 ist die Hessesche Normalform der Ebene gegeben. Damit ist es
recht einfach, alle übrigen Datenpunkte hinsichtlich ihres Abstands von dieser Ebene zu
überprüfen, denn für jeden Punkt p ist der Abstand zur Ebene mit d = |(p–pi)·n0| leicht
zu berechnen. Abhängig von ihrem Abstand d werden nun alle Punkte der Datenmenge
als Inlier oder als Outlier bewertet. Die Prozedur der zufälligen Auswahl dreier Punkte
wird mehrmals wiederholt, wobei letztlich die umfangreichste sich ergebende Inlier-
Menge für die Bestimmung einer Regressionsebene verwendet wird.
2.2 Iterative closest point (ICP)
Wenn wir im Folgenden zwar nicht das eigentliche ICP-Verfahren verwenden, baut un-
ser Ansatz doch darauf auf. Daher werden kurz die Grundideen des ICP-Verfahrens
angesprochen. Der iterative-closest-point Algorithmus von Besl und McKay zielt darauf
ab, eine genaue und effiziente automatische Registrierung von 3D-Objekten zu errei-
chen, seien es Punktmengen, Linienstücke, Kurvensegmente oder parametrisierte
Oberflächen. In unserer Arbeit betrachten wir aufgrund der Beschaffenheit der ALS-
Datensätze nur den Fall der Koregistrierung zweier Punktmengen. In der Regel wird
einer der Datensätze als „Daten“ bezeichnet, der andere als „Modell“. Der Einfachheit
halber bleiben wir im Folgenden bei dieser Terminologie. Es wird angenommen, dass
Daten und Modell sich bereits näherungsweise in Position befinden, was für unsere La-
ser-Messwerte der Fall ist.
Während des ICP-Ablaufs wird der gesamte Datensatz D iterativ zum Modell M hin be-
wegt, sodass beide abschließend bestmöglich zueinander ausgerichtet sind. Jeder ICP-
Iterationsschritt ist in zwei Schritte unterteilt, die nun beschrieben werden. Für die hier
vorliegende Arbeit ist vor allem Schritt 2 relevant.
1. Bilden von Punkt-zu-Punkt Korrespondenzen zwischen D und M. Das klassische
ICP-Vorgehen ist hierbei, einfach die nächstliegenden Punkte einander zuzuord-
nen, daher auch der Name des Algorithmus. Bei unseren unregelmäßig verteilten
3D-Punkten tritt hier insbesondere das Problem auf, Suchoperationen effizient zu
gestalten, z.B. durch ein TIN oder eine Octree-Datenstruktur. Wir umgehen das
Problem der Korrespondenzsuche aber bereits dadurch, dass wir nicht die Origi-
nalpunkte koregistrieren, sondern Schnittpunkte jeweils dreier Ebenen, von de-
nen die jeweils korrespondierenden Ebenen im anderen Datensatz vorab be-
kannt sind, weswegen auch die zugehörigen Schnittpunkte einander zugeordnet
werden können.
2. Bestimmung einer Transformation, bestehend aus Translation und Drehung, die
die einander zugeordneten Punkte bestmöglich in Einklang bringt, im Sinne einer
Minimierung der euklidischen Abstände. Auch hier kann ein RANSAC-Ansatz
sinnvoll sein, wenn die Punkt-zu-Punkt Zuordnungen noch Ausreißer enthalten.
Das Problem der Transformationsbestimmung kann explizit wie im Originalartikel
[2] durch den Einsatz von Quaternionen gelöst werden, für die Implementierung
ist allerdings eine Singulärwertzerlegung geschickter, wie sie in [1] beschrieben
wird. Wenn M‘ = {mi | i=1,…,n} die Teilmenge aller Modellpunkte ist, die Daten-
punkten in D‘ = {di | i=1,…,n} zugeordnet sind, und zwar so, dass mi mit di kor-
respondiert, so werden zunächst die Schwerpunkte dieser Punktmengen be-
rechnet:
1 1
1 1,
n n
m i d i
i in nc m c d (1)
M‘ und D‘ werden anschließend mit ihren Schwerpunkten cm und cd auf den Ur-
sprung des Koordinatensystems bewegt:
| , 1,..., ,
| , 1,...,
i i i m
i i i d
M i n
D i n
m m m c
d d d c (2)
Danach wird eine 3x3 Matrix H definiert als
1
nT
i
Hi i
d m (3)
Die Singulärwertzerlegung dieser Matrix H = UAVT führt zur optimalen Drehung R
und Translation t
,T
m dR VU Rt c c (4)
Der Beweis hierfür wird in [1] gegeben. Nach der Bestimmung von R und t kann
der gesamte Datensatz D nun entsprechend bewegt und gedreht werden.
Da im Standard-ICP-Verfahren korrespondierende Punkte durch nächstliegende Punkte
approximiert werden, müssen dort die Schritte 1 und 2 mehrfach wiederholt werden.
Besl und McKay haben gezeigt, dass ICP dann in ein lokales Minimum der Punktab-
standsfunktion konvergiert, welches bei hinreichend guter Anfangsnäherung auch die
optimale Lösung repräsentiert. In unserem Verfahren werden keine Iterationen durchge-
führt, da von Anfang an nur echte Punktkorrespondenzen gebildet werden. Daher muss
nur einmal obiger Schritt 2 des Verfahrens angewendet werden, um jeweils zwei der
Datensätze aufeinander auszurichten.
3 Verfahrensablauf
3.1 Vorbereitung der Daten für die weitere Verarbeitung
Die Grundidee unseres Ansatzes ist es, zunächst diejenigen Punkte in den Datensätzen
ausfindig zu machen, die am vielversprechendsten für eine erfolgreiche gegenseitige
Koregistrierung sind. Bei der von uns verwendeten Sensorkonfiguration eines luftgetra-
genen Laserscanners in Schrägsicht sind dies am ehesten Messpunkte auf den Dä-
chern von Gebäuden. Der Verfahrensablauf sieht daher zunächst vor, eine Vorklassifi-
zierung der gemessenen 3D-Punkte durchzuführen, und zwar hinsichtlich der Zugehö-
rigkeit zu den Klassen „Fassade“, „Boden“, „Vegetation“ oder „Dach“.
Im ersten Schritt sollen nun alle Datenpunkte der Gebäudefassaden herausgefiltert
werden. Hierzu werden alle Punkte zunächst einem zweidimensionalen, horizontalen
Bodengitter zugeordnet, in dem die Zellgröße in etwa dem durchschnittlichen Punktab-
stand entspricht (bei unseren Daten typischerweise 0.5 m). Trotzdem kann eine Zelle
dieses Arrays auch gar keine oder mehrere 3D-Punkte erfassen, insbesondere dort wo
sich Fassadenpunkte befinden. Zum besseren Verständnis: Bei diesem Vorgehen wird
keine Interpolation oder Projektion der Daten durchgeführt, es handelt sich um eine rei-
ne 2D-Markierung der 3D-Daten, die für die Datensortierung und für Nachbarschafts-
operationen verwendet werden kann.
(a) Flugrichtung: Von Süden nach Norden
(b) Süd-Nord
(c) Nord-Süd
(d) West-Ost
(e) Ost-West
Abb.2: Verteilung der jeweiligen Punktwolke auf ein horizontales Gitter und Vergleich der vier Datensätze
Wenn man nun pro Zelle nur den höchsten vorkommenden Punkt betrachtet, werden
alle Punkte der Fassaden herausgefiltert, aber auch z.B. Äste unter Baumkronen oder
Baumstämme. Man kann das sich ergebende 2D-Gitter mit dem jeweils höchsten pro
Zelle vorkommenden 3D-Punkt nun auch als Bild darstellen. Abbildung 2 zeigt hierzu
ein Beispiel mit gleichzeitiger Verwendung der gemessenen Intensitätswerte der Laser-
pulse zur Grauwertdarstellung der einzelnen Pixel. Die verschiedenen Scanrichtungen
sind in Abbildung 2b-2e anhand der Position der Schatten zu erkennen
3.2 Segmentierung der Bodenfläche
Aufgrund der aktiven Beleuchtung der Szene durch den luftgetragenen Lasersensor
treten die meisten Schatten und Verdeckungen hinter Gebäuden und auf der Bodenflä-
che auf. Als nächster Schritt ist daher eine Segmentierung der Bodenpunkte durchzu-
führen. Hierfür analysieren wir die lokalen Histogramme der Höhenwerte aus der zuvor
generierten Datenmatrix, jeweils in einer geeignet großen Umgebung von z.B. 50 Me-
tern um die jeweilige Position. Jedes lokale Höhenhistogramm zeigt eine multimodale
Verteilung, in der das Bodenniveau als erstes deutliches Maximum auftritt. Direkt über
dem zu diesem Maximum gehörenden Höhenwert kann nun ein lokaler Schwellwert
gesetzt werden, der darüber bestimmt, ob Punkte zum Boden gezählt werden oder
nicht. Abbildung 3a zeigt ein Beispiel für solch eine lokale Verteilung der Höhenwerte,
worin der dort gefundene Schwellwert als rote Linie eingezeichnet ist. Abbildung 3b
zeigt die insgesamt für einen der Datensätze verbleibenden Punkte nach dem Herausfil-
tern der Bodenpunkte.
(a)
(b)
Abb. 3: (a) Bestimmung des Schwellwertes für die „Bodenentscheidung“ anhand eines lokalen Höhenhis-togramms, (b) verbleibende Punkte über dem Bodenniveau
3.3 Segmentierung planarer Strukturen (Dachflächen)
Die verbleibenden Datenpunkte, wie sie als Beispiel in Abbildung 3b zu sehen sind, re-
präsentieren nun hauptsächlich noch die Gebäudedächer, aber auch Vegetation wie
etwa Bäume oder Gebüsch. In der in diesem Abschnitt dargestellten Methode geht es
erstens darum, zwischen diesen zwei Klassen zu unterscheiden und zweitens, im Fall
einer gefundenen Dachfläche auch gleichzeitig die Ebenenparameter zu bestimmen. Es
folgt eine prozedurale Formulierung des Algorithmus und anschließend eine umgangs-
sprachliche Beschreibung der Vorgehensweise:
(1) Wähle eine bislang unmarkierte Position (i, j) zufällig unter den besetzten Zel-
len der Datenmatrix aus.
(2) Fasse Punkte in einer Nachbarschaft dieser Position zur Menge S zusammen.
(3) Setze den Zähler k auf Null.
(4) Falls S genügend Punkte enthält (z.B. mindestens sechs), fahre fort mit (5),
andernfalls markiere diese Position als „Vegetation“ und gehe zu (14).
(5) Erhöhe den Zähler k um Eins.
(6) Führe eine RANSAC-basierte Ebenenanpassung an die Punkte in S durch.
(7) Falls der Anteil der gefundenen Inlier gering ist, markiere diese Position als
„Vegetation“ und gehe zu (14).
(8) Bestimme die Hessesche Normalform (x–pi)·n0 = 0 der Ebene und lege die
Position (i, j) auf einen Stapel (Datenstruktur).
(9) Entnehme das oberste Element (u, v) des Stapels.
(10) Falls der Zähler k das vordefinierte Maximum erreicht hat und S genügend
Elemente enthält, markiere alle zu Punkten aus S gehörenden Positionen als
„prozessiert“ und bestimme die Regressionsebene zu S. Speichere den
Schwerpunkt von S und den Normalenvektor der Regressionsebene.
(11) Prüfe Positionen in einer Nachbarschaft von (u, v), die noch nicht diesbezüg-
lich überprüft wurden, ob sie Datenpunkte enthalten, die die gefundene Ebene
stützen. Falls dies so ist, lege die entsprechenden Positionen auf den Stapel
und füge die zugehörigen Datenpunkte einer neuen Menge S‘ hinzu.
(12) Solange der Stapel nicht leer ist, wiederhole ab (9), andernfalls gehe zu (13).
(13) Falls der Zähler k das vordefinierte Maximum erreicht hat (z.B. drei Durchläu-
fe), setze ihn zurück auf Null und fahre fort mit Schritt (14). Ansonsten gehe
mit der neuen Menge S:=S‘ zurück zu Schritt (4).
(14) Wiederhole ab Schritt (1) bis alle Datenpunkte markiert sind.
Verfahren zur lokalen Anpassung von Ebenen an die Punktmengen
Es handelt sich um ein iteratives Verfahren, bei dem in jedem Durchlauf zunächst eine
noch nicht behandelte Position innerhalb der Datenmatrix zufällig gewählt wird (mit „Da-
tenmatrix“ ist das horizontale Bodengitter mit jeweils zugeordnetem höchsten 3D-Punkt
gemeint, siehe Abschnitt 3.1). Nun wird mittels RANSAC-Verfahren (vgl. Abschnitt 2.1)
versucht, an diese Position mit den Datenpunkten aus der unmittelbaren Nachbarschaft
eine Ebene anzupassen. Falls dies nicht gelingt, also die Anzahl der Outlier immer sehr
hoch ist, wird diese Stelle in der Datenmatrix als „Vegetation“ klassifiziert. Anderenfalls
wird mit einem Füllverfahren (seed fill, Schritte (9) bis (12)) nach weiteren Punkten ge-
sucht, die auf der gefundenen Ebene liegen, d.h. die Ebene wird soweit wie möglich
ausgedehnt. Mit den so gefundenen Punkten kann dann die Ebenenanpassung wieder-
holt werden, um das Ergebnis weiter zu verbessern. Abbildung 4a zeigt beispielhaft die
detektierten Dachflächen für einen der Datensätze, wobei die verbleibenden Daten-
punkte nun entsprechend der Normalenrichtung eingefärbt sind. Das Gesamtergebnis
der automatischen Vorklassifizierung ist in Abbildung 4b zu sehen, wobei hier die Far-
ben die Zugehörigkeit zu „Boden“, „Vegetation“ oder „Dach“ kennzeichnen.
(a)
(b)
Abb. 4: (a) Segmentierte Dachflächen, eingefärbt entsprechend Normalenrichtung, (b) Ergebnis der Vor-klassifikation, Boden (blau), Vegetation (grün), Dächer (rot)
Letztlich ist es das Ziel, die gefundenen Dachflächen für eine gegenseitige Registrie-
rung verschiedener Ansichten des gleichen Stadtgebiets zu verwenden. Man könnte
daher das hier beschriebene Verfahren zunächst unabhängig voneinander auf alle vor-
handenen Datensätze anwenden. Danach würde aber keine Information über korres-
pondierende Flächen in den Punktwolken vorliegen. Die vollständige Abarbeitung aller
Punkte wird daher zunächst nur an einem Datensatz durchgeführt, der hierbei als Refe-
renzdatensatz gilt.
Bei den anderen Datensätzen wird wie folgt vorgegangen: Anstelle der zufälligen Aus-
wahl von Positionen der Datenmatrix in Schritt (1) werden nur die Positionen der
Schwerpunkte von bereits gefundenen Flächen des Referenzdatensatzes als Aus-
gangspunkte für die Ebenenanpassung gewählt. Da man davon ausgehen kann, dass
die einzelnen Datensätze bereits grob registriert sind, ist die Wahrscheinlichkeit groß,
dass man hier gerade die jeweils gleiche Ebene im anderen Datensatz wiederfindet.
Um sicher zu gehen, dass so gefundene Ebenen mit den Referenzflächen korrespon-
dieren, wird zusätzlich der Winkel zwischen den Normalenvektoren überprüft und eben-
so der Flächeninhalt, z.B. der konvexen Hülle der gefundenen Punktmenge. Nur wenn
sich hier gute Übereinstimmungen zeigen, wird dieses Paar korrespondierender Ebe-
nen in eine Liste eingetragen.
3.4 Koregistrierung der Punktwolken
Es bestehen nun mit der Vorarbeit aus Abschnitt 3.3 entsprechend der Anzahl der Da-
tensätze mehrere Listen, in denen jeweils korrespondierende Ebenen des Referenzda-
tensatzes und der dazu auszurichtenden Punktmenge eingetragen sind. Es bezeichne
M den Referenzdatensatz und D die jeweils andere Punktmenge, die zu M koregistriert
werden soll. Es seien Ei, Ej und Ek drei Ebenen die zum Datensatz M gehören, und (Êi,
Êj, Êk) die entsprechenden Ebenen aus D.
Wenn k die Länge der Liste ist, also die Anzahl der paarweise korrespondierenden
Ebenen, so gibt es (1/6)·k·(k-1)·(k-2) Möglichkeiten zur Auswahl eines Ebenentripels.
So groß ist folglich auch die Anzahl möglicher Schnittpunkte dieser. Da die Ebenenkor-
respondenzen bekannt sind, können auch die sich ergebenden Schnittpunkte einander
zugeordnet werden. Um zu vermeiden, dass nahezu parallele Ebenen geschnitten wer-
den, wird vor der jeweiligen Berechnung der Schnittpunkte überprüft, ob die Normalen-
vektoren der Ebenen hinreichend voneinander abweichen. Man kann dazu z.B. das Vo-
lumen des durch die Normalenvektoren ni, nj und nk aufgespannten Spats betrachten,
welches gerade die Determinante der durch die drei Vektoren gebildeten Matrix ist, und
fordern:
i j kdet( , , )n n n s (5)
Der Schwellwert s wird dabei am besten auf einen Wert größer als 0.5 gesetzt, um nur
möglichst schiefstehende sich schneidende Ebenen zuzulassen, was für die numeri-
sche Bestimmung des Schnittpunktes günstiger ist. Auch kann gefordert werden, dass
i j k i j k
ˆ ˆ ˆdet( , , ) det( , , ) , 0n n n n n n (6)
mit einem ε nahe bei Null, um so Ausreißer bei den Schnittpunktpaarungen zu vermei-
den, da man so ein weiteres mal sicherstellen kann, nur jeweils zugehörige Ebenen aus
den beiden Datensätzen für die Schnittpunktbestimmung verwendet zu haben.
Mit den sich ergebenden Mengen MS und DS von Schnittpunkten innerhalb Referenzda-
tensatz M und Datensatz D sowie den gleichzeitig bekannten Zuordnungen zwischen
Elementen aus MS und DS kann nun das in Abschnitt 2.2 (Schritt 2) beschriebene Ver-
fahrung zur Koregistrierung der Punktwolken angewendet werden. Man erhält dann ei-
ne Rotation R und eine Translation t, durch die Datensatz D an M ausgerichtet werden
kann. Dieses Vorgehen wird für alle vorhandenen Punktmengen D wiederholt, um so
das Gesamtresultat einer koregistrierten Punktwolke zu erhalten.
In den nachfolgenden Abbildungen zeigen wir ein Beispiel eines Verarbeitungsergeb-
nisses. Zunächst sind in Abbildung 5 Ausschnitte von vier unabhängig voneinander auf-
gezeichneten Datensätzen des gleichen urbanen Gebietes zu sehen, jeweils in Schräg-
sicht aus einer anderen Richtung kommend abgetastet und zur Unterscheidbarkeit vi-
sualisiert in einer jeweils eigenen Farbe. Datensatz 1 dient als Referenzdatensatz. Die
anderen drei Punktwolken werden durch die jeweils berechnete Rotation R und Transla-
tion t zum Referenzdatensatz hinbewegt. Abbildungen 6a und 6b zeigen den Vergleich
aller Datensätze vor bzw. nach der Koregistrierung. In Abbildung 7 ist das Resultat für
das Testgebiet TUM zu sehen.
(a)
Datensatz 1
Referenzdatensatz
(b)
Datensatz 2
0.9999 0.0003 -0.0059
-0.0003 0.9999 -0.0005
0.0059 0.0005 0.9999
R
2.2132 -2.2679 0.0043T
t =
(c)
Datensatz 3
0.9999 0.001 -0.0023
-0.0011 0.9999 -0.004
0.0023 0.004 0.9999
R
1.9227 -0.2162 0.1492T
t =
(d)
Datensatz 4
0.9999 0.0005 0.0025
-0.0005 0.9999 0.0029
-0.0025 -0.0029 0.9999
R
-2.0898 0.2883 0.2324T
t =
Abb. 5: Koregistrierung von vier Multi-Aspekt-Datensätzen des gleichen urbanen Gebiets
(a)
(b)
Abb. 6: Vergleich der gegenseitigen Ausrichtung, (a) vor der Koregistrierung, (b) nach der Koregistrierung
Abb. 7: Koregistrierung der vier Datensätze des Testgebiets TUM. Punkte des gleichen Datensatzes sind gleich eingefärbt.
4 Abschließende Bemerkungen und Ausblick
Das hier vorgestellte Verfahren zur Koregistrierung mehrerer Punktwolken hat den Vor-
teil, dass es nicht direkt von der Abtastung der Szene abhängt, weil es auf einer Zuord-
nung von Ebenenschnittpunkten beruht. Im Vergleich zum Standard-ICP-Verfahren wird
keine Iterationsschleife benötigt, da die Zuordnung der Ebenenschnittpunkte sich unmit-
telbar aus der Ebenenkorrespondenz ergibt. Aufgrund dieser Verbesserungen bedarf es
auch keiner Suchoperation nach nächstliegenden Punkten („closest point“), welche bei
unregelmäßigen Punktmengen auch nur schwer zu realisieren ist. Insgesamt zeigt sich
daher eine klare Verkürzung der Rechenzeit, wobei man natürlich einen Teil wieder bei
der Ebenensegmentierung einbüßt. Eine Vorklassifizierung der Punkte ist aber ohnehin
in den meisten Fällen für die automatische Auswertung erforderlich. Die Laufzeit des
Verfahrens zur Vorklassifizierung und anschließender Koregistrierung zweier Punktwol-
ken (jeweils 750.000 Punkte) betrug auf einem Standard-PC etwa 8 Minuten. Die der-
zeitige Implementierung ist in MATLAB umgesetzt und nicht bezüglich der Laufzeit op-
timiert. Bei den Tests konnten für mehr als 80% zugeordneter Ebenen erreicht werden,
dass die Flächenschwerpunkte nach der Koregistrierung weniger als 15 cm von der je-
weils anderen Ebene abweichen.
Als Nachteil gegenüber Standard-ICP kann die Abhängigkeit vom Vorhandensein pla-
narer Strukturen angeführt werden. In städtischen Gebieten sind jedoch häufig planare
Flächen zu finden.
In der vorliegenden Arbeit wurde das unabhängige Ausrichten auf jeweils nur zwei Da-
tensätze beschränkt und schrittweise gelöst, anstatt das Problem als Ganzes anzuge-
hen. Zukünftig soll die Registrierung aller Datensätze gleichzeitig durch einen Bündel-
blockausgleich geschehen. Da wir vermuten, dass nach der Korrektur der GPS-
Positionen mit Hilfe von SAPOS-Daten der nun noch verbleibende Versatz der
Punktwolken aus der leichten Dejustierung der Sensoranbringung resultiert, gilt es in
zukünftigen Arbeiten die Ergebnisse der Punktmengen-Koregistrierung auf die Korrektur
zur Sensorposition und Orientierung zu übertragen, um so eine Sensorselbstkalibrie-
rung zu realisieren.
Literaturverzeichnis
[1] Arun, K.S.; Huang, T.S.; Blostein, S.D.: Least square fitting of two 3-d point sets.
IEEE Transactions on Pattern Analysis and Machine Intelligence 9 (5), 1987, pp.
698-700.
[2] Besl, P.J; McKay, N.D.: A method for registration of 3-D shapes. IEEE Transacti-
ons on Pattern Analysis and Machine Intelligence, Vol. 14, No. 2, 1992, pp. 239-
256.
[3] Chen, Y.; Medioni, G.: Object Modelling by Registration of Multiple Range Images.
Image and Vision Computing, Vol. 10, No. 3, 1992, pp. 145-155.
[4] Chetverikov, D.; Stepanov, D.; Krsek, P.: Robust Euclidean alignment of 3D point
sets: the trimmed iterative closest point algorithm. Image and Vision Computing 23
(3), 2005, pp. 299-309.
[5] Fischler, M.A.; Bolles, R.C.: Random sample consensus: a paradigm for model
fitting with applications to image analysis and automated cartography. Communi-
cations of the ACM 24 (6), 1981, pp. 381-395.
[6] Hebel M.; Stilla U.: Automatic registration of laser point clouds of urban areas. In-
ternational Archives of Photogrammetry, Remote Sensing and Spatial Geoinfor-
mation Sciences, Vol. 36 (3/W49A), 2007, pp. 13-18.
[7] Hoover, A.; Jean-Baptiste, G.; Jiang, X.; Flynn, P.J.; Bunke, H.; Goldof, D.B.,
Bowyer, K.; Eggert, D.W.; Fitzgibbon, A.; Fisher, R.B.: An Experimental Compari-
son of Range Image Segmentation Algorithms. IEEE Transactions on Pattern
Analysis and Machine Intelligence 18 (7), 1996, pp. 673-689.
[8] Jutzi, B.; Stilla, U.: Precise range estimation on known surfaces by analysis of full-
waveform laser. Photogrammetric Computer Vision PCV 2006. International Arc-
hives of Photogrammetry and Remote Sensing Vol. 36 (Part 3), 2006.
[9] Kapoutsis, C.A. ; Vavoulidis, C.P. ; Pitas, I. : Morphological techniques in the itera-
tive closest point algorithm. Proceedings of the International Conference on Image
Processing ICIP 1998, Vol. 1 (4-7), pp. 808-812.
[10] Maas, H.-G.: Least-Squares Matching with Airborne Laserscanning Data in a TIN
Structure. International Archives of Photogrammetry and Remote Sensing 33 (3a),
2000, pp. 548-555.
[11] Makadia, A.; Patterson A.; Daniilidis K.: Fully Automatic Registration of 3D Point
Clouds. Proceedings of the IEEE Computer Society Conference on Computer Vi-
sion and Pattern Recognition CVPR 2006, Vol. 1, pp. 1297-1304.
[12] Pottmann, H.; Leopoldseder, S.; Hofer, M.: Registration without ICP. Computer
Vision and Image Understanding 95 (1), 2004, pp. 54-71.
[13] Rabbani, T.; Dijkmann, S.; van den Heuvel, F.; Vosselman, G.: An integrated ap-
proach for modelling and global registration of point clouds. ISPRS Journal of Pho-
togrammetry and Remote Sensing 61 (6), 2007, pp. 355-370.
[14] Rusinkiewicz, S.; Levoy, M.: Efficient Variants of the ICP Algorithm. Proceedings
of 3D Digital Imaging and Modeling 2001, IEEE Computer Society Press, 2001,
pp. 145-152.
[15] Schnabel, R.; Wahl, R.; Klein, R.: Shape Detection in Point Clouds. Technical re-
port No. CG-2006-2, Universität Bonn, ISSN 1610-8892, 2006.
[16] The Digital Michelangelo Project: http://graphics.stanford.edu/projects/mich/
[17] Vosselman, G.; Gorte, B.G.H.; Sithole, G.; Rabbani, T.: Recognising structure in
laser scanner point clouds. International Archives of Photogrammetry, Remote
Sensing and Spatial Information Sciences 46 (8), 2004, pp. 33–38.
[18] Wehr, A.; Lohr, U.: Airborne Laser Scanning - an Introduction and Overview.
ISPRS Journal of Photogrammetry and Remote Sensing, 54, 1999, pp. 68-82.
[19] Zhang, Z.: Iterative Point Matching for Registration of Free-Form Curves and Sur-
faces. International Journal of Computer Vision 13 (2), 1994, pp. 119-152.