Volumenvisualisierung von Georadardaten

Software - Praktikum an der Universität Heidelberg im WS 2007/2008

Jérôme Arnoldy, Jens Fangerau

 

 

Inhalt

 

 

 

Einleitung

Das Praktikum "Volumenvisualisierung von Georadardaten" basiert auf der Idee, aus 2d Daten ein Volumen zu erstellen, indem die Anzahl der Daten die Größe der Z-Dimension bestimmt. Da man aber häufig zu wenig Daten hat und somit auch eine zu kleine Dimension, legt man einen beliebigen Abstand zwischen den Datenscheiben fest und füllt diesen mit geeigneter Interpolation.

Mithilfe des Volumenvisualisierers GVvis ("G" für Georadar) ist die Interpolation und die Darstellung des Volumens möglich und erlaubt das interaktive Arbeiten auf einer benutzerfreundlichen Oberfläche, die mit Qt entwickelt wurde. GVvis wurde hauptsächlich aus dem Grund entwickelt, um Daten aus dem sog. Georadarverfahren zu nutzen, das im Folgenden nun näher erläutert wird.

 

 

Georadarverfahren

Das Georadar ist ein Verfahren der Gesellschaft für Geophysikalische Untersuchungen (http://www.ggukarlsruhe.de) zur Bestimmung von Strukturen oder zur Objektdetektion im Untergrund und im Bauwerk aufgrund elektromagnetischer Wellen.

 

Es werden in sehr dichter Reihenfolge elektromagnetische Wellen aus dem Sender S in ein Medium gesendet, meist in Form von Impulsen mit Dominanzfrequenzen. Der Sender, der im Normalfall neben dem Empfänger ist, befindet sich auf der sog. Geländeoberkante (GOK).

Die Wellenausbreitung hängt von den Materialeigenschaften des Mediums ab. Wesentlich hierfür sind die Dielektrizität Epsilon und die Leitfähigkeit Sigma. Sobald eine Materialänderung auftritt, wird das Signal als Echo zurück zum Empfänger E gesendet.

Die Objektoberkante d berechnet sich aus der Wellengeschwindigkeit v und der Laufzeit t des Signals. Da jedoch die Wellengeschwindigkeit im Medium nie wirklich konstant ist, muss die Bodenbeschaffenheit bekannt sein bzw. angenähert werden, damit eine Umrechnung von Laufzeit zur Tiefe möglich ist. Entscheidend für die Radarreichweite ist die Höhe der Leitfähigkeit, die die Absorption des Signals bestimmt. Registriert wird die Signalamplitude in Abhängigkeit der Laufzeit als sog. Signalspur.

 

 

Durch Aneinanderreihung der Signalspuren erhält man ein Diagramm, in dem die Entfernung entlang der Meßlinie über der Laufzeit aufgetragen ist.

In diesem sog. Radargramm zeigen sich die reflektierenden Strukturen als mehr oder weniger starke Signale.

Die Abbildung rechts zeigt eine Gegenüberstellung von einem Radargramm und einer Untergrundsituation im selben Maßstab.

Das Primärsignal ist mit A gekennzeichnet. B,C und D zeigen die Reflexionsstruktur mit den dazugehörigen Schichtgrenzen b,c und d. E zeigt eine Diffraktion, also eine Brechung an einem Objekt e an.

 

 

Bei flächigen Untersuchungen lassen sich die Messungen rasterförmig anlegen. Wenn man die Radargramme aneinanderreiht, läßt sich die Signalamplitudendarstellung innerhalb der gesamten Meßfläche für eine bestimmte Laufzeit als sog. Zeitscheibe erfassen .

Da sich diese Zeitscheiben einer bestimmten Tiefe zuordnen lassen, entsprechen sie einer grundrißähnlichen Darstellung des Untergrunds mit den entsprechenden Reflexionsstärken.

Nun werden aber die Zeitscheiben aufgrund der Messungenauigkeiten innerhalb eines bestimmten Tiefenbereichs auf eine einzelne Zeitscheibe abgebildet, deren Messung sich aus der Mittelung der Werte der anderen Scheiben ergibt. Letztendlich hat man dann eventuell nur 6 Zeitscheiben, die insgesamt 1.5 m umfassen, wobei jede einen Bereich von 0.25 m abdeckt.

 

Jedoch gehen dadurch Information verloren, die gerade bei der Volumenvisualisierung wichtig sind. Aus diesem Grund versucht man in GVvis mittels geschickter Interpolation den Verlust wieder gut zu machen und ein möglichst realistisches Volumen zu erstellen. Als Arbeitsdaten dienen die Georadarmessungen vom Kloster Lorsch, dessen Georadaruntersuchung vom Januar bis April 2006 stattfand.

Die Abbildung unten zeigt die Messung vom Friedhofsgelände des Klosters. Die roten Linien stellen die Wegmessung entlang der definierten y - Achse dar.

 

Unten sind die dazugehörigen Zeitscheiben mit den entsprechenden Laufzeiten und Tiefen abgebildet. Zu beachten ist, das es sich hier um eine Fläche von 75 m in
X - Richtung und 45 m in Y - Richtung handelt, jedoch die Untersuchung in die Tiefe "nur" 1.5 m reicht. Man sollte deshalb bei der Volumenvisualisierung beachten, das man die Größe der Abstände zwischen den Scheiben geeignet wählt, damit die Z - Dimension insgesamt nicht zu klein ausfällt.

 

Laufzeit in ns 0 - 5   5 - 10   10 - 15  
Tiefe in m 0.0 - 0.25   0.25 - 0.5   0.5 - 0.75  

Laufzeit in ns 15 - 20   20 - 25   25 - 30  
Tiefe in m 0.75 - 1.0   1.0 - 1.25   1.25 - 1.5  

 

Download der Zeitscheiben des Friedhofs:

256 x 128 ~ 399 KB
Original 1318 x 805 ~ 9.38 MB => Bei der Auflösung kam es zu Problemen beim Import in GVvis, die bis jetzt noch nicht behoben wurden.

 

Volumenvisualisierer GVvis

Da GVvis nicht der erste Volumenvisualisierer im Rahmen eines Praktikums ist und damit man grundlegende Dinge nicht neu programmieren musste
(wie Renderer, Transferfunktion etc.), wurde der Code vom Volumenvisualisierer Vvis (http://vvis.sourceforge.net) bzw. dessen Erweiterung QtVvis
(http://pille.iwr.uni-heidelberg.de/~volumen02) genutzt. Um was es sich genau bei Qt handelt, wird im Abschnitt GUI mit Qt erläutert.


Vvis wurde im Rahmen eines Software - Praktikums im Sommersemester 2005 entwickelt. Mit ihm ist es möglich, aus einer Menge von TIFF - Bildern ein Volumen zu erzeugen und dieses im Vrend oder OpenQvis Format abzuspeichern. Vrend und OpenQvis sind selbst Volumenvisualisierer, die ein eigenes Fileformat benutzen, welches auch kompatibel zu GVvis ist:

Um GVvis kompilieren zu können, müssen folgende Voraussetzungen unter Windows erfüllt sein:

  • OpenGl

http://www.opengl.org

  • mindestens Qt Version 4.3.1

http://trolltech.com

  • libtiff

http://libtiff.org

Momentan existiert nur eine Windows Version von GVvis die man hier runterladen kann (~4.4 MB), aber Linux und Mac Versionen sind geplant. Es sei auch erwähnt, das dies noch keine finale Version ist. Es sind noch diverse Erweiterungen geplant und eventuell stößt man bei Benutzung noch auf kleine Fehler.

Im Folgenden werden nun die Funktionen und Erweiterungen von GVvis erläutert.

Zum einen ist die volle Funktionalität von QtVvis erhalten geblieben, d. h. jegliche Operation, die in QtVvis gemacht werden konnte, ist auch in GVvis möglich. Außerdem wurde das sog. StatusWindow eingeführt, mit dem das interaktive Arbeiten übersichtlicher und benutzerfreundlicher wird.

 

  • Renderer

    Hier kann man aus einer ComboBox zwischen den möglichen Renderen wählen und dazugehörige Eigenschaften ändern. Durch Drücken des "Go!" Buttons wird die Darstellung gerendert. Da die Renderer und deren Eigenschaftenanpassung aus dem Code von Vvis übernommen wurde, soll hier auf die Dokumentation von QtVvis hingewiesen werden, in der die Renderer und Einstellungen unter anderem im Detail erläutert wird.

     

  • Interpolation

    Hier läßt sich die Interpolationsart für das gegebene Volumen aus einer ComboBox auswählen. Bei jeder Methode kann man den Abstand zwischen den Scheiben bestimmen, die dann letztendlich die Z - Dimension des Volumens festlegen. Die Interpolationsmethoden werden weiter unten im Detail beschrieben.

  • Terrain

    Der "Edit Terrain" Button öffnet ein neues Fenster, in dem es dann möglich sein soll, das Terrain entlang der x - und y -Achse zu verändern, um eventuelle Unebenheiten, die in der Realität im Gelände auftreten können, nachbauen zu können. Jedoch ist diese Funktion noch in der Entwicklung und daher noch nicht eingebaut.

  • Information

    • Das Informationsfenster gibt Auskunft über die wichtigsten Eigenschaften des aktiven Volumens und wird bei jeder Änderung dieser Werte aktualisiert.

  • Brightness

    • Hier läßt sich die Helligkeit des Volumens mit Hilfe eines Sliders und der LED - Anzeige bestimmen.

  • Screenshot

    • Mit diesem Button läßt sich jederzeit ein Screenshot vom Volumen in seiner jetzigen Position und Darstellung machen, das im TIFF - Format dann abgespeichert werden kann.

  • Transfer Function

    • Durch Drücken dieses Buttons wird die Transferfunktion in einem neuen Fenster geöffnet. Da diese aber auch aus Vvis bzw. QtVvis übernommen wurde, wird hier wieder für nähere Erklärung auf die Dokumentation von QtVvis verwiesen.

     

Nun zum Hauptfenster von GVvis

     

Die Menütitel "File", "Edit", "Renderer", "Display" und "Help" waren auch schon in Vvis bzw. QtVvis vertreten, jedoch wurden sie um ein paar weitere Einträge erweitert. Es werden im Folgenden nur die neu dazugekommenden Einträge erklärt. Die anderen sind wieder in der Dokumentation von QtVvis zu finden.

  • File
    • Save As

    Hiermit besteht die Möglichkeit, das aktuelle Volumen im Vrend oder OpenQvis Format abzuspeichern.

     

    • Import

    Ermöglicht den Import von TIFF - Files und anschließende Interpolation

    Mit "Start at number" und "End at number" legt man die Anfangs- und Endzahl der Tiff - Files fest (z.B. Test1.tif - Test6.tif). "Number of digits" legt die Anzahl der Ziffern fest, die auftreten können. Nun kann man außerdem den Abstand zwischen den Scheiben festlegen und die Interpolationsmethode wählen. Falls die Hermite - Interpolation ausgewählt wurde, werden auch die SpinBoxes "Bias" und "Tension" freigeschaltet. Damit läßt sich die Kurve beeinflussen. Mehr dazu weiter unten.

     

    • Convert

    Das Konvertieren aus Vvis wurde hier um den Befehl "Convert Geo Tiff files" erweitert. Der Unterschied zum "normalen" Konvertieren besteht darin, das nun die Interpolation verfügbar ist. Das nächste Fenster stimmt mit dem des Imports überein. Weiterhin kann man dann das Format wählen, in das das Volumen gespeichert werden soll.

     

  • Edit
    • StatusWindow

    • Hiermit läßt sich das StatusWindow öffnen und auch wieder schließen. Es ist erst aktiviert, sobald man entweder ein Volumen importiert oder konvertiert hat

    .

    • CoSystem

      Dies ermöglicht das Ein- und Ausblenden eines Koordinatensystems im Falle der Darstellung mit dem Texture - Renderer.

 

Interpolationsmethoden

GVvis bietet vier Interpolationsmöglichkeiten an, um zwischen den Zeitscheiben zu interpolieren. Für die nachfolgenden Bilder gilt jeweils immer ein Abstand zwischen den Scheiben von 8 und die selbe Perspektive, in der die längere Seite der X - Achse und die kürzere der Z - Achse entspricht.

  • Lineare Interpolation

    • Bei der linearen Interpolation wird einfach linear zwischen den Stützstellen (hier also den Zeitscheiben) interpoliert. Unten ist der dazugehörige Pseudo - Code abgebildet.

    VoxelType fz; // unsigned char, which holds the grey information of the voxel
    VoxelType fz0,fz1;
    float fInterFactor;
    unsigned int uiIndex;

    // uiIndex gives the lower z-knot value according to the current position
    // uiZ in [uiIndex, uiIndex + n*getDiscDistance()]
    // n <= number of discs
    uiIndex = (int)((uiZ - uiBorder) / getDiscDistance());

    fz0 = m_pOriginalVoxelField->getVoxel(uiX, uiY, uiIndex);
    fz1 = m_pOriginalVoxelField->getVoxel(uiX, uiY, uiIndex + 1);
    fInterFactor = (float)(uiZ - (uiIndex*getDiscDistance()+uiBorder)) / (float)getDiscDistance();

    // linear interpolation between fz0 and fz1 dependent on the fInterFactor
    fz = (VoxelType)((1.0 - fInterFactor) * (float)fz0 + fInterFactor * (float)fz1);

    return fz;

     

  • Kubische Spline Interpolation

Die kubische Interpolation ist ein globales Interpolationsverfahren, d.h. bei jeder Änderung eines Stützpunktes muss die gesamte Kurve neu berechnet werden. Für die Berechnung werden zunächst sogenannte Momente aus einer tridiagonalen Matrix mittels des Iterationsverfahrens bestimmt. Die Koeffizienten des Polynoms werden aus den Momenten berechnet und ausgegeben wird dann der Wert des Polynoms an der jeweiligen Z - Koordinate zwischen den Scheiben.

// get VoxelType belonging to the knots values
for(int i = 0; i < getNumberOfDiscs(); i++)
{
pfz[i] = m_pOriginalVoxelField->getVoxel(uiX, uiY, i);
pMoments[i] = 0;
}

// natural border conditions
pRightSide[0] = pRightSide[getNumberOfDiscs()-1] = 0;
// fill the right side
for(int k = 1; k < getNumberOfDiscs()-1; k++)
{
float factor1 = (pfz[k + 1] - pfz[k]);
float factor2 = (pfz[k] - pfz[k - 1]);
pRightSide[k] = (3./((float)uih * (float)uih)) * (factor1 - factor2);
}

// start the IterationTechnique and write the solution to pMoments
// iMode == JACOBI or GAUSSSEIDEL
IterationTechnique(MatrixA, pMoments, pRightSide, iMode);

// free pointer
delete[] pRightSide;

// distance between the knots, actual different but equidistant here
unsigned int uih = getDiscDistance();

// coefficients of the polynom
float fDelta,fGamma,fBeta,fAlpha;

// gives the index of the knots according to the current position
unsigned int j = (unsigned int)((uiZ - uiBorder) / getDiscDistance());

// set the coefficients
fDelta = pfz[j];
fGamma = ((pfz[j+1] - pfz[j])/(float)uih) - (uih/6.)*(pMoments[j+1] + 2*pMoments[j]);
fBeta = pMoments[j]/2.;
fAlpha = (pMoments[j+1] - pMoments[j])/(6.*(float)uih);

// return the value of the polynom
return (VoxelType)PolynomFunction(fAlpha, fBeta, fGamma, fDelta, (float)uiZ, (float)(j* getDiscDistance() + uiBorder));
}

// polynom function
float CInterpolation::PolynomFunction(float fAlpha, float fBeta, float fGamma, float fDelta, float fx, float fxi)
{
return fAlpha*(fx - fxi) + fBeta*(fx - fxi) + fGamma*(fx - fxi) + fDelta;
}

Das Iterationsverfahren ist eher für dünnbesetzte Matrizen geeignet, jedoch ist es auch hier tauglich, wobei es aber auch einfacher ginge (und vielleicht schneller). Dieses Verfahren wird auf jeden Fall noch mal überarbeitet werden, da auch die optische Ausgabe nicht ganz das ist was man erwartet (man erkennt gewisse "Treppeneffekte").

  • Hermite Interpolation

Die Hemite Interpolation ist ein lokales Interpolationsverfahren, bei dem vier Stützstellen lokal berücksichtigt werden, wobei aber nur zwischen der 2. und 3. Stützstelle interpoliert wird. Außerdem hat man die Möglichkeit mittels "Bias" und "Tension" die Kurve zu beeinflussen. Bias ändert die Ausrichtung der Kurve an den Stützstellen und Tension bestimmt die Spannung der Kurve, also wie scharf sie sich krümmt.

// we need four knots for interpolation
int uiIndex0;
unsigned int uiIndex1;
unsigned int uiIndex2;
unsigned int uiIndex3;

// interpolation between uiIndex1 and uiIndex2
// get uiIndex1 according to current position
uiIndex1 = (unsigned int)((uiZ - uiBorder) / getDiscDistance());
// tu in [0,1]
float tu = (float)(uiZ - (uiIndex1*getDiscDistance() + uiBorder))/(float)getDiscDistance();

//set the other index values depending on the current uIndex1
uiIndex0 = uiIndex1 - 1;
uiIndex2 = uiIndex1 + 1;
uiIndex3 = uiIndex2 + 1;

// in the first interval we need to generate uiIndex0
if(uiIndex0 < 0)
{
uiIndex0 = 0;
}
// in the last interval we need to generate uiIndex3
if(uiIndex3 > getNumberOfDiscs()-1)
{
uiIndex3 = getNumberOfDiscs()-1;
}

float p0,p1,p2,p3;
float t0,t1,tu2,tu3;
float f0,f1,f2,f3;
float result = 0.0;

p0 = m_pOriginalVoxelField->getVoxel(uiX, uiY, uiIndex0);
p1 = m_pOriginalVoxelField->getVoxel(uiX, uiY, uiIndex1);
p2 = m_pOriginalVoxelField->getVoxel(uiX, uiY, uiIndex2);
p3 = m_pOriginalVoxelField->getVoxel(uiX, uiY, uiIndex3);

// ftension : can be used to tighten up the
// curvature (voxel field) at the known points (Spannung)
//
// fbias : is used to twist the curve (voxel field)
// about the known points (Ausrichtung)
tu2 = tu * tu;
tu3 = tu2 * tu;
t0 = (p1 - p0) * (1 + fbias) * (1 - ftension)/2.;
t0 += (p2 - p1) * (1 - fbias) * (1 - ftension)/2.;
t1 = (p2 - p1) * (1 + fbias) * (1 - ftension)/2.;
t1 += (p3 - p2) * (1 - fbias) * (1 - ftension)/2.;

f0 = 2.*tu3 - 3.*tu2 + 1.;
f1 = tu3 - 2*tu2 + tu;
f2 = tu3 - tu2;
f3 = -2*tu3 + 3*tu2;

result = f0*p1 + f1*t0 + f2*t1 + f3*p2;

return (VoxelType)result;

 

  • BSpline Interpolation

 

Die BSpline Interpolation ist wie die Hemite Interpolation ein lokales Interpolationsverfahren, bei dem auch vier Stützstellen für einen Interpolationsvorgang betrachtet werden, wobei aber nur zwischen der 2. und 3. Stützstelle interpoliert wird. Wichtig ist noch, das das BSpline Verfahren nur die erste und vierte Stützstelle interpoliert und sonst approximiert, also sich den inneren Stützstellen nur annähert.

// we need four knots for interpolation
int uiIndex0;
unsigned int uiIndex1;
unsigned int uiIndex2;
unsigned int uiIndex3;

// interpolation between uiIndex1 and uiIndex2
// get uiIndex1 according to current position
uiIndex1 = (unsigned int)((uiZ - uiBorder) / getDiscDistance());
// t in [0,1]
float t = (float)(uiZ - (uiIndex1*getDiscDistance() + uiBorder))/(float)getDiscDistance();

//set the other index values depending on the current uIndex1
uiIndex0 = uiIndex1 - 1;
uiIndex2 = uiIndex1 + 1;
uiIndex3 = uiIndex2 + 1;

// in the first interval we need to generate uiIndex0
if(uiIndex0 < 0)
{
uiIndex0 = 0;
}
// in the last interval we need to generate uiIndex3
if(uiIndex3 > getNumberOfDiscs()-1)
{
uiIndex3 = getNumberOfDiscs()-1;
}

// calculate the BSpline polynom
float p0,p1,p2,p3,f1,f2,f3,f4;
float result = 0.0;
p0 = m_pOriginalVoxelField->getVoxel(uiX, uiY, uiIndex0);
p1 = m_pOriginalVoxelField->getVoxel(uiX, uiY, uiIndex1);
p2 = m_pOriginalVoxelField->getVoxel(uiX, uiY, uiIndex2);
p3 = m_pOriginalVoxelField->getVoxel(uiX, uiY, uiIndex3);

f1 = t*t*t*(-p0 + 3*p1 - 3*p2 + p3);
f2 = t*t*(3*p0 - 6*p1 + 3*p2);
f3 = t*(-3*p0 + 3*p2);
f4 = p0 + 4*p1 + p2;

result = (f1 + f2 + f3 + f4)/6.;

return (VoxelType)result;

 

GUI mit Qt

QT ist eine Klassenbibliothek zur Programmierung graphischer Benutzeroberflächen (Graphical User Interface), die typischerweise als Erweiterung der Sprache C++ verwendet wird (es gibt aber auch Erweiterungen für Java, C#, Python, Perl, Ruby, Pascal, Ada, PHP).

Versionen gibt es u.a. für die Betriebssysteme Windows, Unix/Linux, Mac OS X, deshalb nennt man QT auch plattformübergreifend. Zu den wichtigsten Basisklassen gehören QtCore und QtGUI, es gibt aber auch Schnittstellen für u.a. Datenbank- oder OpenGL-Anwendungen.

 

Zu den bekanntesten Anwendungen, die QT benutzen, gehören u.a. die Desktopumgebung KDE, der Webbrowser Opera, der virtuelle Globus Google Earth und die Telefoniesoftware Skype.

QT ermöglicht die Nutzung von Features, die über den Funktionsumfang der zugrundeliegenden Sprache C++ hinausgehen. Zu nennen wäre zunächst einmal die Introspektion, die es dem Programmierer erlaubt, Informationen über Variablen und Klassen aus dem Quellcode zur Laufzeit des Programms abzufragen.

 

Des Weiteren führt QT das Signals/Slots Konzept ein, welches auf intuitive Weise die Kommunikation zwischen einzelnen Objekten innerhalb eines Programms ermöglicht und somit die in vielen Programmiersprachen üblichen Rückruffunktionen (Callback Functions) ersetzt. Das ausgesandte Signal eines Objekts kann mehrere Slots aktivieren, und umgekehrt können mehrere Signals denselben Slot auslösen.

Für die Umsetzung besagter Features nutzt QT spezielle Tools, von denen der MOC (Meta Object Compiler) am wichtigsten ist. Dieser sammelt sogenannte Meta Informationen über Klassen eines Programms aus dem ursprünglichen Quelltext und erstellt daraus neuen Code, der dann problemlos vom üblichen C++ Compiler übersetzt werden kann. Anhand diesen Prozesses werden die Konzepte der Introspektion und der Signals/Slots erst in die Realität umgesetzt.

Ein weiteres interessantes Tool stellt das eigenständige Programm QT Designer dar, mit dem man ohne vorherige Programmierkenntnisse eigene Elemente einer GUI, d.h. Hauptfenster, Dialoge und graphische Steuerelemente per Drag&Drop erstellen kann. Die daraus resultierenden UI Dateien werden dann typischerweise mit dem UIC (User Interface Compiler) in C++ Code umgewandelt und können somit in ein bestehendes Programm eingebaut werden.

 

 

Quellen

http://vvis.sourceforge.net