Automatische Koregistrierung von Laserdaten aus...

16
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

Transcript of Automatische Koregistrierung von Laserdaten aus...

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.