Thesen zur Dissertation

Simultaneous Activity and Attenuation Reconstruction

in Single Photon Emission Computed Tomography,

a Nonlinear Ill-Posed Problem

eingereicht in der Wissenschaftsdisziplin Numerische Analysis an der

Mathematisch-Naturwissenschaftlichen Fakultät der Universität Potsdam

von Volker Dicken

im November 1997

1. Ziel der Arbeit
   In der nuklearmedizinischen Diagnose wird häufig die Single Photon Emission Computed Tomography (SPECT) eingesetzt. Sie soll Stoffwechselfunktionen sichtbar machen durch indirekte Beobachtung der Verteilung radioaktiv markierter Pharmaka. Aus mehreren Szintigrammen, die aus verschiedenen Richtungen (30-120) mit richtungssensitiven -Kameras aufgenommen werden, ist ein 3D-Modell der Aktivitätsverteilung zu rekonstruieren. Dabei ist die Dämpfung der emitierten Strahlung zu berücksichtigen, diese ist proportional zur Gewebedichte. Ziel ist die Früherkennung krankhafter Veränderungen, bevor jene zu anatomischen Abweichungen führen. Anatomische Veränderungen sind mit klassischen Tomographieverfahren (CT, Ultraschall, NMR/MRT) besser diagnostizierbar.
2. Mathematisches Modell
   Die SPECT-Daten werden für jeden der transaxialen 2D Schnitte zu Sinogrammen angeordnet. Diese können näherungsweise als Meßwerte der gedämpften (engl. attenuated) Radontransformation (ATRT) aufgefaßt werden:

   

   Dabei parametrisiert die Blickrichtung der Kamera und den Abstand vom Zentrum der Rotation. Die Meßdaten sind stark verrauscht, Ursachen sind u.a. das stochastische Rauschen des radioaktiven Zerfallsprozesses und die Vernachlässigung der Streuung im Modell. Weiter führt die faktische Integration über spitze Kegel anstelle der Linien im Modell zu systematischen Fehlern in der Messung.

3. ATRT: Ein schlecht gestelltes nichtlineares Problem

   Die Rekonstruktionsaufgabe der simultanen Approximation von und aus gegebenen verrauschten Meßwerten der ATRT-Daten ist nichtlinear in bezug auf und ein im mathematischen Sinne schlecht gestelltes Problem. Da die ATRT hochfrequente Anteile der Signale und stark dämpft, müßte ein (möglicherweise gar nicht existierender) inverser ATRT Operator unbeschränkt, also unstetig sein, um die hochfrequenten Anteile zu rekonstruieren.

4. Tikhonov Regularisierung
   Schlecht gestellte, nichtlineare Probleme können durch Regularisierung näherungsweise gelöst werden. Ein im linearen Fall gut untersuchtes Verfahren sucht ein Minimum des Tikhonov-Funktionals

   

   als Approximation der Lösung (mit ungestörten Daten ). Engl, Kunisch und Neubauer studierten das Verfahren für nichtlineare Parameteridentifikationsprobleme bei partiellen Differentialgleichungen. Sie zeigten (Linz, 1989), daß unter gewissen Bedingungen jede hinreichend gute Approximation des das Tikhonov-Funktional minimierenden Vektors eine ordnungsoptimale Approximation der wahren Lösung ergibt. Die Bedingungen richten sich an den Operator und die Hilberträume sowie an das Verhältnis des Regularisierungsparameters zum Fehlerniveau mit . Ferner muß der Abstand der wahren Lösung zur Startlösung im Strafterm ein hinreichend kleines Urbild unter dem Adjungierten der Ableitung des Operators haben.

   
   Die wichtigsten theoretischen Ergebnisse der Arbeit

5. Umdefinierter ATRT Operator
   Der originale ATRT Operator erfüllt die Bedingungen der Linzer Arbeit nur unter unrealistischen Annahmen an die Regularität der Funktionen . Nach Umdefinition der ATRT für aus physikalischen Gründen ausgeschlossene Werte von (Integration über endliche Strecken, Umdefinition der Exp-Abbildung für negative Argumente) erfüllt die ATRT die Voraussetzungen der genannten Arbeit unter realistischen Voraussetzungen an . Hinreichende Bedingungen dafür sind, daß über dem gewichteten Produkt-Sobolevraum als Operator nach aufgefaßt wird und gilt. Darin bezeichnet ein Kreisgebiet mit Radius in und eine beschränkte Gewichtsfunktion darüber. Schränkt man auf beschränkte Funktionen, genauer ein, so sind die Bedingungen bereits für und erfüllt.

   Daher gilt das zentrale theoretische Resultat: Für konvergieren die Tikhonov-Lösungen des nichtlinearen ATRT Problems mit ordnungsoptimaler Geschwindigkeit gegen eine verallgemeinerte Lösung, falls die Startlösung hinreichend gut ist und geeignet gewählt wird. Eine notwendige Bedingung (zumindest für ) an die Startlösung ist dabei . Der optimale Regularisierungsparameter ist dann durch bestimmt. Analoge Ergebnisse konnten für zwei semidiskrete Varianten der ATRT erzielt werden.

   Der Beweis stützt sich wesentlich auf die beiden folgenden funktionalanalytischen Resultate und bekannte Einbettungssätze über Sobolevräume und einen neuen über gewisse Quasisobolevräume über .

6. Glättungseigenschaften der Fan-beam-Transformation
   Die Fan-beam-Transformation (Halbstrahl-Radontransformation)

   

   hat Glättungseigenschaften als Operator zwischen und dem Schnitt zweier gewisser Quasisobolevräume der Ordnung bzw. über . Der Schnittraum, der das Bild von enthält, ist stetig in den Raum eingebettet, falls .

7. Fréchet-Differenzierbarkeit der Exponentialabbildung
   Eine für Argumente mit negativen Spitzen geeignet umdefinierte Exponentialabbildung mit ist ein n-fach Fréchet-differenzierbarer Operator zwischen und , falls und ein endliches Volumen hat. Die (operatorwertige) n-te Abbleitung ist Lipschitz-stetig, falls und hinreichend regulär ist.

   
   Numerische Realisierung

8. Minimierung des Tikhonov-Funktionals
   Numerisch zu lösen bleibt die Minimierung des Tikhonov-Funktionals. Die Struktur legt einen Gauß-Newton-Ansatz nahe. Iterativ werden quadratische Formen

   

   minimiert und gesetzt, wobei abwechselnd Suchrichtungen der Form und zugelassen werden. Dies halbiert die Größe der jeweils zu lösenden Gleichungssysteme. Konvergenzaussagen sind (lokal) beweisbar, doch der implementierte Algorithmus besitzt einen großen numerischen Konvergenzbereich.

9. Numerische Realisierung mittels konjugierter Gradienten
   Die Minimierung von verlangt das Lösen der linearen Gradientengleichung im Produkt-Sobolevraum . Dies ist äquivalent zur Lösung von in . Darin bezeichnet die beiden Komponenten der gedämpften Radon-Rückprojektion, einen mit dem Sobolevraum assoziierten Pseudodifferentialoperator (PDO) und den Multiplikationsoperator zum Gewicht . Das letztere lineare Gleichungssystem ist mit einem konjugierten Gradientenverfahren (CG) effizient lösbar. Ein früher Abbruch der CG Iteration ergibt dabei einen zusätzlichen regularisierenden Effekt.

10. Implementierung der Operatoren
    1. ATRT, Ableitungen und Rückprojektionen
       Die Diskretisierung der Operatoren ist mit ausreichender Genauigkeit und effizient durch bilineare Interpolation auf gegen das Pixelgitter rotierten Hilfgittern mittels vorberechneter quantisierter Interpolationsgewichte möglich.

    2. Realisierung des Sobolev-PDO
       Der Operator kann sehr effizient durch eine Faltung mit kleiner Maske (3x3) approximiert werden, da nur qualitative Eigenschaften von wesentlich sind.
    3. Bestimmung der Gewichtsfunktion
       Eine gute Wahl des Gewichts für den Operator ergibt sich aus der nichlinearen Radon-Rückprojektion (für ein geeignetes ) als (mit ). Solch ein Gewicht erzwingt eine stärkere Annäherung an die Startlösung in jenen Bereichen, über die die gemessenen Daten wenig Informationen enthalten.

11. Beschleunigung der Implementation
    1. Die mehrfache Auswertung der Exponentialfunktion von ähnlichen Halbstrahlintegralen kann durch eine rekursive Implementation mit Taylorapproximationen beschleunigt werden.
    2. Eine schnelle Auswertung von und ist durch Speicherung einiger größerer Datensätze und Zwischenergebnisse zu Begin der CG Iterationen für möglich.
    3. Einige CG-Iterationen lassen sich einsparen, wenn dem CG Verfahren ein guter Startwert vorgegeben wird. Ein solcher ist für leicht durch 2-3 Iterationen einer neuen dämpfungskorrigierten Variante der gefilterten Rückprojektion erhältlich.
    4. Ein guter Startwert für die Dämpfung ergibt sich aus der inversen Gewichtsfuntion zu (mit z.B. der Dämpfungkonstanten von Wasser ), da die Funktion unter den medizinischen Bedingungen eine gute Approximation des Trägers von ist.

12. Numerische Resultate
    In Testrechnungen mit Phantomdaten einer Herzuntersuchung mit 20%multiplikativem Rauschen zeigt der Algorithmus folgende Eigenschaften.
    1. Für festes wird sehr gut aus verrauschten Daten rekonstruiert, falls eine gute Approximation von bekannt ist.
    2. Das nichtlineare Rekonstruktionsproblem für bei gegebenem und ist gut lösbar, wenn das wahre bekannt ist. Mit einer glatten Approximation von anstelle des wahren tritt ein Negativ von im Fehlerbild für bereits für ungestörte Daten auf.
    3. Werden iterativ und rekonstruiert, so wird in guter Qualität gefunden. Die Rekonstruktion von zeigt das erwähnte Negativ von . Dieser Effekt trat in zwei neueren Publikationen zum Thema sehr deutlich auf (Young und Krol u.a., Syracuse NY, '95,'96). Diese Arbeiten setzten gänzlich verschiedene numerische Verfahren (ART-IntraSPECT, EM-IntraSPECT) ein. Durch eine größere Zahl von Stabilisierungsparametern im vorgestellten Verfahren sollte es möglich sein, den Effekt durch geschickte Parameterwahl weiter einzuschränken.

13. Ergebnisse mit realen Daten

    Der Algorithmus erlaubt die Rekonstruktion von aus realen Daten (Bergman Klinikum Potsdam, Uniklinik Göttingen, Siemens Erlangen) in einer dem klinischen Standard mindestens vergleichbaren Qualität. Ziel der Arbeit ist die Korrektur von Artefakten, die aufgrund fehlerhafter Dämpfungskorrekur entstehen. Aufgrund des starken stochastischen Rauschens in den zur Verfügung stehenden Daten verschiedener Myocard Untersuchungen werden diejenigen Artefakte, welche zu unberücksichtigten Dichtevariationen korreliert sind, durch das hohe Rauschniveau in den Rekonstruktionen aus realen Daten überdeckt. Ein Nachweis der Wirksamkeit des Verfahrens ist daher mit den vorhandenen Daten schwierig. Die zum Teil fehlende Antisymmetrieinformation bei Datensätzen aus älteren Geräten verhindert eine Anwendung des Algorithmus zur Approximation von in einigen Fällen. Deutliche Fortschritte für die Rekonstruktion von ergeben sich aber bereits schon dort durch den Einsatz der dämpfungskorrigierten gefilterten Rückprojektion mit . Eine erfolgreichere Anwendung des Algorithmus auf reale Daten erfordert noch eine Vorverarbeitung der Daten zur Verbesserung des Signal-Rausch-Verhältnisses oder bessere Meßdaten mit höheren Zählraten, so daß Artefakte welche von der Diskrepanz zwischen und herrühren bzw. deren Unterdrückung, überhaupt in Erscheinung treten können. Für herkömmliche Daten sind die Resultate bislang kaum von denjenigen einfacher moderner dämpfungskorrigierter Rekonstruktionsalgorithmen zu unterscheiden.