Partielle Differentialgleichungen in der Bildverarbeitung

Digitale Bildverarbeitung ist ein sehr vielseitiges Gebiet, für das es neben einer Vielzahl von Anwendungen beinahe ebenso viele Lösungsvorschläge gibt. In früheren Jahren war die digitale Bildverarbeitung eine Domäne von Ingenieuren und Physikern, die mit Hilfe linearer und nichtlinearer Filtertechniken, Spektralanalysen oder einfachen Wahrscheinlichkeitsmodellen Bilddaten analysieren oder spezielle Informationen aus diesen extrahieren konnten. Seit etwa 20 Jahren entwickelt sich die digitale Bildverarbeitung dynamisch zu einem eigenen Forschungsgebiet der Angewandten Mathematik. Abhängig von der Antwort auf die für den Mathematiker fundamentale Frage "Was ist ein (digitales) Bild?" haben sich verschiedene Ansätze zur mathematischen Modellierung und Analyse digitaler Bilddaten als sinnvoll erwiesen. Die wichtigsten basieren auf stochastischen Modellen, Wavelets und auf partiellen Differentialgleichungen¹, so genannten PDEs (Partial Differential Equations), bzw. variationellen Ansätzen. PDE-Ansätze werden heute erfolgreich in vielen Bereichen der Bildverarbeitung, wie etwa der Bildverbesserung, der Segmentierung oder auch zur Berechnung des so genannten optischen Flusses, eingesetzt.

Die numerische Lösung partieller Differentialgleichungen bildet einen Forschungsschwerpunkt innerhalb der Angewandten Mathematik an der Heinrich-Heine-Universität Düsseldorf. Die Notwendigkeit der Entwicklung effizienter Lösungsverfahren kam ursprünglich aus den Anforderungen der Physik und der Mechanik, wie etwa im Bereich der Strömungsmechanik. In den letzten Jahren wurden durch neuere Entwicklungen, beispielsweise in der Biologie, der Medizin, der Finanzmathematik sowie in der mathematischen Bildverarbeitung, die existierenden Lösungsverfahren auf die jeweiligen Anwendungen adaptiert. Hierbei werden die zu untersuchenden Problemstellungen zunächst durch geeignete partielle Differentialgleichungen modelliert und anschließend auf dem durch die Bildauflösung vorgegebenen (meist äquidistanten) Gitter diskret betrachtet.

In diesem Beitrag wird ein Überblick über PDE-basierte Ansätze zur Bildregistrierung und deren Anwendungen in der Medizin, der Biologie und der Meteorologie vorgestellt. Die Untersuchung solcher Ansätze bildet einen weiteren Schwerpunkt innerhalb der Angewandten Mathematik an der Heinrich-Heine-Universität Düsseldorf. Im Vordergrund stehen hierbei die Untersuchung der mathematischen Grundlagen, geeignete Modellierungen der Probleme aus den Anwendungsgebieten, die Entwicklung leistungsfähiger Algorithmen sowie deren effiziente Implementierung.

Das Registrierungsproblem

Unter der Registrierung zweier Bilder versteht man die räumliche Anpassung eines Modellbildes an ein Referenzbild, so dass korrespondierende Bildelemente nach der Registrierung auch räumlich überlagert sind (vgl. Abb. 1). Registrierungsalgorithmen sind leistungsfähige Verfahren zur Anpassung digitaler Bilddaten. Die Registrierung (auch matching, mapping oder fusion genannt) zwei- und dreidimensionaler Bilddatensätze ist ein aktueller Forschungsschwerpunkt in der mathematischen Bildverarbeitung, gehört zu den anspruchvollsten Problemen der digitalen Bildverarbeitung und findet Anwendung in vielen Bereichen der Naturwissenschaften, der Medizin und der Technik. Gerade in der Medizin besteht ein hoher Bedarf an Registrierungsalgorithmen, beispielsweise, wenn

• Bilder eines Objektes zu unterschiedlichen Zeiten,
• Bilder mehrerer unterschiedlicher Objekte, von unterschiedlichen Perspektiven oder
• durch unterschiedliche bildgebende Verfahren

erzeugt werden und korrespondierende Strukturen einander zugeordnet werden sollen. Mathematisch fällt das Registrierungsproblem in die Klasse inverser Probleme. Ziel hierbei ist die Identifikation von Größen (hier: Verschiebungen), auf die der Mediziner etwa durch direkte Messungen keinen Zugang hat, sondern auf die man durch Beobachtung anderer verfügbarer Größen (hier: digitale Bilddatensätze) indirekt schließen muss. Das inverse Problem ist nichtlinear und im Sinne von Hadamard schlecht bzw. inkorrekt gestellt. Die berühmte These von Hadamard² aus dem Jahr 1923 besagt, dass ein Problem korrekt (gut, sachgemäß) gestellt ist, falls

1. zu jedem Datensatz eine Lösung existiert,
2. die Lösung eindeutig bestimmt ist und
3. die Lösung stetig von den Daten abhängt.

Probleme, die eine oder mehrere der Bedingungen (1) bis (3) nicht erfüllen, nennt man inkorrekt oder schlecht gestellt. In diesem Zusammenhang sei betont, dass ein inkorrekt gestelltes Problem im Allgemeinen nicht bedeutet, dass die mathematische Modellierung des physikalischen Problems falsch ist, das heißt schlecht oder nicht exakt vorgenommen wurde, sondern dass die Schlechtgestelltheit eine natürliche Eigenschaft des Problems ist. In diesem Sinne ist das Registrierungsproblem schlecht gestellt, weil im Allgemeinen weder die Existenz, die Eindeutigkeit noch die stetige Abhängigkeit der Lösung von den Daten gewährleistet ist. Durch diese Aufgabenstellung ergibt sich auf natürliche Weise das folgende Optimierungsproblem: Bestimme eine Koordinatentransformation u aus einem Parameterraum X, so dass das durch u modellierte mathematische Problem möglichst gut mit dem beobachteten Verhalten des realen (physikalischen) Systems übereinstimmt.

---

---

Die mathematische Modellierung des Problems erfolgt durch die kontinuierlichen Signale (Bilder) T,R : R, mit dem so genannten Referenzbild R(x), das keine Deformation enthält (vgl. beispielsweise Abb. 1 links), und dem als deformiert angenommenen so genannten Templatebild T(x) (vgl. beispielsweise Abb. 1 rechts). Die Deformation eines Objekts wird durch die Abbildung

beschrieben. Sie transformiert ein Bildelement in eine eindeutig bestimmte Position (x,u). Die gesuchten Verschiebungen werden durch das folgende Optimierungsproblem bestimmt:

Finde u X, so dass D[R,T,u(x)] minimal wird.

Hierbei wird durch das Funktional D ein Maß für die Ungleichheit der Bilder definiert, aus dessen Minimierung auf geeignete Weise auf die gesuchte Transformation u(x) geschlossen werden soll. Die geeignete Wahl von D[R,T,u(x)] hängt wesentlich von der betrachteten Aufgabenstellung ab (vgl. Abschnitt Ähnlichkeitsmaße).

In der Praxis ist jedoch die Erfassung dieses physikalischen Gesamtzustandes nicht möglich. Deshalb werden die Bilder digitalisiert, so dass T^h und R^h Abbildungen aus dem Bildbereich _h in die Menge der möglichen Grauwerte [0,g_max] mit maximalem Grauwert g_max (typisch g_max = 255 bei Grauwertbildern) sind. Für numerische Berechnungen wird die Grauwertquantisierung häufig in ein Intervall [0,1] R übertragen. Das Diskretisierungsgitter _h wird durch das durch ein Gitter überdeckte Einheitsquadrat [0,1]² bzw. den Einheitswürfel [0,1]³ für zwei- bzw. dreidimensionale Probleme ersetzt. Abhängig von der Problemstellung gehören bei der Entwicklung von Registrierungsalgorithmen

• die Definition geeigneter Ähnlichkeitsmaße (Zielfunktionale) D,
• die Bestimmung eines geeigneten Parameterraums X,
• eine sachgerechte Modellierung der Verschiebungsfelder,
• geeignete Einbindung von a priori Informationen sowie
• die robuste und effiziente numerische Berechnung geeigneter Verschiebungsfelder

zu den wichtigsten Aufgaben.

Ähnlichkeitsmaße

Abhängig von der betrachteten Problemstellung unterscheiden sich sinnvolle Maße für die Ungleichheit der Bilder. Bei der Registrierung von Bilddaten, die mit dem selben bildgebenden Verfahren aufgenommen wurden, erhält man durch die Betrachtung des Referenzbildes R(x) und des deformierten Templatebildes T(x) die Zustandsgleichung:

(Z)

Sie führt auf ein weit verbreitetes Ähnlichkeitsmaß für zwei Bilddatensätze (mit vergleichbarer Grauwertinformation)

mit einem Abweichungsmaß g : R R⁺.

Die Summe der Abweichungen aller Bildpunkte ist immer positiv, wenn die Zustandsgleichung (Z) nicht erfüllt ist, und null, wenn die Bilder gleich sind. Hierbei wird das Abweichungsmaß der Fehlerquadrate g : x x² in vielen Fällen bevorzugt. Eine anschauliche Rechtfertigung hierfür beschrieb Gauß:³ Die Bestimmung einer Größe durch eine einem größeren oder kleineren Fehler unterworfene Beobachtung wird nicht unpassend mit einem Glücksspiel verglichen, in welchem man nur verlieren, aber nicht gewinnen kann, wobei also jeder zu befürchtende Fehler einem Verlust entspricht.

Das Risiko eines solchen Spiels wird nach dem wahrscheinlichen Verlust geschätzt, d. h. nach der Summe der Produkte der einzelnen möglichen Verluste in die zugehörigen Wahrscheinlichkeiten.

Welchem Verlust man aber jeden einzelnen Beobachtungsfehler gleichsetzen soll, ist keineswegs an sich klar; hängt doch vielmehr diese Bestimmung zum Teil von unserem Ermessen ab. Den Verlust dem Fehler selbst gleichzusetzen, ist offenbar nicht erlaubt; würden nämlich positive Fehler wie Verluste behandelt, so müssten negative als Gewinne gelten. Die Größe des Verlustes muss vielmehr durch eine solche Funktion des Fehlers ausgedrückt werden, die ihrer Natur nach immer positiv ist. Bei der unendlichen Mannigfaltigkeit derartiger Funktionen scheint die einfachste, welche diese Eigenschaft besitzt, vor den übrigen den Vorzug zu verdienen, und diese ist unstreitig das Quadrat. Dem Mediziner stehen heute unterschiedliche Bildaufnahmegerätetechniken (z. B. T1 und T2 gewichtete MR-Sequenzen) zur Verfügung, die ihm ermöglichen, verschiedenste Informationen über das Körperinnere zu gewinnen. Jedes Aufnahmeverfahren veranschaulicht jedoch immer nur einen Teil der tatsächlich vorhandenen Informationen. Durch das Zusammenführen verschiedener Bilddaten (z. B. Verknüpfung von histologischen und MR-Bildern) wird der Informationsgehalt erhöht. Abhängig von der Art der Bildaufnahme unterscheiden sich die Grauwerte gleicher Strukturen und der Kontrast zu umliegenden Strukturen grundlegend. In diesem Fall ist D^LS kein Maß mehr für die Ähnlichkeit der Datensätze. Ein besser geeignetes Distanzmaß⁴ ist die aus der Informationstheorie⁵ bekannte Mutual Information (oder: gemeinsame Information):

zwischen den Bildern T und R. Hierbei bezeichnet H(I) die mittlere Information eines Bildes I, seine Entropie. Die Entropie ist ein Maß für die Unordnung eines Bildes I und kann (mit der relativen Häufigkeit p) als Summe des Informationsgehalts aller Pixel definiert werden:

Die gemeinsame Entropie der beiden Bilder I und J

(mit der relativen Häufigkeit q) ist ein Maß für den gemeinsamen Informationsgehalt. Die Minimierung der gemeinsamen Entropie führt zu einer Maximierung überlappender Bildregionen und somit zu einer guten Registrierung. Deshalb wird in diesem Fall das Funktional D^MI[R, T, u(x)] maximiert oder -D^MI[R, T, u(x)] minimiert. Weitere geeignete Distanzmaße, wie etwa die Normalised Mutual Information,

werden in Studholme (1997) vorgestellt und analysiert. In den Arbeiten von Wells et al. (1996) und Viola und Wells (1997) wurde die Mutual Information benutzt, um Bilder verschiedener Modalitäten (z. B. MR-CT, CT-PET) global aneinander anzupassen. Das Verfahren kann für die Berechnung so genannter affin linearer Transformationen

(Drehungen und Verschiebungen) verwendet werden und liefert gute Resultate für unterschiedliche Bilddaten. Ein Verfahren zur Berechnung nichtlinearer Verschiebungen wurde von Henn und Witsch (demnächst) vorgestellt.

Modellierung der Verschiebungsvektoren

Bei der Lösung inkorrekt gestellter Probleme kann es schon bei kleinsten Mess- oder Modellfehlern zu Resultaten kommen, die weit von einem exakten Ergebnis entfernt liegen. Dieser Instabilität begegnet man in der Registrierung, wie bei vielen schlecht gestellten Problemen, durch die Methode der Regularisierung. Beispielsweise wird bei dem klassischen Regularisierungsverfahren von Tikhonov⁶ das Problem durch einen Regularisierungs- oder Strafterm G[u] mit geeigneten Eigenschaften ergänzt und im Folgenden eine Minimallösung des regularisierten Problems (dem so genannten Tikhonov-Funktional)

mit einem so genannten Regularisierungsparameter > 0 berechnet. Ist G eine symmetrisch positiv definite Bilinearform B[u,v], dann ergibt sich die notwendige Bedingung (die so genannten Euler-Lagrange Gleichungen) für ein Minimum als Lösung der Variationsgleichung

Das ist die schwache Form der partiellen Differentialgleichung

mit einem Differentialoperator L. Durch die Wahl der Bilinearform B wird die Menge der möglichen Minimallösungen mitbestimmt. Hierdurch werden Lösungen, die gewünschte Eigenschaften besitzen, bevorzugt, andere Lösungen hingegen werden durch den Term bestraft. Häufig sollen die Lösungen sehr "glatt" sein, so dass G[u] auch als Glättungsterm bezeichnet wird. Im Folgenden werden drei Ansätze vorgestellt.

Das Modell der eingespannten elastischen Membran

Diese Vorgehensweise lässt sich durch ein einfaches Membranmodell⁷ veranschaulichen. Hierzu werden die Bilddaten auf einer elastischen Membran dargestellt. Dort, wo die Verschiebungsvektoren bekannt sind, wird die Membran in die entsprechende Richtung verschoben.

Die aus der Glattheitsbedingung resultierenden inneren elastischen Kräfte sorgen dafür, dass sich die Verzerrungen möglichst gleichmäßig über das ganze auf der Membran dargestellte Bild ausbreiten. Die äußeren Kräfte greifen über Federn die Membran an. Die Ausgangslage des Federsystems wird um den berechneten Verschiebungsvektor versetzt. Die inneren Kräfte wirken dabei den äußeren Kräften entgegen und versuchen, die unterschiedlichen Verzerrungen auszugleichen.

Die Membran ist eingespannt, deshalb können die Bildelemente am Rand nicht verschoben werden; d. h., es gilt u(x) = 0 für alle Pixel x am Bildrand. Durch diese so genannten Dirichlet Randbedingungen werden gleichmäßige Bewegungen der Membran bestraft, d. h., Scherungen, Rotationen und Translationen der Bildpunkte sind durch dieses Modell nicht berücksichtigt.

Das Modell der elastischen Membran mit freiem Rand

In manchen Anwendungen ist es sinnvoll, eine gleichmäßige Verschiebung der Bildpunkte in eine Richtung zuzulassen. Dies führt auf das Modell einer (elastischen) Membran mit freiem Rand⁸ mit dem folgenden Strafterm

und den dimensionslosen Konstanten und sowie den linearen Verzerrungen _i,j(u) = 1/2(u_i/x_j + u_j/x_i) für 1 < i,j < n für die Verschiebung u. Die Euler-Lagrange Gleichungen führen in diesem Fall auf eine partielle Differentialgleichung mit Neumann Randbedingungen, die eine Verschiebung der Bilddaten über den Rand hinaus zulassen.

Das Modell der frei schwebenden Platte

Modelle, die gleichmäßige Verschiebungen jeder Art erlauben, lassen sich durch das Modell der frei schwebenden Platte veranschaulichen. Die Vorstellung hierbei ist, dass das Templatebild T auf einer frei durch den Raum schwebenden Platte dargestellt ist, die ohne jede Randeinwirkung auf das Referenzbild geführt wird. Das Modell führt auf einen Glättungsterm der Form

In diesem Modell führen die Euler-Lagrange Gleichungen auf eine partielle Differentialgleichung vierter Ordnung

und Randbedingungen höherer Ordnung:

für alle Bildpunkte (x,y) mit Tangentialrichtung s und den gegebenen Funktionen g₁ bzw. g₂.

Unstetige Verschiebungsvektoren

Ein aktuelles Forschungsziel ist die Berechnung von Verschiebungsfeldern, die gewisse Unstetigkeiten zulassen. Die oben vorgestellten Registrierungsalgorithmen modellieren die gesuchten Verschiebungen so, dass Glattheitsanforderungen erfüllt sind. Dies geschieht, indem man Nebenbedingungen an die Verschiebungen stellt und damit ungewünschtes Verhalten (wie etwa zu große Ableitungen, Verzerrungen oder eine zu große Norm der Verschiebungen) ausschließt. In der Praxis treten sehr häufig Unstetigkeiten in den Verschiebungsfeldern auf, wie zum Beispiel bei den folgenden Fragestellungen.

Wichtige Größen bei der kurz- und mittelfristigen Wettervorhersage sind die Windrichtung und die Windgeschwindigkeit. Durch die Assimilation von zu verschiedenen Zeitpunkten von einem Satelliten gesendeten Wolkenbildern kann auf die Verschiebung der Wolken und damit auch auf die gesuchten Winddaten geschlossen werden. In Zusammenarbeit mit Dr. Holmlund von der Firma Eumetsat (Europe’s Meteorological Satellite Organisation) in Darmstadt wurden die vorhandenen Algorithmen auf das beschriebene Problem übertragen. Eine direkte Übertragung der vorhandenen Algorithmen war hierbei, bedingt durch die unterschiedlichen Problemstellungen, nicht möglich.

Zwar sind die Windgeschwindigkeiten sicherlich stetig, jedoch sind bei Satellitenbildsequenzen, bedingt durch die Projektion auf die zweidimensionale Bildebene, Phänomene wie

• mehrfache Bildobjekte (Wolken), die sich in unterschiedlichen Höhen oft auch in unterschiedliche Richtungen fortbewegen,
• teilweise verdeckte Bildobjekte (Wolken) oder
• Bildobjekte, die sich über einen ähnlichen Hintergrund bewegen (z. B. eine weiße Wolke, die sich über einen schneebedeckten Berggipfel bewegt),

zu beobachten, die sich nur durch unstetige Verschiebungen modellieren lassen.

Beim so genannten Tumor-Screening wird in der Medizin über einen längeren Zeitraum ein Tumor beobachtet, indem er in regelmäßigen Abständen aufgenommen wird. Durch das Wachstum des Tumors kann es hierdurch zu erheblichen Unstetigkeiten innerhalb dieser Bildsequenzen kommen. Denkt man beispielsweise zur Operationsplanung an die Anpassung eines MR-Bildes mit Hirntumor an ein Referenzgehirn, so sind die resultierenden Verschiebungen im Tumorbereich unstetig. Zudem ändert das Gehirn während einer Operation seine Form und Größe, etwa durch den Verlust von Flüssigkeit nach der Öffnung der Schädeldecke (so genanntes Brainshift).

Beim computergestützten Operieren verwenden Neurochirurgen zur Navigation ein vor der Operation entstandenes Gehirnbild (das präoperative Bild). Die durch die Registrierung berechneten Verschiebungen sind unstetig. Zur Modellierung solcher Verschiebungsfelder wird ein Modellierungsterm

mit einer Funktion : R⁺ R⁺ für die gesuchte Verschiebung u(x) benutzt. Beispielsweise erhält man für (x) = x die Seminorm des Raumes der Funktionen mit beschränkter Total-Variation und für die Funktion (x) = x² die Norm des Sobolev-Raumes H¹.

Numerische Berechnung der Verschiebungsvektoren

Eine wichtige Aufgabe der Angewandten Mathematik besteht in der stabilen und effizienten Berechnung der Verschiebungsvektoren. Dies führt, wie oben dargestellt, auf die Minimierung der Ähnlichkeitsfunktionale. Obwohl durch die Einführung eines Regularisierungsterms das Registrierungsproblem "gutartiger" wird, müssen bei der Wahl des Minimierungsverfahrens mehrere Aspekte berücksichtigt werden. Zum einen sollen die iterativ berechneten Näherungslösungen gegen eine geeignete Lösung des Registrierungsproblems konvergieren. Andererseits spielen praktische Belange eine erhebliche Rolle. Die Anzahl der Daten ist in praktischen Anwendungen immens. Hierdurch wird der Einsatz schneller Algorithmen erforderlich. Wichtige Hilfsmittel hierzu sind

• Optimierungsverfahren und
• die effiziente numerische Lösung partieller Differentialgleichungen.

Optimierungsverfahren

Die berechnete Minimallösung des Tikhonov-Funktionals hängt stark vom Regularisierungsparameter ab. Es zeigt sich, dass für große eine Minimallösung numerisch stabil berechnet werden kann. Die berechneten Verschiebungen bewirken jedoch kaum Verschiebung im Templatebild T. Anderseits wird für kleine Werte von die numerische Lösung instabil, weil der Einfluss des Regularisierungsterms im Tikhonov-Funktional abgeschwächt wird. Deshalb wird eine Folge von Minimallösungen⁹

betrachtet. Die Iteration startet mit einem relativ großen Wert ₀, so dass zunächst sehr glatte Lösungen stabil berechnet werden können und im weiteren Minimierungsprozess als Startnäherung eingehen. Im Laufe der Iteration wird durch die Wahl von _k+1 =

_k der Einfluss des Glättungsterms immer weiter reduziert. Diese Vorgehensweise wird als Iterative Tikhonov-Regularisierung bezeichnet. Andere Ansätze, wie die Landweber-Iteration¹⁰, basieren auf der Linearisierung des Ähnlichkeitsfunktionals D. Eine Interpretation dieser Methoden als Gradienten-Fluss-Verfahren findet man in der Arbeit von Clarenz et al. (2002). Einen weitreichenden Überblick über moderne Optimierungsverfahren und deren Anwendungen findet man im Buch von Jarre und Stoer (2003).

Die Mehrgitteridee

Durch die Diskretisierung der partiellen Differentialgleichungen entstehen abhängig vom Optimierungsverfahren sehr große, dünn besetzte lineare oder nichtlineare Gleichungssysteme, die erst durch moderne Hochleistungsrechner behandelt werden können. Typischerweise enthalten in der Medizin zweidimensionale Schnittbilder bis zu 512 × 512 Bildpunkte und dreidimensionale Volumendatensätze 256 × 256 × 128 Bildpunkte. In der Biologie werden zweidimensionale Protein-Gelbilder mit bis zu 1024 × 1024 Bildpunkten aufgelöst, und in der Meteorologie enthalten Aufnahmen der Erdoberfläche bis zu 4096 × 4096 Bildpunkte.

Klassische direkte Verfahren zur Lösung solcher Probleme, wie etwa die Gauß-Elimination, können wegen zu großer Speicherplatzanforderungen und des damit verbundenen Rechenaufwandes nicht verwendet werden. Traditionelle iterative Methoden, wie etwa das Gauß-Seidel-Verfahren oder die effizientere SOR-Iteration, konvergieren in den ersten Iterationsschritten sehr schnell, jedoch reduziert sich im Laufe der Iteration die Konvergenzgeschwindigkeit deutlich, so dass sich eine schlechte asymptotische Konvergenzgeschwindigkeit einstellt. Das führt dazu, dass die Zahl der Iterationen, die nötig sind, um eine Lösung mit vorgegebener Genauigkeit zu bestimmen, mit der Anzahl der Bildpunkte steigt. Das macht solche Iterationsverfahren für die oben aufgeführten Probleme unpraktikabel.

Zerlegt man den Startfehler in seine Frequenzanteile, so stellt man fest, dass durch die Iterationsverfahren nur dessen hochfrequente Anteile vernünftig schnell geglättet, seine niedrigfrequenten Anteile jedoch nur sehr wenig verkleinert werden.

---

---

Der Grundgedanke so genannter Mehrgitterverfahren ist einfach. Sie basieren darauf, dass die niedrigfrequenten Anteile des Fehlers fast ohne Informationsverlust auch auf gröberen Gitterebenen repräsentiert werden können, hochfrequente Anteile jedoch nicht. Das Mehrgitterprinzip beruht auf der Kombination aus traditionellen Iterationsverfahren auf feinen und groben Gittern sowie Gittertransferoperatoren auf einer Folge immer gröber werdender Gitter (vgl. Abb. 2).

Multiskalenansatz

Bei genauer Analyse der iterativen Tikhonov-Regularisierung zur Bildregistrierung erkennt man, dass hierdurch kleine (hochfrequente) Bildstrukturen bevorzugt angepasst und große vernachlässigt werden. Das führt dazu, dass die Verfahren in vielen Fällen unbrauchbare lokale Minimalstellen berechnen. Zudem benötigen die Verfahren hierfür sehr viele Iterationsschritte, was gerade bei Bildern mit sehr feiner Auflösung sehr lange dauert.

Analog zu der oben beschriebenen Eigenschaft des Fehlers bei der iterativen Lösung von partiellen Differentialgleichungen besagt das Abtasttheorem der digitalen Bildverarbeitung, dass eine periodische Bildstruktur aus den Abtastwerten nur dann richtig rekonstruiert werden kann, wenn die kleinste Wellenlänge mindestens zweimal abgetastet wird. Das bedeutet, dass bei Bildern mit einer Abtastschrittweite h nur Bildstrukturen mit einer Größe von mindestens 2h präsent sind. Diese Tatsache wird bei der Bildregistrierung ausgenutzt, indem die Bilddaten zunächst auf gröbere Auflösungen transformiert werden. Hier sind als Folge des Abtasttheorems große Strukturen der Originalbilder hochfrequent und können mit wenigen Schritten eines der beschriebenen Optimierungsverfahren angepasst werden. Die anschließende Überführung der Verschiebungen auf die nächstfeinere Gitterebene sorgt für eine gute Startnäherung des Registrierungsproblems.

Beispiele

Wie leistungsfähig die präsentierten Registrierungsverfahren sind, wird an dem nachfolgenden synthetischen Beispiel deutlich.

---

---

Das Quadrat (Abb. 3 links oben) soll so deformiert werden, dass die Grauwertdifferenz D^LS aller Bildpunkte mit dem Kreis (Abb. 3 rechts unten) möglichst minimal wird.

---

---

Affin lineare Transformationen führen nicht zu einer Verringerung von D^LS, so dass das Modell der eingespannten elastischen Membran sinnvoll ist. Nach wenigen Iterationen eines linearisierten Verfahrens ist das deformierte Templatebild nicht mehr von dem Referenzbild zu unterscheiden. Abbildung 4 zeigt die resultierenden Verschiebungsvektoren der wichtigsten Bereiche des in Abbildung 1 gestellten Registrierungsproblems. Die Vektoren sind der Deutlichkeit halber dem Templatebild überlagert.

Bibliographie

CATTERMOLE, Kenneth. Statistische Analyse und Struktur von Information. New York 1988.

CLARENZ, Ulrich, Stefan HENN, Martin RUMPF und Kristian WITSCH. "Relations between optimization and gradient flow methods with application to image registration", Proceedings of the 18^th GAMM-Seminar Leipzig (2002), 11-30.

HADAMARD, Jacques. Lectures on the Cauchy Problem in Linear Partial Differential Equations. New York 1923.

HENN, Stefan, Thorsten SCHORMANN, Knut ENGLER, Karl ZILLES und Kristian WITSCH. "Elastische Anpassung in der digitalen Bildverarbeitung auf mehreren Auflösungsstufen mit Hilfe von Mehrgitterverfahren", Informatik aktuell, Mustererkennung. 19. DAGM-Symposium (1997).

HENN, Stefan und Kristian WITSCH. "A Multigrid Approach for Minimizing a Nonlinear Functional for Digital Image Matching", Computing 64(4) (1999), 339-348.

HENN, Stefan und Kristian WITSCH. "Iterative Multigrid Regularization Techniques For Image Matching", SIAM Journal on Scientific Computing (SISC), 23(4) (2001), 1077-1093.

HENN, Stefan und Kristian WITSCH. "Multi-Modal Image Registration Using A Variational Approach", SIAM Journal on Scientific Computing (SISC), (demnächst).

JARRE, Florian und Josef STOER Optimierung, Berlin, 2003.

SCHMIDT, Günter. "Gaußsche Fehlerquadratmethode", Numerisches Rechnen, Algorithmen und Computer. Mainz 1985, 117-147.

STUDHOLME, Colin. Measures of 3D Medical Image Alignment. Dissertation. London 1997.

TIKHONOV, Andrei Nikolaevich. "Solution of incorrectly formulated problems and the regularization method", Soviet Math Dokladi (1963a), 1035-1038.

TIKHONOV, Andrei Nikolaevich. "Regularization of incorrectly posed problems", Soviet Math Dokladi (1963b), 1624-1627.

VIOLA, Paul und William WELLS. "Alignment by maximization of Mutual Information", International Journal of Computer Vision, 24(2) (1997), 137-154.

WELLS, William, Paul VIOLA, Hideki ATSUMI, Shin NAKAJIMA und Ron KIKINIS. "Multi-modal volume registration by maximization of mutual information", Medical Image Analysis 1 (1996), 35-51.