ANGEWANDTE MATHEMATIK

Tomographische Bildrekonstruktion mit direkten algebraischen Verfahren

Inauguraldissertation zur Erlangung des Doktorgrades

der Naturwissenschaften im Fachbereich Mathematik der Westfälischen Wilhelms-Universität Münster

vorgelegt von

Willehad Klaverkamp

aus Münster Dekan: Prof.’in Dr. Elisabeth Peveling Erster Referent: Prof. Dr. Frank Natterer Zweiter Referent: Prof. Dr. Colin W. Cryer Mündliche Prüfung: Tag der Promotion:

Inhaltsverzeichnis

Einleitung 2

1 Mathematische Grundlagen und Bezeichnungen 6 1.1 Die Radon-Transformation und ihre Singulärwertzerlegung . . . . . . . . 6 1.2 Die Fouriertransformation . . . . . . . . . . . . . . . . . . . . . . . . . . 11 1.3 Bandbeschränkung und Prolate Spheroidal Wave Functions (Pswf’s) . . 12 1.4 Besselfunktionen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

2 Mathematisches Modell 16 2.1 Meßgeometrien und das Inversionsproblem . . . . . . . . . . . . . . . . . 16 2.2 Zurückführung auf Punktauswertungsfunktionale . . . . . . . . . . . . . 21 2.3 Problematisierung des in 2.2 skizzierten Vorgehens . . . . . . . . . . . . 22

3 Auflösungsuntersuchungen 25

3.1 Die Standard-Geometrie . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.1.1 Abschätzungen für die Fehler r_{1} und r_{2} . . . . . . . . . . . . . . 28 3.1.2 Stabilität . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

3.2 Die Interlaced-Geometrie . . . . . . . . . . . . . . . . . . . . . . . . . . 39

4 Der Direkte Algebraische Algorithmus (DAA) 52 4.1 Aufstellung des linearen Gleichungssystems . . . . . . . . . . . . . . . . 52 4.2 Notwendige Grundlagen aus der Darstellungstheorie . . . . . . . . . . . 54 4.3 Anwendung der Gruppentheoretischen Methode auf die diskrete Radon-

Transformation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 4.4 Verallgemeinerte Lösung der Gleichungssysteme. . . . . . . . . . . . . . 67 5 Numerische Ergebnisse 73

5.1 Testobjekte . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 5.2 Rekonstruktionen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75

Einleitung

Eine klassische und gut untersuchte Aufgabenstellung der Integralgeometrie ist das Radonsche Inversionsproblem der Rekonstruktion einer Funktion f : IR^{n} → IR aus ihren Integralen über alle Hyperebenen des IR^{n}. Die Behandlung und Lösung des Problems durch J. Radon im Jahre 1917 trägt noch ganz die Grundzüge der klassischen Analysis; Radon konnte durch direkte Berechnungen eine konstruktive Inversionsformel angeben. In den fünfziger und sechziger Jahren ist seine Fragestellung unter zum Teil wesentlich abstrakteren Gesichtspunkten von Mathematikern wie I. M. Gelfand und S. Helgason untersucht und erweitert worden. Daraus resultierte eine Fülle schöner Ergebnisse und Einsichten, und das Radonsche Problem schien eine gänzlich der Reinen Mathematik vorbehaltene Fragestellung zu sein. Dies änderte sich jedoch im Jahre 1963, als A. Cormack im Zusammenhang mit einer damals neuen Aufgabenstellung der medizinischen Röntgen-Diagnostik sich genau diesem Problem (für den Fall n = 2) gegenübergestellt sah. Er leitete (in Unkenntnis der Arbeit Radons) eine eigene Inversionsformel her, die er auch praktisch anwendete und deren Unbrauchbarkeit er recht bald einsehen mußte.

Mit dem Auftauchen des Radonschen Problems in einer konkreten praktischen Anwendung (heutzutage als Computertomographie jedem bekannt) begannen sich sofort auch Mathematiker der angewandten Richtung für das Thema zu interessieren. Dabei zeigte sich, daß die bereits gewonnenen Antworten für praktische Belange nicht ausreichten, so ist z.B. die Radonsche Formel für eine Implementation ungeeignet. Die praktischen Probleme resultieren aus dem Zusammenwirken zweier Tatsachen. Zum einen ist stets eine Diskretisierung erforderlich: es können nicht die Integrale über alle Geraden des IR^{2} gemessen werden. Die analytischen Inversionsformeln sind also stets durch gewisse Näherungen zu ersetzen, die zu Fehlern führen. Zum zweiten ist das Problem im Sinne Hadamards schlecht gestellt, so daß mit großen Verstärkungen von Daten- und systematischen Fehlern gerechnet werden muß. Desweiteren spielen für praktische Belange Effizienzfragen eine nicht unerhebliche Rolle. Die Anzahl der Daten ist, zumindest in medizinischen Anwendungen immens, so daß der Einsatz schneller Algorithmen erforderlich ist. Erst im Jahre 1971 konnte von G. N. Ramachandran und A. V. Lakshminarayanan ein Algorithmus vorgeschlagen werden, der den Anforderungen der Praxis weithin gerecht wurde. Er firmiert unter dem Namen gefilterte Rückprojektion und ist heute auf wohl allen kommerziellen Computertomographen implementiert. Obgleich dieser Algorithmus die Möglichkeiten tomographischer Rekonstruktion sicherlich weitgehend ausnutzt, ging und geht die Suche nach besseren Algorithmen natürlich weiter. Eine mögliche Reaktion auf die oben geschilderten Schwierigkeiten ist es, das Problem schon auf einer früheren Stufe zu diskretisieren. Dies kann geschehen, indem man sich gleich auf die Suche einer approximativen Lösung aus einem endlichdimensionalen Teilraum beschränkt. Die gesuchten Koeffizienten der Lösungsfunktion werden sich dann aus einem linearen Gleichungssystem ergeben, welches infolge der zu erwartenden schlechten Kondition mit Regularisierungstechniken zu lösen ist. Das Radonsche Problem ist somit in ein Problem der linearen Algebra transformiert worden, dessen Behandlung im Prinzip ohne irgendwelche Kenntnisse der theoretischen Grundlagen möglich ist. Ein ganz aus diesem ”numerischen Geist“ geborener Algorithmus ist der von F. Natterer im Jahre 1982 vorgestellte Direkte Algebraische Algorithmus (DAA). Er berechnet eine L_{2}-Bestapproximation der gesuchten Funktion f aus einem Teilraum V, der durch charakteristische Funktionen schmaler Streifengebiete aufgespannt wird. Gerade darin liegt aber auch schon ein gewisses Problem: Anders als z. B. im Fall der Methode der Finiten Elemente (FEM) kennen wir die Approximationsgüte dieses Unterraumes nicht! Bei der FEM haben die Ansatzfunktionen charakteristischerweise einen kleinen Träger. Diese Eigenschaft ist wesentlich für den Beweis der Approximationsgüte des FEM-Raumes. Verfeinert man hingegen beim DAA die Diskretisierung, so werden die Streifen nur schmaler, behalten aber ihre Länge, so daß FEM-Abschätzungen nicht zugänglich sind. Die eigentlich sehr befriedigende Eigenschaft der Lösung, L_{2}-Proximum an f in V zu sein, hängt aber ohne jeden Nachweis der Approximationseigenschaften von V förmlich etwas in der Luft. Die Aussage, daß die von der gefilterten Rückprojektion gelieferte Näherungslösung auch in V liegt, stimmt nur cum grano salis. Sie ist nur richtig, falls im letzten Schritt der gefilterten Rückprojektion die sogenannte ”nearest neighbour interpolation“ verwendet wird, und genau dies ist als unzureichend erkannt worden und wird folglich in der Praxis vermieden. Die Wahl des Unterraumes V beim DAA folgt einer anderen Intension: Um das entstehende Gleichungssystem (man hat mit ca. 10^{5} Unbekannten zu rechnen) noch mit direkten Methoden lösen zu können, ist folgende Invarianzeigenschaft erforderlich: mit f ∈ V muß auch f ∘ D_{j} ∈ V gelten, wobei D_{j} die Drehung um ein Winkelinkrement beschreibt. Die Effizienz des DAA leitet sich aus einem Zerfall des sehr großen Gesamtsystems in viele nur noch mäßig große Einzelsysteme ab. Dieser Zerfall in Einzelsysteme ist Konsequenz der Invarianzeigenschaft. Ein weiterer Vorteil des DAA besteht in der realistischeren Modellierung des Meßvorgangs. Wir müssen nicht mehr von der unzutreffenden Annahme ausgehen, Linienintegrale zu messen. Vielmehr ist es möglich, die Abmessungen und Empfindlichkeiten der Detektoren mit in das Modell einzubauen. Wir stellen in dieser Arbeit der oben skizzierten Philosophie des DAA eine ganz neue und andere gegenüber, indem wir ihn als eine diskrete Form der Invertierung der Radon-Transformation über ihre Singulärwertzerlegung begreifen. Um dies zu können, haben wir allerdings die oben erwähnten (symmetrischen) Einzelsysteme über die Entwicklung nach Eigenvektoren zu lösen. In Kapitel 1 und 2 werden zunächst die mathematischen Grundlagen bzw. die Modellierung des Meßvorgangs vorgestellt. Das genaue Studium der Singulärwertzerlegung in Kapitel 3.1 ermöglicht dann die folgende erstaunliche Erkenntnis: Die Eigenvektoren der Einzelsysteme zerfallen nahezu eindeutig in zwei Klassen, von denen die zweite zum Aufbau einer guten Rekonstruktion von f überhaupt nicht benötigt wird. Voraussetzung dazu ist allerdings eine gewisse Regularität von f, nämlich eine näherungsweise Bandbeschränktheit. Das hocherfreuliche Ergebnis ist nun, daß die zweite Klasse gerade aus den Eigenvektoren zu extrem kleinen Eigenwerten zusammengesetzt ist. Läßt man diese Eigenvektoren bei der Entwicklung der Lösung einfach weg, so erreicht man automatisch eine ”optimale“ Stabilität der Gleichungssysteme! Ein weiterer wesentlicher Vorteil unseres Standpunktes ist es, den DAA erstmalig auch für den Interlaced-Datensatz modifizieren zu können. Dabei handelt es sich um ein Abtastschema, bei dem, gegenüber der Standard-Abtastung, nur die Hälfte der Daten gemessen wird. Es wurde von A. Cormack im Jahre 1978 durch eine heuristische Überlegung gefunden, kann aber durch ein verallgemeinertes Abtasttheorem auch mathematisch begründet werden. Der Nachweis, daß die gefilterte Rückprojektion aus diesem Datensatz gute Rekonstruktionen liefert, ist 1986 von H. Kruse geliefert worden. Das dortige Vorgehen lieferte aber überhaupt keinen Hinweis darauf, wie im Fall des DAA zu verfahren ist. Wir können in Kapitel 3.2 folgendes zeigen: Benutzt man den Interlaced-Datensatz, so wird die rechte Seite des Gleichungssystems zwar erheblich verfälscht, der Fehler hat jedoch eine angenehme Eigenschaft. Er steht nahezu senkrecht auf den Eigenvektoren, die zum Aufbau des Bildes benötigt werden. Der DAA besticht hier durch seine Einfachheit: Alles was wir bei seiner Anwendung auf Interlaced-Daten zu tun haben, ist auf die Entwicklung nach den ohnehin irrelevanten Eigenvektoren zu verzichten!

Bei der gefilterten Rückprojektion ist dagegen recht erheblicher zusätzlicher Rechenaufwand erforderlich. Das kuriose Verhalten des Fehlers beruht auf besonderen Eigenschaften der Singulärfunktionen bzw. ihrer Fouriertransformierten. Im 4. Kapitel wird der DAA genau beschrieben, wobei insbesondere auf die Rolle eingegangen wird, welche die Invarianzeigenschaft der Meßgeometrie für die Effizienz des Algorithmus spielt. Im 5. Kapitel werden dann Rekonstruktionen vorgestellt, welche die theoretischen Überlegungen bestätigen. Abschließend kann festgehalten werden, daß auch für das Verständnis des einfachen, auf rein numerischen Prinzipien beruhenden DAA, eine genaue Kenntnis der Abbildungseigenschaften der Radontransformation und ihrer Singulärwertzerlegung vonnöten ist. Dies verleitet vielleicht dazu, unter Ausbeutung eines Aperçus von S. Beckett, die Hoffnung auszusprechen, daß unsere Behandlung des DAA ein Beispiel für die beruhigende Vorstellung von Reiner Mathematik und Numerik als einem singenden Mohrenduo aus dem Teatro dei Piccoli ist, ebenso beruhigend wie der Anblick eines sorgfältig zusammengeklappten Schinkenbrots.

An dieser Stelle möchte ich mich herzlich für die ”TEX-nische“ Hilfe bei der Erstellung dieser Arbeit bedanken. Es haben getippt die Diplom-Mathematikerinnen A. Fuhrmeister, A. Mondry, D. Weßels, sowie die männlichen (Diplom-) Kollegen R. Holtrode und F. Wübbeling. Ebenso möchte ich mich bei Dr. A. Faridani für einige Korrekturhinweise und bei Prof. Dr. F. Natterer für die Anregung zu dieser Arbeit bedanken.

1 Mathematische Grundlagen und Bezeichnungen

Wir geben einen Überblick über die von uns verwendeten Bezeichnungen, Hilfsmittel und mathematischen Grundlagen.

1.1 Die Radon-Transformation und ihre Singulärwertzerlegung

Von zentraler Bedeutung ist folgende Integraltransformation, die nach dem österreichischen Mathematiker Johann Radon benannt ist. Wir definieren sie zunächst auf einer kleinen Funktionenklasse, dem Schwartzschen Raum S(IR^{n}) der schnell fallenden C^{∞}- Funktionen.

Definition 1.1.1 Sei f ∈ S(IR^{n}), ω ∈ S^{n−1} = {x ∈ IR^{n}, ∥x∥ = 1} , s ∈ IR, dann heißt die Abbildung

f ↦ Rf mit Rf(ω, s) = ∫_{<x,ω>=s} f(x)dx  (1.1)

Radon-Transformation. Hierbei steht dx für das Lebesgue-Maß auf der (n-1)-dimensio-nalen Hyperebene, senkrecht zu ω mit euklidischem Abstand s vom Ursprung.

Bemerkung 1.1.1 Die vorausgesetzte Regularität von f sichert die Existenz des Integrals. Die Funktion Rf ist auf Z := S^{n−1} × IR definiert, und man sieht leicht ein, daß Rf ∈ S(Z) ist, worunter der Vektorraum der C^{∞}-Funktionen auf Z, welche in der letzten Variable schnell fallen, zu verstehen ist. Um die Abhängigkeit von der zweiten Variablen bei festgehaltenem ω ∈ S^{n−1} zu betonen, werden wir gelegentlich die Schreibweise

R_{ω}f := Rf(ω, ·)  (1.2)

benutzen, R_{ω} bildet dann also von S(IR^{n}) in den S(IR) ab.

Wir werden ausschließlich den Fall n = 2 betrachten, in welchem wir ω = (cos φ, sin φ)^{⊤} ∈ S^{1} mit dem Winkel φ ∈ [0, 2π[ identifizieren können, so daß wir auch R_{φ} für R_{ω} schreiben werden.

Im Hinblick auf Anwendungen im Zusammenhang mit der Computertomographie ist es sicherlich wünschenswert, den obigen Definitionsbereich von R abändern zu können. Es dürfte kaum sinnvoll sein, zweidimensionale Bilddichten als Elemente des metrischen Raumes S(IR^{2}) zu interpretieren. Zwar liegt S(IR^{2}) dicht in wesentlich größeren Funktionenräumen, doch ist die Topologie natürlich gänzlich unpassend: zwei Bilder, die sich augenscheinlich (und das ist nun einmal das einzig sinnvolle Kriterium) sehr gleichen, können als Elemente von S(IR^{2}) gedeutet, einen sehr großen Abstand haben. Einen normierten Raum zu finden, dessen Topologie zu der des Augenscheins ”äquivalent“ ist, ist nicht einfach. Am besten dürfte der Sobolevraum H^{β} mit β ”etwas unter 1/2“ geeignet sein (siehe Natterer [19]).

Wir wollen zunächst einmal nur verlangen, daß f eine L_{2}-Funktion mit kompaktem Träger sein soll. Sei also D = D_{1}(0) der Kreis mit Radius 1 um Null, L_{2}(D) die Menge der quadratintegrablen Funktionen mit Träger in D, dann können wir Rf auch für f ∈ L_{2}(D) definieren.

Satz 1.1.2 Die Operatoren R_{ω} und R aus (1.1) bzw. (1.2) lassen sich eindeutig zu stetigen Operatoren von L_{2}(D) nach L_{2}([−1, 1], w^{−1}) bzw. L_{2}(D) nach L_{2}(Z, w^{−1}) fortsetzen. Die Fortsetzungen bezeichnen wir ebenfalls mit R_{ω} bzw. R. Hierbei ist w die Gewichtsfunktion w(s) := \sqrt{1 − s^{2}}.

Beweis: Der Fortsetzungsprozeß ist kanonisch; da S dicht in L_{2}(D) liegt, haben wir nur zu zeigen, daß

R_{ω} : L_{2}(D) ∩ S → L_{2}([−1, 1], w^{−1})

stetig ist, bekanntlich folgen dann Existenz und Eindeutigkeit der Fortsetzung, entsprechendes gilt für R. Sei also f ∈ L_{2}(D) ∩ S, so liefert die Cauchy-Schwarzsche Ungleichung:

∫_{−1}^{1} |R_{ω}f(s)|^{2} w^{−1}(s) ds = ∫_{−1}^{1} | ∫_{−w(s)}^{w(s)} f(sω + tω^{⊥})dt|^{2} w^{−1}(s)ds

≤ 2 · ∫_{−1}^{1} ∫_{−w(s)}^{w(s)} |f(sω + tω^{⊥})|^{2} dtds
= 2 · ∥f∥_{L_{2}(D)}^{2} .

Demnach auch:

∥Rf∥_{L_{2}(Z,w^{−1})}^{2} = ∫_{0}^{2π}∫_{−1}^{1} |R_{ω}f(s)|^{2}w^{−1}(s) ds dφ ≤ 4π · ∥f∥_{L_{2}(D)}^{2} .

□

Bemerkung 1.1.3 Die in den obigen Abschätzungen angegebenen Konstanten sind nicht zu verbessern, es gilt also ∥R_{ω}∥ = \sqrt{2}, ∥R∥ = 2\sqrt{π}.

Nun lassen sich auch leicht die stetigen Adjungierten R_{ω}^{∗}, R^{∗} berechnen, es gilt

R_{ω}^{∗}h = (w^{−1} h) (< ·, ω >) ,  (1.3)
R^{∗}g = ∫_{0}^{2π}g(ω, < ·, ω >) w^{−1}(< ·, ω >)dφ .  (1.4)

Besonders übersichtlich werden die Abbildungseigenschaften eines Operators, wenn man ihn durch geeignete Basiswahl in eine Diagonalform überführt. Dies ist bei der Singulärwertzerlegung der Fall; genauer wird das Tupel ((f_{k})_{k}, (g_{k})_{k}, (σ_{k})_{k}) singuläres System eines stetigen Operators A : H_{1} → H_{2} (H_{1}, H_{2} seien hierbei separable Hilberträume) genannt, falls

1. {f_{k}}_{k} ein ONS in H_{1} ist,

2. {g_{k}}_{k} ein ONS in H_{2} ist,

3. σ_{k} ∈ IR^{+}, ∃M ∈ IR, so daß σ_{k} < M ∀k gilt,

4. ∀f ∈ H_{1} gilt: Af = ∑_{k=1}^{∞} σ_{k}(f, f_{k})_{H_{1}} g_{k}.

Die Darstellung von A in der Form 4. wird Singulärwertzerlegung von A genannt.

Leicht zu ziehende Schlußfolgerungen sind die folgenden Darstellungen:

A^{∗}g = ∑_{k=1}^{∞} σ_{k}(g, g_{k})_{H_{2}} f_{k} ,  (1.5)
AA^{∗}g = ∑_{k=1}^{∞} σ_{k}^{2}(g, g_{k})_{H_{2}} g_{k} ,  (1.6)
A^{∗}Af = ∑_{k=1}^{∞} σ_{k}^{2}(f, f_{k})_{H_{1}} f_{k} .  (1.7)

Die ONSe {g_{k}}_{k}, {f_{k}}_{k} lassen sich also als Eigenfunktionen der selbstadjungierten Operatoren AA^{∗} bzw. A^{∗}A (zu von Null verschiedenen Eigenwerten) berechnen; es besteht der Zusammenhang

A^{∗}g_{k} = σ_{k}f_{k} ,  (1.8)
Af_{k} = σ_{k}g_{k} .  (1.9)

Die verallgemeinerte Inverse des Operators A läßt sich nun besonders bequem mit Hilfe der Singulärwertzerlegung angeben. Dabei versteht man unter der verallgemeinerten Inversen A^{+} von A diejenige lineare Abbildung, welche jedem g ∈ range(A)+(range(A))^{⊥} die Minimum-Norm-Lösung der Gleichung Af = g zuordnet. D.h. A^{+}g minimiert den Defekt Af − g in der Norm von H_{2} und hat unter allen Lösungen dieser Minimierungsaufgabe die kleinste Norm in H_{1}. Es ist nicht schwer, die Wohldefiniertheit und die Äquivalenz zu der Charakterisierung von A^{+}g als (eindeutige) Lösung der Normalgleichungen

A^{∗}Af = A^{∗}g

in \overline{range(A^{∗})} zu zeigen. Für A^{+} existiert die folgende elegante Darstellung:

Satz 1.1.4 Für g ∈ range (A) + (range (A))^{⊥} gilt in H_{1}:

A^{+}g = ∑_{k=1}^{∞} σ_{k}^{−1}(g, g_{k})_{H_{2}} f_{k} .

Beweis: (siehe Natterer [20], S. 87)

Bemerkung 1.1.5 Wir können die ONSe stets so umordnen, daß die zugehörigen Singulärwerte {σ_{k}}_{k∈IN} eine schwach monoton fallende Folge bilden. Gilt dann σ_{k} ↘ 0, so ist A^{+} unstetig, und wir werden i.a. noch einmal Abstriche an der angestrebten ”Lösung“ der Gleichung Af = g machen müssen, indem wir die Reihenentwicklung in Satz 1.1.4 vorzeitig abbrechen oder auch andere Regularisierungstechniken, wie die Tychonoff-Phillips-Regularisierung verwenden. Hinter diesem Konzept verbirgt sich natürlich stets die Hoffnung, daß die sehr kleinen Singulärwerten entsprechenden Singulärfunktionen f_{k} in irgendeinem Sinne unwichtig für die gesuchte ”Lösung“ von Af = g sind.

Ist A injektiv (dies ist bei der Radon-Transformation der Fall), so gilt A^{−1} = A^{+}, und wir erhalten für alle f ∈ H_{1} mit g = Af:

f = A^{−1} (Af) = ∑_{k=1}^{∞} σ_{k}^{−1} (Af, g_{k})_{H_{2}} f_{k} = ∑_{k=1}^{∞} (f, σ_{k}^{−1}A^{∗}g_{k})_{H_{1}} f_{k}
= ∑_{k=1}^{∞} (f, f_{k})_{H_{1}} f_{k} ,

d.h. {f_{k}}_{k} ist ein vollständiges ONS in H_{1}.

Die Singulärwertzerlegung der Radon-Transformation wurde von Davison [1] und Louis [13] gefunden.

Es gilt folgendes: Die Funktionen

g_{ml}(φ, s) = \{ \begin{matrix}\frac{1}{π} u_{m}(s)e^{ilφ}, & |s| ≤ 1 \\ 0, & sonst \end{matrix} ,  (1.10)

(m ∈ IN, |l| ≤ m, m + l gerade) sind die (in L_{2}(Z, w^{−1})) orthonormierten Eigenfunktionen des Operators RR^{∗} zu den Eigenwerten

σ_{ml}^{2} = \frac{4π}{m + 1} .  (1.11)

Hierbei ist u_{m} = w U_{m}, wobei U_{m} wiederum für das m-te Tschebyscheff-Polynom 2. Art steht. Nach (1.8) bilden dann die Funktionen

f_{ml} = \frac{1}{σ_{ml}} R^{∗} g_{ml}  (1.12)

ein ONS in L_{2}(D), und die Singulärwertzerlegung von R stellt sich demnach wie folgt dar (σ_{m} := σ_{ml}):

Rf(φ, s) = \frac{1}{π} ∑_{m=0}^{∞} σ_{m} ∑_{|l|≤m}′ (f, f_{ml})_{L_{2}(D)} u_{m}(s) e^{ilφ} .  (1.13)

Hierbei soll der Strich an dem Summenzeichen andeuten, daß nur über die Werte von l zu summieren ist, für welche m + l gerade ist, d.h. l und m haben stets gleiche Parität. Das Gleichheitszeichen ist selbstverständlich im Sinne der L_{2}-Konvergenz zu verstehen.

Die Funktionen f_{ml} lassen sich auch geschlossen angeben (der Radialanteil besteht aus den sogenannten Zernike-Polynomen), sie sind jedoch weniger ”handlich“ als die g_{ml}, alles was wir über sie wissen müssen, ist die Identität (1.12).

Bemerkung 1.1.6 Wir ziehen aus dem Verhalten der Singulärwerte (1.11) noch zwei Folgerungen. Da die Folge der σ_{m} monoton gegen 0 fällt, konvergieren die abgeschnittenen Singulärwertzerlegungen in der starken Operatortopologie gegen R. Die Radon-Transformation ist also als gleichmäßiger Grenzwert von Operatoren endlichen Ranges kompakt. Die Asymptotik σ_{m}^{−1} = O (m^{\frac{1}{2}} ) bedeutet, daß es sich bei der Invertierung der Radon-Transformation um ein mäßig schlecht gestelltes Problem handelt.

Ein weiteres wesentliches Hilfsmittel beim Studium der Radon-Transformation ist der Projektionssatz (s. Natterer [20], S. 11).

Satz 1.1.7 Für f ∈ L_{2}(D) und beliebiges ω ∈ S^{1} gilt:

R̂_{ω}f(σ) = (2π)^{\frac{1}{2}} f̂(σω) ∀σ ∈ IR^{1} .  (1.14)

Bemerkung 1.1.8 Nach unserer Voraussetzung ist f ∈ L_{2}(D), dies impliziert einerseits f ∈ L_{1}(D), andererseits R_{ω}f ∈ L_{2}([−1, 1]) und demnach auch R_{ω}f ∈ L_{1}([−1, 1]). Die Fouriertransformation auf beiden Seiten von Gleichung (1.14) kann also im klassischen Sinn verstanden werden. Die transformierten Funktionen sind demnach stetig, so daß (1.14) punktweise gelesen werden darf.

1.2 Die Fouriertransformation

Wir benutzen die ”symmetrische“ Definition der Fouriertransformation, d.h. wir schreiben für f ∈ L_{1}(IR^{n})

f̂(ξ) = (2π)^{− \frac{n}{2}} ∫_{IR_{n}} e^{−i<x,ξ>}f(x)dx ,  (1.15)
f̃(x) = (2π)^{− \frac{n}{2}} ∫_{IR_{n}} e^{i<x,ξ>}f(ξ)dξ .  (1.16)

Gilt sowohl f ∈ L_{1} als auch f̂ ∈ L_{1}, so besteht die Identität (nach eventueller Abänderung auf einer Nullmenge)

˜f̂ = ˆ f̃ = f .

Natürlich werden wir die Symbole ˆ und ˜ auch für die bekannten Erweiterungen der Fouriertransformation auf L_{2}(IR^{n}) und S′(IR^{n}) benutzen.

Eine hilfreiche Identität ist die Poissonsche Summenformel (s. E. M. Stein, G. Weiss [27]):

Satz 1.2.1 Sei f̂(ξ) = (2π)^{− \frac{n}{2}} ∫_{IR}_{n} e^{−i<x,ξ>}f(x)dx, f(x) = (2π)^{− \frac{n}{2}} ∫_{IR}_{n} e^{i<x,ξ>}f̂(ξ)dξ, f bzw. f̂ genüge den Wachstumsbedingungen:

|f(x)| ≤ A (1 + ∥x∥)^{−n−δ} , |f̂(ξ)| ≤ A (1 + ∥ξ∥)^{−n−δ} ,

mit Konstanten A, δ > 0, dann gilt für alle h > 0

∑_{m∈ZZ^{n}} f̂(ξ − \frac{2π}{h} m) = (2π)^{\frac{−n}{2}} h^{n} ∑_{m∈ZZ^{n}} f(hm)e^{−ih<ξ,m>} ,

insbesondere also

∑_{m∈ZZ^{n}} f̂(\frac{2π}{h} m) = (2π)^{\frac{−n}{2}} h^{n} ∑_{m∈ZZ^{n}} f(hm) ,  (1.17)

wobei alle obigen Reihen absolut konvergieren.

1.3 Bandbeschränkung und Prolate Spheroidal Wave Functions (Pswf’s)

Von fundamentaler Bedeutung in der Signal- und Bildverarbeitung ist der Begriff der Bandbeschränktheit; wir nennen eine Funktion f b-bandbeschränkt, falls gilt

supp(f̂) ⊂ Db(0) = {ξ ∈ IR^{n}, ∥ξ∥ < b}.

Eine einfache Anwendung der Differentiationsregeln für parameterabhängige Integrale zeigt, daß eine b-bandbeschränkte Funktion holomorph ist. Der Identitätssatz lehrt dann, daß eine solche Funktion (mit Ausnahme der Nullfunktion) nicht auch ortsbeschränkt sein kann, d.h. keinen kompakten Träger haben kann.

Da die Forderung nach (strikter) Bandbeschränktheit offensichtlich zu einschränkend ist (wir wollen Bilddichten stets als ortsbeschränkt ansehen können), müssen wir diesen Begriff geeignet abschwächen. Dazu dient

Definition 1.3.1 Mit Bb : L_{2}(IR^{n}) → L_{2}(IR^{n}) bezeichnen wir den (Projektions-)Ope-rator der Bandbeschränkung:

Bbf = (f̂χ_{\frac{1}{b}} )˜ ,  (1.18)

mit χ_{c}(x) = χ(cx) und χ(x) = \{ \begin{matrix}1, & ∥x∥ ≤ 1 \\ 0, & sonst \end{matrix} .

Damit vereinbaren wir nun:

Definition 1.3.2 Die Funktion f ∈ L_{2}(IR^{n}) heißt annähernd b-bandbeschränkt, falls für ein hinreichend kleines ϵ > 0 gilt:

∥f − Bbf∥_{L_{2}} ≤ ϵ ∥f∥_{L_{2}} .  (1.19)

Wegen ∥f − Bbf∥_{L_{2}}^{2} = ∥f∥_{L_{2}}^{2} − ∥Bbf∥_{L_{2}}^{2} ist Def. 1.3.2 äquivalent zu der Forderung

∥Bbf∥_{L_{2}}^{2} ≥ (1 − ϵ^{2}) ∥f∥_{L_{2}}^{2} .  (1.20)

Eine weitere triviale, nur die Isometrie-Eigenschaft der Fouriertransformation benutzende, Umschreibung von Def. 1.3.2 ist

∫_{∥ξ∥>b} |f̂(ξ)|^{2}dξ ≤ ϵ^{2} ∥f∥_{L_{2}}^{2} .  (1.21)

Von großem Interesse ist nun das Produkt der Projektoren DrBb, wobei Dr : L_{2}(IR^{n}) → L_{2}(IR^{n}) durch

Drf(x) = \{ \begin{matrix}f(x), & ∥x∥ ≤ r \\ 0, & sonst \end{matrix}  (1.22)

definiert sei. Das Studium dieses Operators (bzw. des Operators BbDr) ist Kernstück der berühmten und grundlegenden Artikelserie [23], [9], [10], [24], [25]. Von Slepian, Landau und Pollak wird darin untersucht, inwieweit eine Funktion gleichzeitig annähernd band-und ortsbeschränkt sein kann. Für den Fall n = 1 wird das in [23] und [9] beantwortet, wobei u.a. die Abhängigkeit der Antwort von dem Produkt b·r genau behandelt wird. Wesentliche Hilfsmittel sind dabei die Eigenfunktionen von BbDr, die sogenannten ”prolate spheroidal wave functions“, bzw. die Berechnung der zugehörigen Eigenwerte.

Die wichtigsten Eigenschaften der Eigenfunktion ψ_{j} und Eigenwerte λ_{j} seien kurz genannt.

1. (ψ_{i}, ψ_{j})_{L_{2}(IR)} = (ψ̂_{i}, ψ̂_{j})_{L_{2}([−b,b])} = δ_{ij}

2. (ψ_{i}, ψ_{j})_{L_{2}([−r,r])} = λ_{j}δ_{ij}

3. Das System {ψ_{j}}_{j=0,1,...} ist vollständig in L_{2} ((−r, r)) und im Raum der b-band-beschränkten Funktionen.

4. λ_{j} > λ_{j+1} ∀j ∈ IN

5. ψ_{j}(x) hat die Parität von j.

6. ψ_{j} hat genau j Nullstellen in (−r, r).

7. Die ψ_{j} sind Eigenfunktionen eines Differentialoperators zweiter Ordnung, welcher bei der Separation der (reduzierten) Wellengleichung in gestreckt-rotationsellip-tischen Koordinaten auftritt. Dieser Eigenschaft verdanken sie ihren Namen.

Von Landau und Widom wird in [11] bewiesen, daß die Eigenwertverteilung für große Werte von b · r die in Abbildung 1.1 dargestellte dichotomische Struktur hat.

Abbildung 1.1:

Dabei sind (etwas ungenau aber suggestiv) die diskreten Eigenwerte λ_{n} durch eine stetige Kurve verbunden worden. Es gibt demnach, abgesehen von einer schmalen Übergangszone, ungefähr \frac{2}{π} · b · r Eigenwerte λ_{n} ≈ 1, die restlichen Eigenwerte sind sämtlich ≈ 0 (tatsächlich gilt λ_{n} ∈]0, 1[ ∀n). Die Größe der Übergangszone wächst dabei logarithmisch mit b · r. Abb. 1.1 enthält fast schon die ”ganze Wahrheit“ über den berühmten, von Ingenieuren erfolgreich angewandten, von Mathematikern mit skeptischem Kopfschütteln bedachten ”Lehrsatz“, wonach die Dimension des Raumes der gleichzeitig b-bandbeschränkten und r-ortsbeschränkten Funktionen für große Produkte b · r nahezu gleich \frac{2}{π} · b · r ist, siehe dazu auch [10] und den Übersichtsartikel von D. Slepian [26]. Die Arbeiten [24] und [25] beschäftigen sich mit der Übertragung der o.a. Ergebnisse auf den Fall mehrerer Dimensionen bzw. auf den Fall diskreter Folgen.

1.4 Besselfunktionen

Bei der Fouriertransformation der u_{m} aus (1.10) treten Besselfunktionen erster Art auf, die wir mit J_{k} bezeichnen und über das Integral

J_{k}(x) = \frac{i^{−k}}{2π} ∫_{0}^{2π}e^{ix cos φ−ikφ}dφ  (1.23)

definieren wollen. Den Index k bezeichnet man als die Ordnung der Besselfunktion.

Sehr wichtig wird es für uns sein, die folgende Eigenschaft der Besselfunktion auszunutzen. Es gilt (s. Natterer [], S. 198) für ϑ ∈]0, 1[:

0 ≤ J_{k}(ϑk) ≤ \begin{matrix}e & − \frac{k}{3}(1−ϑ^{2})^{\frac{3}{2}} \\ \sqrt{2πk} & (1 − ϑ^{2}) & \frac{1}{4} \end{matrix} .  (1.24)

Die Abschätzung (1.24) ist eine (für den Zweck der leichteren Handhabbarkeit) geringfügig abgeschwächte Form der Debyeschen Formel. Sie besagt, grob gesprochen, daß die Besselfunktionen für Argumente, welche kleiner als die Ordnung sind, fast verschwinden.

2 Mathematisches Modell

2.1 Meßgeometrien und das Inversionsproblem

Wir wollen in diesem Abschnitt das Rekonstruktionsproblem mathematisch genau formulieren und vereinbaren zu diesem Zweck einige Bezeichnungen.

Es sei d eine positive, symmetrische, quadratintegrable Funktion, deren Träger im Intervall [\frac{1}{2}, −\frac{1}{2}] liegen möge. Als Detektorfunktion mit der Detektorweite h definieren wir nun die Funktion

d_{h}(x) = \frac{1}{h}d(\frac{x}{h}).

Als unmittelbare Konsequenz dieser Definition läßt sich festhalten:

supp(d_{h}) ⊂ [−\frac{h}{2}, \frac{h}{2}] ,  (2.1)
∥ d_{h} ∥_{L_{1}} = \frac{1}{h} ∫− \frac{h}{2}^{\frac{h}{2}} |d(\frac{x}{h})|dx = ∥ d ∥_{L_{1}} ,  (2.2)
∥ d_{h} ∥_{L_{2}} = (\frac{1}{h^{2}} ∫− \frac{h}{2}^{\frac{h}{2}} |d(\frac{x}{h})|^{2}dx)^{\frac{1}{2}} = h^{− \frac{1}{2}} ∥ d ∥_{L_{2}} ,  (2.3)
d̂_{h}(ξ) = d̂(hξ) und demnach
∥ d̂_{h} ∥_{L_{1}} = \frac{1}{h} ∥ d̂∥_{L_{1}} , sowie  (2.4)
∥ d̂_{h} ∥_{L_{2}} = ∥ d_{h} ∥_{L_{2}}= h^{− \frac{1}{2}} ∥ d ∥_{L_{2}} .  (2.5)

Wir werden zwei verschiedene (Meß-) Geometrien unterscheiden:

1. Die parallele Geometrie

Im allgemeinen Fall gehen wir hierbei davon aus, daß die Meßwerte

g(φ_{j}, s_{ν,j}) := ∫_{−1}^{1} ∫_{−\sqrt{1−s^{2}}}^{√\overline{1−s^{2}}} f(sω_{j} + tω_{j}^{⊥})dt d_{h}(s − s_{ν,j})ds  (2.6)
(ω_{j} = ( \begin{matrix}cos φ_{j} \\ sin φ_{j} \end{matrix} ))  (2.7)

für gewisse Werte φ_{j} ∈ [0, 2π[, j = 0,…, 2p − 1 und s_{ν,j} ∈ [−1, 1], ν = 1,…, q(j), vorliegen.

D.h. die Daten sind Faltungen der Radon-Transformierten von f mit der Detektorfunktion d_{h}:

g(φ_{j}, s_{ν,j}) = Rf(φ_{j}, ·) ∗ d_{h}(s_{ν,j}) .  (2.8)

Wir wollen im folgenden stets eine gewisse Regelmäßigkeit der Geometrie voraussetzen und spezifizieren:

φ_{j} = \frac{jπ}{p} , j = 0,…, 2p − 1 ,  (2.9)

d.h. wir betrachten nur gleichförmiges Abtasten im Winkel sowie auch (fast ausschließlich) äquidistantes, vom Winkel unabhängiges Abtasten in lateraler Richtung, wir setzen

s_{ν,j} = s_{ν} = −1 − \frac{h}{2} + νh, ν = 1,…, 2q, h = \frac{1}{q} ,  (2.10)

oder auch

s_{ν,j} = s_{ν} = −1 − h + νh, ν = 1,…, 2q + 1, h = \frac{1}{q} ,  (2.11)

also damit insbesondere s_{q+1} = 0.

Sprachlich unterscheiden wir beide Fälle, indem wir von der Geometrie ohne (2.10) bzw. mit (2.11) Zentralstrahl sprechen.

Eine Ausnahme von der Winkelunabhängigkeit der lateralen Abtastung bildet die sogenannte Interlaced-Geometrie, bei der in Winkelrichtung abwechselnd jeweils die Messungen mit geradem bzw. ungeradem Index übersprungen werden. Arbeiten wir ohne Zentralstrahl, so setzen wir

s_{ν,j} = −1 − \frac{h}{2} + 2νh, ν = 1,…, q, j mod 2 = 0 ,  (2.12)

sowie

s_{ν,j} = −1 − \frac{h}{2} + (2ν − 1)h, ν = 1,…, q, j mod 2 = 1 ,  (2.13)

im anderen Fall

s_{ν,j} = −1 + (2ν − 1)h, ν = 1,…, q , j mod 2 = 0 ,  (2.14)
s_{ν,j} = −1 + 2(ν − 1)h, ν = 1,…, q + 1 , j mod 2 = 1 .  (2.15)

Wir wollen nun noch eine andere Interpretation der Daten herleiten. Substituieren wir in (2.6) x = sω_{j} + tω_{j}^{⊥}, so erhalten wir

g(φ_{j}, s_{ν,j}) = ∫_{D} f(x)d_{h}(< x, ω_{j} > −s_{ν,j})dx
= (f, L_{h,j,ν})_{L_{2}}  (2.16)

mit L_{h,j,ν}(x) = \{ \begin{matrix}d_{h}(< x, ω_{j} > −s_{ν,j}), & x ∈ D_{1}(0) \\ 0, & sonst \end{matrix} .

Eine kleine Rechnung zeigt die Abschätzung ∥ L_{h,j,ν} ∥_{L_{2}}≤ \sqrt{\frac{2}{h}} ∥ d ∥_{L_{2}}, wir können also die Daten auch, unter Vermeidung des Begriffs der Radon-Transformation, direkt als L_{2}-stetige lineare Funktionale der gesuchten Funktion f deuten.

2. Die Fächer-Geometrie

Bei dieser, vor allem in medizinischen Anwendungen (und dort fast ausschließlich) gebräuchlichen Geometrie liegen die folgenden Meßwerte vor:

g(β_{j}, α_{ν,j}) = ∫ ᾱ ; −ᾱ Df(r ( \begin{matrix}cos β_{j} \\ sin β_{j} \end{matrix} ) , β_{j} + α − π)d_{h}(α − α_{ν,j})dα,
mit Df(a, ψ) = ∫_{0}^{∞} f(a + t ( \begin{matrix}cos ψ \\ sin ψ \end{matrix} ))dt ,  (2.17)
r > 1 ,  (2.18)
ᾱ = arcsin(\frac{1}{r}) ,  (2.19)
β_{j} = j \frac{2π}{p} ,  (2.20)
α_{ν,j} ∈ [−ᾱ, ᾱ] .  (2.21)

Setzen wir (2.17) in (2.17) ein, substituieren

x = r ( \begin{matrix}cos β_{j} \\ sin β_{j} \end{matrix} ) + t ( \begin{matrix}cos(β_{j} + α − π) \\ sin(β_{j} + α − π) \end{matrix} )  (2.22)

und schreiben a_{j} = r ( \begin{matrix}cos β_{j} \\ sin β_{j} \end{matrix} ), so erhalten wir

g(β_{j}, α_{ν,j}) = ∫_{D} f(x) ∥ x − a_{j} ∥^{−1} ·  (2.23)

d_{h}(sign(< x, a_{j}^{⊥} >) arccos(< \frac{x − a_{j}}{∥ x − a_{j} ∥}, −\frac{a_{j}}{r} >) − α_{ν,j})dx.

(2.24)

Die Daten sind also auch in diesem Fall als stetige lineare Funktionale von f zu interpretieren:

g(β_{j}, α_{ν,j}) = (f, C_{h,j,ν})_{L_{2}}  (2.25)
mit C_{h,j,ν}(x) = \{ \begin{matrix}^{\frac{d_{h}(± arccos(< \frac{x−a_{j}}{∥x−a_{j} ∥} ,− \frac{a_{j}}{r} >)−α_{ν,j})}{∥x−a_{j}∥}}, & x ∈ D \\ 0, & sonst \end{matrix} .

(2.26)

wobei das Vorzeichen mit dem von < x, a_{j}^{⊥} > übereinstimmt. Wie im parallelen Fall wählt man die α_{ν,j} meist unabhängig von j und äquidistant im Intervall [−ᾱ, ᾱ], wobei man sich dafür entscheiden kann, den Zentralstrahl durch den Ursprung gehen zu lassen (d.h. α_{q+1} = 0, bei 2q + 1 Messungen pro Richtung) oder ihn aber um eine halbe oder auch viertel Detektorbreite zu versetzen.

Da wir uns intensiv mit der Abbildung

f ↦ g(φ_{j}, s_{ν,j}) , j = 0,…, 2p − 1 , ν = 1,…, q(j) ,  (2.27)

mit g wie in (2.8) bzw. äquivalent dazu (2.16), beschäftigen wollen, geben wir ihr einen eigenen Namen und ein eigenes Symbol.

Definition 2.1.1 Wir nennen die Abbildung (2.27) die diskrete Radon-Transformation und schreiben dafür

g = Rf .  (2.28)

Den Vektor g ∈ IR^{N}, N = ∑_{j=0}^{2p−1} q(j), werden wir gelegentlich in

g = ( g^{0} ; ⋮ ; g^{2p−1} )  (2.29)

aufteilen. Die Vektoren g^{j} ∈ IR^{q(j)} sind dann über

g_{ν}^{j} = g(φ_{j}, s_{ν,j}), ν = 1,…, q(j)  (2.30)

definiert. Für die zugehörige Abbildung schreiben wir

g^{j} = R_{j}f ,  (2.31)

womit R also die Aufteilung

R = ( R_{0} ; ⋮ ; R_{2p−1} )  (2.32)

erfährt.

Eine kleine Rechnung liefert uns die Adjungierten von R bzw. R_{j}, es gilt:

(R_{j}^{∗}y)(x) = ∑_{ν=1}^{q(j)} y_{ν}L_{h,j,ν}(x) ,  (2.33)
R^{∗}g(x) = ∑_{j=0}^{2p−1}(R_{j}^{∗}g^{j})(x)
= ∑_{j=0}^{2p−1} ∑_{ν=1}^{q(j)} g_{ν}^{j}L_{h,j,ν} .  (2.34)

Natürlich ist im Fall äquidistanter Abtastung in lateraler Richtung die Hälfte der Daten redundant, so daß tatsächlich nur g^{0}, ⋯, g^{p−1} gemessen wird. Die formale Ergänzung der Daten auf dem Vollkreis erweist sich jedoch für viele Anwendungen als nützlich. Um diese Redundanz auch für die Interlaced-Geometrie zu ”sichern“, ist jedoch eine Paritätsforderung an p zu stellen (s. Kapitel 3, S. 40).

Wir können das Inversionsproblem jetzt wie folgt formulieren: Gegeben sei

g = Rf mit f ∈ L_{2}(D) , f(x) ∈ IR ,  (2.35)

gesucht ist f.

Da R selbstverständlich nicht injektiv ist, suchen wir natürlich eigentlich eine möglichst gute Approximation an f.

In [14] wird nachgewiesen, daß der für endliche Mengen A nichttriviale Unterraum

{f ∈ L_{2}(D) : Rf(ω, ·) = 0 , ω ∈ A}  (2.36)

ausschließlich aus hochfrequenten Funktionen besteht. Wir werden daher stets, im Sinne einer regularisierenden a priori-Information, von f verlangen, annähernd b-bandbe-schränkt zu sein.

Das heißt wir nehmen an: ∃ ε_{f} << 1 , mit

∫_{∥ξ∥>b} |f̂(ξ)|^{2}dξ ≤ ε_{f}^{2} ∥ f ∥_{L_{2}}^{2} .  (2.37)

2.2 Zurückführung auf Punktauswertungsfunktionale

In der bestehenden Literatur über Auflösungsfragen bei der Inversion der Radon-Trans-formation geht man davon aus, daß die Daten als Punktauswertungen der Radon-Transformierten der gesuchten Funktion f vorliegen.

Wir wollen zeigen, daß auch in unserem Fall ein ähnlicher Standpunkt eingenommen werden kann, jedoch auch auf die Schwierigkeiten hinweisen, die mit der dann möglichen Übernahme der bisherigen Theorie einhergehen.

Neben den beiden bisher gewonnenen Darstellungen

g(φ_{j}, s_{ν}) = R_{φ_{j}}f ∗ d_{h}(s_{ν}), bzw.  (2.38)
g(φ_{j}, s_{ν}) = (f, L_{h,j,ν})_{L_{2}}  (2.39)

zeigen wir nun eine dritte Interpretationsmöglichkeit auf. Ausgangspunkt ist dabei die Identität ([20], S.13):

R_{φ}(f ∗ g) = R_{φ}f ∗ R_{φ}g.  (2.40)

Ein Blick auf (2.38) und (2.40) zeigt: Finden wir eine Funktion D_{h}, so daß d_{h} = R_{φ}D_{h} ∀φ gilt, so besteht die Identität

g(φ_{j}, s_{ν}) = R_{φ_{j}}f ∗ d_{h}(s_{ν}) = R_{φ_{j}}(f ∗ D_{h})(s_{ν}) ,  (2.41)

d.h. die Daten lassen sich als Punktauswertungen der Radon-Transformierten von f ∗D_{h} deuten. Wir wollen nun herleiten, wie D_{h} mit der Detektorfunktion d_{h} zusammenhängt und benutzen dazu zweckmäßig die Cormacksche Inversionsformel für den Fourierkoeffizienten der Ordnung 0 ([20], S. 28). Da R_{φ}D_{h} = d_{h} unabhängig von φ ist, gilt dies auch für D_{h}; schreiben wir D_{h}(x) = e_{h}(∥ x ∥), so besteht die Beziehung

e_{h}(r) = −\frac{1}{π} ∫_{r}^{∞} \frac{d′_{h}(s)}{\sqrt{s^{2} − r^{2}}} ds .  (2.42)

Insbesondere gilt also e_{h}(r) = 0 ∀r > \frac{h}{2}.

Wir wollen dieses Ergebnis im folgenden Satz festhalten:

Satz 2.2.1 Mit e_{h}(r) = −\frac{1}{π} ∫_{r}^{∞} \frac{d′_{h}(s)}{\sqrt{s^{2}−r^{2}}} ds, sowie D_{h}(x) = e_{h}(∥ x ∥) gilt:

R_{φ}f ∗ d_{h} = R_{φ}(f ∗ D_{h}) .  (2.43)

Beispiel:

Für das Rechteckprofil

d_{h}(x) = \{ \begin{matrix}\frac{1}{h}, & |x| < \frac{h}{2} \\ 0, & sonst \end{matrix}  (2.44)

ergibt sich mit d′_{h} = \frac{1}{h}(δ_{− \frac{h}{2}} − δ_{\frac{h}{2}} ):

D_{h}(x) = \{ \begin{matrix}\frac{1}{πh} (\frac{h^{2}}{4} − ∥ x ∥^{2})^{− \frac{1}{2}} & , ∥ x ∥< \frac{h}{2} \\ 0 & sonst \end{matrix} .  (2.45)

Für das glattere Parabelprofil

d_{h}(x) = \{ \begin{matrix}\frac{6}{h^{3}} (\frac{h^{2}}{4} − x^{2}) & , |x| < \frac{h}{2} \\ 0 & sonst \end{matrix}  (2.46)

berechnen wir:

D_{h}(x) = \{ \begin{matrix}\sqrt[\frac{12}{πh^{3}}]{\frac{h^{2}}{4} − ∥ x ∥^{2}} & , ∥ x ∥< \frac{h}{2} \\ 0 & , sonst \end{matrix} .  (2.47)

Wir haben damit die von uns vorgestellte Meßsituation formal auf Punktauswertungen der Radon-Transformierten von f ∗D_{h} mit bekannten Faltungskernen D_{h} zurückgeführt. Da die D_{h} für sinnvolle Detektorfunktionen d_{h} hinreichend vernünftig sind, kann man also hoffen, die nun zugänglichen Aussagen über die mögliche Auflösung von f ∗ D_{h} für ebensolche über f verwenden zu können.

2.3 Problematisierung des in 2.2 skizzierten Vorgehens

Es liegt also nahe, die von Natterer ([20]) entwickelten und von Kruse ([8]) und Faridani ([4]) ausgebauten Ideen für die Diskussion von Auflösungsfragen auch hier aufzugreifen.

Die in den o.a. Arbeiten verwendete Funktionenklasse der ”wesentlich b-bandbeschränk-ten“ Funktionen wird über die Größe

ϵ_{0}(f, b) = ∫_{∥ξ∥>b} |f̂(ξ)|dξ  (2.48)

charakterisiert; man nennt eine Funktion wesentlich bandbeschränkt, wenn ϵ_{0} numerisch ”sehr klein“ ist.

Betrachten wir also

ϵ_{0}(f ∗ D_{h}, b) = (2π)^{\frac{1}{2}} ∫_{∥ξ∥>b} |f̂(ξ)| |D̂_{h}(ξ)|dξ  (2.49)
= ∫_{∥ξ∥>b} |f̂(ξ)| |d̂_{h}(∥ ξ ∥)|dξ .  (2.50)

Unter unseren bisherigen Annahmen macht dieser Ausdruck nicht notwendig einen Sinn, von f wollen wir nur Quadratintegrierbarkeit verlangen, und für das Rechteckprofil ist D_{h} ∉ L_{2}. Obiges Integral wird also i.a. nicht existieren.

Natürlich lassen sich auch glattere Detektorfunktionen (wie z.B. das Parabelprofil in (2.46)) proklamieren, so daß D̂_{h} schneller fällt und zumindest aus L_{2} ist. Es bleibt dann aber immer noch die folgende, spezifische Schwierigkeit: Wir wollen natürlich die zu verarbeitende Datenmenge, repräsentiert durch \frac{1}{h}, nicht unabhängig von der gewünschten Auflösung, repräsentiert durch b, wählen. Vielmehr ist von einer linearen Relation b = \frac{c}{h} auszugehen. Die oben zitierten Arbeiten lassen dabei, durch das Abtasttheorem begründet, an eine Konstante von c ≈ π denken. Die Cauchy-Schwarzsche Ungleichung ergibt nun:

ϵ_{0}(f ∗ D_{h}, \frac{c}{h}) ≤ ∥ f̂ ∥_{L_{2}(DC ; \frac{c}{h} (0))} · ∥ d̂_{h}(∥ · ∥) ∥_{L_{2}(DC ; \frac{c}{h} (0))} .  (2.51)

Für den zweiten Faktor berechnen wir:

∫_{∥ξ∥> \frac{c}{h}} |d_{h}(∥ ξ ∥)|^{2}dξ = 2π ∫_{ρ> \frac{c}{h}} |d̂_{h}(ρ)|^{2}ρ dρ  (2.52)
= \frac{2π}{h^{2}} ∫_{σ>c} |d̂(σ)|^{2}σdσ .  (2.53)

Zwar existiert das Integral für hinreichend glatte Detektorfunktionen (für das parabelförmige Profil ergibt sich z.B. d̂(σ) = const · σ^{− \frac{3}{2}} J_{\frac{3}{2}} (\frac{1}{2}σ)), wir verbleiben jedoch mit dem störenden Faktor \frac{1}{h} vor ∥ f̂ ∥_{L_{2}(D_{c/h}^{C}(0))}, der Größe, die wir als klein im Sinne einer regularisierenden Annahme ansehen wollen.

Auch sind dem Versuch, das Problem auf die Wahl einer ”guten“ Detektorfunktion zu verschieben, theoretische Grenzen gesetzt. Wünschenswert wäre ein d_{h}, welches einerseits gänzlich auf [−\frac{h}{2}, \frac{h}{2}] konzentriert ist (also supp(d_{h}) ⊂ [−\frac{h}{2}, \frac{h}{2}]), dessen Fouriertransformierte andererseits weitestgehend auf [−\frac{c}{h}, \frac{c}{h}] konzentriert ist.

In der bereits zitierten Arbeit [9] wird gezeigt, daß der höchstmögliche Konzentrationsgrad im Sinne der L_{2}-Norm vom Parameter c abhängt und für c = π nur ungefähr \sqrt{0.5} beträgt. Wir werden deshalb den in Kapitel 2.2 beschriebenen Weg nicht weiter verfolgen und statt dessen versuchen die Daten weiterhin als stetige Funktionale auf den Hilberträumen L_{2}(D) bzw. L_{2}([−1, 1], w^{−1}) zu behandeln.

3 Auflösungsuntersuchungen

Wir werden in diesem Kapitel einen (hypothetischen) Rekonstruktionsalgorithmus für unser in (2.35) formuliertes Problem untersuchen, dessen theoretische Grundlage die Singulärwertzerlegung (1.13) der Radon-Transformation bildet. Ziel wird es sein, hinreichend viele Entwicklungskoeffizienten der gesuchten, annähernd b-bandbeschränkten Funktion aus den Meßdaten zu berechnen. Wir wollen dabei aber die in Kapitel 2.1 beschriebenen Meßanordnungen genau berücksichtigen. Dies bedeutet u.a. auch, daß Satz 1.1.4 für das Rekonstruktionsproblem nicht direkt angewendet werden kann, da die Funktion g = Rf selbst gar nicht bekannt ist, somit auch das Skalarprodukt

(g, g_{ml})_{L_{2}(Z,w^{−1})}  (3.1)

nicht exakt berechnet werden kann. Wird (3.1) durch Näherungsformeln, wie etwa die Trapezregel ersetzt, so treten Fehler auf, die zu berücksichtigen sind. Die Idee, die Singulärwertzerlegung als Hilfsmittel bei der Diskussion von Auflösungsfragen zu verwenden, ist nicht neu, siehe [12]. Dort wird allerdings ein ganz anderer Ansatz verfolgt.

Zunächst werden wir nachweisen, daß sich das Problem (2.35) ohne allzu großen Genauigkeitsverlust auf die verallgemeinerte Lösung endlich vieler endlicher linearer Gleichungssysteme zurückführen läßt. Die Unbekannten x_{m,l} sind dabei die Entwicklungskoeffizienten von f nach den Singulärfunktionen, also x_{m,l} = (f, f_{ml})_{L_{2}(D)}. Wir werden allerdings zunächst versuchen, die Gleichungssysteme nach den Größen

z_{m,l} = σ_{m} x_{m,l}

aufzulösen. Daß wir danach in einem zweiten Schritt beim Übergang von z_{m,l} zu den eigentlich gesuchten x_{m,l} durch die (eine Nullfolge bildenden) Singulärwerte σ_{m} = 2\sqrt{\frac{π}{m+1}} zu dividieren haben, ist Ausdruck der Schlechtgestelltheit des Problems. Die damit einhergehende Fehlerverstärkung läßt sich nicht vermeiden.

3.1 Die Standard-Geometrie

Zunächst gilt nach (1.13) im Sinne der L_{2}-Konvergenz:

Rf(φ, s) = \frac{1}{π} ∑_{m=0}^{∞} σ_{m} ∑_{m=|l|}′(f, f_{ml})_{L_{2}(D)} u_{m}(s)e^{ilφ}  (3.2)
= ∑_{|l|=0}^{∞} e^{ilφ} \underbrace{\frac{1}{π} ∑_{m=|l|}^{∞}′(f, f_{ml})_{L_{2}(D)} σ_{m}u_{m}(s)}_{=: h_{l}(s)} .  (3.3)

Leicht einzusehen sind die folgenden Eigenschaften von h_{l}:

h_{l}(−s) = (−1)^{l}h_{l}(s) ,  (3.4)
h_{−l} = \overline{h_{l}} .  (3.5)

Die Gleichung (3.4) ist Konsequenz aus Rf(φ+π, s) = Rf(φ, −s) bzw. gleichbedeutend mit der Paritätsbedingung an l und m in (1.13), (3.5) folgt aus f(x) ∈ IR.

Die Reihe konvergiert für festes φ gegen eine Funktion R_{φ}f ∈ L_{2}([−1, 1], w^{−1}), so daß wir (3.3) in Winkelrichtung auch punktweise lesen dürfen, d.h. für φ = φ_{j} = \frac{jπ}{p} gilt:

Rf(φ_{j}, ·) = ∑_{l∈ZZ} e^{ilj \frac{π}{p}} h_{l}  , j = 0,…, p − 1 .  (3.6)

Für |l′| ≤ p berechnen wir die folgenden Koeffizienten (s. (3.7)) der diskreten Fouriertransformation der Länge 2p. Hierfür gibt es gleich zwei Gründe. Zum einen können wir die Fourierkoeffizienten als Diskretisierung des Integrals ∫_{0}^{2π} e^{−il′φ}g(φ, s) dφ durch die Trapezregel zu den Punkten φ_{j} auffassen. Dieses Integral ist ein ”Teil“ des zu berechnenden Skalarproduktes (3.1). Zum anderen, und dies ist für uns von größerem Interesse, entspricht dies exakt dem Vorgehen bei dem in Kapitel 4 beschriebenen und von uns untersuchten Algorithmus. Die Begründung für die Berechnung der Fourierkoeffizienten ist dort allerdings eine (vordergründig) vollkommen andere, sie ist notwendig, um ein sehr großes lineares Gleichungssystem auf Blockdiagonalgestalt zu transformieren.

g_{l′}(s) := \frac{1}{2p} ∑_{j=0}^{2p−1} e^{−il′j \frac{π}{p}} Rf(φ_{j}, s)  (3.7)
= ∑_{|l|=0}^{∞} \underbrace{\frac{1}{2p} ∑_{j=0}^{2p−1} e^{i(l−l′)j \frac{π}{p}}}_{=\{ \begin{matrix}1, & l − l′ = 0 mod 2p \\ 0, & sonst \end{matrix}} h_{l}(s)

= h_{l′}(s) + ∑_{k≠0} h_{l′+2pk}(s) .  (3.8)

Diese Gleichung wird Ausgangspunkt für unsere Untersuchungen werden: Nach erfolgter diskreter (Block-) Fouriertransformation der Datenvektoren (g(φ_{j}, ·) ∗ d_{h}(s_{ν}))_{ν=1,...,2q} bezüglich des Winkelindex j haben wir Zugang zu den Größen

g_{l′} ∗ d_{h}(s_{ν})  , \begin{matrix}l′ = −p,…, p \\ ν = 1,…, 2q \end{matrix} ,  (3.9)

und (3.8) wird zu einem unterbestimmten System für die gesuchten Koeffizienten (f, f_{ml′})_{L_{2}(D)}:

g_{l′} ∗ d_{h}(s_{ν}) = h_{l′} ∗ d_{h}(s_{ν}) + (∑_{k≠0} h_{l′+2pk}) ∗ d_{h}(s_{ν}) , \begin{matrix}|l′| ≤ p , \\ ν = 1,…, 2q . \end{matrix}  (3.10)

Um daraus ein endliches System zu machen, schneiden wir die Reihenentwicklung der h_{l′} jeweils nach dem K-ten Koeffizienten ab, definieren demnach

h_{l′,K} := \frac{1}{π} ∑_{m=|l′|}^{K}′(f, f_{ml′})σ_{m}u_{m} ,  (3.11)
sowie  r_{1} := ∑_{k≠0} h_{l′+2pk} ,  (3.12)
und  r_{2} := \frac{1}{π} ∑_{m=K+1}^{∞}′(f, f_{ml′})σ_{m}u_{m} ,  (3.13)

womit (3.10) zu dem nun endlichen (aber durch den Fehler der Trapezregel r_{1} und den Abbruchfehler r_{2} gestörten) System

\begin{matrix}\frac{1}{π} ∑_{m=|l′|}^{K}′σ_{m}(f, f_{ml′})d_{h} ∗ u_{m}(s_{ν}) = g_{l′} ∗ d_{h}(s_{ν}) + (r_{1} + r_{2}) ∗ d_{h}(s_{ν}) \\ |l′| ≤ p, ν = 1,…, 2q \end{matrix}  (3.14)

für die Unbekannten z_{m,l′} = σ_{m}(f, f_{ml′}); m = |l′|,…, K; |l′| ≤ p; m + l′ gerade, wird. Dieses zerfällt offensichtlich in 2p+1 Einzelsysteme mit jeweils 2q Gleichungen. Da h_{l′,K} (bis auf den Faktor w) ein Polynom mit der Parität von l′ ist, sind in (3.14) tatsächlich aber jeweils nur q Gleichungen zu berücksichtigen: für ν = 1,…, q ergibt sich die Gleichung für −s_{ν} = s_{2q+1−ν} aus derjenigen für s_{ν} durch Multiplikation mit (−1)^{l′}. Die Anzahl der Unbekannten beträgt [\frac{K−|l′|}{2}] + 1. Desweiteren ist offensichtlich, daß man durch Addition bzw. Subtraktion der Gleichungen für l′ und −l′ in (3.14) zu zwei getrennten Gleichungssystemen für den Real- bzw. Imaginärteil von z_{m,l′} kommt. Soll (3.14) für ein Rekonstruktionsverfahren verwendet werden, so sind vorab zwei Punkte zu klären:

1.) Die Fehler r_{1}, r_{2} – die Abhängigkeit der Fehler von l′ und p bzw. K wird zwar durch diese Notation unterschlagen, dies sollte jedoch nicht zu Verständnisschwierigkeiten führen – müssen ∀l′ für eine geeignete Wahl von p, K klein sein.

2.) Das Gleichungssystem muß hinreichend stabil auflösbar sein.

Wir wollen diese Punkte nun der Reihe nach behandeln.

3.1.1 Abschätzungen für die Fehler r_{1} und r_{2}

Lemma 3.1.1 Es gilt mit p = ϑ^{−1}b , 0 < ϑ < 1

∥r_{1}∥_{L_{2}([−1,1],w^{−1})} ≤ (\frac{2}{\sqrt{p}} ε_{f} + \frac{c}{p}q0^{\frac{p+1}{2}} ) ∥f∥

mit q_{0} = e^{− \frac{2}{3} (1−ϑ^{2})^{\frac{3}{2}}} und einer nur von ϑ abhängigen Konstanten c.

Beweis:

∥r_{1}∥_{L_{2}([−1,1],w^{−1})}^{2} = ∥\frac{1}{π} ∑_{k≠0} ∑_{m≥|l′+2kp|}′(f, f_{m,l′+2kp})σ_{m}u_{m}∥_{L_{2}([−1,1],w^{−1})}^{2}
= ∥\frac{1}{π} ∑_{m=c_{l′}}^{∞} σ_{m} u_{m} ∑_{{|l|≤m \atop l=l′ mod2p}}′(f, f_{ml}) δ_{l,l′}∥_{L_{2}([−1,1],w^{−1})}^{2} , (3.15)
mit c_{l′} := min(|l′ + 2p|, |l′ − 2p|) ≥ p, da |l′| ≤ p .

Mithin folgt:

∥r_{1}∥_{L_{2}([−1,1],w^{−1})}^{2} = ∑_{m=c_{l′}}^{∞} \frac{σ_{m}^{2}}{π^{2}} ∫_{−1}^{1} U_{m}^{2}(s)w(s)ds | ∑_{{|l|≤m \atop l=l′ mod2p}}′ δ_{l,l′}(f, f_{ml})|^{2}
= 2 ∑_{m=c_{l′}}^{∞} \frac{1}{m + 1} | ∑_{{|l|≤m \atop l=l′ mod2p}}′ δ_{l,l′}(f, f_{ml})|^{2}
≤ 2 ∑_{m=c_{l′}}^{∞} \frac{\frac{m}{p} + 1}{m + 1} ∑_{{|l|≤m \atop l=l′ mod2p}}′ |(f, f_{ml})|^{2}
≤ \frac{4}{p} ∑_{m=c_{l′}}^{∞} ∑_{{|l|≤m \atop l=l′ mod2p}}′ |(f, f_{ml})|^{2} .  (3.16)

Dabei haben wir im vorletzten Schritt die Anzahl der Summanden in der inneren Summe abgezählt und die Cauchy-Schwarzsche Ungleichung angewendet.

Nun nutzen wir die a priori-Voraussetzung (2.37) an f aus:

∥r_{1}∥_{L_{2}([−1,1],w^{−1})} ≤ \frac{2}{\sqrt{p}} ( ∑_{m=c_{l′}}^{∞} ∑_{{|l|≤m \atop l=l′ mod2p}}′ |(f − Bbf, f_{ml}) + (Bbf, f_{ml})|^{2})^{\frac{1}{2}}

≤ \frac{2}{\sqrt{p}} ( ∑_{m=c_{l′}}^{∞} ∑_{{|l|≤m \atop l=l′ mod2p}}′ |(f − Bbf, f_{ml})|^{2})^{\frac{1}{2}}

+\frac{2}{\sqrt{p}} ( ∑_{m=c_{l′}}^{∞} ∑_{{|l|≤m \atop l=l′ mod2p}}′ |(f, Bbf_{ml})|^{2})^{\frac{1}{2}} =: S_{1} + S_{2} . (3.17)

Die erste Summe schätzen wir mit der Besselschen Ungleichung ab:

S_{1} ≤ \frac{2}{\sqrt{p}}∥f − Bbf∥ ≤ \frac{2}{\sqrt{p}} ε_{f} ∥f∥ .  (3.18)

Für die zweite Summe nutzen wir zunächst die Isomorphie-Eigenschaft der Fouriertransformation aus:

S_{2} = \frac{2}{\sqrt{p}} ( ∑_{m=c_{l′}}^{∞} ∑_{{|l|≤m \atop l=l′ mod2p}}′ |(f̂, B̂bf_{ml})|^{2})^{\frac{1}{2}}

≤ \frac{2}{\sqrt{p}} ∥f∥ ( ∑_{m=c_{l′}}^{∞} ∑_{{|l|≤m \atop l=l′ mod2p}}′ ∥B̂bf_{ml}∥^{2})^{\frac{1}{2}} ,  (3.19)

mit B̂bf_{ml} = f̂_{ml} · χ_{\frac{1}{b}} . Wir müssen also zunächst f̂_{ml} berechnen, was am einfachsten über den Projektionssatz geschieht:

f̂_{ml}(σ · w) = (2π)^{− \frac{1}{2}} (R_{φ}f_{ml}) (σ)
= (2π)^{− \frac{1}{2}} (RR^{∗}g_{ml} \frac{1}{σ_{m}}) (φ, σ)
= \frac{σ_{m}}{\sqrt{2π}} ĝ_{m,l}(φ, σ)
= \frac{σ_{m}}{2π^{2}} ∫_{−1}^{1} e^{−isσ}(1 − s^{2})^{\frac{1}{2}} U_{m}(s)ds e^{ilφ}
= \sqrt{m + 1}\frac{i^{m}}{\sqrt{π}} \frac{J_{m+1}(σ)}{σ}e^{ilφ} .  (3.20)

Dabei haben wir von Formel 7.321 in [6] Gebrauch gemacht.

⟹ ∥B̂f_{ml}∥^{2} = \frac{m + 1}{π} ∫_{0}^{b}∫_{0}^{2π} dφ \frac{J_{m+1}^{2}(σ)}{σ}dσ

= 2(m + 1) ∫_{0}^{b} \frac{J_{m+1}^{2}(σ)}{σ}dσ .  (3.21)

Aus ([15], S. 256) entnehmen wir nun

∫_{0}^{b} \frac{J_{m+1}^{2}(σ)}{σ}dσ = \frac{1}{m + 1} ∑_{k=m+2}^{∞} J_{k}^{2}(b) + \frac{1}{2(m + 1)}J_{m+1}^{2}(b) .  (3.22)

Wir erinnern uns nun daran, daß die Reihenentwicklung in S_{2} mit c_{l′} ≥ p beginnt, wir also die Debyesche Beziehung für die gewünschte Abschätzung heranziehen können, wenn wir sicherstellen, daß m ≥ c_{l′} ≥ p ≥ ϑ^{−1} · b gilt. Setzen wir also p = \frac{1}{ϑ} b , dann gilt (s. (1.24))

J_{k}^{2} = J_{k}^{2} (\frac{p}{k} ϑ · k) ≤ \frac{1}{2πk} \frac{1}{(1 − (\frac{pϑ}{k})^{2})^{\frac{1}{2}}} e^{−k \frac{2}{3}(1−( \frac{pϑ}{k})^{2})^{\frac{3}{2}}}  (3.23)

und somit

∑_{k=m+2}^{∞} J_{k}^{2}(b) ≤ \frac{1}{2π} ∑_{k=m+2}^{∞} \frac{1}{k} \frac{q^{k}(k)}{(1 − (\frac{pϑ}{k})^{2})^{\frac{1}{2}}}  (3.24)
mit q(k) = e^{− \frac{2}{3} (1−( \frac{pϑ}{k})^{2})^{\frac{3}{2}}} .

Wegen m ≥ p gilt q(k) < q_{0} = e^{− \frac{2}{3} (1−ϑ^{2})^{\frac{3}{2}}} , k ≥ m + 2 , was zur Abschätzung

∑_{k=m+2}^{∞} J_{k}^{2}(b) ≤ \frac{1}{2π(m + 1)(1 − ϑ^{2})^{\frac{1}{2}}} · \frac{q_{0}^{m+2}}{1 − q_{0}}  (3.25)

führt. Dies verwenden wir in (3.22) und erhalten

∫_{0}^{b} \frac{J_{m+1}^{2}(σ)}{σ}dσ ≤ \frac{1}{4π(1 − ϑ^{2})^{\frac{1}{2}}} \frac{1}{(m + 1)^{2}} \frac{1 + q_{0}}{1 − q_{0}} q_{0}^{m+1} ,  (3.26)

um damit schließlich für S_{2} ( s. (3.19) u. (3.21)) die Abschätzung

S_{2} ≤ \frac{2}{\sqrt{p}} ∥f∥ \sqrt{\frac{1 + q_{0}}{1 − q_{0}}} (\sqrt{2π}(1 − ϑ^{2})^{\frac{1}{4}} )^{−1} ( ∑_{m=c_{l′}}^{∞} \frac{1}{m + 1} ∑_{{|l|≤m \atop l=l′ mod2p}}′ q_{0}^{m+1})^{\frac{1}{2}}

≤ \frac{2}{p}\sqrt{2} ∥f∥ \sqrt{\frac{1 + q_{0}}{1 − q_{0}}} (\sqrt{2π}(1 − ϑ^{2})^{\frac{1}{4}} )^{−1} ( ∑_{m=c_{l′}}^{∞} q_{0}^{m+1})^{\frac{1}{2}}
≤ \frac{∥f∥}{p} c q0^{\frac{p+1}{2}}  (3.27)

zu gewinnen. Die Kombination von (3.18) und (3.27) beendet den Beweis.

□

Den Fehler r_{2} können wir nun mit ganz ähnlichen Techniken abschätzen.

Lemma 3.1.2 Es gilt mit K = p = \frac{1}{ϑ}b , q_{0} wie in Lemma 3.1.1 und einer nur von ϑ abhängigen Konstanten c

∥r_{2}∥_{L_{2}([−1,1],w^{−1})} ≤ (\sqrt{\frac{2}{p}} ε_{f} + \frac{c}{p}q0^{\frac{p+1}{2}} ) ∥f∥_{L_{2}} .

Beweis:

∥r_{2}∥_{L_{2}([−1,1],w^{−1})} = \frac{1}{π} ( ∑_{m=K+1}^{∞}′ |(f, f_{ml′})|^{2}σ_{m}^{2} ∫_{−1}^{1} U_{m}^{2}(s) w(s) ds)^{\frac{1}{2}}
= \sqrt{2} ( ∑_{m=K+1}^{∞}′ |(f, f_{ml′})|^{2} \frac{1}{m + 1})^{\frac{1}{2}}
≤ \sqrt{\frac{2}{p}} ( ∑_{m=K+1}^{∞}′ |(f − Bbf, f_{ml′})|^{2})^{\frac{1}{2}} + \sqrt{\frac{2}{p}} ( ∑_{m=K+1}^{∞}′ |(f, Bbf_{ml′})|^{2})^{\frac{1}{2}}
=: S_{1} + S_{2}  (3.28)

Für S_{1} erzielen wir wie im Lemma 3.1.1 die Abschätzung

S_{1} ≤ \sqrt{\frac{2}{p}} ε_{f} ∥f∥ .

Analog zu (3.19) erhalten wir hier für S_{2}:

S_{2} ≤ \begin{matrix}√ & \overline{2} \\ √ & \overline{p} \end{matrix}∥f∥ ( ∑_{m=K+1}^{∞} ∥Bf_{ml′}∥^{2})^{\frac{1}{2}}  (3.29)
≤ \begin{matrix}√ & \overline{2} \\ √ & \overline{p} \end{matrix}∥f∥\sqrt{\frac{1 + q_{0}}{1 − q_{0}}} (\sqrt{2π}(1 − ϑ^{2})^{\frac{1}{4}} )^{−1} ( ∑_{m=K+1}^{∞} \frac{1}{m + 1}q_{0}^{m+1})^{\frac{1}{2}}  (vgl. (3.27))
≤ c \frac{∥f∥}{p} q0^{\frac{p+1}{2}} .  (3.30)
□

Ein numerisches Beispiel zeigt, daß bei realistischen Werten von b und ϑ die Fehler r_{1} und r_{2} mit hinreichend kleinen Faktoren gegen ∥f∥ abgeschätzt werden können.

Beispiel: Sei b = 400, ϑ = 0.90, sowie K = p = 445. Dann ergibt die explizite Berechnung der Konstanten in (3.27)

S_{2} ≤ 5.9 · 10^{−7} .

Wir sind also der Klärung des ersten Punktes ein gutes Stück nähergekommen.

Wir berechnen noch die l_{2}-Vektornorm des Residuenvektors r_{ν}^{D} := (r_{1} + r_{2}) ∗ d_{h}(s_{ν}), ν = 1,…, 2q; für j = 1, 2 gilt ( mit d_{h_{ν}} = d_{h}(· − s_{ν}), I_{ν} = [s_{ν} − \frac{h}{2}, s_{ν} + \frac{h}{2}]):

∑_{ν=1}^{2q} |(r_{j}, d_{h_{ν}})|^{2} = ∑_{ν=1}^{2q} |(r_{j} · w^{− \frac{1}{2}} , w^{\frac{1}{2}} d_{h_{ν}})|^{2}
≤ ∑_{ν=1}^{2q} ∥r_{j_{|_{Iν}}} ∥_{L_{2}(w^{−1})}^{2} · ∥d_{h}∥_{L_{2}(w)}^{2}
= ∥r_{j}∥_{L_{2}(w^{−1})}^{2} · ∥d_{h}∥_{L_{2}(w)}^{2}
≤ h^{−1}∥r_{j}∥_{L_{2}(w^{−1})}^{2} · ∥d∥_{L_{2}(w)}^{2} .  (3.31)

Demnach, unter Verwendung von Lemma 3.1.1 und Lemma 3.1.2:

∥r^{D}∥ ≤ \sqrt{h^{−1}} (∥r_{1}∥_{L_{2}(w^{−1})}^{2} + ∥r_{2}∥_{L_{2}(w^{−1})}^{2})^{\frac{1}{2}} ∥d∥_{L_{2}}
≤ ∥f∥ (\frac{c̃}{\sqrt{hp}} ε_{f} + \frac{˜c̃}{p\sqrt{h}} q0^{\frac{p+1}{2}} ) .  (3.32)

Da hp ≈ π zu wählen ist, zeigt (3.32), daß bei allen Gleichungssystemen in (3.14) die Residuenvektoren als sehr klein angenommen werden können.

3.1.2 Stabilität

Wir kommen nun zum zweiten Punkt, der Stabilität der Gleichungssysteme in (3.14). Wir werden dabei einsehen müssen, daß die Systeme (l′ = −p,…, p)

\begin{matrix}\frac{1}{π} ∑_{m=|l′|}^{p}′σ_{m} (f, f_{ml′}) u_{m} ∗ d_{h}(s_{ν}) = g_{l′} ∗ d_{h}(s_{ν}) + r_{ν}^{D} \\ ν = 1,…, 2q \end{matrix}  (3.33)

für vernünftige Werte von q hoffnungslos instabil sind, d.h. obwohl r_{ν}^{D} verschwindend klein ist, können wir prinzipiell nicht damit rechnen, die Unbekannten z_{m,l′}, und damit erst recht nicht die Koeffizienten (f, f_{ml′}), mit genügender Genauigkeit berechnen zu können. Auf den ersten Blick scheint dies im eklatanten Widerspruch zur Praxis zu stehen, die ja die Möglichkeit(en) tomographischer Bildrekonstruktionen täglich neu aufzeigt. Wir können diesen scheinbaren Widerspruch auflösen, indem wir erkennen, daß es sich bei h_{l′} ∀l′ um eine nahezu b-bandbeschränkte Funktion handelt. Die ”approximative Dimension“ (s. Abbildung 1.1) dieser Funktionenmenge ist jedoch \frac{2b}{π} (für l′ = 0). Dies erzwingt, daß einige Linearkombinationen der u_{m} (nahezu) verschwinden müssen. Wir werden sehen, daß es gerade diese Linearkombinationen sind, welche den Nullraum bzw. die Singulärvektoren zu extrem kleinen Singulärwerten von (3.33) aufspannen. Alles was wir also bei einer ”Lösung“ von (3.33) zu tun haben, ist diese Linearkombinationen in Form einer abgeschnittenen Singulärwertzerlegung des Gleichungssystems wegzulassen.

Wir wollen zunächst einmal zeigen, daß (3.33) für genügend großes q eindeutig lösbar ist.

Lemma 3.1.3 Für 2q > p ist die einzige Lösung von

( ∑_{m=|l′|}^{p} α_{m} u_{m}) ∗ d_{h}(s_{ν}) = 0 ν = 1,…, 2q

der Nullvektor α_{m} = 0, m = |l′|,…, p .

Beweis: Mit h_{l′,p}(s) = ∑_{m=|l′|}^{p} α_{m} u_{m} gilt h_{l′,p} = w · q_{p}, wobei q_{p} ein Polynom vom Grade p ist.

Aus ∫_{I_{ν}} d_{h,ν}(s) h_{l′,p}(s) ds = 0 für ν = 1,…, 2q, folgt nach dem Mittelwertsatz der Integralrechnung die Existenz von 2q Nullstellen von h_{l′,p} und damit von q_{p}, mithin ist q_{p} = 0, und es folgt die Behauptung.

□

Wir wollen nun durch den Vergleich mit der klassischen Aufgabenstellung der Polynominterpolation zeigen, daß weit mehr als p Messungen pro Richtung erforderlich sind, um (3.33) stabil zu lösen. Dies gilt allerdings nur für den Fall äquidistant verteilter Knoten, wesentlich günstiger im Hinblick auf Stabilität ist es, die Knoten in die Nullstellen der Tschebyscheff-Polynome zu legen. Dies zeigt eine Möglichkeit auf, der drohenden Instabilität zu entgehen. Für den Fall, daß wir h_{l′,p} (punktweise) in den Tschebyscheff-Knoten kennen, lassen sich die Gleichungssysteme (3.33) stabil über die diskrete Orthogonalitätsrelation (s. [16], S. 312-313) der U_{m} lösen. Für die Tschebyscheff-Polynome 2. Art kann eine den Lemmata 1 u. 2 in ([16], S. 307 u. 311) analoge Formel bewiesen werden. Dies erlaubt dann sogar die Behandlung der Interlaced-Geometrie, die in diesem Fall mit der PET-Geometrie zusammenfällt ([20], S. 76). Da wir dafür immer noch \frac{p}{2} Messungen in lateraler Richtung benötigen,verfolgen wir dies hier und bei der Behandlung der Interlaced-Geometrie im nächsten Abschnitt nicht weiter.

Betrachten wir zunächst die Spezialfälle l′ = 0 und d_{h} = δ (d.h. Punktauswertung), dann ist nach Division durch w(s_{ν}) das Problem (3.33) die klassische Interpolationsaufgabe. A. Schönhage betrachtet in [21] die folgende Fragestellung: Seien

M = \{p ∈ P_{n}, |p(x_{ν})| ≤ 1, ν = 0,…, m\} und M_{n,m} = max_{{x_{0}≤ξ≤xm \atop p∈M}} |p(ξ)| ,

wobei m ≥ n gelte. Welche Schranke läßt sich für M_{n,m} angeben, wenn wir die x_{ν} äquidistant verteilen? Für n = m erhält er das deprimierende Ergebnis

Satz 3.1.4 Mit M_{n} := M_{n,n} und der Eulerschen Konstanten C gilt:

M_{n} ∼ \frac{2^{n+1}}{e · n · (lg n + C)} .

Bemerkung 3.1.5 M_{n} wird üblicherweise als Lebesgue-Konstante bezeichnet. Da die U_{m} auf [−1, 1] durch m + 1 beschränkt sind, bedeutet Satz 3.1.4 automatisch, daß auch die Entwicklungskoeffizienten des Polynoms nach den U_{m} sehr schlecht durch die Interpolationsaufgabe bestimmt sind.

Die Abschätzung für den allgemeinen Fall sind später von H. Ehlich und K. Zeller in [2] verbessert worden. Dort wird folgendes gezeigt:

Satz 3.1.6 Bei äquidistant verteilten Knoten (x_{ν} = −1 + \frac{2ν}{m}) besteht für n ≥ 2 die Abschätzung:

M_{n,m} ≤ (1 − ϱ)^{−1} ,

wenn immer m so groß gewählt wird, daß

ϱ := \frac{n^{2}(n^{2} − 1)}{6m^{2}} < 1 ist.

Wir sehen also, daß m ≥ \frac{n^{2}}{6} gewählt werden muß, um diesen Satz überhaupt anwenden zu können. Für m = \frac{1}{2}n^{2} ergibt der Satz eine Schranke von M_{n, \frac{n^{2}}{2}} < 3. Für die Interpolation in Tschebyscheff-Punkten (x_{ν} = \frac{cos 2ν−1}{2m}, ν = 1,…, 2m ) können wesentlich günstigere Resultate erzielt werden (s. [2]); für explizite Darstellungen der Lebesgue-Konstanten (im Fall m = n) siehe auch [3].

Die o.a. negativen Resultate sind für den Grenzfall d_{h} = δ(∉ L_{2}) hergeleitet worden. Daß sie sich bei Verwendung von mittelnden Detektorfunktionen d_{h} nicht verbessern lassen, zeigen Computerberechnungen. Nimmt man als sozusagen entgegengesetzten Extremfall die Mittelung über die Intervallänge h, also für d_{h} das Rechteckprofil, und berechnet die Singulärwerte der Systemmatrix von (3.33), also (u_{m} ∗ d_{h}(s_{ν}))_{{m=l′,l′+2,...,p \atop ν=1,...,q}} für verschiedene l′, p, q , so zeigt sich folgendes. Selbst im deutlich überbestimmten Fall ergeben sich (numerisch) nichttriviale ”Nullvektoren“, die erst dann verschwinden, wenn q in der Größenordnung von p^{2} (bei l′ = 0) gewählt wird, d.h. die enorme Schlechtgestelltheit der Interpolationsaufgabe spiegelt sich auch hier wieder. Wir wollen nun zeigen, daß diese ”Nullvektoren“ in den h_{l′,p} nicht vorkommen können und bei einer verallgemeinerten Lösung von (3.33) einfach weggelassen werden dürfen. Dazu dient als erster Schritt folgender

Satz 3.1.7 Mit den Bezeichnungen von Lemma 3.1.2 gilt:

∥ĥ_{l′,p}∥_{L_{2}([−b,b]^{c})} ≤ ((\sqrt{\frac{2}{p}} + \sqrt{\frac{2}{b}}) ε_{f} + \frac{c}{p} q0^{\frac{p+1}{2}} ) ∥f∥_{L_{2}} .

Beweis: Unmittelbar aus der Definition ergibt sich

h_{l} = \frac{1}{2π} ∫_{0}^{2π} Rf(φ, s) e^{−ilφ} dφ
⟹  ĥ_{l}(σ) = \frac{1}{2π} ∫_{0}^{2π} R̂_{φ}f(σ) e^{−ilφ} dφ
= (2π)^{− \frac{1}{2}} ∫_{0}^{2π} f̂(σω) e^{−ilφ} dφ ,  ω = \binom{cos φ}{sin φ}
⟹  ∫_{b}^{∞} |ĥ_{l}(σ)|^{2}dσ = \frac{1}{2π} ∫_{b}^{∞}|∫_{0}^{2π} f̂(σω) e^{−ilφ} dφ|^{2} dσ
≤ ∫_{b}^{∞}∫_{0}^{2π} |f̂(σω)|^{2} dφdσ
= ∫_{∥ξ∥>b} |f̂(ξ)|^{2} ∥ξ∥^{−1} dξ
≤ b^{−1} ε_{f}^{2} ∥f∥^{2}

⟹ ∥ĥ_{l}∥_{L_{2}([−b,b]^{c})} ≤ \sqrt{2} b^{− \frac{1}{2}} ε_{f} ∥f∥ .

Und nun mit r_{2} nach (3.13) und unter Verwendung von Lemma 3.1.2:

∥ĥ_{l′,p}∥_{L_{2}([−b,b]^{c})} = ∥ĥ_{l′} − r̂_{2}∥_{L_{2}([−b,b]^{c})}
≤ ∥ĥ_{l′}∥_{L_{2}([−b,b]^{c})} + ∥r_{2}∥_{L_{2}}
≤ ((\sqrt{\frac{2}{p}} + \sqrt{\frac{2}{b}}) ε_{f} + \frac{c}{p} q0^{\frac{p+1}{2}} ) ∥f∥_{L_{2}} .
□

Als nächstes bemerken wir, daß es sich bei h_{l′,p} um eine endliche Linearkombination der u_{m} handelt, wegen û_{m}(ξ) = O(ξ^{− \frac{3}{2}} ) ist es sinnvoll, folgende Annahme zu treffen: ∃ ε_{1} << 1, so daß

|ĥ_{l′,p}(σ)| ≤ \frac{ε_{1}}{1 + σ^{\frac{3}{2}}} ∀|σ| ≥ b .  (3.34)

Die Abschätzung (3.34) gilt erst recht für die Funktion

f_{1} := h_{l′,p} ∗ d_{h} ,  (3.35)

da lim_{σ→∞} |d̂_{h}(σ)| = 0 und wir sogar d̂_{h}(σ) = O(σ^{−1}) erwarten dürfen (s. auch (3.67)). Aus technischen Gründen wollen wir hier d̂_{h}(σ) = O(σ^{− \frac{1}{2}} ) verlangen, so daß für f_{1} aus (3.35) die Ungleichung:

|f̂_{1}(σ)| ≤ \frac{ε_{1}}{1 + σ^{2}} ∀|σ| ≥ b  (3.36)

besteht.

Der Satz 1.2.1 liefert nun auf f_{1}^{2} angewandt:

∑_{ν=1}^{2q} f_{1}^{2}(s_{ν}) = ∑_{ν∈ZZ} f_{1}^{2}(νh + \frac{h}{2})
= (2π)^{1/2} h^{−1} ∑_{m∈ZZ}(−1)^{m} f̂_{1}^{2}(\frac{2π}{h} m)
= h^{−1} (2π)^{1/2} f̂_{1}^{2}(0) + h^{−1} ∑_{m≠0}(−1)^{m} f̂_{1} ∗ f̂_{1}(\frac{2π}{h} m)
= h^{−1} ∥f_{1}∥_{L_{2}}^{2} + h^{−1} ∑_{m≠0}(−1)^{m} f̂_{1} ∗ f̂_{1}(\frac{2π}{h} m) .

Damit sind wir in der Lage, den folgenden Satz zu beweisen.

Satz 3.1.8 Für f_{1} aus (3.35) bestehe die Ungleichung

∑_{ν=1}^{2q} f_{1}^{2}(s_{ν}) < ε^{2} ,

dann gilt unter der Annahme (3.36) und falls h = \frac{π}{b} gewählt wird, die Abschätzung:

∥f_{1}∥_{L_{2}} ≤ c_{1} b^{− \frac{3}{2}} ε_{1} + \sqrt{ε_{1}^{2} (c_{2}b^{−3} + c_{3}b^{−2}) + ε^{2} \frac{π}{b}} ,  (3.37)

mit positiven Konstanten c_{i}, i = 1,…, 3. Beweis:

∑_{ν=1}^{2q} f_{1}^{2}(s_{ν}) = \frac{b}{π} ∥f_{1}∥_{L_{2}}^{2} + \frac{b}{π} ∑_{m≠0}(−1)^{m} f̂_{1} ∗ f̂_{1}(2bm)
⇒ ∥f_{1}∥_{L_{2}}^{2} ≤ \frac{π}{b} ε^{2} − ∑_{m≠0}(−1)^{m} f̂_{1} ∗ f̂_{1}(2bm)
≤ \frac{π}{b} ε^{2} + ∑_{m≠0} ∫ ; IR |f̂_{1}(x − 2bm) f̂_{1}(x)| dx
≤ \frac{π}{b} ε^{2} + ∑_{m≠0} ( b ; ∫ ; −b |f̂_{1}(x − 2bm) f̂_{1}(x)| dx

+ b ; ∫ ; −b |f̂_{1}(x + 2bm) f̂_{1}(x)| dx + ε_{1}^{2} ∫ ; IR \frac{1}{(1 + (x − 2bm)^{2} ) (1 + x^{2})} dx)

≤ \frac{π}{b} ε^{2} + (∥f̂_{1}∥_{L_{2}([−b,b])} c ε_{1} b^{− \frac{3}{2}} ∑_{m≠0} m^{−2}

+ ε_{1}^{2} 2π ∑_{m≠0}(4b^{2}m^{2} + 4)^{−1}) ,

wobei wir zweimal die Cauchy-Schwarzsche Ungleichung und die Abschätzung (3.36) angewendet, sowie das dritte Integral direkt berechnet haben.

⇒ ∥f_{1}∥_{L_{2}}^{2} − \frac{c_{1}ε_{1}}{b^{3/2}} ∥f_{1}∥_{L_{2}} ≤ \frac{π}{b} ε^{2} + c_{2} ε_{1}^{2} b^{−2}
⇒ ∥f_{1}∥_{L_{2}} ≤ \frac{c_{1}ε_{1}}{2b^{3/2}} + \sqrt{ε_{1}^{2} (\frac{c_{1}^{2}}{4b^{3}} + \frac{c_{2}}{b^{2}} ) + \frac{π}{b} ε^{2}}
□

Bemerkung 3.1.9 Der Satz 3.1.8 besagt also, daß diejenigen Funktionen f_{1} welche der Bedingung (3.36)genügen und deren Abtastwerte sämtlich klein sind, selbst auch nahezu verschwinden müssen. Hierzu war die Koppelung von h an b in der Form h = \frac{π}{b} erforderlich.

Wir müssen nun noch die Abschätzung für f_{1} = h_{l′,p} ∗ d_{h} auf h_{l′,p} übertragen:

∥h_{l′,p}∥_{L_{2}}^{2} = ∥ĥ_{l′,p}∥_{L_{2}([−b,b])}^{2} + ∥ĥ_{l′,p}∥_{L_{2}([−b,b]^{C})}^{2}
≤ sup_{σ∈[−b,b]} |d̂_{h}^{−2}(σ)| · ∥ĥ_{l′,p} · d̂_{h}∥_{L_{2}([−b,b])}^{2} + ∥ĥ_{l′,p}∥_{L_{2}([−b,b]^{C})}^{2}
= (2π)^{−1} sup_{σ∈[−b,b]} |d̂_{h}^{−2}(σ)| · ∥f̂_{1}∥_{L_{2}([−b,b])}^{2} + ∥ĥ_{l′,p}∥_{L_{2}([−b,b]^{C})}^{2}
≤ (2π)^{−1} sup_{σ∈[−b,b]} |d̂_{h}^{−2}(σ)| · ∥f_{1}∥_{L_{2}}^{2} + ∥ĥ_{l′,p}∥_{L_{2}([−b,b]^{C})}^{2} .

Hierbei sind wir davon ausgegangen, daß d̂_{h} im Intervall [−b, b] streng positiv ist. Für das Rechteckprofil ist dies erfüllt, und es gilt

sup_{σ∈[−b,b]} |d̂_{h}^{−2}(σ)| = \frac{π^{3}}{2} ,

so daß jetzt die Kombination von Satz 3.1.7 und Satz 3.1.8 eine Abschätzung für h_{l′,p} liefert.

Als Quintessenz können wir festhalten, daß die Funktion h_{l′,p}, deren Fouriertransformierte fast gänzlich auf [−b, b] konzentriert ist, durch die Abtastungen mit der Detektorfunktion d_{h} stabil bestimmt ist, wenn h = \frac{π}{b} gewählt wird. Dies war im Hinblick auf das Abtasttheorem zu erwarten. Bei der verallgemeinerten Lösung der Gleichungssysteme in (3.33), dürfen wir also die Singulärvektoren zu sehr kleinen Singulärwerten ohne großen Genauigkeitsverlust unberücksichtigt lassen. Unsere Untersuchungen haben gezeigt, daß es möglich ist, alle wesentlich von Null verschiedenen Entwicklungskoeffizienten (f, f_{ml}) hinreichend genau aus den Daten zu berechnen. Dazu war im Prinzip nur die Annahme der annähernden b-Bandbeschränktheit von f erforderlich. Funktionen dieses Typs können also bei richtiger Wahl der Parameter von dem Verfahren ”aufgelöst“werden.

Wir werden die gewonnenen Erkenntnisse im 4. Kapitel auf den Direkten Algebraischen Algorithmus anwenden. Dieser unterscheidet sich von dem in Kapitel 3 untersuchten Vorgehen durch eine weitere Diskretisierungsstufe. Wir betrachten dort die Eigenvektoren von RR^{∗} anstelle der Eigenfunktionen von RR^{∗}.

3.2 Die Interlaced-Geometrie

Wir wollen nun eine analoge Untersuchung für die Interlaced-Geometrie durchführen. Dabei zeigt sich zunächst eine unüberwindlich scheinende Schwierigkeit. Der zu (3.12) analoge Fehler läßt sich nicht mehr ∀l′ ∈ {−p,…, p} abschätzen; er nimmt typischerweise mit wachsendem Betrag von l′ zu, so daß die (3.14) entsprechenden Gleichungssysteme z.T. extrem gestört werden. Es läßt sich jedoch zeigen, daß der Fehler eine bestimmte Struktur hat, die es erlaubt, die Gleichungssysteme trotzdem zu ”lösen“, vorausgesetzt man hält gewisse Bedingungen an p und q genau ein.

Wir starten wieder mit

Rf(φ_{j}, s) = ∑_{l∈ZZ} e^{ilj \frac{π}{p}} h_{l}(s), mit  h_{l}(s) = \frac{1}{π} ∑_{m≥|l|}′σ_{m}(f, f_{ml})u_{m}(s) (3.38)

und behandeln nun aber die geraden und ungeraden Winkelindizes getrennt:

g_{l′}^{I} := \frac{1}{p} ∑_{j=0}^{p−1} e^{−il′2j \frac{π}{p}} Rf(φ_{2j}, s) = ∑_{l∈ZZ} \underbrace{\frac{1}{p} ∑_{j=0}^{p−1} e^{−i(l−l′)2j \frac{π}{p}}}_{=\{ \begin{matrix}1, & l − l′ = 0 mod p \\ 0, & sonst \end{matrix}} h_{l}(s)

= h_{l′}(s) + ∑_{k≠0} h_{l′+pk} .  (3.39)

Ebenso:

g_{l′}^{II} := \frac{1}{p} ∑_{j=0}^{p−1} e^{−il′(2j+1) \frac{π}{p}} Rf(φ_{2j+1}, s)
= ∑_{l∈ZZ} e^{i(l−l′) \frac{π}{p}} \frac{1}{p} ∑_{j=0}^{p−1} e^{−i(l−l′)2j \frac{π}{p}}
= h_{l′}(s) + ∑_{k≠0} (−1)^{k} h_{l′+pk} .  (3.40)

Wir numerieren (entsprechend (2.10) bzw. (2.11) ) die Abtastpunkte s_{ν} von 1 bis 2q(+1) durch und können die Interlaced-Geometrie dann auch wie folgt beschreiben. Für gerade Winkelindizes j werden nur die Werte

Rf(φ_{j}, ·) ∗ d_{h}(s_{2ν}), ν = 1,…, q  (3.41)

und für ungerade j nur die Werte

Rf(φ_{j}, ·) ∗ d_{h}(s_{2ν−1}), ν = 1,…, q(+1)  (3.42)

gemessen, wobei insgesamt die Winkel φ_{0},…, φ_{p−1} durchlaufen werden. Bei der formalen Ergänzung dieser Daten für die Winkel φ_{p} bis φ_{2p−1}, wie sie für (3.39) und (3.40) benötigt wird, haben wir darauf zu achten, daß der Datensatz für φ_{j} zu dem für φ_{j+p} äquivalent ist. Um dies zu erreichen, muß im Fall der Abtastung mit Zentralstrahl p gerade, im anderen Fall p ungerade sein. Im ersten Fall gehen nämlich die Abtastschemata (3.41) und (3.42) durch Drehung um 180^{∘} jeweils in sich, im zweiten Fall in das jeweils komplementäre über. Würden wir die o.a. Bedingung an p verletzen, so würden sich die Datensätze zu φ_{j} und φ_{j+p} zum vollen Datensatz ergänzen, und es läge Standard-Geometrie vor.

Um nicht durch viele unübersichtliche Fallunterscheidungen aufgehalten zu werden, führen wir die weitere Rechnung nur noch für die Geometrie ohne Zentralstrahl durch. Die Behandlung des anderen Falls ergibt sich durch einfache und unwesentliche Modifikationen.

Wir fassen (3.39) und (3.40) zusammen und erhalten mit der Setzung

g_{l′} ∗ d_{h}(s_{ν}) = \{ \begin{matrix}g_{l′}^{I} ∗ d_{h}(s_{ν}), & ν gerade \\ g_{l′}^{II} ∗ d_{h}(s_{ν}), & ν ungerade \end{matrix}  (3.43)

für l′ = −p,…, p die zu (3.10) analoge Gleichung

g_{l′} ∗ d_{h}(s_{ν}) = h_{l′} ∗ d_{h}(s_{ν}) + \{ \begin{matrix}∑_{k≠0} h_{l′+pk} ∗ d_{h}(s_{ν}), & ν gerade \\ ∑_{k≠0} (−1)^{k} h_{l′+pk} ∗ d_{h}(s_{ν}), & ν ungerade \end{matrix} . (3.44)

Man erkennt nun sofort den deutlichen Unterschied zur Standard-Geometrie. Der erste Störterm in der Summe auf der rechten Seite von (3.44) trägt nun den Index l′−sign(l′)p, im Gegensatz zu l′ − sign(l′)2p in (3.10). D.h. der Betrag des Index ist in jedem Fall gleich p−|l′|. Da wir h_{l} aber nur für |l| ≥ p als sehr klein ansehen können, ist der Fehler in (3.44), mit Ausnahme von l′ = 0, keinesfalls vernachlässigbar klein. Wir müssen vielmehr damit rechnen, daß die Norm der Fehlerkomponenten

r_{1,l′}(s_{ν}) := \{ \begin{matrix}∑_{k≠0} h_{l′+pk} ∗ d_{h}(s_{ν}), & ν gerade \\ ∑_{k≠0} (−1)^{k} h_{l′+pk} ∗ d_{h}(s_{ν}), & ν ungerade \end{matrix} , ν = 1,…, 2q, (3.45)

für l′ > p/2 wesentlich größer ist als die Norm von (h_{l′} ∗ d_{h}(s_{ν}))_{ν=1,...,2q}. Wir werden im folgenden sehen, daß die Abhängigkeit von r_{1,l′} vom Index l′ nun wesentlich wird, weshalb wir dies in der Notation berücksichtigen. Desweiteren verifiziert man leicht, daß der Störterm r_{1,l′} die Parität von l′ besitzt, d.h. es gilt (beachte s_{2q+1−ν} = −s_{ν}, p ungerade):

r_{1,l′}(s_{2q+1−ν}) = (−1)^{l′} r_{1,l′}(s_{ν}), ν = 1,…, q .  (3.46)

Da das Gleiche für h_{l′} ∗ d_{h} gilt, können wir r_{1,l′} nicht (wie man vielleicht zunächst vermuten würde) durch (anti-)symmetrische Mittelbildung eliminieren. Der Grund für die erfolgreiche Anwendbarkeit der Interlaced-Geometrie liegt hier viel versteckter. Wir werden zeigen können, daß der Vektor (r_{1,l′}(s_{ν}))_{ν=1,...,2q} nahezu senkrecht auf (h_{l′} ∗ d_{h}(s_{ν}))_{ν=1,...,2q} steht.

Als Vorbereitung dazu dient folgendes

Lemma 3.2.1 Es seien f_{i}, i = 1, 2 stetige Funktionen mit

f_{i}(x) = 0 ∀ x ∈ ]−1 − \frac{h}{2}, 1 + \frac{h}{2}[^{C}

und

f̂_{1}f_{2} = O (∥ξ∥^{−1−δ}) , mit δ > 0 .

Dann gilt für s_{ν} = −1 − \frac{h}{2} + νh, ν = 0, 1,…, 2q, q = 1/h:

∑_{ν=1}^{q} f_{1} · f_{2}(s_{2ν−1}) − f_{1} · f_{2}(s_{2ν}) = h^{−1} (−1)^{q} ∑_{l ungerade} i^{l} f̂_{1} ∗ f̂_{2} (\frac{π}{h}l) .  (3.47)

Bemerkung 3.2.2 Insbesondere fehlt auf der rechten Seite von (3.47) der Term für l = 0.

Beweis:

∑_{ν=1}^{q} f_{1}f_{2}(s_{2ν−1}) − f_{1}f_{2}(s_{2ν})

= ∑_{ν=1}^{q}(f_{1}f_{2})_{−1− \frac{3h}{2}} (2νh) − ∑_{ν=1}^{q}(f_{1}f_{2})_{−1− \frac{h}{2}} (2νh) ,

mit (f_{1}f_{2})_{c}(x) = f_{1}f_{2}(x + c) .

Die Poissonsche Summenformel (1.17) liefert nun

∑_{ν=1}^{q} f_{1}f_{2}(s_{2ν−1}) − f_{1}f_{2}(s_{2ν})

= (2π)^{1/2} (2h)^{−1} ∑_{l∈ZZ}((f_{1}f_{2})_{−1− \frac{3h}{2}} )^{∧}(\frac{π}{h}l) − ((f_{1}f_{2})_{−1− \frac{h}{2}} )^{∧}(\frac{π}{h}l)
= (2h)^{−1} ∑_{l∈ZZ} (e^{−i (1+ \frac{3h}{2} ) \frac{π}{h} l} − e^{−i (1+ \frac{h}{2} ) \frac{π}{h} l}) f̂_{1} ∗ f̂_{2}(\frac{π}{h}l)
= (2h)^{−1} ∑_{l∈ZZ}(−1)^{ql} i^{l} (1 − (−1)^{l}) f̂_{1} ∗ f̂_{2}(\frac{π}{h}l)
= h^{−1} (−1)^{q} ∑_{l ungerade} i^{l} f̂_{1} ∗ f̂_{2}(\frac{π}{h}l) .
□

Wir wollen Lemma 3.2.1 auf die durch r_{1,l′} gestörten Systeme (3.44) anwenden. Dazu ist zunächst, analog zu (3.11), h_{l′} durch h_{l′,K} mit K = p zu ersetzen, was wiederum zu einem (natürlich von l′ und p abhängigen) Fehler r_{2} führt. Somit erhalten wir, vollkommen analog zu dem Vorgehen im Fall der Standard-Geometrie, die Gleichungssysteme (ν = 1,…, 2q):

\frac{1}{π} ∑_{m=|l′|}^{p}′ (f, f_{ml′}) σ_{m}u_{m} ∗ d_{h}(s_{ν}) = g_{l′} ∗ d_{h}(s_{ν}) + r_{1,l′} ∗ d_{h}(s_{ν}) + r_{2} ∗ d_{h}(s_{ν}). (3.48)

Der Fehler r_{2} ist identisch mit dem aus dem Standardfall (3.13), so daß er durch Lemma 3.1.2 abgeschätzt ist.

Desweiteren ist r_{1,l′} noch einmal in zwei Komponenten zu zerlegen:

r_{1,l′} ∗ d_{h}(s_{ν}) = \underbrace{∑_{k≠0} h_{l′+2kp} ∗ d_{h}(s_{ν})}_{= r_{1,l′}^{s} ∗ d_{h}(s_{ν})} + (−1)^{ν} \underbrace{∑_{k≠0} h_{l′+(2k+1)p} ∗ d_{h}(s_{ν})}_{= r_{1,l′}^{i} ∗ d_{h}(s_{ν})} .  (3.49)

Der Anteil r_{1,l′}^{s} ∗ d_{h}(s_{ν}) ist identisch mit dem Fehler r_{1}^{D} im Fall der Standard-Geometrie und muß nach Lemma 3.1.1 bzw. (3.31) nicht weiter verfolgt werden. Für den diskreten Interlaced-Fehler führen wir die Abkürzung

r_{ν}^{iD} = (−1)^{ν} r_{1,l′}^{i} ∗ d_{h}(s_{ν}) , ν = 1,…, 2q  (3.50)

ein.

Wir wollen nun Lemma 3.2.1 anwenden, um zu zeigen, daß die inneren Produkte der Vektoren (r_{1,l′} ∗ d_{h})(s_{ν}))_{ν=1,...,2q} mit den Vektoren der Form (^{∑′}_{m=|l′|}^{p}α_{m}u_{m} ∗ d_{h}(s_{ν}))_{ν=1,...,2q} für solche Werte von α_{m}, m = |l|′, |l′| + 2,…, p klein sind, für die ^{∑′}_{m=|l′|}^{p}α_{m}u_{m} eine näherungsweise b-bandbeschränkte Funktion ist. Satz 3.1.7 zeigt, daß es sinnvoll ist, sich bei einer verallgemeinerten Lösung von (3.48) auf genau solche (α_{m})_{m} zu beschränken.

Der Beweis dazu wird leider sehr technisch werden, eine vollkommene Einsicht in die wirkenden Prinzipien bekommen wir jedoch schon, wenn wir grob argumentieren und einige Größen in der Rechnung gänzlich vernachlässigen, obgleich sie eigentlich nur numerisch ”klein“ sind. In diesem Sinn wollen wir zunächst f_{1} := h_{l′,p} ∗ d_{h} und f_{2} := ∑_{k≠0} h_{l′+(2k+1)p} ∗ d_{h} als strikt b-bandbeschränkt annehmen. Nun beginnt die Entwicklung nach den u_{m} für h_{l′,p} mit dem Index |l′| und für ∑_{k≠0} h_{l′+(2k+1)p} mit p − |l′|; die Rechnung in (3.20) im Verein mit der Debyeschen Beziehung zeigt nun, daß ĥ_{l′p}(σ) für |σ| < ϑl′ und ∑_{k≠0} ĥ_{l′+(2k+1)p}(σ) für |σ| < ϑ (p − |l′|) sehr klein ist (wobei ϑ ”etwas kleiner als 1“ angenommen sei).

Wir vergröbern stark und nehmen an

supp \binom{}{f_{1}} ⊂ [−b, −|l′|] ∪ [|l′|, b] ,

supp \binom{}{f_{2}} ⊂ [−b, −(b − |l′|)] ∪ [b − |l′|, b] .

Wir wenden Lemma 3.2.1 an und erhalten

|∑_{ν=1}^{2q} h_{l′,p} ∗ d_{h}(s_{ν}) · r_{ν}^{iD}| = |∑_{ν=1}^{q} f_{1}f_{2}(s_{2ν−1}) − f_{1}f_{2}(s_{2ν})|

≤ h^{−1} ∑_{l ungerade} |f̂_{1} ∗ f̂_{2} (\frac{π}{h}l)| = 0, falls h = \frac{π}{b} und p =^{(>)} b . (3.51)

Denn bei dieser Setzung der Parameter p, h überschneiden sich die bei den Faltungen gemäß (3.51) zu verschiebenden Träger von f̂_{1} und f̂_{2} höchstens in Randpunkten.

Um dies schnell einzusehen, mag ein Blick auf die Abbildung 3.1 hilfreich sein.

Für die Berechnung von (3.51) sind die Träger von f̂_{1} und f̂_{2} um \frac{π}{h}(2k+1) = b(2k+1) gegeneinander zu verschieben. Wir sehen, daß für die Fälle k = 0, bzw. k = −1, supp \binom{}{f_{2}} genau in die ”Lücke“ von supp \binom{}{f_{1}} fällt. Für k ≠ 0, −1 ist f̂_{1} ∗ f̂_{2} (b(2k + 1)) trivialerweise gleich Null. Der Beweis von Lemma 3.2.1 zeigt, daß das alternierende Vorzeichen in r_{ν}^{iD} dafür verantwortlich ist, daß in der Summe in (3.51) der Term für l = 0 fehlt. Dies stellt sich jetzt als wesentlich heraus. Unverzichtbar für die obige Argumentation ist weiterhin die Abhängigkeit des Fehlers

Abbildung 3.1:

r_{1,l′}vom Index l′, genauer: die Berücksichtigung der Tatsache, daß die Reihenentwicklung von r_{1,l′} bzw. r_{1,l′}^{i} erst mit dem l′-ten Term beginnt und die daraus resultierende Einschränkung des Trägers von r̂_{1,l′}^{i}. Zum dritten ist das Einhalten der Relationen \frac{π}{h} = b, p =^{(>)} b zwingend notwendig. Dies ist das genaue Analogon zum Vorgehen im Beweis zu Theorem 3.1 in [20].

Eine genaue Abschätzung von ⟨(r_{1,l′} ∗ d_{h}(s_{ν}))_{ν} , ( ∑_{m=|l′|}^{p}′α_{m}u_{m} ∗ d_{h}(s_{ν}))_{ν}⟩_{IR_{2q}} ist nun

reine ”Epsilontik“, für die wir zunächst Vorbereitungen schaffen wollen.

Im folgenden seien f̂_{1}, f̂_{2} zunächst beliebige stetige Funktionen, die (für große |ξ| ) der Bedingung |f̂_{i}(ξ)| = O (|ξ|^{−1−δ}), δ > 0, i = 1, 2 genügen. Wir setzen für ϑ ∈]0, 1[

I_{1} := ]−b, −ϑ|l′|[ ∪ ]ϑ|l′|, b[ , I_{2} := ]−b, − (b − ϑ|l′|)[ ∪ ]b − ϑ|l′|, b[ . (3.52)

Dann zerlegen wir f̂_{1} und f̂_{2} wie folgt:

f̂_{1} = χ_{I_{1}} · f̂_{1} + (1 − χ_{I_{1}}) f̂_{1} ,  (3.53)
f̂_{2} = χ_{I_{2}} · f̂_{2} + (1 − χ_{I_{2}}) f̂_{2} ,  (3.54)

wobei χ_{I_{i}}(x) = \{ \begin{matrix}1, & x ∈ I_{i} \\ 0, & sonst \end{matrix} , (i = 1, 2) ist und berechnen damit

∑_{l ungerade} f̂_{1} ∗ f̂_{2}(l b) = ∑_{l ungerade} (f̂_{1}χ_{I_{1}}) ∗ (f̂_{2}χ_{I_{2}}) (l b)

+ ∑_{l ungerade} [((1 − χ_{I_{1}}) f̂_{1}) ∗ (f̂_{2}χ_{I_{2}}) + (χ_{I_{1}} f̂_{1}) ∗ ((1 − χ_{I_{2}}) f̂_{2})] (l b)

+ ∑_{l ungerade} ((1 − χ_{I_{1}}) f̂_{1}) ∗ ((1 − χ_{I_{2}}) f̂_{2}) (l b)

= ∑_{1} + ∑_{2} + ∑_{3} .  (3.55)

Wir erhalten damit

∑_{1} = ∑_{l ungerade} ∫_{I_{2}} f̂_{1}(l b − y)f̂_{2}(y) χ_{I_{1}}(l b − y)dy ,  (3.56)

wegen

χ_{I_{1}}(l b − ·)|_{I_{2}} ≡ 0 ∀ l ≠ 0 ,  (3.57)

gilt also:

∑_{1} = 0 .  (3.58)

Dabei macht man sich die Identität (3.57) an einer zu Abbildung 3.1 analogen Skizze klar (|l′| ist dabei durch ϑ|l′| zu ersetzen), oder rechnet sie mithilfe von Fallunterscheidungen direkt nach.

Um ∑_{2} abzuschätzen, behandeln wir zunächst

∑_{l ungerade} |∫_{I_{2}} (1 − χ_{I_{1}}) f̂_{1}(l b − y)f̂_{2}(y)dy|
≤ ∑_{l ungerade} ∥(1 − χ_{I_{1}}) f̂_{1}∥_{L_{2}(I_{2}+l b)} ∥f̂_{2}∥_{L_{2}(I_{2})}

≤ ∥f̂_{2}∥_{L_{2}(I_{2})} (∥(1 − χ_{I_{1}}) f̂_{1}∥_{L_{2}(I_{2}±b)} + ∑_{|l|≥3 l ungerade} ∥f̂_{1}∥_{L_{2}([−b,b]+l b)})

= ∥f̂_{2}∥_{L_{2}(I_{2})} (∥f̂_{1}∥_{L_{2}([0,ϑ|l′|])} + ∥f̂_{1}∥_{L_{2}([2b−ϑ|l′|,2b])} + ∥f̂_{1}∥_{L_{2}([−ϑ|l′|,0])}
+ ∥f̂_{1}∥_{L_{2}([−2b,−2b+ϑ|l′|])} + ∑_{|l|≥3 l ungerade} ∥f̂_{1}∥_{L_{2}([−b,b]+lb)}) .  (3.59)

Aus ”Symmetriegründen“ gilt dann

∑_{l ungerade} |∫_{I_{1}} f̂_{1}(y) (1 − χ_{I_{2}}) f̂_{2}(l b − y)dy|

≤ ∥f̂_{1}∥_{L_{2}(I_{1})} (∥f̂_{2}∥_{L_{2}([b+ϑ|l′|,2b])} + ∥f̂_{2}∥_{L_{2}([−2b,−b−ϑ|l′|])} + ∥f̂_{2}∥_{L_{2}([ϑ|l′|−b,b−ϑ|l‘|])}
+ ∑_{|l|≥3} ∥f̂_{2}∥_{L_{2}([−b,b]+lb)}) .  (3.60)

Eine Abschätzung für ∑_{2} ergibt sich also durch Kombination von (3.59) und (3.60). Wir behandeln nun zunächst noch kurz ∑_{3}.

∑_{3} = ∑_{l ungerade} ∫_{I_{2}^{C}} (1 − χ_{I_{1}}) f̂_{1}(lb − y)f̂_{2}(y)dy

≤ ∑_{l ungerade} ∫_{I_{2}^{C}∩(I_{1}+lb)^{C}} |f̂_{1}(lb − y)f̂_{2}(y)| dy ,  (3.61)

wobei wir später I_{2}^{C} ∩ (I_{1} + lb)^{C} in

([−b, b]^{C} ∩ [(l − 1)b, (l + 1)b]^{C})

∪ ([−b + ϑ|l′|, b − ϑ|l′|] ∩ [(l − 1)b, (l + 1)b]^{C})
∪ ([−b, b]^{C} ∩ [−ϑ|l′| + bl, ϑ|l′| + bl])
∪ ([−b + ϑ|l′|, b − ϑ|l′|] ∩ [−ϑ|l′| + bl, ϑ|l′| + bl])

zerlegen werden. Die Herleitung von (3.53)–(3.61) geschah zunächst für weitgehend beliebige Funktionen; wir wollen die gewonnenen Formeln jetzt konkret auf die Fouriertransformierten von

f_{1} = h_{l′,p} ∗ d_{h}  (3.62)
f_{2} = ∑_{k≠0} h_{l′+(2k+1)p} ∗ d_{h} = r_{1,l′}^{i} ∗ d_{h}  (3.63)

anwenden. Dazu sind einige Zusatzannahmen erforderlich. Aus unseren bisherigen Bedingungen an f ließ sich nur R_{φ}f ∈ L_{1} und damit R̂_{φ}f ∈ C([−1, 1]) mit lim_{|σ|→∞} R̂_{φ}f(σ) = 0 schließen. Wir wollen nun eine weitere Wachstumsbedingung an R_{φ}f stellen, die dann ebensolche für f_{1} und f_{2} impliziert. Wir wollen verlangen, daß gilt

R̂_{φ}f(σ) = O (|σ|^{−1}) .  (3.64)

Dies dürfte für praktische Belange nicht die geringste Restriktion bedeuten, denn falls f die charakteristische Funktion eines n-Ecks ist, so gilt bekanntermaßen

R̂_{φ}f(σ) = O (|σ|^{−2})

für alle φ bis auf die Kantenrichtungen φ_{1},…, φ_{n}, für die

R̂_{φ_{j}}f(σ) = O (|σ|^{−1}) (j = 1,…, n)  (3.65)

gilt. Ist f die charakteristische Funktion einer Ellipse, so rechnet man leicht

R̂_{φ}f(σ) = O (|σ|^{− \frac{3}{2}} )  (3.66)

nach. Diese beiden Fälle (bzw. die Linearkombinationen davon) dürften aber in der Praxis bereits zu den am schwierigsten zu rekonstruierenden Bilddichten gehören.

Zusätzlich soll für die Detektorfunktion

d̂_{h}(σ) = O (|σ|^{−1})  (3.67)

gelten. Auch dies dürfte in der Praxis keine Einschränkung bedeuten, da die Bedingung vom Rechteckprofil erfüllt wird, und i.a. sogar etwas ”glattere“ Detektorfunktionen Verwendung finden dürften.

Die Kombination von (3.64) und (3.67) führt zur Annahme

(R_{φ}f ∗ d_{h})^{∧} (σ) = O (|σ|^{−2})  (3.68)

für alle φ.

Es ist daher sinnvoll, für h_{l′,p} ∗ d_{h} und r_{1,l′}^{i} ∗ d_{h} folgende Zusatzannahmen zu treffen

|ĥ_{l′,p} · d̂_{h}(σ)| ≤ c\frac{ϵ}{1 + σ^{2}} ∀ |σ| ≥ b ,  (3.69)
|r̂_{1,l′}^{i} · d̂_{h}(σ)| ≤ c\frac{ϵ}{1 + σ^{2}} ∀ |σ| ≥ b .  (3.70)

Aus der Debyeschen Beziehung (1.24) folgt nach einer kleinen Rechnung

∥ĥ_{l′,p}∥_{L_{2}([−ϑ|l′|,ϑ|l′|])} ≤ c e^{− \frac{|l′|}{3} (1−ϑ^{2})^{3/2}}, ∀ |l′| ≥ 1 ,  (3.71)

bzw.

∥r̂_{1,l′}^{i}∥_{L_{2}([−b+ϑ|l′|,b−ϑ|l′|])} ≤ c e^{− \frac{(b−ϑ|l′|)}{3}(1−ϑ^{2})^{3/2}}, ∀ |l′| < p ,  (3.72)

mit von ϑ abhängigen Konstanten c.

Wegen

∥f̂· d̂_{h}∥_{L_{2}} ≤ ∥d̂_{h}∥_{L_{∞}} ∥f̂∥_{L_{2}} ≤ ∥d_{h}∥_{L_{1}} ∥f̂∥_{L_{2}} = ∥d ∥_{L_{1}} ∥f̂∥_{L_{2}}

für beliebige f ∈ L_{2}, gelten die Abschätzungen (3.71), (3.72) (natürlich mit anderen Konstanten) auch für f̂_{1} = (h_{l′,p} ∗ d_{h})^{∧}, bzw. f̂_{2} = (r_{1,l′}^{i} ∗ d_{h})^{∧}.

Die Zusatzannahmen (3.69), (3.70) erlauben die Abschätzungen

∥f̂_{1}∥_{L_{2}([−2b,−2b+ϑ|l′|])} ≤ cϵ\begin{matrix}(ϑ|l′|) & \frac{1}{2} \\ b & 2 \end{matrix} ,  (3.73)
∥f̂_{1}∥_{L_{2}([2b−ϑ|l′|,2b])} ≤ cϵ\begin{matrix}(ϑ|l′|) & \frac{1}{2} \\ b & 2 \end{matrix} ,  (3.74)

bzw.

∥f̂_{2}∥_{L_{2}([b+ϑ|l′|,2b])} ≤ cϵb^{− \frac{3}{2}} (1 − \frac{ϑ|l′|}{b})^{\frac{1}{2}} ,  (3.75)
∥f̂_{2}∥_{L_{2}([−2b,−b+ϑ|l′|])} ≤ cϵb^{− \frac{3}{2}} (1 − \frac{ϑ|l′|}{b})^{\frac{1}{2}} ,  (3.76)

die man einfach durch direktes Nachrechnen erhält.

Wir wollen nun folgendes festhalten:

Lemma 3.2.3 Mit f̂_{1}, f̂_{2} wie in (3.62), (3.63) unter der Zusatzannahme (3.69), (3.70) gilt mit positiven Konstanten c_{i} > 0, (i = 1,…, 5) und ∑_{2} wie in (3.55):

|b∑_{2}| ≤ b (c_{1}e^{− \frac{|l′|}{3} (1−ϑ^{2})^{3/2}} + c_{2}e^{− \frac{(b−ϑ|l′|)}{3} (1−ϑ^{2})^{3/2}})

+ c_{3} \frac{ϵ (ϑ|l′|)^{\frac{1}{2}}}{b} + c_{4}b^{− \frac{1}{2}} ϵ (1 − \frac{ϑ|l′|}{b})^{\frac{1}{2}} + c_{5}ϵb^{− \frac{1}{2}} .

Beweis: Wir setzen in (3.59) und (3.60) die Ungleichungen (3.71)–(3.76) ein und sehen, daß nur noch die Ausdrücke

∑_{|l|≥3 l ungerade} ∥f̂_{i}∥_{L_{2}([−b,b]+lb)} , (i = 1, 2)

abzuschätzen sind. Dies aber ist wiederum eine einfache Anwendung von (3.69), (3.70):

∑_{|l|≥3 l ungerade} ∥f̂_{i}∥_{L_{2}([−b,b]+lb)} = ∑_{|l|≥3 l ungerade} ( ∫_{(l−1)b}^{(l+1)b}|f̂_{i}(x)|^{2} dx)^{\frac{1}{2}}

≤ ϵ ∑_{|l|≥3 l ungerade} ( ∫_{(l−1)b}^{(l+1)b}x^{−4}dx)^{\frac{1}{2}} ≤ cϵ ∑_{|l|≥3 l ungerade} l^{−2}b^{− \frac{3}{2}}
≤ cϵb^{− \frac{3}{2}} .

□

Wir wollen nun nur noch sehen, daß auch ∑_{3} aus (3.55) mit (3.69), (3.70) weiter abgeschätzt werden kann.

Lemma 3.2.4 Mit f̂_{1}, f̂_{2} wie in (3.62), (3.63) unter der Zusatzannahme (3.69), (3.70) gilt mit positiven Konstanten c_{i} > 0, (i = 1,…, 3) und ∑_{3} wie in (3.55):

|b ∑_{3}| ≤ \frac{ϵ}{b} · (c_{1}ϵ + c_{2} \sqrt{b − ϑ|l′|} e^{− \frac{(b−ϑ|l′|)}{3} ·(1−ϑ^{2})^{3/2}} + c_{3} \sqrt{ϑ|l′|} e^{− \frac{|l′|}{3} ·(1−ϑ^{2})^{3/2}}) . (3.77)

Beweis: Wir verwenden (3.61) und die dort angesprochene Zerlegung der Menge I_{2}^{C} ∩ [I_{1} + lb]^{C}. Wir beginnen mit M_{1} := [−b, b]^{C} ∩ [(l − 1)b, (l + 1)b]^{C}:

b ∑_{l ungerade} ∫_{M_{1}} |f̂_{1}(lb − y) f̂_{2}(y)|dy

≤ b ∑_{l ungerade} ∫_{M_{1}} \frac{ϵ^{2}}{(1 + (lb − y)^{2}) (1 + y^{2})} dy
≤ ϵ^{2} b ∑_{l ungerade} ∫_{IR} \frac{1}{(1 + (lb − y)^{2}) (1 + y^{2})} dy
= 2π ϵ^{2} b ∑_{l ungerade} \frac{1}{(lb)^{2} + 4}
≤ 2π ϵ^{2} b^{−1} ∑_{l ungerade} l^{−2} ≤ c ϵ^{2} b^{−1} .

Nun wird M_{2} := [−b + ϑ|l′|, b − ϑ|l′|] ∩ [(l − 1)b, (l + 1)b]^{C} behandelt:

b ∑_{l ungerade} ∫_{M_{2}} |f̂_{1}(lb − y) f̂_{2}(y)|dy

≤ b ∥f̂_{2}∥_{L_{2}([−b+ϑ|l′|,b−ϑ|l′|])} ∑_{l ungerade} ∥f̂_{1}(lb − ·)∥_{L_{2}(M_{2})}

≤ b ∥f̂_{2}∥_{L_{2}([−b+ϑ|l′|,b−ϑ|l′|])} [( ∫_{−b+ϑ|l′|}^{0} \frac{ϵ^{2}}{(1 + (b − x)^{2})^{2}} dx)^{1/2}

+ ( ∫_{0}^{b−ϑ|l′|}\frac{ϵ^{2}}{(1 + (−b − x)^{2})^{2}} dx)^{1/2}

+ ∑_{|l|≥3 l ungerade} ( ∫_{−b+ϑ|l′|}^{b−ϑ|l′|} \frac{ϵ^{2}}{(1 + (lb − x)^{2})^{2}} dx)^{1/2}]

≤ b ∥f̂_{2}∥_{L_{2}([−b+ϑ|l′|,b−ϑ|l′|])} ·

[2ϵ ( ∫_{b}^{2b−ϑ|l′|}x^{−4})^{−1/2} + ∑_{|l|≥3 l ungerade} ϵ ( ∫_{(l−1) b+ϑ|l′|}^{(l+1) b−ϑ|l′|} x^{−4} dx)^{1/2}]

≤ c ∥f̂_{2}∥_{L_{2}([−b+ϑ|l′|,b−ϑ|l′|])} (\frac{ϵ}{b}\sqrt{b − ϑ|l′|} + \frac{ϵ}{b} ∑_{|l|≥3 l ungerade} l^{−2} \sqrt{b − ϑ|l′|})
≤ c \frac{ϵ}{b} \sqrt{b − ϑ|l′|} · e^{− \frac{(b−ϑ|l′|)}{3} (1−ϑ^{2})^{3/2}} ,

wobei im letzten Schritt (3.72) benutzt wurde. Nun zu M_{3} := [−b, b]^{C} ∩ [−ϑ|l′| + bl, ϑ|l′| + bl]:

b ∑_{l ungerade} ∫_{M_{3}} |f̂_{1}(lb − y) f̂_{2}(y)|dy

≤ b ∥f̂_{1}∥_{L_{2}([−ϑ|l′|,ϑ|l′|])} ∑_{l ungerade} ∥f̂_{2}∥_{L_{2}(M_{3})}
= b ∥f̂_{1}∥_{L_{2}([−ϑ|l′|,ϑ|l′|])} (( ∫_{b}^{b+ϑ|l′|} |f̂_{2}(x)|^{2} dx)^{1/2}

+ ( ∫_{−ϑ|l′|−b}^{−b} |f̂_{2}(x)|^{2} dx)^{1/2} + ∑_{|l|≥3 l ungerade} ∥f̂_{2}∥_{L_{2}([−ϑ|l′|+bl,ϑ|l′|+bl])}) .

≤ c \frac{ϵ}{b} \sqrt{ϑ|l′|} e^{− \frac{|l′|}{3} (1−ϑ^{2})^{3/2}} .

Zum Schluß erfreulich einfach ist der Fall der Mengen M_{4} := [−ϑ|l′| + bl, ϑ|l′| + bl] ∩ [−b + ϑ|l′|, b − ϑ|l′|] .

Es gilt:

b ∑_{l ungerade} ∫_{M_{4}} |f̂_{1}(lb − y) f̂_{2}(y)|dy = 0 ,

denn M_{4} = ∅ für |l| > 1 und M_{4} = {±(b − ϑ|l′|} für l = ±1. □ Wir haben nun das angestrebte Ziel, für < r_{ν}^{iD}, h_{l′,p} ∗ d_{h}(s_{ν}) >_{IR}_{2q} eine Abschätzung zu finden, erreicht.

Satz 3.2.5 Mit f̂_{1}, f̂_{2} wie in (3.62), (3.63) p = \frac{1}{ϑ} b, \frac{1}{q} = h = \frac{π}{b}, und den Bedingungen (3.69), (3.70) gilt für |l′| ≥ 1 und ϑ |l′| < b mit positiven Konstanten c_{i}, i = 1,…, 8:

|∑_{ν=1}^{2q} h_{l′,p} ∗ d_{h}(s_{ν}) r_{1,l′}^{i}(s_{ν})|

≤ b (c_{1} e^{− \frac{|l′|}{3} (1−ϑ^{2})^{3/2}} + c_{2} e^{− \frac{(b−ϑ|l′|)}{3} (1−ϑ^{2})^{3/2}})

+ c_{3} (ϵ |l′|)^{1/2} + c_{4} ϵ b^{−1/2} (1 − \frac{ϑ|l′|}{b})^{1/2} + c_{5} ϵ b^{−1/2}

+ \frac{ϵ}{b} (c_{6} ϵ + c_{7} \sqrt{b − ϑ|l′|} e^{− \frac{(b−ϑ|l′|)}{3} (1−ϑ^{2})^{3/2}} + c_{8} \sqrt{ϑ|l′|} e^{− \frac{|l′|}{3} (1−ϑ^{2})^{3/2}}) .

Beweis: Die Lemmata 3.2.3, 3.2.4 im Verein mit der Zerlegung (3.55) sowie (3.58). □

Bemerkung 3.2.6 i) Die unerfreuliche Abschätzungsarbeit in den obigen Lemmata sollte nicht den Blick darauf versperren, daß (3.58) das eigentliche Kernstück des Beweises ist. ii) Die Terme vor c_{1}, c_{2} sind für kleine Werte von l′ bzw. b − ϑ|l′| nicht notwendig klein. Das Problem läßt sich durch geringfügige Vergrößerung von p (über \frac{b}{ϑ} hinaus) leicht umgehen; für kleine Werte von |l′| ist dann der Index des ersten Störterms von r_{1,l′}^{i} (vgl. (3.49)) betraglich gleich p − |l′| > \frac{b}{ϑ}, so daß ∥r_{1,l′}^{i}∥_{L_{2}(w^{−1})} klein bleibt. Ähnliches gilt für die Indizes mit |l′| ≈ \frac{b}{ϑ}. Wir können uns auf die Lösung von (3.44) für Indizes |l′| ≤ \frac{b}{ϑ} beschränken, da h_{l′} für |l′| > \frac{b}{ϑ} verschwindend klein ist. Wählt man nun p > \frac{b}{ϑ}, so bekommt man für den zu c_{2} gehörigen Exponenten:

−\frac{1}{3} (pϑ − ϑ|l′|) (1 − ϑ^{2})^{3/2} < −\frac{ϑ}{3} (p − \frac{b}{ϑ}) (1 − ϑ^{2})^{3/2} .

iii) Von dem Gesamtfehler im Fall der Interlaced-Geometrie

r = r_{1,l′}^{i} + r_{1,l′}^{s} + r_{2}

hatten wir nur den Anteil r_{1,l′}^{i} zu berücksichtigen, da die Fehler r_{1,l′}^{s} und r_{2} schon in Kapitel 3.1 abgeschätzt wurden. Der Satz 3.2.5 zusammen mit (3.32) zeigt, daß auch der Gesamtfehler nahezu senkrecht auf (h_{l′,p} ∗ d_{h}(s_{ν}))_{ν=1,...,2q} steht.

4 Der Direkte Algebraische Algorithmus (DAA)

Wir wollen an dieser Stelle einen Rekonstruktionsalgorithmus vorstellen, welcher unmittelbaren Bezug auf die in Kapitel 2.1 beschriebene Modellierung des Meßvorganges nimmt. Das Rekonstruktionsproblem wird dabei umgesetzt in ein Problem der linearen Algebra – das Lösen eines (allerdings sehr großen) linearen Gleichungssystems. Aufgrund von Symmetrieüberlegungen ist es möglich, dieses System trotz seiner immensen Größe mit direkten Methoden (im Gegensatz zu Iterationsverfahren) zu lösen, womit auch schon der etwas befremdliche Name erklärt ist. Wir werden uns bei der Beschreibung im wesentlichen auf die parallele Geometrie beschränken, gelegentlich jedoch auch den Fan-Beam-Fall beachten.

4.1 Aufstellung des linearen Gleichungssystems

Der einfachste Zugang dürfte der folgende sein: Wir können natürlich nicht erwarten, daß endlich viele Riesz-Funktionale

g_{j} = (f, l_{j})_{L_{2}(D)} j = 1,…, N

eine unbekannte Funktion f aus dem unendlichdimensionalen Vektorraum L_{2}(D) eindeutig bestimmen. Wir entscheiden uns deshalb für eine L_{2}-Bestapproximation an f aus einem endlichdimensionalen Teilraum V, den wir, und das ist das Charakteristische am DAA, von den Riesz-Funktionen selber aufspannen lassen wollen, d.h. wir suchen die (eindeutige) Lösung f_{V} des Problems

(f − f_{V} , l_{j})_{L_{2}(D)} = 0 , j = 1,…, N
f_{V} ∈ span(l_{1},…, l_{N}) =: V .

Setzen wir also f_{V} = ∑_{i=1}^{N} α_{i}l_{i}, so ist das lineare Gleichungssystem

g_{j} = (f, l_{j})_{L_{2}(D)} = ∑_{i=1}^{N} α_{i}(l_{i}, l_{j})_{L_{2}(D)} , j = 1,…, N  (4.1)

zu lösen. Daß im Fall der linearen Abhängigkeit der l_{j} die Gramsche Matrix ( (l_{i}, l_{j})_{L_{2}})_{i,j} nur positiv semidefinit, das lineare Gleichungssystem (4.1) also nicht eindeutig lösbar ist, kann nicht weiter stören: mit zwei Lösungen (α_{i})_{i} und (β_{i})_{i} von (4.1) gilt für f_{V}^{I} = ∑ α_{i}l_{i} , f_{V}^{II} = ∑ β_{i}l_{i} natürlich

f_{V}^{I} − f_{V}^{II} ∈ V ∩ V^{⊥} = {0} .

Wir können also z.B. stets die Minimum-Norm-Lösung von (4.1) nehmen. Schwerer hingegen wiegt die Tatsache, daß es sich bei (4.1) um die Diskretisierung eines schlechtgestellten Problems handelt, so daß bei dem System (4.1) mit einer großen Konditionszahl zu rechnen ist. Wir müssen demzufolge regularisierte, d.h. insbesondere auch nur noch näherungsweise Lösungen von (4.1) ins Auge fassen.

Bemerkung 4.1.1 Formal können wir obige Vorgehensweise in die Klasse der Momentenverfahren einordnen: Es wird < e_{i}, Rf_{V} − g >= 0 , i = 1,…, N gelöst, wobei e_{i} den i-ten kanonischen Einheitsvektor bezeichnet.

Wir paraphrasieren (4.1) noch einmal, wobei wir auf die Notation in Kapitel 2 zurückgreifen. Gesucht ist f_{V} ∈ V = range(R^{∗}) mit R(f_{V} − f) = 0. Dies ist äquivalent zu:

löse  RR^{∗}α = Rf ,
setze  f_{V} = R^{∗}α .

(4.2)

Nach Konstruktion ist f_{V} sowohl Lösung der Minimierungsaufgabe:

min_{α∈IR^{N}} ∥f − R^{∗}α∥_{L_{2}} ,

als auch Minimum-Norm-Lösung von

Rf = g

in L_{2}(D).

Wir wollen uns nun näher mit der Systemmatrix A := RR^{∗} beschäftigen. Es ist zweckmäßig

Rf = ( g^{0} ; ⋮ ; g^{2p−1} ) , mit g^{j} = R_{j}f ,

sowie

R^{∗} ( x^{0} ; ⋮ ; x^{2p−1} ) = ∑_{j=0}^{2p−1} R_{j}^{∗}x^{j}

zu schreiben, um eine natürliche Aufteilung von A in Blockmatrizen R_{j}R_{l}^{∗} zu erhalten:

A = (R_{j}R_{l}^{∗})_{j,l=0,...,2p−1} .

Lösen wir das Gleichungssystem in (4.2) über die Entwicklung von g = Rf nach Eigenvektoren von RR^{∗} zu von Null verschiedenen Eigenwerten, so besteht der DAA in der Berechnung der Minimum-Norm-Lösung der Gleichung

Rf = g

über die Singulärwertzerlegung von R (Satz 1.1.4).

Ein standardmäßiges Vorgehen wäre es nun, auszunutzen, daß für die Blöcke die Beziehung:

R_{j}R_{l}^{∗} = R_{(j−l)mod2p}R_{0}^{∗} ,

besteht, daß es sich bei A also um eine sogenannte blockzyklische Faltung handelt. Dies bedeutet dann, daß wir (4.2) mithilfe von Block-Fouriertransformationen effizient lösen können. In dieser Version ist der DAA erstmals von Natterer [18] vorgestellt worden. Wir wollen das Problem jedoch in einen allgemeineren Rahmen stellen und es als ”System mit Symmetrien“ behandeln. Diese generellere Methode ist in dem einführenden Lehrbuch von Stiefel und Fässler [28] beschrieben und benutzt einige einfache Tatsachen aus der Darstellungstheorie (endlicher) Gruppen.

4.2 Notwendige Grundlagen aus der Darstellungstheorie

Wir geben hier eine kurze Erläuterung der verwendeten Begriffe, wobei wir uns auf das Allereinfachste beschränken können, und geben (ohne Beweis) die für die Theorie notwendigen Aussagen an. Die Beweise können sämtlich in [28] nachgelesen werden.

Definition 4.2.1 Es sei (G, ∘) eine Gruppe und V ein n-dimensionaler Vektorraum über C, mit n ∈ IN. Dann nennen wir einen Homomorphismus ϑ von G in die Gruppe GL(V) eine (lineare) n-dimensionale Darstellung von G; V wird auch Darstellungsraum genannt.

Bemerkung 4.2.1 Eine Darstellung ϑ ordnet also jedem Element s ∈ G eine invertierbare lineare Transformation P_{s} = ϑ(s) auf V zu, und diese Zuordnung respektiert die Gruppenstruktur, d.h. P_{s∘t} = P_{s}P_{t} ∀s, t ∈ G, also insbesondere P_{e} = 1l und P_{s^{−1}} = (P_{s})^{−1} ∀s ∈ G.

Definition 4.2.2 Die Darstellung ϑ heißt treu, falls sie injektiv ist; zwei Darstellungen ϑ_{1}, ϑ_{2} heißen äquivalent, falls ein T ∈ GL(V) existiert, mit

ϑ_{1}(s) = T^{−1}ϑ_{2}(s)T ∀s ∈ G,

d.h. dem Übergang von einer Darstellung ϑ_{1} zu einer äquivalenten ϑ_{2} entspricht eine Koordinatentransformation in V.

Definition 4.2.3 Eine Darstellung heißt unitär bezüglich eines Skalarproduktes [., .], falls gilt

[P_{s}x, P_{s}y] = [x, y] ∀x, y ∈ V , ∀s ∈ G .

Man sagt dann auch, daß das Skalarprodukt invariant gegenüber der Darstellung ist.

Definition 4.2.4 Die Darstellung ϑ heißt reduzibel, falls ein echter Teilraum V_{1} existiert, welcher invariant unter allen P_{s} = ϑ(s), s ∈ G ist. Es gilt dann also

P_{s}(V_{1}) ⊂ V_{1} ∀s ∈ G .

Existiert kein solcher Teilraum, so nennt man ϑ irreduzibel.

Bemerkung 4.2.2 Es mag vielleicht verwundern, warum man in der obigen Definition, scheinbar an der sprachlichen Logik vorbei, von dem Teilraum V_{1} nicht gleich verlangt, reduzierend zu sein. Man zeigt jedoch leicht folgendes:

a) für jede Darstellung einer endlichen Gruppe existiert ein Skalarprodukt bezüglich dessen sie unitär ist,

b) ist V_{1} invariant unter der ( bezüglich irgendeines Skalarproduktes [., .]) unitären Darstellung ϑ, so auch das ( bezüglich [., .] gebildete) orthogonale Komplement V_{1}^{⊥}, d.h. das Paar (V_{1}, V_{1}^{⊥}) reduziert ϑ(s) ∀s ∈ G.

Definition 4.2.5 Ist (V_{1}, V_{2}) ein reduzierendes Paar für die Darstellung ϑ, gilt also V = V_{1} ⊕ V_{2} mit jeweils invarianten Unterräumen, so spricht man von einem Zerfall von ϑ in die Teildarstellungen ϑ_{1}(s) = ϑ(s)_{|_{V_{1}}} , ϑ_{2}(s) = ϑ(s)_{|_{V_{2}}} ∀s ∈ G. Man nennt eine Darstellung vollreduzibel, wenn sie entweder irreduzibel ist, oder in lauter irreduzible Darstellungen zerfällt. Die Zerlegung in irreduzible Darstellungen nennt man Ausreduktion.

Bemerkung 4.2.3 Man erhält nun mühelos, daß jede Darstellung einer endlichen Gruppe vollreduzibel ist, d.h. jede Darstellung ist direkte Summe irreduzibler Darstellungen.

Es ist an dieser Stelle vielleicht angebracht, etwas über den von uns benötigten Rahmen der endlichdimensionalen Darstellungen endlicher Gruppen hinauszugehen. Die oben angesprochene Konstruktion eines invarianten Skalarproduktes läßt sich auf kompakte und sogar auf lokalkompakte Gruppen übertragen. Der wesentliche (von Hurwitz stammende) Trick besteht bei den endlichen Gruppen in einer gewissen Summenbildung über alle Gruppenelemente, die im allgemeinen Fall durch eine Integration bzgl. eines (rechts-) invarianten Maßes zu ersetzen ist. Die Existenz eines solchen Maßes wurde 1933 von A. Haar für (separable) lokalkompakte Gruppen gezeigt. Damit war es möglich, die Theorie der unitären Darstellung auch auf solche Gruppen zu erweitern. Z.B. gilt der Satz über den Zerfall unitärer Darstellungen in endlichdimensionale, irreduzible Teildarstellungen auch im Fall kompakter topologischer Gruppen. Für einen historischen Überblick zur Harmonischen Analyse und ihrer vielfachen Anwendung in Physik und weiten Teilen der Mathematik siehe G.W. Mackey [17].

Definition 4.2.6 Wir sagen eine Matrix A habe die Symmetrien der Darstellungen ϑ, falls gilt

P_{s} A = A P_{s} ∀P_{s} = ϑ(s) , s ∈ G .

Bemerkung 4.2.4 Diese Relation ist natürlich unabhängig von der Wahl einer Basis in V.

Nun können wir das für unsere Belange zentrale Ergebnis der Darstellungstheorie formulieren:

Satz 4.2.5 Es sei ϑ eine vollreduzible Darstellung einer Gruppe G in den n-dimensionalen Vektorraum V mit der ausreduzierten Form

ϑ = c_{1}ϑ_{1} + c_{2}ϑ_{2} + … + c_{l}ϑ_{l} ,

d.h. die ϑ_{j} sind irreduzibel und paarweise inäquivalent, die c_{j} ∈ IN sind die Vielfachheiten der ϑ_{j}, mit n_{j} sei die Dimension der ϑ_{j} bezeichnet. Die Matrix A habe die Symmetrien von ϑ, also

P_{s} A = A P_{s} ∀P_{s} = ϑ(s) , s ∈ G .

Dann existiert eine Basis von V, die sogenannte symmetriegerechte Basis von V, so daß A, bezüglich dieser Basis dargestellt, die folgende Blockdiagonalstruktur hat:

n_{1} (c_{1}, c_{1})-Matrizen D_{1}
n_{2} (c_{2}, c_{2})-Matrizen D_{2}
⋮ ⋮ ⋮
n_{l} (c_{l}, c_{l})-Matrizen D_{l} .

(4.3)

Anders ausgedrückt: Sei A bezüglich irgendeiner Basis dargestellt (in aller Regel bzgl.der kanonischen), dann existiert eine Koordinatentransformation S derart, daß gilt

S^{−1}AS = ( \begin{matrix}D_{1} \\ & D_{1} & & & & & & 0 \\ & & . & . & . \\ & & & & & D_{1} \\ & & & & & & D_{2} \\ . & . & . \\ & 0 & & & & & & & ... \\ & & & & & & & & & D_{l} \\ & & & & & & & & & & . & . & . \\ & & & & & & & & & & & & & D_{l} \end{matrix} ) .

Bemerkung 4.2.6 Einen Vorteil wird man sich von dieser Darstellung nur dann erhoffen dürfen, wenn die Matrizen S, S^{−1} ”einfach“ sind, d.h. sie sollten mit geringem Aufwand abspeicherbar und schnell anzuwenden sein.

Wir werden im nächsten Abschnitt ausnutzen, daß unser Problem die Symmetrien der Diedergruppe aufweist. Dabei tritt als Transformationsmatrix im wesentlichen die unitäre Matrix der Blockfouriertransformation auf. Diese besitzt die oben genannte Eigenschaft der ”Einfachheit“ in hohem Maße.

4.3 Anwendung der Gruppentheoretischen Methode auf die diskrete Radon-Transformation

Wir wollen hier zeigen, daß die Matrix

A = RR^{∗}

die Symmetrien der Diedergruppe besitzt. Voraussetzung dafür ist allerdings gleichmäßiges Abtasten in Winkelrichtung und winkelunabhängiges Abtasten in lateraler Richtung. Um nicht durch Fallunterscheidungen unnötig aufgehalten zu werden, beschränken wir uns auf die Geometrie ohne Zentralstrahl und gehen von 2p gleichverteilten Richtungen (bzw. p Quellenpositionen) und dann jeweils 2q gleichverteilten Messungen pro Richtung aus. Für die parallele Geometrie sind damit die Hälfte der Daten redundant, was wir später berücksichtigen werden. Zunächst können jedoch Parallel- und Fan-Beam-Geometrie simultan behandelt werden.

Wir wollen die in Kapitel 3 von [28] angegebenen Tabelle symmetriegerechter Basisfunktionen verwenden, wozu wir den Datenvektor g als Funktion auf einer (2p) · (2q)- elementigen Punktmenge zu interpretieren haben. Wir identifizieren daher einen Vektor g ∈ IR^{2p·2q} mit der Wertebelegung von jeweils 2q gleichverteilten Punkten auf den Seiten eines regelmäßigen 2p-Ecks. Dies sei am Beispiel p = 2, q = 3 illustriert:

Die Operationen der Diedergruppe D_{2p} setzen sich bekanntlich aus der Drehung d = d_{φ} mit φ = \frac{π}{p} und der Spiegelung s an irgendeiner Symmetrieachse zusammen, genauer läßt sich jedes Element t ∈ D_{2p} schreiben als

t = s^{j} d^{k} k = 0,…, 2p − 1 , j = 0, 1 .

Diese Operationen überführen das reguläre 2p-Eck in sich und permutieren dabei die 2p · 2q Randpunkte, d.h. wir haben so eine Darstellung der Diedergruppe im Raum IR^{2p·2q} durch Permutationsmatrizen gewonnen, die mit ϑ_{per} bezeichnet sei. Die Angabe der Permutationen ist denkbar einfach: (wieder sei ϑ_{per}(t) = P_{t})

P_{d^{k}} ( g^{0} ; g^{1} ; ⋮ ; g^{2p−1} ) = ( g^{k} ; g^{k+1} ; ⋮ ; g^{2p−1} ; g^{0} ; ⋮ ; g^{k−1} ) , P_{s} ( g^{0} ; g^{1} ; ⋮ ; g^{2p−1} ) = ( T g^{0} ; T g^{2p−1} ; T g^{2p−2} ; ⋮ ; T g^{1} )

mit g^{j} ∈ IR^{2q}, T = T_{2q} = \binom{0 1}{1 0} ∈ Mat(2q, 2q).

Bemerkung 4.3.1 Der Index an der Matrix T soll ihre Größe bezeichnen, wir werden ihn im folgenden immer dann angeben, wenn wir es für klärend halten. Ähnlich halten wir es mit der Einheitsmatrix 1l und der expliziten Angabe der Räume bei den inneren Produkten.

Nun sind die Vertauschungsrelationen

P_{s^{j}d^{k}}RR^{∗} = RR^{∗}P_{s^{j}d^{k}} , k = 0,…, 2p − 1 , j = 0, 1
⟺ P_{s^{j}}P_{d^{k}}RR^{∗} = RR^{∗}P_{s^{j}}P_{d^{k}} , k = 0,…, 2p − 1 , j = 0, 1

zu beweisen, und dies ist natürlich äquivalent zu dem

Lemma 4.3.2 Es gilt P_{s}RR^{∗} = RR^{∗}P_{s}, sowie P_{d}RR^{∗} = RR^{∗}P_{d} .

Beweis:

< P_{s}Rf, x >IR_{4pq} = < Rf, P_{s}x >IR_{4pq} = ∑_{j=0}^{2p−1} ∑_{l=1}^{2q}(f, L_{h,j,l})_{L_{2}}(T x^{2p−j})_{l}
= ∑_{j=0}^{2p−1} ∑_{l=1}^{2q}(f, L_{h,2p−j,2q+1−l})_{L_{2}} x_{l}^{j} = ∑_{j,l}(f, L_{h,j,l} ∘ S)_{L_{2}}x_{l}^{j}
= ∑_{j,l}(f ∘ S, L_{h,j,l})_{L_{2}}x_{j}^{l} = < R(f ∘ S), x >IR_{4pq}

mit S : IR^{2} ⟶ IR^{2}, \binom{x}{y} ↦ \binom{−x}{y} (Spiegelung an der y-Achse). Also gilt:

P_{s}R = RQ_{S} ,
mit Q_{S} : L_{2} → L_{2} , f ↦ f ∘ S .
⇒ R^{∗}P_{s} = Q_{S}R^{∗}
⇒ P_{s}RR^{∗}P_{s} = RQ_{S}^{2}R^{∗}
⇒ P_{s}RR^{∗} = RR^{∗}P_{s}

Noch einfacher ist:

< P_{d}Rf, x > = ∑_{j,l}(f, L_{h,j+1,l})x_{l}^{j} = ∑_{j,l}(f, L_{h,j,l} ∘ D^{⊤})x_{l}^{j}
= ∑_{j,l}(f ∘ D, L_{h,j,l})x_{l}^{j} = < R(f ∘ D), x > ,

mit D = \binom{cos φ −sin φ}{sin φ cos φ} , φ = \frac{π}{p} .

Also analog zu obigem Vorgehen:

P_{d}R = RQ_{D} ,
mit Q_{D} : L_{2} → L_{2} , f ↦ f ∘ D .
⇒ R^{∗}P_{d}^{∗} = Q_{D}^{∗}R^{∗}
⇒ P_{d}RR^{∗}P_{d}^{∗} = RQ_{D}Q_{D}^{∗}R^{∗}
⇒ P_{d}RR^{∗} = RR^{∗}P_{d}
□

Bemerkung 4.3.3 i.) Die Symmetrien des Abtastschemas haben wir durch die Gleichungen

L_{h,j,l} ∘ S = L_{h,2p−j,2q+1−l}  (4.4)
L_{h,j,l} ∘ D^{⊤} = L_{h,j+1,l}  (4.5)

ausgenutzt. Sie lassen sich direkt aus der Definition von L_{h,j,l} ablesen. ii.) Der Beweis kann mit ganz geringfügigen Abänderungen auch auf die Fan-Beam-

Geometrie übertragen werden. Dort ist

C_{h,j,l} ∘ S = C_{h,p−j,2q+1−l}  (4.6)
C_{h,j,l} ∘ D^{⊤} = C_{h,j+1,l}  (4.7)

auszunutzen, wobei S nun die Spiegelung an der x-Achse bedeutet, und D für die Drehmatrix zum Winkel \frac{2π}{p} steht.

Nachdem wir die Voraussetzungen geklärt haben, können wir nun direkten Gebrauch von der in [28] bereitgestellten Tabelle symmetriegerechter Basislösungen machen (s. S. 56, bzw. S. 58 in [28] für den reellen Fall.) In der dortigen Bezeichnungsweise haben wir es mit q Koronen vom Typ t_{1} zu tun. Fassen wir nun die Basisfunktionen als Spaltenvektoren der (orthogonalen) Transformationsmatrix S zusammen, so liest man für S aus der Tabelle ab:

S = ( \begin{matrix}S^{I} & +S^{II} & S_{0,1}^{III} & S_{0,1}^{IV} & S_{0,2}^{III} & S_{0,2}^{IV} & ⋯ & S_{0,p−1}^{III} & S_{0,p−1}^{IV} \\ S^{I} & −S^{II} & S_{1,1}^{III} & S_{1,1}^{IV} & S_{1,2}^{III} & S_{1,2}^{IV} & ⋯ & S_{1,p−1}^{III} & S_{1,p−1}^{IV} \\ S^{I} & +S^{II} & S_{2,1}^{III} & S_{2,1}^{IV} & S_{2,2}^{III} & S_{2,2}^{IV} & ⋯ & S_{2,p−1}^{III} & S_{2,p−1}^{IV} \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋯ & ⋮ & ⋮ \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋯ & ⋮ & ⋮ \\ S^{I} & −S^{II} & S_{2p−1,1}^{III} & S_{2p−1,1}^{IV} & S_{2p−1,2}^{III} & S_{2p−1,2}^{IV} & ⋯ & S_{2p−1,p−1}^{III} & S_{2p−1,p−1}^{IV} \end{matrix} ) ,

(4.8)

mit S^{I} = \frac{1}{2 \sqrt{p}} \binom{1l_{q} 1l_{q}}{T_{q} −T_{q}} , S^{II} = \frac{1}{2 \sqrt{p}} \binom{1l_{q} −1l_{q}}{T_{q} T_{q}} ,

sowie S_{i,j}^{III} = \frac{1}{\sqrt{2p}} \binom{c_{ij}1l_{q} −s_{ij}1l_{q}}{c_{ij}T_{q} s_{ij}T_{q}} , S_{i,j}^{IV} = \frac{1}{\sqrt{2p}} \binom{s_{ij}1l_{q} −c_{ij}1l_{q}}{s_{ij}T_{q} −c_{ij}T_{q}} ,

wobei c_{l} = cos \frac{πl}{p} , s_{l} = sin \frac{πl}{p} bedeutet.

Wir lesen weiterhin ab (S. 58 in [28]), daß in ϑ_{per} die irreduziblen Darstellungen der Diedergruppe mit den folgenden Vielfachheiten auftreten:

Darstellung  Dimension Vielfachheit
ϑ_{1}^{1}  1  q
ϑ_{2}^{1}  1  q
ϑ_{3}^{1}  1  q
ϑ_{4}^{1}  1  q  .
ϑ_{l}^{2}, l = 1,…, p − 1  2  2q

Der Satz 4.2.5 besagt nun, daß RR^{∗} = A nach Übergang auf die symmetriegerechte Basis wie folgt zerfällt:

D := S^{⊤}AS = ( \begin{matrix}D_{1}^{I} \\ & DI_{2} & & & & 0 \\ & & D_{3}^{I} \\ & & & D_{4}^{I} \\ & & & & D_{1}^{II} \\ & & & & & D_{1}^{II} \\ & 0 & & & & & D_{2}^{II} \\ & & & & & & & D_{2}^{II} \\ & & & & & & & & . & . & . \\ & & & & & & & & & & & D_{p−1}^{II} \\ & & & & & & & & & & & & D_{p−1}^{II} \end{matrix} )

(4.9)

mit den vier (q, q)-Matrizen D_{l}^{I}, l = 1,…, 4 und den p−1 jeweils doppelt auftretenden (2q, 2q)-Matrizen D_{l}^{II}, l = 1,…, p − 1.

Für die Fan-Beam-Geometrie haben wir damit schon alle Vereinfachungsmöglichkeiten ausgeschöpft, wir wollen jedoch nun zeigen, daß die Redundanzen im Falle der parallelen Geometrie noch zu weiteren, erheblichen Reduktionen in den Gleichungen führen.

Wir berücksichtigen nun die Tatsache, daß es sich bei Rf im Fall paralleler Daten um eine gerade Funktion im folgenden Sinn handelt:

(Rf)_{p+j} = T (Rf)_{j} ∀f ∈ l_{2}(D) .

Dies ist eine triviale Konsequenz aus:

Rf(ω_{p+j}, s) = Rf(−ω_{j}, s) = Rf(ω_{j}, −s) .

Wir schreiben das System (A = RR^{∗}):

Ax = g äquivalent um zu

S^{⊤}ASy = S^{⊤}g , x = Sy .

Um die Wirkung von S^{⊤} auf Vektoren aus range(R) zu berechnen, ist also (vgl. (4.8))

S^{I^{⊤}} + S^{I^{⊤}}T,  S^{II^{⊤}} + (−1)^{p}S^{II^{⊤}}T,
sowie  S_{i,j}^{III^{⊤}} + S_{p+i,j}^{III^{⊤}}T  und S_{i,j}^{IV^{⊤}} + S_{p+i,j}^{IV^{⊤}}T

für i = 0,…, p − 1; j = 0,…, p − 1 auszuwerten. Berücksichtigung von c_{(p+i)j} = (−1)^{j}c_{ij}, s_{(p+i)j} = (−1)^{j}s_{ij} ergibt:

S_{i,j}^{III^{⊤}} + S_{p+i,j}^{III^{⊤}}T = \{ \begin{matrix}\frac{2}{\sqrt{2p}}c_{ij} ( \begin{matrix}1l & T \\ 0 & 0 \end{matrix} ) , & j gerade \\ \frac{2}{\sqrt{2p}}s_{ij} ( \begin{matrix}0 & 0 \\ −1l & T \end{matrix} ) , & j ungerade \end{matrix} ,

sowie

S_{i,j}^{IV^{⊤}} + S_{p+i,j}^{IV^{⊤}}T = \{ \begin{matrix}\frac{2}{\sqrt{2p}}s_{ij} ( \begin{matrix}1l & T \\ 0 & 0 \end{matrix} ) , & j gerade \\ \frac{2}{\sqrt{2p}}c_{ij} ( \begin{matrix}0 & 0 \\ 1l & −T \end{matrix} ) , & j ungerade \end{matrix} .

Weiterhin rechnet man nach

S^{I^{⊤}} + S^{I^{⊤}}T = \frac{1}{\sqrt{p}} ( \begin{matrix}1l & T \\ 0 & 0 \end{matrix} ) ,

S^{II^{⊤}} + (−1)^{p}S^{II^{⊤}}T = \frac{1}{\sqrt{p}} \{ \begin{matrix}( \begin{matrix}1l & T \\ 0 & 0 \end{matrix} ) , & p gerade \\ ( \begin{matrix}0 & 0 \\ −1l & T \end{matrix} ) , & p ungerade \end{matrix} .

Fassen wir nun alle Indizes, auf denen S^{⊤}R nach obiger Rechnung verschwindet, zu der Indexmenge I zusammen, so gilt:

(S^{⊤}Rf)_{i} = 0 ∀i ∈ I , ∀f ∈ L_{2} ,

also auch

< \underbrace{S^{⊤}RR^{∗}S}_{=D} x, e_{j} > = 0  ∀ x ∈ IR^{4pq}, j ∈ I
Sym.^{⟹}  D_{ij} = 0 ,  falls i oder j ∈ I .

Also gilt:

D_{2}^{I} = 0 ,

sowie

D_{4}^{I} = 0 ,  falls p gerade ,
D_{3}^{I} = 0 ,  falls p ungerade ,

und D_{l}^{II} hat die Form : \{ \begin{matrix}( \begin{matrix}Dl^{′II} & 0 \\ 0 & 0 \end{matrix} ) , & für l gerade \\ ( \begin{matrix}0 & 0 \\ 0 & Dl^{′II} \end{matrix} ) , & für l ungerade \end{matrix} .

Nun wird demonstriert, wie die Blockdiagonalmatrizen schnell berechnet werden können; es gilt

DS^{⊤} = S^{⊤}A  , also insbesondere
(DS^{⊤})_{erste Spalte} = (S^{⊤}A)_{erste Spalte} .  (4.10)

Wird die erste (Block-)Spalte von A mit A_{0},…, A_{2p−1} bezeichnet, so folgt aus

(Rf)_{p+j} = T (Rf)_{j}

sofort

A_{k+p} = T A_{k} .

Ferner haben alle A_{j} die Form

A_{j} = \binom{a_{j} b_{j}T_{q}}{T_{q}b_{j} T_{q}a_{j}T_{q}}

mit gewissen symmetrischen (q, q)-Untermatrizen a_{j}, b_{j}, dies ist gleichbedeutend mit

A_{j} = A_{j}^{⊤} ∧ T A_{j} = A_{j}T ,

d.h. die A_{j} sind sowohl symmetrisch als auch persymmetrisch, was unmittelbar anhand der Definitionen von A_{j} nachgewiesen werden kann.

Wir zeigen am Beispiel D_{1}^{I}, wie beide Relationen auszunutzen sind, und geben dann für die anderen Matrizen gleich das Ergebnis an.

Aus (4.10) folgt:

( \begin{matrix}D_{1}^{I} & 0 \\ 0 & 0 \end{matrix} ) · ( \begin{matrix}1l & T_{q} \\ 1l & −T_{q} \end{matrix} ) = 2 \sqrt{p} S^{I^{⊤}} ∑_{j=0}^{2p−1} A_{j}

= 2 \sqrt{p} S^{I^{⊤}}(1l + T_{2q})∑_{j=0}^{p−1} A_{j}

= 2 ( \begin{matrix}1l & T_{q} \\ 0 & 0 \end{matrix} ) ∑_{j=0}^{p−1} A_{j}

= 2 ∑_{j=0}^{p−1} ( \begin{matrix}a_{j} + b_{j} & (a_{j} + b_{j})T_{q} \\ 0 & 0 \end{matrix} ) ,

also

D_{1}^{I} = 2∑_{j=0}^{p−1}(a_{j} + b_{j}) .  (4.11)

Analoge Berechnungen ergeben:

D_{3}^{I} = 2 ∑_{j=0}^{p−1}(−1)^{j}(a_{j} + b_{j})  , falls p gerade ,
bzw. D_{4}^{I} = 2 ∑_{j=0}^{p−1}(−1)^{j}(a_{j} − b_{j})  , falls p ungerade ,

(4.12)

Dl^{′II} = 2 ∑_{j=0}^{p−1} cos(\frac{jlπ}{p})(a_{j} + b_{j}) , falls l gerade ,
Dl^{′II} = 2 ∑_{j=0}^{p−1} cos(\frac{jlπ}{p})(a_{j} − b_{j}) , falls l ungerade ,

(4.13)

wobei l von 1 bis p − 1 läuft.

Wir sehen also, daß wir es mit einer Block-Kosinustransformation der (q, q)-Matrizen (a_{j} ±b_{j}) zu tun haben. Rechnet man standardmäßig, so ist der Aufwand von der Größenordnung q^{2}p^{2}, es gibt aber auch schnelle Versionen speziell für die Kosinustransformation (siehe [7]), bei deren Verwendung sich der Aufwand auf O(q^{2}p log(p)) mit günstigen Vorfaktoren herabsetzen läßt. Man beachte jedoch, daß die Matrizen nur von der Meßgeometrie abhängen und somit auch nur einmal berechnet werden müssen.

Die Berechnung von S^{⊤}g mit g ∈ range(R) kann, wenn die auftretenden Nullen weggelassen werden, nach dem oben beschriebenen wie folgt zusammengefaßt werden: 1.) Berechnung von

(S^{⊤}g)^{1} := ∑_{j=0}^{p−1}(1l|T_{q})g^{j}

2.) Berechnung von

(S^{⊤}g)^{2} := \{ \begin{matrix}∑_{j=0}^{p−1}(−1)^{j}(1l|T_{q})g^{j} & , für p gerade, \\ ∑_{j=0}^{p−1}(−1)^{j}(−1l|T_{q})g^{j} & , für p ungerade. \end{matrix}

3.) Für l = 1,…, p − 1 ist dann zu berechnen:

\begin{matrix}(S^{⊤}g)^{l′,1} := & \begin{matrix}√ & \overline{2} \\ √ & \overline{p} \end{matrix} ∑_{j=0}^{p−1} c_{jl′}(1l|T_{q})g^{j} \\ (S^{⊤}g)^{l′,2} := & \begin{matrix}√ & \overline{2} \\ √ & \overline{p} \end{matrix} ∑_{j=0}^{p−1} s_{jl′}(1l|T_{q})g^{j} \end{matrix} \}  falls l′ gerade,

bzw.

\begin{matrix}(S^{⊤}g)^{l′,1} := & \begin{matrix}√ & \overline{2} \\ √ & \overline{p} \end{matrix} ∑_{j=0}^{p−1} s_{jl′}(−1l|T_{q})g^{j} \\ (S^{⊤}g)^{l′,2} := & \begin{matrix}√ & \overline{2} \\ √ & \overline{p} \end{matrix} ∑_{j=0}^{p−1} c_{jl′}(1l| − T_{q})g^{j} \end{matrix} \}  falls l′ ungerade.

Die ”Datenfunktion“ wird also ∀ j zunächst in ihren geraden ((1l|T)g^{j}) und ihren ungeraden ((1l| − T)g^{j}) Anteil zerlegt, und dann sind abwechselnd die Fourierkoeffizienten dieser Vektoren zu berechnen. Man sieht auch, daß der p-te Fourierkoeffizient durch die von uns getroffene Wahl der Basisvektoren an die zweite Stelle gerückt worden ist. Der Aufwand für die Berechnung von S^{⊤}g ist bei Verwendung der schnellen Fouriertransformation von der Ordnung qp log(p).

Wir wollen nun noch einmal zusammenfassen, was bei der numerischen Invertierung der Radon-Transformation mittels des DAA zu tun ist: 0.) Einmalige (Vor-)Berechnung der Blockdiagonalmatrix (4.9). Dabei sind die Ver-

einfachungen ((4.11), (4.12), (4.13)) zu berücksichtigen.

Die Matrizen D_{1}^{I}, D_{3}^{I} (bzw. D_{4}^{I}), Dl^{′II}, l = 1,…, p − 1 können in O(p log(p)q^{2}) Rechenoperationen berechnet werden. Die anschließende Berechnung der Eigenvektoren und Eigenwerte erfordert O(pq^{3}) Rechenoperationen.

1.) Berechnung der transformierten rechten Seite nach 1.), 2.), 3.) auf Seite 65. Hierzu sind O(p log(p)q) Rechenoperationen erforderlich.

2.) Verallgemeinerte Lösung der 2p Gleichungssysteme:

\begin{matrix}D_{1}^{I} y^{1} & = & (S^{⊤}g)^{1} \\ D_{3}^{I} y^{2} & = & (S^{⊤}g)^{2} & (bzw. D_{4}^{I}y^{2} = (S^{⊤})^{2}) \\ Dl^{′II} y^{l′,1} & = & (S^{⊤}g)^{l′,1} \\ Dl^{′II} y^{l′,2} & = & (S^{⊤}g)^{l′,2} & l′ = 1,…, p − 1 . \end{matrix} \}  (4.14)

Wir setzen dann y = (y^{1^{⊤}}, y^{2^{⊤}}, y^{1,1^{⊤}}, y^{1,2^{⊤}},…, y^{p−1,2^{⊤}})^{⊤}.

Die Gleichungssysteme haben die Ordnung q, so daß hierfür O(pq^{2}) Rechenoperationen erforderlich sind.

3.) Rücktransformation des Vektors y in O(p log(p)q) Rechenoperationen. 4.) Berechnung von R^{∗}(Sy) auf einem hinreichend feinem Gitter. Es ist naheliegend,

ein 2q×2q-Gitter zu verwenden, so daß dies in O(q^{2}p) Rechenoperationen bewerkstelligt werden kann. Unter Berücksichtigung der Symmetrien der Diedergruppe D_{4} kann man dabei zu günstigen Vorfaktoren kommen.

Der einzig problematische Punkt ist 2.). Die bestehende Literatur ([18], [5]) schlägt eine Regularisierung nach Tychonoff-Phillips vor, wobei in ([5]) Vorschläge zur Bestimmung des Regularisierungsparameters (der für alle Systeme in (4.14) gleich gewählt wird) gemacht werden. Dabei konnte allerdings nicht berücksichtigt werden, daß die Systeme in (4.14) eine unterschiedliche innere Struktur besitzen. Verantwortlich für diese Unterschiede ist jedoch der Umstand, daß die Reihenentwicklung der Funktionen h_{l} in der Singulärwertzerlegung (3.3) mit dem l-ten Tschebyscheff-Polynom beginnt. Wir wollen diese Eigenschaft der Radon-Transformation im folgenden ausnutzen.

4.4 Verallgemeinerte Lösung der Gleichungssysteme.

Wir wollen jetzt unsere Untersuchungen aus dem Kapitel 3 für die Gewinnung einer geeigneten Strategie bei der Lösung von (4.14) ausnutzen.

Dazu ist zu allererst zu bemerken, daß der erste Schritt beim DAA, nämlich die Anwendung der Transformationsmatrix S^{⊤} äquivalent mit der diskreten Fouriertransformation in (3.7) ist. Die beiden einzigen Unterschiede bestehen darin, daß in S^{⊤} die Zerlegung in Real- und Imaginärteil schon eingerechnet ist und daß (bedingt durch die Symmetrie der Diedergruppe) die Parität auf beiden Seiten der transformierten Gleichung bereits berücksichtigt wird.

Läge nun keine weitere Störung durch den Übergang von RR^{∗} zu RR^{∗} vor, so wären nach (3.14) bzw. (3.33) die Eigenvektoren (v^{m}) ∈ IR^{q}, m = 1,…, q von D_{1}^{I} von der Form (mit u_{m} und σ_{m} aus (1.10) und (1.11)):

v_{ν}^{m} = u_{2(m−1)}(s_{ν}) , ν = 1,…, q ,  (4.15)

und die Eigenwerte λ_{m}, m = 1,…, q wären durch

λ_{m} = σ_{2(m−1)}^{2}  (4.16)

gegeben. Dies alles natürlich nur bis auf einen gemeinsamen Proportionalitätsfaktor.

Für Dl^{′II} ergäben sich die Eigenvektoren v^{m} bzw. Eigenwerte λ_{m} jeweils für l = 1, 2,…, p − 1 durch (m = 1,…, q):

v_{ν}^{m} = u_{2(m−1)+l}(s_{ν})  (4.17)
λ_{m} = σ_{2(m−1)+l}^{2} .  (4.18)

Berücksichtigt man jedoch, daß die Fouriertransformierte û_{m}(ξ) proportional zu ξ^{−1} J_{m+1}(ξ) ist, so wird klar, daß die Übereinstimmung der diskreten Eigenvektoren mit den Eigenfunktionen für große Indizes nicht mehr gut sein kann. Die Diskretisierung ist mit einer Abtastrate von h vorgenommen worden, so daß die (für große m) ausschließlich hochfrequenten Funktionen u_{m} nicht angemessen dargestellt werden können. Für niedrige Indexwerte ist die Übereinstimmung jedoch beeindruckend gut. Insbesondere ”wiederholen“ sich die Eigenvektoren und Eigenwerte beinah exakt nach dem in (4.15) und (4.16) bzw. (4.17) und (4.18) vorausgesagten (zyklischen) Schema. D. h. der zweite Eigenvektor von D_{1}^{I} stimmt fast genau mit dem ersten Eigenvektor von D2^{′II} überein, der zweite Eigenvektor von D2^{′II} ist fast exakt der erste von D4^{′II} usw. . Die Übereinstimmung bei den Eigenwerten ist ähnlich überzeugend. Dieses Wiederholungsmuster bricht aber dann bei größeren Indexwerten abrupt ab, für die Matrizen Dl^{′II} mit l > 1 sehen wir einen regelrechten Einbruch des Spektrums bei den kleinen Eigenwerten (siehe die Abbildungen 4.1 und 4.2). Dies können wir nach dem oben Gesagten nur damit erkären, daß die entsprechenden Eigenvektoren zu ausschließlich hochfrequentierten Bilddichten gehören. Dies bestätigt auch eine explizite Berechnung von R^{∗}v für Eigenvektoren v, die zu diesem Teil des Spektrums gehören. Unsere Lösungsstrategie benutzt die theoretischen Erkenntnisse aus Kapitel 3 und den offensichtlichen Zerfall des Spektrums in zwei disjunkte Mengen. Die Eigenvektoren der Matrizen Dl^{′II}, welche zu den Eigenwerten λ gehören, die wie \frac{1}{m} fallen, sind bei einer ”Lösung“ der Gleichungssysteme zu berücksichtigen. Sie gehören zu näherungsweise bandbeschränkten Eigenfunktionen von R^{∗}R. Die Eigenvektoren, welche zu Eigenwerten in der zweiten Zone gehören, entsprechen ausschließlich hochfrequenten Eigenfunktionen von R^{∗}R und sind für die Rekonstruktion näherungsweise bandbeschränkter Funktionen unwichtig. Sie können und sollten aus Stabilitätsgründen bei einer näherungsweisen Lösung der Gleichungssysteme unberücksichtigt bleiben. Die Suche nach ”optimalen“ Regularisierungsparametern, wie sie in [5] für den DAA beschrieben wurde, erweist sich also (zumindest für moderate Rauschanteile in den Daten) als unnötig! Die Untersuchungen in Kapitel 3.2 über die Interlaced-Geometrie besagen, daß wir in diesem Fall bei einer ”Lösung“ des Gleichungssystems die Entwicklung nach Eigenvektoren aus der zweiten Menge unbedingt vermeiden müssen. Eine bloße Verwendung der Tychonoff-Phillips-Regularisierung muß in diesem Fall scheitern! Wir können die Richtigkeit unserer Interpretation des Spektrums eindrucksvoll durch die Rekonstruktionen des DAA bestätigen. Wir wollen zum Schluß nicht verschweigen, daß unsere Implementation etwas von der Darstellung in Abschnitt 4.1 abweicht. Um die Situation im Kontinuierlichen möglichst genau nachzubilden, statten wir auch den IR^{N} mit einer gewichteten Norm aus. Wir betrachten also

R : L_{2}(D) → IR^{2p·2q} ,

wobei der IR^{2p·2q} über:

∥(x^{0}^{⊤},…, x^{2p−1}^{⊤})^{⊤}∥_{W} = ∑_{j=0}^{2p−1}(x^{j}, W x^{j})

normiert wird. Hierbei ist

W = diag (∥L_{h,j,1}∥_{L_{1}}^{−1},…, ∥L_{h,j,q}∥_{L_{1}}^{−1}) (j beliebig).

Dies führt dazu, daß wir auch die Matrix A = RR^{∗} bzw. deren Eigenvektoren transformieren müssen. Die Details sind elementar, und wir verzichten auf ihre Darstellung.

Abbildung 4.1: 71 Abbildung 4.2:

72 Die Abbildungen 4.1 und 4.2 zeigen Eigenwertverläufe von verschiedenen Blockdiagonalmatrizen. Hierbei haben wir die jeweils 128 Eigenwerte durch Polygonzüge verbunden. Man achte auf die zum Teil erheblich voneinander abweichenden Skalierungen! Die Abbildung 4.1 zeigt die Eigenwerte der folgenden Matrizen: in der ersten Spalte von oben nach unten D_{1}^{I}, D18^{′II}, D38^{′II}, und in der zweiten Spalte (von oben nach unten) D58^{′II},

D78^{′II}, D98^{′II}.

In Abbildung 4.2 sind die Eigenwerte von D148^{′II}, D198^{′II}, D248^{′II} (erste Spalte, von oben nach unten) und diejenigen von D298^{′II}, D348^{′II}, D398^{′II} (zweite Spalte, von oben nach unten) dargestellt. Die Parameter sind: p = 403, q = 128. Man erkennt deutlich den abrupten Abfall der Eigenwerte an der ”Bruchstelle“. Dieses markante Verhalten ist nur bei den Matrizen mit hohen Indizes etwas ”verschmiert“. Bei den Dl^{′II} mit l ≥ 390 erreichen die Eigenwerte die Bruchstelle überhaupt nicht mehr.

5 Numerische Ergebnisse

Wir wollen in diesem Kapitel den praktischen Nachweis der Rekonstruktionsqualitäten des diskreten algebraischen Algorithmus erbringen. Unser besonderes Interesse richtet sich dabei auf die Fähigkeit dieses Algorithmus, auch aus den Daten der Interlaced-Geometrie gute Rekonstruktionen berechnen zu können. Den theoretischen Untersuchungen zufolge kann dies nur gelingen, wenn wir die Relation p ≈ πq einhalten. Wir werden die Notwendigkeit dieser Bedingung durch Experimente nachweisen. Die wesentliche Philosophie bei der verallgemeinerten Lösung der Gleichungssysteme besteht darin, die Entwicklung nach Eigenvektoren von RR^{∗} auf dem richtigen Niveau abzuschneiden.

5.1 Testobjekte

Als erstes Testobjekt verwenden wir den in der Literatur über Rekonstruktionsalgorithmen allgegenwärtigen sogenannten ”Shepp-Logan-Kopf“. Dieses, von Shepp und Logan [22] entwickelte Querschnittsmodell des menschlichen Kopfes besteht aus Linearkombinationen der charakteristischen Funktionen von 11 Ellipsen. Wir verwenden dieselben Parameter wie Kruse [8], bzw. Faridani [4].

Das zweite Testobjekt besteht aus den charakteristischen Funktionen von 20 Kreisen mit unterschiedlichen Radien. An ihm lassen sich sehr schön Auflösungsgrenzen aufzeigen. Die Kreise sind durch die folgende Tabelle beschrieben.

\begin{matrix}x & y & r & f \\ 0.6 & 0.6 & 0.1 & 900 \\ 0.3 & ” & 0.09 & ” \\ 0.0 & ” & 0.08 & ” \\ -0.3 & ” & 0.07 & ” \\ -0.6 & ” & 0.06 & ” \\ 0.6 & 0.3 & 0.05 & ” \\ 0.3 & ” & 0.045 & ” \\ 0.0 & ” & 0.040 & ” \\ -0.3 & ” & 0.035 & ” \\ -0.6 & ” & 0.030 & ” \end{matrix}  \begin{matrix}x & y & r & f \\ 0.6 & 0.0 & 0.025 & 900 \\ 0.3 & ” & 0.020 & ” \\ 0.0 & ” & 0.015 & ” \\ -0.3 & ” & 0.014 & ” \\ -0.6 & ” & 0.013 & ” \\ 0.6 & -0.3 & 0.012 & ” \\ 0.3 & ” & 0.011 & ” \\ 0.0 & ” & 0.010 & ” \\ -0.3 & ” & 0.009 & ” \\ -0.6 & ” & 0.008 & ” \end{matrix}

Hierbei sind x, y die Koordinaten des Mittelpunktes, r der Radius und f der Funktionswert. Der kleinste Kreis stellt ein Detail der Größe 0.016 dar, sollte also bei Verwendung der Schrittweite h = 0.008 = \frac{1}{125} gerade noch auflösbar sein.

Bei allen Rekonstruktionen wurde das Rechteckprofil, sowie die Geometrie ohne Zentralstrahl verwendet. Wir arbeiteten mit den Datensätzen

\begin{matrix}p = 403 & p = 303 & p = 303 \\ q = 128 & q = 96 & q = 128 \end{matrix}

wobei die letzte Setzung deutlich die Relation p ≈ π q verletzt.

5.2 Rekonstruktionen

Bild 1

Das erste Bild zeigt oben links das Original des Kopfphantoms von Shepp und Logan, rechts oben die Rekonstruktion aus Standard-Daten, rechts unten diejenige aus Interlaced-Daten. Beide Rekonstruktionen sind mit den Parametern p = 403, q = 128, sowie der Abschneidegrenze σ_{ab} = 1.99 erstellt worden. Letztere ist genau die ”Bruchstelle“ in den Eigenwertdarstellungen der Diagonalmatrizen. Die Bilder sind mit einer Skala von 0 bis 90 (dies entspricht genau 10% des Maximalwertes der Originalfunktion) dargestellt, um die Detailunterschiede besser sichtbar zu machen. Das Hauptaugenmerk der Fachwelt konzentriert sich bei diesem Phantom auf die Trennung der drei kleinen Objekte im unteren Drittel des Bildes. Der horizontale Querschnitt unten links zeigt, wie gut dies dem DAA gelingt. Standard- und Interlaced-Rekonstruktion unterscheiden sich kaum.

Bild 2

Um noch einmal zu demonstrieren, wie wenig sich Standard- und Interlaced-Rekonstruk-tionen unterscheiden, zeigt Bild 2 noch einmal oben links die Standard, oben rechts die Interlaced-Rekonstruktion wie in Bild 1. Unten links sieht man dann das Differenzbild, wobei hier die Skala zwischen -5 und 5 liegt, dies sind bereits die maximalen Abweichungen! Dies ist schon verblüffend, wenn man bedenkt, daß die Interlaced-Rekonstruktion über die verallgemeinerte Lösung eines sehr schlecht konditionierten Gleichungssystems berechnet wurde, dessen rechte Seite um 100% gestört war! Unten rechts sieht man einen vertikalen Querschnitt durch die Bildmitte beider Rekonstruktionen.

In Bild 3 demonstrieren wir die Wirkung der verschiedenen Abschneidegrenzen σ_{ab} für Standard- (obere Reihe) und Interlaced- (untere Reihe) Rekonstruktionen. Die Abschneidegrenzen sind für beide Fälle gleich und sind von links nach rechts

σ_{ab} = 1.99; 1.90; 1.00; 0.10 .

Hier wird ganz augenfällig, daß das Abschneiden der Entwicklung nach Eigenvektoren für den Standard-Fall ”nur“ regularisiernde Bedeutung hat, für den Interlaced-Datensatz jedoch zwingend erforderlich ist! Die anderen Parameter sind wieder p = 403, q = 128. Der Anteil der Eigenvektoren, welche in die Entwicklung mit aufgenommen wurden, beträgt für die verschiedenen Werte von σ_{ab} (von 1.99 bis 0.10):

39.43% 39.48% 42.06% 77.19% .

Das Bild 4 zeigt deutlich die Notwendigkeit der Relation p ≈ π q für die Interlaced-Rekonstruktion. Wir haben hier p = 303, q = 128 gewählt. Die ”Bruchstelle“ im Spektrum liegt für diese Parameter bei σ = 1.57 . Die obere Reihe zeigt die Rekonstruktionen von Standard-Daten zu Abschneideparametern (von links nach rechts)

σ_{ab} = 1.57; 1.00; 2.00 .

Die Rekonstruktion ist, wie zu erwarten, unempfindlich gegenüber der Wahl des Parameters σ_{ab}. Für die Interlaced-Rekonstruktion (untere Reihe, gleiche Abschneideparameter) läßt sich jedoch keine sinnvolle Wahl von σ_{ab} treffen. Dies ist nur zu verständlich: die rechte Seite des Gleichunssystems weist einen unstrukturierten, erheblichen Fehler auf, der durch Regularisieren nicht zu korrigieren ist.

Das Bild 5 zeigt das zweite Phantom, oben links das Original, oben rechts einen horizontalen Querschnitt durch die untere Reihe der Kreise. Die mittlere Bilderreihe zeigt Rekonstruktionen vom Standard- (links) und Interlaced- (mitte) Datensatz zu den Parametern p = 403, q = 128, σ_{ab} = 2.0, sowie rechts den horizontalen Querschnitt durch die untere Kreisreihe. Die untere Bildzeile zeigt dieselbe Konstellation, nur daß wir jetzt die Parameter p = 303, q = 96, σ_{ab} = 2.01 (”Bruchstelle“) gewählt haben. Man sieht deutlich, daß der kleinste Kreis bei dem größeren Datensatz noch in voller Größe rekonstruiert wird, in der unteren Reihe jedoch schon Einbußen zu konstatieren sind.

Bild 6 dient noch einmal der Demonstration von Auflösungsgrenzen. Der DAA berechnet Rekonstruktionen der Form

f_{r}(x) = R^{∗} α(x), α ∈ IR^{2p 2q} ,

mit α = A_{r}(Rf), wobei A_{r} irgendeine Regularisierung von (RR^{∗})^{+} ist. Also gilt (A_{r} = A_{r}^{∗}):

f_{r}(x) = < α_{j,ν}, L_{h,j,ν}(x) >IR_{2p 2q}
= (f, R^{∗}A_{r}(L_{h,j,ν}(x) )_{j,ν})_{L_{2}}
= (f, K_{x}) .

Im Idealfall wäre dabei K_{x} = δ_{x}, wir können jedoch von K_{x} sinnvollerweise nur erwarten, eine stark auf den Punkt x ∈ D konzentrierte Funktion zu sein. Die ”mittlere Breite“ von K_{x} ist ein Maß für die Auflösung. Bild 6 zeigt K_{x} für drei verschiedene Punkte x und horizontale Querschnitte durch diese Funktionen. Dabei ist p = 303, q = 96 und σ_{ab} = 2.01 gewählt worden. Das ”Ausmessen“ der Kernfunktion von K_{x} ergibt auch hier, in Übereinstimmung zu Bild 5, eine Auflösungsgrenze von ≈ \frac{2}{q}. Die erreichbare Auflösung ist dabei im wesentlichen unabhängig von x.

Literaturverzeichnis

[1] Davison, M. E. (1981): A singular value decomposition for the Radon transform in n-dimensional euclidean space. Numer. Func. Anal. and Opt. 3, 321-340.

[2] Ehlich, H. und Zeller, K. (1964): Schwankungen von Polynomen zwischen Gitterpunkten. Math. Zeitschr. 86, 41-44.

[3] Ehlich, H. und Zeller, K. (1966): Auswertung der Normen von Interpolationsoperatoren. Math. Annalen 164, 105-112.

[4] Faridani, A. (1988): Abtastbedingungen und Auflösung in der Beugungstomographie. Dissertation, Westf. Wilhelms-Universität Münster.

[5] Girard, D. A. (1987): Optimal regularized reconstruction in computerized tomography. SIAM J. Sci. Stat. Comput. 8, 934-950.

[6] Gradshteyn, I. S. und Ryzhik, I. M. (1965): Table of Integrals, Series and Products. Academic Press.

[7] Hockney, R. W. (1970): The potential calculation and some applications. Meth. Comput. Physics 9, 135-211.

[8] Kruse, H. (1986): Die Auflösung von Rekonstruktionsverfahren in der Computer-Tomographie. Dissertation, Westf. Wilhelms-Universität Münster.

[9] Landau, H. J. und Pollak, H. O. (1961): Prolate spheroidal wave functions, Fourier analysis and uncertainty - II. Bell Syst. Tech. J. 40, 65-84.

[10] Landau, H. J. und Pollak, H. O. (1962): Prolate spheriodal wave functions, Fourier analysis and uncertainty - III. The dimension of the space of essentially time and band-limited signals. Bell Syst. Tech. J. 41, 1295-1336.

[11] Landau, H. J. und Widom, H. (1980): Eigenvalue distribution of time and frequency limiting. J. Math. Anal. Appl. 77, 469-481.

[12] Louis, A. K. (1981): Ghosts in tomography - The null space of the Radon transform. Math. Meth. in the Appl. Sci. 3, 1-10.

[13] Louis, A. K. (1984a): Orthogonal function series expansions and the null space of the Radon transform. SIAM J. Math. Anal. 15, 621-633.

[14] Louis, A.K. (1984b): Nonuniqueness in inverse Radon problems: The frequency distribution of the ghosts. Math. Zeitschr. 185, 429-440.

[15] Luke, Y. (1962): Integrals of Bessel Functions. McGraw-Hill.

[16] Luke, Y. (1969): The Special Functions and their Approximations, Vol. 1,. Academic Press.

[17] Mackey, G. W. (1980): Harmonic analysis as the exploitation of symmetry - A historical survey. Bull. AMS 3, 543-698.

[18] Natterer, F. (1980a): A Sobolev space analysis of picture reconstruction. SIAM J. Appl. Math. 39, 402-411.

[19] Natterer, F. (1980b): Efficient implementation of ’optimal’ algorithms in computerized tomography. Math. Meth. in Appl. Sci. 2, 545-555.

[20] Natterer, F. (1986): The Mathematics of Computerized Tomography. J. Wiley/B.G. Teubner.

[21] Schönhage, A. (1961): Fehlerfortpflanzung bei Interpolation. Num. Math. 3, 62-71.

[22] Shepp, L. A. und Logan, B. F. (1974): The Fourier reconstruction of a head section. IEEE Trans. Nucl. Sci., NS-21, 21-43.

[23] Slepian, D. und Pollak, H.O. (1961): Prolate spheroidal wave functions, Fourier analysis and uncertainty - I. Bell Syst. Tech. J. 40, 43-63.

[24] Slepian, D. (1964): Prolate spheroidal wave functions, Fourier analysis and uncertainty - IV: Extensions to many dimensions; Generalized prolate spheroidal functions. Bell Syst. Tech. J. 43, 3009-3057.

[25] Slepian, D. (1978): Prolate spheroidal wave functions, Fourier analysis and uncertainty - V: The discrete case. Bell Syst. Tech. J. 57, 1371-1430.

[26] Slepian, D. (1983): Some comments on Fourier analysis, uncertainty and modeling. SIAM Review 25, 379-393.

[27] Stein, E. M. und Weiss, G. (1971): Introduction to Fourier Analysis on Euclidean Spaces. Princeton University Press.

[28] Stiefel, E. und Fässler, A. (1979): Gruppentheoretische Methoden und ihre Anwendung. B.G. Teubner.

Lebenslauf

1. Persönliche Daten:

Name: Willehad Klaverkamp Adresse: Dechaneistr. 11, 4400 Münster Geburtstag und -ort: 11.9.1956 in Münster

2. Ausbildung:

1963-1966 Besuch der Grundschule in Münster 1966-1972 Besuch der Realschule 1972-1974 Fachoberschule Münster, Fachrichtung Technik 8.7.1974 Fachhochschulreife

WS 1974/75 - Studium an der Fachhochschule Münster, WS 1977/78 Fachbereich Bauingenieurwesen 9.2.1978 Abschluß als Dipl.-Ingenieur für Bauingenieurwesen WS 1978/79 Aufnahme des Studiums der Mathematik mit Nebenfach Physik

an der Westf. Wilhelms-Universität Münster

22.9.1980 Vordiplom 20.8.1984 Abschluß des Studiums als Diplom-Mathematiker

Thema der Diplomarbeit:

Das ”Beschränkte-Winkel-Problem“ in der Computertomographie

3. Berufliche Tätigkeit:

WS 1981/82- Unterrichtsauftrag für Mathematik an der Fachhochschule Münster WS 1986/87 SS 1981 - SS 1984 Studentische Hilfskraft im Fachbereich Mathematik

der Westf. Wilhelms-Universität Münster

1.9.1984-31.12.1985 Zivildienst im Elektroniklabor der Thorax-, Herz- und Gefäßchirugie

der Universitäts-Klinik Münster

2.1.1986-31.12.1990 Wissenschaftlicher Mitarbeiter am Institut für Numerische und instrumentelle Mathematik der Westf. Wilhelms-Universität Münster. Seit 1.1.1991 wissenschaftliche Hilfskraft

2.1.1986 Beginn der Dissertation