DE60313742T2

DE60313742T2 - Verfahren und geräte für schnelle fächerstrahl-tomographie

Info

Publication number: DE60313742T2
Application number: DE60313742T
Authority: DE
Inventors: Shu San Jose XIAO; Yoram Urbana BRESLER; David C. Jr. Dexter MUNSON
Original assignee: University of Illinois; University of Illinois at Urbana Champaign
Current assignee: University of Illinois; University of Illinois at Urbana Champaign
Priority date: 2002-02-28
Filing date: 2003-02-28
Publication date: 2008-01-24
Anticipated expiration: 2023-03-01
Also published as: AU2003230578A8; US20030161443A1; AU2003230578A1; EP1478274A4; US6771732B2; CN100475146C; JP2005518892A; WO2003075036A3; EP1478274B1; EP1478274A2; ATE362152T1; DE60313742D1; WO2003075036A2; JP4550426B2; CN1658795A

Description

TECHNISCHES GEBIET
Die Erfindung betrifft Verfahren und Vorrichtungen für die Divergentstrahl-Tomographie und insbesondere Verfahren und Vorrichtungen für die Divergentstrahl-Tomographie, welche die Rechenkomplexität verringern ohne eine mit den Augen wahrnehmbare Verschlechterung oder signifikante Zahlenungenauigkeit.
TECHNISCHER HINTERGRUND
Tomographische Rekonstruktion ist eine bekannte Technik, die praktisch allen diagnostischen Bildgebungsmodalitäten zugrunde liegt, einschließlich Röntgen-Computertomographie (CT), Positronenemissionstomographie (PET), Einzelphotonenemissionstomographie (SPECT) und bestimmter Aufnahmeverfahren für das Magnetresonanz-Imaging (MRI). Sie wird zudem weit verbreitet bei der Herstellung für die nichtdestruktive Untersuchung (non-destructive evaluation, NDE) und seit neuestem zur Gepäckinspektion auf Flughäfen eingesetzt. Der vorherrschende Datenerfassungsmodus bei derzeitigen CT-Scannern ist die Fächerstrahl-Geometrie, oft in Kombination mit einer helixförmigen oder spiralförmigen Abtastung.
Im Interesse eine schnelleren Bilderfassung bewegt sich die CT-Scanner-Technologie auf eine Kegelstrahl-Erfassung unter Verwendung von Flächendetektoren zu. Viele der hergestellten CT-Scanner haben bereits einen Mehrreihen-Detektor (in der Regel 16 Reihen), was zu einer (Flachwinkel-)Kegelstrahl-Erfassung führt. Flächendetektoren für die Weitwinkel-Kegelstrahl-Erfassung werden in den nächsten Jahren für medizinische CT-Scanner eingeführt. Zurzeit besteht ein noch ungestillter Bedarf an tomographischer Hochgeschwindigkeits-3D-Bildgebung zur Ermittlung versteckter Waffen und Explosionsvorrichtungen. Bei anderen Modalitäten, wie PET oder SPECT, ist die Kegelstrahl-Erfassung bereits das vorherrschende Format.
Das Rekonstruktionsproblem bei der Fächerstrahl- und Kegelstrahl-CT besteht darin, ein 2D- oder 3D-Bild von einem Objekt zurückzugewinnen aus einem Satz von dessen Kurvenintegral-Projektionen entlang Kurven, die in Form eines Fächers bzw. Kegels von einer Source divergieren, die sich auf einem festgelegten Erfassungsrichtstrahl oder Orbit bewegt. Daher sind die Fächerstrahl- und die Kegelstrahl-Geometrie Fälle der Divergentstrahl-Geometrie.
Das Verfahren der Wahl für die tomographische Rekonstruktion ist die gefilterte Rückprojektion (filtered backprojection, FBP) oder die Rückprojektion mit Faltung (Konvolution) (convolution backprojection, CBP), die beide einen gewichteten Rückprojektionsschritt verwenden. Dieser Schritt ist der Rechen-Flaschenhals bei der Technik, wobei die Rechenanforderungen ein Ausmaß von N³ für ein N×N-Pixel-Bild in 2D und ein Ausmaß von mindestens N⁴ für ein N×N×N-Voxel-Bild in 3D annehmen. Somit führt eine Verdopplung der Bildauflösung von N auf 2N zu einer ungefähr 8-fachen (oder in 3D 16-fachen) Zunahme bei der Berechnung. Zwar sind die Computer viel schneller geworden, aber mit dem Aufkommen neuer Technologie, die immer größere Mengen an Daten in Echtzeit sammeln können (z.B. Herzbildgebung mit Mehrreihen-Detektoren, Interventionsbildgebung), und der Hinwendung zu einer 3D-Kegelstrahl-Geometrie gibt es immer noch einen steigenden Bedarf an schnellen Rekonstruktionstechniken. Die schnelle Rekonstruktion kann den Bilderstellungsprozess beschleunigen, die Kosten für einen speziellen Bildrekonstruktionscomputer (siehe unten) verringern oder beides.
Der Tandemschritt zur Rückprojektion ist die Reprojektion, also der Vorgang der Berechnung der Projektionen eines elektronisch gespeicherten Bildes. Auch dieser Vorgang spielt bei der tomographischen Rekonstruktion eine grundlegende Rolle. Eine Kombination von Rückprojektion und Reprojektion kann auch zur Konstruktion schneller Rekonstruktionsalgorithmen für das Problem des langen Objektes bei der helixförmigen Kegelstrahl-Geometrie verwendet werden, das der Schlüssel für die praktische 3D-Bildgebung von Menschen darstellt. Zudem ist es bei verschiedenen Anwendungen vorteilhaft oder sogar erforderlich, dass iterative Rekonstruktionsalgorithmen verwendet werden, in denen Rückprojektions- und Reprojektionsschritte zur Rekonstruktion eines einzigen Bildes mehrmals durchgeführt werden. Eine Beschleunigung der Rückprojektions- und Reprojektionsschritte bestimmt die ökonomische Durchführbarkeit solcher iterativen Verfahren.
Bei der Erläuterung schneller Verfahren muss man zwischen den beiden hauptsächlichen Formaten tomographischer Daten unterscheiden: (i) Parallelstrahl und (ii) Divergentstrahl und seine Spezialfälle. Das Parallelstrahl-Format ist für theoretische und Rechenmanipulationen am besten geeignet. Divergentstrahl ist dagegen das Format, das man am häufigsten in kommerziellen medizinischen und industriellen Scannern findet. Die Fächerstrahl-Rekonstruktion war ursprünglich ein 2D-Verfahren und ist die Schlüsselkomponente bei Verfahren des Standes der Technik für die helixförmige und die Mehrschicht-helixförmige-3D-(volumetrische) Rekonstruktion. Es wird erwartet, dass neue und zukünftige Bauweisen zu Kegelstrahl-Geometrie übergehen. Somit ist bei 2D und 3D die Divergentstrahl-Geometrie das vorherrschende Bildgebungsformat und bleibt dies wahrscheinlich auch in der näheren Zukunft.
Begrenzt auf einen großen Source-Abstand von dem Objekt, reduziert sich die Divergentstrahl-Geometrie auf die Parallelstrahl-Geometrie. Deshalb sind alle Datenver arbeitungsverfahren für die Divergentstrahl-Geometrie (einschließlich der erfindungsgemäßen Verfahren) auch auf die Parallelstrahl-Geometrie anwendbar. Das Gegenteil gilt jedoch nicht: Bei den kleinen oder mäßigen Source-Abständen, die man in der Praxis findet, sind die beiden Geometrien genügend verschieden, dass Parallelstrahl-Verarbeitungsverfahren schlechtere oder nicht zufrieden stellende Ergebnisse liefern, werden sie direkt auf Divergentstrahl-Daten angewendet. Somit erfordern Divergentstrahl-Geometrien eine andere Verarbeitung als Parallelstrahl-Geometrien.
Die Erfindung betrifft die Divergentstrahl-Tomographie. Einige Rekonstruktionsverfahren beruhen auf einem als Rebinning bezeichneten Verfahren, wobei Divergentstrahl-Daten zu Parallelstrahl-Daten umgeordnet (oder interpoliert) werden, die dann durch einen Parallelstrahl-Algorithmus verarbeitet werden. Andere Verfahren für 3D-Divergentstrahl-Geometrien beruhen auf der Transformation der Kegelstrahl-Daten in 3D-Radontransformationsdaten (oder ein Derivat davon). Stattdessen verarbeiten native Divergentstrahl-Verfahren (von denen die Erfindung eine schnelle Version darstellt) definitionsgemäß die Divergentstrahl-Daten direkt ohne vorheriges Rebinning zu Parallelstrahl-Daten oder Transformieren der Daten in die 3D-Radontransformationsdomäne. Das Rekonstruktionsverfahren bei nativen Divergentstrahl-Verfahren besteht aus der Vorverarbeitung der Projektionen (z.B. Gewichten und Filtern) und einem anschließenden gewichteten Divergentstrahl-Rückprojektionsschritt und eventuell einem Divergentstrahl-Reprojektionsschritt oder einer Reihe solcher Schritte. Man nennt solche Verfahren daher gefilterte Divergentstrahl-Rückprojektions-(DB-FBP-)Algorithmen. Ihre Berechnung wird wie im Fall des Parallelstrahls gewöhnlich von den Rückprojektions- und Reprojektionsschritten beherrscht.
Eine Divergentstrahl-Reprojektion kann durchgeführt werden, indem zuerst eine Parallelstrahl-Reprojektion und dann ein Rebinning der Ergebnisse durchgeführt wird. Auch hier ist es jedoch vorteilhaft, wenn ein nativer Divergentstrahl-Reprojektionsalgorithmus verwendet wird, der kein Rebinning einsetzt.
Die Nachteile des Rebinning sind u.a. mögliche Artefakte und zusätzliche Berechnung, wodurch die Beschleunigung beschränkt wird. Rebinning erfordert zudem die Erfassung und Manipulation eines großen Teils der Daten, bevor die Verarbeitung beginnen kann, wodurch wiederum die Geschwindigkeit beschränkt wird. Verfahren, die auf einer Transformation in die 3D-Radon-Domäne beruhen, haben ähnliche Nachteile. Somit besteht ein Bedarf an Verfahren und Vorrichtungen für die Divergentstrahl-Tomographie, die diese Nachteile überwinden, hochflexibel sind und große Leistungsgewinne verglichen mit Verfahren bereitstellen, die eine herkömmliche Divergentstrahl-Rückprojektion einsetzen.
Spezial-Hardware ist der herkömmliche Weg zur Beschleunigung des Rückprojektionsverfahrens gewesen. Spezielle Computertypen unter Verwendung anwendungsspezifischer Chips oder mehrerer Prozessoren oder Kombinationen davon sind derart gebaut, dass die die notwendigen Berechnungen durchführen.
Nachteile dieses Ansatzes sind u.a. die Kosten für mehrere Prozessoren oder anwendungsspezifische Hardware und die Tatsache, dass mit der starken Steigerung in der Leistung von Allzweckcomputern die erforderlichen Spezialarchitekturen schnell obsolet werden. Somit besteht ein Bedarf an schnellen Verfahren für die Divergentstrahl-Geometrie, die keine Spezial-Hardware erfordern und leicht auf Standard-Seriell- oder -Parallelarchitekturen zu implementieren sind, so dass sie kostengünstiger werden.
AUFGABEN DER ERFINDUNG
Eine Aufgabe der Erfindung ist also die Bereitstellung neuer und besserer Verfahren und Vorrichtungen für die Divergentstrahl-Tomographie.
Eine weitere Aufgabe ist die Bereitstellung neuer und besserer Verfahren und Vorrichtungen für die Divergentstrahl-Tomographie, die die Rechenkomplexität ohne eine mit dem Auge wahrnehmbare Verschlechterung oder signifikante Zahlenungenauigkeit verringern.
OFFENBARUNG DER ERFINDUNG
Unter einem Aspekt der Erfindung wird ein Verfahren vorgeschlagen zum Erzeugen eines elektronischen Bildes aus einem vorverarbeiteten Divergentstrahl-Sinogramm, das einer Rückprojektion unterworfen werden kann, wobei das Sinogramm eine Sammlung von Divergentstrahl-Projektionen ist. Das Verfahren kann die folgenden Schritte umfassen: Unterteilen des Sinogramms in eine Mehrzahl Unter-Sinogramme; Durchführen einer gewichteten Rückprojektion der Unter-Sinogramme in das an dem Objekt fixierte globale Koordinatensystem, wodurch eine Mehrzahl entsprechender Unter-Bilder an den richtigen Stellen im globalen Koordinatensystem erzeugt wird; und Aggregieren der Unter-Bilder, wodurch das elektronische Bild erzeugt wird. Die Unterteilung des Sinogramms in Unter-Sinogramme kann rekursiv erfolgen, bis jedes Unter-Sinogramm ein Unter-Bild mit gewünschter Größe darstellt, wobei die Unterteilungsschritte eine gewünschte Anzahl genauer Unterteilungen und eine gewünschte Anzahl annähernder Unterteilungen umfassen.
Unter einem anderen Aspekt der Erfindung wird ein Verfahren vorgeschlagen zum Reprojizieren eines elektronischen Bildes, d.h. zum Erzeugen eines Divergentstrahl- Sinogramms von dem Bild. Das Verfahren kann die folgenden Schritte enthalten: Unterteilen des Bildes in eine Mehrzahl Unter-Bilder; Berechnen von Unter-Sinogrammen von jedem Unter-Bild in dem globalen Koordinatensystem; und Aggregieren der Unter-Sinogramme, wodurch das Sinogramm erzeugt wird. Die Unterteilung des Bildes kann rekursiv erfolgen, bis jedes Unter-Bild eine gewünschte Größe hat. Die Berechnung der Unter-Sinogramme kann in einer gewünschten Anzahl an Ebenen der Rekursion annähernd und in den verbleibenden Ebenen der Rekursion genau sein.
Verglichen mit herkömmlichen Rückprojektions- und Reprojektionsalgorithmen liefern diese Verfahren eine ähnliche Beschleunigung wie die FFT. Man erzielt erhebliche Recheneinsparungen ohne wahrnehmbare Verschlechterung oder signifikanten Verlust an Zahlengenauigkeit.
KURZE BESCHREIBUNG DER ZEICHNUNGEN
Die oben erwähnten und weitere Merkmale der Erfindung und die Weise, wie sie erzielt werden, werden besser ersichtlich und die Erfindung selbst wird am besten verständlich anhand der folgenden Beschreibung einer Ausführungsform der Erfindung in Verbindung mit den beigefügten Zeichnungen. Es zeigt/zeigen:
1 ein Blockdiagramm von einer Vorrichtung, die mit der Erfindung verwendet wird;
2 ein Schema einer planaren, gleichmäßig beabstandeten Kegelstrahl-Geometrie;
3 ein Schema einer planaren, gleichmäßig beabstandeten, helixförmigen Kegelstrahl-Geometrie;
4 ein Schema zur Veranschaulichung einer Kegelstrahl-Tomographie mit kreisförmigem Source-Orbit und einem planaren, gleichmäßig beabstandeten Detektor;
5 ein Schema einer colinearen, gleichmäßig beabstandeten Fächerstrahl-Geometrie;
6 ein Schema zur Veranschaulichung einer Kegelstrahl-Tomographie mit kreisförmigem Source-Orbit und einem zylindrischen Detektor;
7 ein Schema einer gleichwinkelingen Fächerstrahl-Geometrie;
8 ein Schema einer colinearen, gleichmäßig beabstandeten Fächerstrahl-Projektionsgeometrie für ein M×M-Unter-Bild;
9 ein Schema einer Unterteilung eines Bildes in Unter-Bilder;
10 ein Schema zur Veranschaulichung einer genauen Zerlegung einer Fächerstrahl-Rückprojektion;
11 ein Fließdiagramm von dem Verfahren, das durch den Sinogramm-Zerlegungsoperator O[L,M,δ ⇀] durchgeführt wird;
12 ein Schema von einer annähernden Zerlegung einer Fächerstrahl-Rückprojektion;
13 ein Schema einer dyadischen Zwei-Ebenen-Bildunterteilung;
14 ein Schema einer annähernden Zwei-Ebenen-Zerlegung einer Fächerstrahl-Rückprojektion;
15 ein Schema einer gemischten (einer genauen, einer annähernden) Zwei-Ebenen-Zerlegung einer Fächerstrahl-Rückprojektion;
16 einen Pseudo-Code für eine hierarchische schnelle Fächerstrahl-Rückprojektion;
17 ein Schema zur Veranschaulichung einer genauen Zerlegung einer Fächerstrahl-Rückprojektion;
18 ein Fließdiagramm von dem Verfahren, das durch den Projektionenvermehrungsoperator ν[L,M,δ ⇀] durchgeführt wird;
19 ein Schema von einer annähernden Zerlegung einer Fächerstrahl-Reprojektion;
20 ein Schema einer annähernden Zwei-Ebenen-Zerlegung einer Fächerstrahl-Reprojektion;
21 ein Schema einer gemischten (einer genauen, einer annähernden) Zwei-Ebenen-Zerlegung einer Fächerstrahl-Reprojektion;
22 ein Schema einer gleichwinkeligen Fächerstrahl-Projektionsgeometrie für ein Unter-Bild;
23 ein Schema einer planaren, gleichmäßig beabstandeten Kegelstrahl-Projektionsgeometrie für ein Unter-Bild;
24 ein Schema von einer Unterteilung eines 3D-Bildes in Unter-Bilder;
25 ein Schema einer nicht-isotropen Teilung eines 3D-Bildes in Unter-Bilder;
26 ein Schema einer nichtgleichmäßigen rekursiven Zwei-Ebenen-Unterteilung eines 3D-Bildes;
27(a)-27(d) Bilder, in denen genaue und schnelle Fächerstrahl-Rückprojektionen für colineare, gleichmäßig beabstandete Detektoren verglichen werden;
28(a)-28(d) Bilder, in denen genaue und schnelle Fächerstrahl-Rückprojektionen für gleichwinkelige Detektoren verglichen werden;
29(a) und 29(b) Diagramme durch Schnitte der Rückprojektionen der 28(a)-28(d);
30(a) und 30(b) einen Vergleich von Punktverbreiterungsfunktionen;
31(a) und 31(b) Vergleiche von Diagrammen von Schnitten durch Punktverbreiterungsfunktionen;
32 eine Reihe von Bildern von 3D-Kegelstrahl-Rekonstruktionen;
33 eine weitere Reihe von Bildern von 3D-Kegelstrahl-Rekonstruktionen;
34 eine Reihe von Diagrammen von Schnitten durch die Bilder der 32 und
35 eine Reihe von Diagrammen von Schnitten durch die Bilder der 33.
EINGEHENDE BESCHREIBUNG EINSCHLIEßLICH DER BESTEN ART UND WEISE
Die Erfindung findet Anwendung bei einer Reihe von Bildgebungsvorrichtungen, einschließlich CT-Scannern. Eine übliche Bildgebungsvorrichtung 10 (1) umfasst einen Scanner 12, der Daten von einem Objekt, wie einem Kopf, erfasst und Rohdaten, die Divergentstrahl-Projektionen 14 entsprechen, an einen Projektionsvorprozessor 16 sendet. Der Projektionsvorprozessor 16 wendet auf die Daten verschiedene Umwandlungen, Normalisierungen und Korrekturen sowie Gewichtung und Filtern an, die verschiebungsvariierend (shift varying) sein können. Die Ausgabe des Projektionsvorprozessors 16 ist ein Divergentstrahl-Sinogramm₁ 18, das in einen Sinogramm-Update-Prozessor 20 eingespeist wird. Der Sinogramm-Update-Prozessor 20 modifiziert möglicherweise das eingegebene Sinogramm, 18 unter Verwendung von Information aus dem Divergentstrahl-Sinogramm₂ 34, wobei zum Beispiel verschiedene Artefakte korrigiert werden, einschließlich Strahlaufhärtung, oder als Teil eines Mehrschritt- oder iterativen Rekonstruktionsverfahrens.
Die Ausgabe des Sinogramm-Update-Prozessors 20 ist ein Divergentstrahl-Sinogramm₃ 22, das in einen schnellen Rückprojektionsprozessor 24 eingegeben wird. Der schnelle Rückprojektionsprozessor 24 ist gewöhnlich ein Computer oder Spezial-Hardware oder eine Kombination davon eines beliebigen geeigneten Typs, der/die derart programmiert ist, dass er/sie die hier beschriebenen genauen und annähernden Rückprojektionsalgorithmen durchführt.
Die Ausgabe des schnellen Rückprojektionsprozessors 24 ist ein elektronisches Bild₁ 26, das in einen Bildkonditionierungsprozessor 28 eingegeben wird. Der Bildkonditionierungsprozessor 28 führt eine notwendige Nachbearbeitung des elektronischen Bildes, möglicherweise einschließlich Identifizierung und Entfernung von Artefaktbildern, oder der Bilder für die weitere Verarbeitung in einem Mehrschritt- oder iterativen Rekonstruktionsverfahren durch.
Wenn gewünscht, kann der Bildkonditionierungsprozessor 28 ein elektronisches Bilde 30 erzeugen, das in einen schnellen Reprojektionsprozessor 32 eingespeist wird. Der schnelle Reprojektionsprozessor 32 ist in der Regel ein Computer oder Spezial-Hardware oder eine Kombination davon eines beliebigen geeigneten Typs, der/die derart programmiert ist, dass er/sie die hier beschriebenen genauen und annähernden Reprojektionsalgorithmen durchführt. Wenn gewünscht, kann dieser Prozessor den gleichen Computer und die gleiche Hardware verwenden, die vom Rückprojektionsprozessor 24 eingesetzt werden.
Die Ausgabe des schnellen Reprojektionsprozessors 32 ist ein Divergentstrahl-Sinogramm₂ 34, das in den Sinogramm-Update-Prozessor 20 zurück gespeist wird. Das Rückprojektions-/Reprojektionsverfahren kann weiterlaufen, bis geeignete Ergebnisse erhalten worden sind. Eine Reprojektion wird nicht immer benötigt, kann jedoch in vielen Situationen hilfreich sein.
Ist das elektronische Bild₁ 26 geeignet, erzeugt der Bildkonditionierungsprozessor 28 ein elektronisches Bild₃ 36, das in eine Speicher-/Analyse-/An-zeigevorrichtung 38 eingespeist wird. Es wird erwogen, dass das elektronische Bild₃ 36 in einem Computerspeicher gespeichert und/oder elektronisch zum Beispiel auf Anomalien oder gefährliche Materialien analysiert und/oder angezeigt und/oder in einer beliebigen sichtbaren Form gedruckt werden kann.
Die erfindungsgemäßen schnellen Divergentstrahl-Rückprojektionsalgorithmen beruhen auf dem schnellen hierarchischen Rückprojektionsalgorithmus (FHBP) für die Parallelstrahl-Tomographie, der im U.S.-Patent 6 282 257 beschrieben ist. Ebenso beruhen die erfindungsgemäßen schnellen Divergentstrahl-Reprojektionsalgorithmen auf den schnellen hierarchischen Reprojektionsalgorithmen (FHRP) für die Parallelstrahl-Tomographie, die in den U.S.-Patenten 6 263 096 und 6 351 548 beschrieben sind. Rückprojektion und Reprojektion hängen eng zusammen, und ihre schnellen Algorithmen rühren von denselben Prinzipien her. Die am häufigsten verwendete schnelle Rückprojektion wird eingehender beschrieben.
Zum Verständnis des FHBP-Algorithmus betrachte man ein N×N-Bild mit einem daran gebundenen globalen Koordinatensystem, wobei der Ursprung sich im Mittelpunkt des Bildes befindet. Man stelle sich vor, dass dieses große N×N-Bild in vier kleinere Bilder zerlegt wird mit der Größe N/2×N/2, die man als Unter-Bilder bezeichnet und die sich jeweils in einem anderen Quadranten des globalen Koordinatensystems befinden. Zwei Eigenschaften der Radontransformation werden bei der Ableitung des FHBP-Algorithmus verwendet. (i) Die "Fliegen"-(bow-tie-)Eigenschaft besagt, dass für ein am Ursprung des globalen Koordinatensystems zentriertes Bild mit der halben Größe der spektrale Träger der Projektionen in Bezug auf die Winkelvariable ebenfalls um die Hälfte abnimmt. Dies impliziert, dass ein zentriertes Bild mit der halben Größe aus einem Sinogramm rekonstruiert werden kann, das die Hälfte der Anzahl an Projektionen umfasst. (ii) Die Verschiebungseigenschaft besagt, dass die Parallelprojektionen eines verschobenen Bildes verschobenen Parallelprojektionen des ursprünglichen Bildes entsprechen.
Man nehme an, dass ein Sinogramm mit P gefilterten Projektionen für die Rekonstruktion des gesamten Bildes f verfügbar ist. Zur Rekonstruktion von einem der Unter-Bilder f_i werden die P Projektionen zunächst auf den Träger der Projektion des Unter-Bildes abgeschnitten, dann verschoben, so dass die einer zentrierten Version des Unter-Bildes entsprechen, und in Bezug auf die Winkelrichtung schließlich auf P/2 Projektionen dezimiert. Die zentrierte Version des Unter-Bildes f_i wird dann aus diesen P/2 Projektionen rekonstruiert, und das Unter-Bild wird an seine korrekte Stelle im globalen Koordinatensystem zurück verschoben. Dies wird für jedes der vier Unter-Bilder durchgeführt, die, wenn sie miteinander aggregiert werden, eine Rekonstruktion des gesamten Bildes f liefern. Insgesamt 4c(P/2)(N/2)² = cPN²/2 Arbeitsschritte sind zum Rekonstruieren von vier Unter-Bildern notwendig, wobei c eine Konstante ist. Dies verringert die ursprünglichen Rekonstruktionskosten auf die Hälfte. Durch rekursive Anwendung der Zerlegung, wodurch die Bildgröße in jedem Schritt um die Hälfte verringert wird und Unter-Bilder mit der halben Größe unter Verwendung der Hälfte der Anzahl an Projektionen hergestellt werden können, verringern sich die Gesamtberechnungskosten auf o(N²logN).
Das Tandem dieses Zerlegungs- und Aggregationsansatzes wird bei dem schnellen hierarchischen Parallelstrahl-Reprojektionsalgorithmus FHRP verwendet, der in den U.S.-Patenten 6 263 096 und 6 351 548 beschrieben ist.
Die Erfindung für die Divergentstrahl-Rückprojektion und -Reprojektion basiert auf der Idee einer hierarchischen Zerlegung des Bildes und der anschließenden Rekonstruktion (oder Reprojektion) von Unter-Bildern unter Verwendung eines Divergentstrahl-Unter-Sinogramms, das weniger Divergentstrahl-Projektionen umfasst. Versucht man jedoch, die Ideen der Parallelstrahl-FHBP oder -FRBP direkt auf das Divergentstrahl-Szenario anzuwenden, tritt sofort die Schwierigkeit auf, dass die obige Verschiebungseigenschaft (ii) nicht mehr gilt. D.h. die Divergentstrahl-Projektionen eines verschobenen Unter-Bildes entsprechen nicht den verschobenen Divergentstrahl-Projektionen des ursprünglichen Unter-Bildes. Weil Parallelstrahl-FHBP und -FRBP diese Eigenschaft nutzen, muss der Parallelstrahl-Ansatz für das Divergentstrahl-Szenario modifiziert werden.
Die erfindungsgemäßen Algorithmen verwenden ein anderes Verfahren für die Zerlegung, bei dem es nicht zum Verschieben der Unter-Bilder in Bezug auf den Mittelpunkt des globalen Koordinatensystems kommt, das am Objekt, d.h. am Bild der vollen Größe, fixiert ist. Wie bei der Parallelstrahl-FHBP, haben die neuen Divergentstrahl-Rückprojektionsalgorithmen zwei Arten der Zerlegungen, eine genaue, aber langsame, und eine annähernde, aber schnelle. Die genaue Divergentstrahl-Zerlegung der neuen Algorithmen braucht überhaupt keine Datenverschiebung zu umfassen; sie umfasst nur das Abschneiden des Sinogramms in Unter-Sinogramme, die den Unter-Bildern entsprechen. Diese Unter-Sinogramme werden dann in das globale Koordinatensystem rückprojiziert unter Erzeugung der entsprechenden Unter-Bilder an den richtigen Stellen im globalen Koordinatensystem. Die annähernde Divergentstrahl-Zerlegung verwendet ein anderes Verfahren als es bei der Parallelstrahl-FHBP eingesetzt wird zum Verringern der Anzahl an Projektionen, die zur Rekonstruktion eines kleineren Unter-Bildes verwendet werden. Bei dem neuen Verfahren werden nur die Projektionen, nicht aber die Unter-Bilder abgeschnitten und verschoben. Das neue Verfahren kann Gewichtung umfassen, die bei der FHBP nicht verwendet wird. Analoge Aussagen gelten auch für das neue Verfahren für die Divergentstrahl-Reprojektion.
Zur Rekonstruktion eines Unter-Bildes der halben Größe werden annähernd abgeschnittene (gewichtete und gefilterte) Divergentstrahl-Projektionen um einen Abstand verschoben, der bestimmt wird durch die Position der Projektion des Unter-Bild-Mittelpunkts im globalen Koordinatensystem. Dann werden sie in den Source-Positionsparametern durch zwei dezimiert und um den gleichen Betrag zurück verschoben. Das Unter-Bild kann dann aus diesem reduzierten Salz an Projektionen in guter Annäherung rekonstruiert werden. Für die exakte und die annähernde gewichtete Rückprojektion auf ein Unter-Bild erfolgen alle Berechnungen, die Bildkoordinaten umfassen, einschließlich der Gewichtung, in den absoluten oder globalen Bildkoordinaten, wobei man sich auf den Ursprung des Bildes mit der vollen Größe bezieht. Eine Kombination der genauen und der annähernden Zerlegungen in rekursiver Form führt zu schnellen hierarchischen Rückprojektionsalgorithmen für Divergentstrahl-Geometrien. Wie ihre Parallelstrahl-FHBP-Gegenstücke bieten diese Algorithmen eine Ausgewogenheit zwischen Berechnung und Genauigkeit, indem die Anzahl der genauen gegenüber den annähernden Zerlegungen ausgewählt wird und verschiedene Algorithmusparameter eingestellt werden.
Divergentstrahl-Bildgebungsgeometrien
Die allgemeine Divergentstrahl-Erfassungsgeometrie für die 3D-Tomographie umfasst einen Scheitel (in der Regel eine Source) divergenter Strahlen, der sich auf einem Richtstrahl im 3D-Raum um ein abgebildetes Objekt bewegt, und eine Detektoroberfläche, auf der Kurvenintegrale entlang der Source-Strahlen durch das abgebildete Objekt gemessen werden. Weitere Bildgebungsmodalitäten mit Konvergent- anstelle von Divergentstrahl-Geometrien, wie sie bei der Emissionstomographie (SPECT und PET) auftreten, können in diese Geometrie transformiert werden.
Die Erfassungsgeometrie wird gewöhnlich nach der Gestalt des Detektors und dem damit einhergehenden Strahlabtastschema sowie nach dem Source-Richtstrahl oder -Orbit klassifiziert. Die beiden üblichsten Detektorklassen sind (i) eine planare Detektoroberfläche mit gleichmäßig beabstandeten Detektoren und (ii) eine zylindrische Detektoroberfläche, wobei die Strahlen im Azimuth-Winkel gleichmäßig beabstandet und in der Höhe linear gleichmäßig beabstandet abgetastet werden. Bei der 3D-Bildgebung mit einem Oberflächendetektor nennt man Divergentstrahl-Geometrien auch Kegelstrahl-Geometrien.
Bei dem Spezialfall der 2D-Divergentstrahl-Erfassung ist die Detektoroberfläche auf eine einzige Linie oder einen Bogen von Detektoren reduziert, wodurch man die so genannten colinearen gleichmäßig beabstandeten bzw. gleichwinkeligen Fächerstrahl-Geometrien erhält. Diese Fächerstrahl-Geometrien werden in gegenwärtigen kommerziellen Systemen weit verbreitet eingesetzt.
Unter den verschiedenen vorgeschlagenen Source-Orbits sind die üblichsten der einzelne kreisförmige Orbit und der helixförmige (oder Spiral-)Orbit. Von praktischem Interesse sind jedoch auch anderen Orbits, beispielsweise zwei Kreise oder ein Kreis und eine Linie. In der Praxis treten auch willkürlich gestörte Versionen dieser verschiedenen Orbits auf.
Das erfindungsgemäße Verfahren, die erfindungsgemäßen Algorithmen und die erfindungsgemäße Vorrichtung sind auf die allgemeine Divergentstrahl-Geometrie anwendbar. Sie können auf jede der zuvor erwähnten besonderen Geometrien oder auf andere Fälle der Divergentstrahl-Geometrie zugeschnitten werden. Ebenso sind die erfindungsgemäßen Verfahren auf jeden mathematisch gleichwertigen Datensatz anwendbar, einschließlich Konvergentstrahl-Geometrien, aber nicht darauf beschränkt. Zu Veranschaulichungszwecken werden die Verfahren für die folgenden Fälle ausführlicher dargestellt:

1. Einzelner kreisförmiger Source-Orbit (a) colineare gleichmäßig beabstandete Fächerstrahl-Geometrie (b) gleichwinkelige Fächerstrahl-Geometrie (c) planare Kegelstrahl-Geometrie
2. Allgemeine Kegelstrahl-Geometrie mit planarem Detektor
3. Helixförmiger Source-Orbit

Es wird gezeigt, wie das allgemeine Verfahren auf diese Fälle zugeschnitten wird. Daraus wird ersichtlich, dass der Algorithmus für eine Detektorform und einen Source-Orbit sich leicht für eine andere Detektorform und einen anderen Source-Orbit modifizieren lässt.
Divergentstrahl-Bildgebungsgeometrien mit willkürlichem Orbit
In 2 ist die planare gleichmäßig beabstandete Kegelstrahl-Geometrie dargestellt, bei der die Detektoren auf einer planaren Oberfläche gleichmäßig beabstandet sind. In diesem Fall lässt man die Detektorebene für jede Source-Position eine willkürliche Position und Orientierung einnehmen.
Die Source S von divergenten Strahlen bewegt sich um das 3D-Objekt auf einem Orbit, der durch den Source-Positionsvektor s →(θ) = [x(θ),y(θ),z(θ)]^T relativ zum Mittelpunkt o eines an dem Objekt fixierten globalen Koordinatensystems definiert ist. wobei die Source-Position auf diesem Orbit durch den Parameter θ genauer angegeben wird, wobei θ_min ≤ θ ≤ θ_max. Für jedes θ wird ein an den Detektor gebundenes 3D-Koordinatensystem definiert, so dass der Ursprung o_θ die orthogonale Projektion des Orbitpunktes s →(θ) auf die Detektorebene ist. Die orthonormalen Vektoren t ^1 und t ^2 werden innerhalb der Detektorebene gewählt, wohingegen ŵ der Einheitsvektor ist, der orthogonal zur Detektorebene ist und von o_θ zur Source hin zeigt. Ein Punkt t → = t₁t ^₁ + t₂t ^₂ in der Detektorebene ist anhand seiner Koordinaten (t₁,t₂) definiert. In diesem System wird der Source-Orbit-Punkt definiert als D(θ)ŵ, wobei D(θ) der Abstand zwischen der Source und dem Detektor ist. (Bei SPECT ist D(θ) die fokale Länge des Kegelstrahl-Kollimators).
Die Position in dem 3D-Bild (d.h. dem Objekt) sei gegeben durch r → = [x,y,z]^T, und die t₁,t₂-Position des Schnittpunktes mit der Detektorebene des Source-Strahls, der durch den Punkt r → geht, sei τ →(θ,r →) = [τ₁(θ,r →),τ(θ,r →)]^T (siehe 2). Die Divergentstrahl-Projektion (Pf)(θ,t₁,t₂) von dem Objekt f an der Source-Orbit-Position θ und der Detektorposition (t₁,t₂) ist das Kurvenintegral entlang des Source-Strahls mit den Parametern (θ,t₁,t₂) und ist gegeben durch
wobei δ(r →) die Dirac'sche Deltafunktion ist. Die gleiche Erfassungsgeometrie kann zur Beschreibung einer Situation verwendet werden, wobei die Source-Strahlen zur Source hin konvergieren und nicht von ihr divergieren. Daher sind die erfindungsgemäßen Verfahren ebenso auf Konvergentstrahl-Geometrien anwendbar.
Projektionen werden erfasst an P getrennten Source-Positionen, die definiert sind durch einen Parameter entlang des Source-Orbit, θ_p = pΔθ, p = 0, ... P – 1, und in der Regel, aber nicht notwendigerweise, einen gleichmäßigen Abstand Δθ = (θ_max – θ_min)/P haben. Der Ausdruck (Pf)(θ_p,·) (für alle Werte von t →) wird als Divergentstrahl-Projektion an der Source-Position θ_p bezeichnet. Die Detektorebene wird gewöhnlich auf einem gleichmäßigen rechteckigen Gitter mit möglicherweise unterschiedlichen Intervallen T₁ und T₂ auf den Achsen t₁ und t₂ abgetastet.
Eine Sammlung von Divergentstrahl-Projektionen an unterschiedlichen Source-Positionen wird als Divergentstrahl-Sinogramm bezeichnet. Der gleiche Begriff wird auch für einen Satz von vorverarbeiteten Divergentstrahl-Projektionen oder jeden anderen, mathematisch gleichwertigen Datensatz verwendet. Die Sammlung von Projektionen in einem Sinogramm kann unvollständig sein und bis zu nur einer einzigen Projektion enthalten.
Im Folgenden werden einige Spezialfälle der allgemeinen Divergentstrahl-Geometrie beschrieben.
Helixförmige Kegelstrahl-Tomographie mit planarem Detektor
Siehe 3: In diesem Fall ist der Source-Orbit eine Helix oder Spirale um die Achse, definiert durch s →(θ) = [Dcos(θ),Dsin(θ),hθ/2π], wobei der Radius D konstant ist und h als Kippung der Helix bekannt ist.
Der Parameter θ gibt in diesem Fall den Source-Winkel in der (x,y)-Ebene relativ zur x-Achse an. Ohne die Allgemeingültigkeit zu verlieren wird angenommen, dass die Detektorebene die z-Achse enthält und mit der Source-Bewegung rotiert, so dass sie für jedes θ senkrecht zur Projektion des Source-Vektors s → auf die (x,y)-Ebene ist. Der Ursprung o_θ in der Detektorebene fällt in diesem Fall auf die z-Achse. Die t₂-Achse fällt mit der z-Achse zusammen, und die t₁-Achse ist parallel zur (x,y)-Ebene. (Weitere übliche Auswahlen, bei denen z.B. t ^₁ parallel zur Tangente des Orbits an der Position θ und t ^₂ orthonormal dazu ist, können durch Koordinatenrotation auf den hier beschriebenen Fall reduziert werden.)
In der Praxis kann die Detektorebene zwar anders positioniert sein, aber eine einfache lineare Koordinatentransformation wandelt Projektionen, die auf der tatsächlichen Detektorebene erfasst werden, in Projektionen um, die auf dem oben angenommenen virtuellen Detektor erfasst werden, bei gleichmäßigem Abstand entlang der Koordinaten (t₁,t₂).
Kegelstrahl-Bildgebung mit einem einzigen kreisförmigen Orbit und planarem Detektor
Siehe 4: In diesem Fall ist der Source-Orbit ein Spezialfall des helixförmigen Orbits mit der Neigung Null, h = 0, d.h. er liegt auf einem einzigen Kreis s →(θ) = [Dcos(θ),Dsin(θ),0] mit dem am Ursprung o beginnenden Radius D in der x,y-Ebene. Ohne die Allgemeingültigkeit zu verlieren wird angenommen, dass die Detektorebene die z-Achse enthält und senkrecht zu der Geraden SO zwischen Source und Mittelpunkt verläuft.
Man betrachte einen Punkt r → = [x,y,z]^T in dem Objekt. Die Position τ →(θ,r →) seiner Kegelstrahl-Projektion auf den Detektor ist für den Detektorwinkel θ gegeben durch
Diese expliziten Ausdrücke sind später für die Erläuterung der schnellen Algorithmen nützlich.
Fächerstrahl-Tomographie mit kreisförmigem Orbit und colinearem gleichmäßig beabstandeten Detektor
Siehe 5: In diesem Fall wird ein einziger Querschnitt des Objekts abgebildet und ein 2D-Bild rekonstruiert. Auch hier ist der Source-Orbit ein Kreis mit dem am Ursprung o zentrierten Radius D in der x,y-Ebene, wobei die Source-Position durch den Source-Winkel θ angegeben ist. Der Detektor ist jedoch auf eine einzige Gerade (die t-Achse) in der x,y-Ebene beschränkt. Ohne die Allgemeingültigkeit zu verlieren wird angenommen, dass sie durch den Rotationsmittelpunkt o der Source verläuft und senkrecht zu der Geraden SO zwischen Source und Mittelpunkt ist. Die Detektorelemente sind auf der t-Achse gleichmäßig beabstandet.
Die Position des 2D-Bildes sei r → = [x,y]^T und τ(θ,r →) die t-Position des Schnittpunktes von dem Source-Strahl, der durch den Punkt r → geht, mit der t-Achse (siehe 5). Die Fächerstrahl-Projektion (Pf)(θ,t) des Objekts f bei dem Source-Winkel θ und der Detektorposition t ist das Kurvenintegral entlang des Source-Strahls mit den Parametern (θ,t) und ist gegeben durch
Projektionen erfasst werden bei P verschiedenen Source-Winkeln θ_p = pΔθ, p = 0, ... P – 1 mit dem gleichmäßigen Abstand Δθ = θ_max/P. Der Ausdruck (Pf)(θ_p,·) (für alle Werte von t) ist eine Projektion bei dem Source-Winkel θ_p, und eine Sammlung von Projektionen für verschiedene p ist ein Sinogramm. Der maximale Source-Winkel θ_max für eine vollständige Abtastung ist gleich 2π, kann aber im Fall einer kurzen Abtastung oder einer Über-Abtastung andere Werte annehmen. Partielle Datensätze können Projektionen bei jedem beliebigen Unter-Satz von Winkeln umfassen.
In jedem Fall von Divergentstrahl-Geometrie kann eine gebogene anstelle einer planaren Detektoroberfläche verwendet werden, oft mit einem modifizierten Strahlenabtastschema. Bei der Kegelstrahl-Bildgebung mit helixförmigem oder einem einzigen kreisförmigen Orbit wird die Detektoroberflächen zum Beispiel oft so gewählt, dass sie auf der Oberfläche eines Zylinders liegt, dessen Mittelpunkt sich auf der Source befindet. Die Strahlenabtastung erfolgt gleichwinkelig im Azimuth und gleichmäßig beabstandet in der Höhe z. Dies ist in 6 für den einzelnen kreisförmigen Source-Orbit dargestellt. Dabei ist der gebogene Detektor derart gewählt, dass er die z-Achse enthält.
Die Beziehung zwischen den zylindrischen und planaren Detektorformen im Allgemeinen lässt sich durch Untersuchung der Fächerstrahl-Geometrie verstehen. Siehe 7: Bei der Fächerstrahl-Geometrie mit gleichwinkeliger Strahlenabtastung sind die Detektoren d in Winkeln mit dem Abstand Δt auf einem Bogen gleichmäßig beabstandet, dessen Mittelpunkt auf der Source S liegt. Somit ersetzt bei dieser Bildgebungsgeometrie der Fächerwinkel t die Detektorposition t, die bei der colinearen gleichmäßig beabstandeten Fächerstrahl-Geometrie auftritt. Nun bezeichne τ(θ,r →) die Fächerwinkelposition von dem Strahl, der durch den Punkt r → geht, siehe 7. Die Projektionen sind dann wiederum durch Formel (4) gegeben.
Es wird in Betracht gezogen, dass es andere Fälle der Fächerstrahl-Geometrie gibt, die ebenfalls Spezialfälle der Divergentstrahl-Geometrie sind, auf welche die Erfindung angewendet werden kann.
Gefilterte Divergentstrahl-Rückprojektionsrekonstruktion und Reprojektion
Eine native (direkte) Inversion von Divergentstrahldaten, ob angenähert oder genau, kann als gewichtete gefilterte Rückprojektion formuliert werden. Zunächst wird das Divergentstrahl-Sinogramm vorverarbeitet, zum Beispiel werden die Projektionen, die das Sinogramm umfassen, gewichtet und gefiltert, wodurch das modifizierte Divergentstrahl-Sinogramm erzeugt wird, das die Projektionen g(p,t →) umfasst, die Source-Positionen θ_p entsprechen. Das Filtern kann einfaches verschiebungsinvariantes Ramp-Filtern sein, wie im Fall von Fächerstrahl- oder annähernder Kegelstrahl-Rekonstruktion, oder viel stärker eingehendes verschiebungsvariierendes Filtern, wie im Fall der genauen Kegelstrahl-Rekonstruktion. Das Gewichten und Filtern kann sogar noch allgemeiner sein und zwei oder mehr Sätze von Projektionen hervorbringen, deren Beiträge summiert werden müssen, wie im Fall genauer oder quasi-genauer Verfahren zur Rekonstruktion aus helixförmigen Kegelstrahl-Daten mit einem teilweise abgetasteten langen Objekt. Aus Gründen der Kürze und Allgemeingültigkeit werden die Daten, die rückprojiziert werden sollen, gleich welche Vorverarbeitung verwendet wurde, um sie aus den ursprünglichen Divergentstrahl-Projektionsdaten zu erhalten, als "Divergentstrahl-Projektionsdaten" oder "Divergentstrahl-Sinogramm" bezeichnet. Für ein festes p werden sie "Divergentstrahl-Projektion" (an der Source-Position θ_p) genannt und mit der Funktion g(p,·) bezeichnet.
Das 3-D-Bild wird dann durch die herkömmliche Divergentstrahl-Rückprojektion bei diskreten θ rekonstruiert
wobei W(θ,r →) eine annähernde Gewichtungsfunktion ist. Diese diskrete Rückprojektionsformel approximiert den Integralausdruck für die Rückprojektion, der für alle θ gemessene Projektionen verwendet. In der Praxis wird aufgrund der Verwendung getrennter Detektoren g(p,t →) auch in t → abgetastet. Dies erfordert die Interpolation von g in der Variablen t → zur Implementierung der Rückprojektion, weil τ →(pΔθ,r →) gewöhnlich nicht einer verfügbaren Abtastposition entspricht. Die Rückprojektionsformel für den 2D-Fächerstrahl-Fall hat eine identische Form, ausgenommen dass τ skalar und r → zweidimensional ist.
Die Rechenkosten für eine 3D-Kegelstrahl-Rückprojektion für ein N×N×N-Bild betragen cN³P, wobei die Konstante c von Einzelheiten der Implementierung abhängt, wie der Komplexität der Interpolation. Dagegen betragen die Rechenkosten für das Gewichten und Ramp-Filtern nur O(PN²logN), wenn das Filtern verschiebungsinvariant erfolgt und die Konvolution unter Verwendung von FFT durchgeführt wird. Noch ausgefeiltere Formen des Filterns können oft zu gleichen Kosten durchgeführt werden. Daher beherrschen die Kosten für die Rückprojektion die Kosten der herkömmlichen Kegelstrahl-Rekonstruktion, welche die Kosten O(PN³) oder O(N⁴) hat, wenn P = O(N), wie es oft der Fall ist. Die Situation bei der 2D-Fächerstrahl-Rekonstruktion ist ähnlich, bei der die Komplexitäten der Filter- und Rückprojektionsschritte O(N²logN) bzw. O(N³) betragen.
Die native (direkte) Divergentstrahl-Reprojektion eines elektronisch gespeicherten Bildes f ist eine Implementierung der Gleichung (1). In der Praxis ist f nur in diskreter Form verfügbar, und nur Abtastungen in t von (Pf)(θ_p,t →) werden berechnet. Zum Beispiel kann f dargestellt werden durch die Reihe:
wobei die Koeffizienten β_ijk die diskrete Darstellung von f sind und ⌀_i,j,k(r →) lokalisierte Basisfunktionen sind, wie ein Tensorprodukt von Splines, oder die Standard-Pixel-Basisfunktion, d.h. die Indikatorfunktion des ijk-ten Voxels, wie z.B. in den U.S.-Patenten Nr. 6 263 096 und 6 287 257 für 2D und im U.S.-Patent Nr. 6 332 035 für 3D beschrieben. Die Projektion kann dann berechnet werden durch
wobei ∅ ~(θ_p,t →) = (P⌀_ijk)(θ_p,tt →) die Divergentstrahl-Projektion der Basisfunktion ⌀_ijk(r →) ist, wie durch Gleichung (1) definiert. Diese Projektionen der Basisfunktionen können in der Regel unter Verwendung eines analytischen Ausdrucks berechnet oder in einer Aufruftabelle gespeichert werden. Wie bei der Rückprojektion betragen die Kosten für die native Reprojektion O(N⁴) bei 3D oder O(N³) bei 2D.
Im Fall eines gekrümmten Detektors, wobei t₁,t₂ geeignete Koordinaten auf dem Detektor sind, soll τ →(θ,r →) jetzt die (t₁,t₂)-Position auf der Detektorebene von dem Source-Strahl bezeichnen, der durch den Punkt r → geht. Die Projektionen, die wiederum durch Formel (1) gegeben sind, können unter Verwendung von Gleichung (7) berechnet werden, wobei ⌀ ~_i,j,k(θ_p,t →) = (P⌀_i,j,k)(θ_p,t →) die Divergentstrahl-Projektion der Basisfunktion ⌀_ijk(r →) auf dem gekrümmten Detektor angibt. Ebenso kann die Rekonstruktion wiederum als eine gewichtete Rückprojektion geeignet modifizierter Projektionen g(p,t →) ausgedrückt werden, wobei t₁,t₂ geeignete Koordinaten auf dem gekrümmten Detektor sind. Die Rückprojektion wird durch die gleiche Formel (5) wie oben angegeben, aber mit einer anderen Gewichtungsfunktion W(θ,r →). Zum Beispiel ersetzt bei der gleichwinkligen 2D-Fächerstrahl-Geometrie der Fächerwinkel t die Versetzung t auf dem geraden Detektor bei der gleichmäßig beabstandeten colinearen Geometrie. Die aus der Diskretheit von t → hervorgehenden Überlegungen sind ebenso wie die Überlegungen zur Berechnungskomplexität die gleichen wie bei den Algorithmen für planare Detektoren.
Wie bereits erläutert ist, ob in 2D oder 3D, der Rückprojektionsschritt (und die Reprojektion, wenn verwendet) gewöhnlich hinsichtlich der Berechnung am teuersten und daher das Ziel der Erfindung. Es wird ein neues allgemeines Verfahren zur schnellen Rückprojektion beschrieben, das auf alle Divergentstrahl-Geometrien anwendbar ist. Zur Vereinfachung der Erläuterung werden jedoch zunächst die schnelle Rückprojektion und Reprojektion für die 2D-Fächerstrahl-Geometrien beschrieben.
Schnelle native Rückprojektions- und Reprojektionsalgorithmen für die colineare gleichmäßig beabstandete Fächerstrahl-Geometrie
Siehe 8: Man betrachte den Rückprojektionsschritt für ein Unter-Bild f'(r →) von f. K_M[δ →] sei ein Bildabschneideoperator, der definiert ist durch
wobei ∥r →∥ = max{x,y}. Somit ist f' = K_M[δ →]f ein Unter-Bild von f der Größe M×M mit dem Mittelpunkt bei O' = δ → ∊ R². Die Rückprojektion in das globale Koordinatensystem auf das Unter-Bild f' (das sich an seiner korrekten Position im globalen Koordinatensystem befindet) unter Verwendung von Q Projektionen bei den Source-Winkeln pΔθ,p = 0...Q – 1 ist gegeben durch
Der Ausdruck B_M,Q[δ →] bezeichnet den zugehörigen Rückprojektionsoperator, also genauer f = B_N,P[0 →]g.
Aufgrund der Stellung der Rückprojektion trägt nur ein Teil der Projektion g(p,·) zur Rückprojektion auf f' bei. Dieser Teil wird als K ^_M[δ →]g bezeichnet, wobei K ^_M[δ →] der Operator ist, der g(p,t) in der Variablen t für jedes θp auf den Träger Ξ = AB der Projektion des Trägers von dem Unter-Bild f' = K_M[δ →]f abschneidet (siehe 8). Die Abschneideintervalle, die K ^_M[δ →] festlegen, können für alle Source-Winkel θp sowie Unter-Bild-Größen M und Stellungen δ → von Interesse zuvor berechnet werden. Ersatzweise kann eine einfachere konservative Annäherung verwendet werden zu leicht höheren Rechenkosten in dem hierarchischen Algorithmus aufgrund einer Überschätzung des benötigten Trägers von Projektionen, wozu Verschiebung und Dezimierung nötig sind.
Aus der Lokalisierung der Rückprojektion folgt, wenn f = B_N,P[0 →]g, dass dann f' = KM[δ →]f = BM,P[δ →]K ^M[δ →]g, (10)d.h. die Rückprojektion auf f' kann durch eine Rückprojektion der Größe (M,P) der geeignet abgeschnittenen Projektionen erhalten werden.
Man betrachte jetzt eine Unterteilung des Bildes f in J = (N/M)² nichtüberlappende Unter-Bilder, jeweils mit der Größe M×M.
wobei die Vektoren δ →_j = 1, ... J die Mittelpunkte der Unter-Bilder im globalen Koordinatensystem sind, wie in 9 gezeigt. Durch Anwenden der Gleichung (10) erhalten wir die folgende genaue Zerlegung für die Rückprojektion in Rückprojektionen auf die Unter-Bilder
Diese Zerlegung ist in 10 für N/M = 2 dargestellt. Das Divergentstrahl-Sinogramm₃ 22, das in Gleichung (12) auch mit g bezeichnet ist, wird in eine Mehrzahl an Unter-Sinogrammen geteilt, die den gerade beschriebenen Unter-Bildern entsprechen. Die Unter-Sinogramme können jede gewünschte Größe haben und nur ein Pixel oder ein Voxel ausmachen. Die durch die Abschneideoperatoren 1010, 1012, 1014 und 1016 erzeugten Unter-Sinogramme werden in 1018, 1020, 1022 bzw. 1024 einzeln in das globale Koordinatensystem rückprojiziert. Die rückprojizierten Unter-Sinogramme erzeugen die Unter-Bilder f₁, f₂, f₃ und f_4, die bei 1026 aggregiert werden, wodurch das elektronische Bild₁ 26 in 1 erzeugt wird. Dieses Verfahren erzeugt eine genaue Zerlegung der Fächerstrahl-Rückprojektion. Für N/M = 2 erfolgt die Zerlegung in vier Unter-Sinogramme und Unter-Bilder, wie 10 zeigt.
Siehe wiederum 10: P unter der Linie 22 gibt die Anzahl der Projektionen in einem Sinogramm oder Unter-Sinogramm an. Diese Zerlegung liefert selbst keine Beschleunigung verglichen mit der Ein-Schritt-Rückprojektion B_N,P. Während nämlich die Rechenkosten cM²P für die Rückprojektion eines einzelnen Unter-Bildes B_M,P(N/M)² = J-mal kleiner sind als bei der Ein-Schritt-B_N,P für das ganze Bild, bleiben die Gesamtkosten für die erforderlichen J solcher Unter-Bild-Rückprojektionen gleich (möglicherweise mit einem kleinen zusätzlichen Überhang für die Buchführung). Die genaue Zerlegung kann immer noch dazu verwendet werden, dass man eine Beschleunigung durch Untertteilen der Rückprojektion zwischen J parallelen Prozessoren bereitgestellt oder eine wiederholte Verwendung eines einzigen Prozessors zugelassen wird, was J-mal weniger Speicher benötigt.
Der schnelle Rückprojektionsalgorithmus verwendet die Zerlegungsidee mit der folgenden zusätzlichen Eigenschaft. Für eine feste Bildauflösung (Bandbreite) ist die Anzahl der für die Rekonstruktion des Bildes benötigten Projektionen proportional zur Größe des Bildes. Wenn also P Projektionen zur Rekonstruktion eines N×N-Bildes benötigt werden, dann reichen P' = (M/N)P Projektion zur Rekonstruktion eines M×M-Bildes mit der gleichen Auflösung.
Daher basiert der schnelle Rückprojektionsalgorithmus auf der Rekonstruktion eines Unter-Bildes f' = K_M[δ →]f aus einer kleineren Anzahl an Projektionen. Dieses Reduktionsverfahren wird durchgeführt mit einem Operator O[L,M,δ →], der selbst aus mehreren Operatoren besteht, die im Folgenden definiert werden.
Projektionenverschiebeoperatoren M ^_-[δ →] und M ^₊[δ →], die jede Projektion um das Ausmaß ±r(θ,δ →) entlang der t-Achse verschieben, sind definiert als M ^_[δ →]g(θ,t) = g[θ,t – τ(θ,δ →)] M ^+[δ →]g(θ,t) = g[θ,t + τ(θ,δ →)] (13)
Ein L-facher Winkel-Dezimierungsoperator D_↓L verringert durch Winkelfiltern und anschließende Winkel-Unterabtastung die Anzahl der Ansichten in einem Satz von P Projektionen auf P/L Projektionen. Das Filtern kann als einfache Durchschnittsbildung von Projektionen bei benachbarten Winkeln erfolgen oder ausgefeilter sein und wird angegeben durch g ~(p,t) = hD(p,t)*g(p,t) (14),wobei * die Konvolution angibt und h_D der Filterkern ist. Das Filtern kann nur in p (d.h. bei diskreten θ) für jedes feste t erfolgen oder kombiniertes Filtern in θ und t sein, die trennbar oder nicht trennbar sein können. Das Winkel-Unterabtasten bewahrt eine aus jeweils L gefilterten Projektionen g ~(p,t) auf.
Der Operator O[L,M,δ →], eine Zusammensetzung aus den Abschneide-, Verschiebungs- und Winkeldezimierungsschritten, ist dann definiert durch g'(q,t) = O[L,M,δ →]g(q,t) = M ^_[δ →]W2D↓LW1M ^+[δ →]K ^M[δ →]g(q,t), q = 0, ..., P' – 1, (15) wobei P' = P/L. Der Operator O[L,M,δ →] beinhaltet eine multiplikative Gewichtung mit den Gewichtungsfunktionen W_k(θ,t),k = 1,2 vor und nach der Dezimierung. Eine mögliche Auswahl könnte W₁(θ,t) = W[θ,τ(θ,δ →)] und W₂(θ,t) = 1/W₁(θ,t) sein. Allgemeiner gesagt, kann das Gewichten derart gewählt werden, dass die Genauigkeit des Algorithmus optimiert wird, oder völlig weggelassen werden.
11 zeigt ausführlicher die Art und Weise, wie die Unter-Sinogramme verarbeitet werden können. Zum Beispiel kann das Abschneiden K ^_M[δ →] zuerst im Schritt 1110 erfolgen. Darauf folgen der Verschiebeoperator M ^₊[δ →] im Schritt 1112, ein Gewichtungsschritt 1114, darauf ein Dezimierungsschritt 1116, ein weiterer Gewichtungsschritt 1118 und ein zweiter Verschiebeschritt 1120. Dieses Verfahren beschreibt den Sinogramm-Zerlegungsoperator O[L,M,δ →].
Mit diesen Definitionen wird die genaue Formel (10) für die Rückprojektion auf das Unter-Bild f' durch folgende Annäherung ersetzt f' = KM[δ →]f = BM,P/L[δ →]g' (16) = BM,P/L[δ →]O[L,M,δ →]g, (17)wobei B_M,P/L eine Rückprojektion in das globale Koordinatensystem auf ein M×M-Unter-Bild unter Verwendung von P/L Projektionen ist, wie durch Gleichung (9) definiert. Dies führt zu einer annähernden Zerlegung des Rückprojektionsschrittes für ein unterteiltes Bild, die analog zu Gleichung (12) ist
Diese Zerlegung ist in 12 dargestellt. Gezeigt ist eine annähernde Zerlegung für eine Fächerstrahl-Rückprojektion mit einer kleineren Anzahl an Projektionen, wobei L = N/M = 2. Zur Ausführung der annähernden Zerlegung wird der schnelle Rückprojektionsprozessor 24 derart programmiert, dass er das Divergentstrahl-Sinogramm₃, das in 12 mit g bezeichnet ist und P Projektionen umfasst, in vier Unter-Sinogramme unterteilt, die durch 1210, 1212, 1214 und 1216 erzeugt werden und mit g₁, g₂, g₃ und g₄ bezeichnet sind und jeweils P / 2 Projektionen umfassen. Die Unter-Sinogramme g₁, g₂, g₃ und g₄ werden durch 1218, 1220, 1222 und 1224 in das globale Koordinatensystem rück projiziert, wobei die Unter-Bilder f₁, f₂, f₃ und f₄ an den korrekten Stellen im globalen Koordinatensystem erzeugt werden. Die Unter-Bilder werden bei 1226 aggregiert, so dass das elektronische Bild 26 erzeugt wird, das in 12 mit f bezeichnet ist.
Siehe wiederum 12: Die Anzahl der Schritte für die Rückprojektion B_M,P/L auf jedes der J = (N/M)² Unter-Bilder beträgt cM²P/L, was insgesamt cN²P/L für alle Rückprojektionen ergibt. Dies ist L-mal weniger als bei der Ein-Schritt-Rückprojektion B_N,P. Auf Basis der früheren Beobachtung hinsichtlich der Größenveränderung der Anzahl an Projektionen mit der Bildgröße kann L bis zu einer Höhe von L = (N/M) gewählt werden.
Werden die Rechenkosten von O[L,M,δ →] hinzugezählt, kann eine optimale Anzahl an Unter-Bildern für eine Ein-Ebenen-Zerlegung bestimmt werden, wie im U.S.-Patent 6 263 096 dargelegt.
Bei einer bevorzugten Ausführungsform des schnellen Rückprojektionsalgorithmus werden die Zerlegungen der Gleichungen (12) und (18) rekursiv angewendet, wobei auf jeder Ebene die Unter-Bilder gemäß Gleichung (11) weiter in kleinere Unter-Bilder zerlegt werden. Ein Beispiel für eine dyadische (N/M = 2) Zwei-Ebenen-Unterteilung des Bildes ist in 13 dargestellt. Ebenfalls dargestellt sind die Mittelpunktpositionen der Unter-Bilder auf den verschiedenen Ebenen. Man beachte, dass sie sich alle auf den Bildursprung O der nullten Ebene, d.h. auf das globale Koordinatensystem, beziehen. Die entsprechenden Beispiele für eine annähernde Zwei-Ebenen-Zerlegung und eine Zwei-Ebenen-Zerlegung, die aus einer Kaskade von einer genauen Zerlegung und einer annähernden Zerlegung besteht, sind in den 14 bzw. 15 dargestellt. Siehe 14: das Divergentstrahl-Unter-Sinogramm₃ 22, das P Projektionen umfasst, wird im schnellen Rückprojektionsprozessor 24 verarbeitet, indem das Sinogramm in die Unter-Sinogramme 1410, 1412, 1414 und 1416 unterteilt wird, die jeweils P/2 Projektionen umfassen. Diese annähernde Unterteilung erfolgt durch den Operator für die annähernde Zerlegung, der von dem Fließdiagramm der 11 beschrieben wird. Weitere annähernde Zerlegungen des Unter-Sinogramms 1410 werden bei 1418, 1420, 1422 und 1424 erzeugt, wobei wiederum der von dem Fließdiagramm der 11 beschriebene Operator verwendet wird, diesmal jedoch mit anderen Parametern, die den kleineren Unter-Bildern in 13 entsprechen. Jedes dieser Unter-Sinogramme umfasst P/4 Projektionen und wird bei 1426, 1428, 1430 und 1432 in das globale Koordinatensystem zurückprojiziert, so dass eine Mehrzahl an Unter-Bildern f_1,1, f_1,2, f_1,3, f_1,4 an den richtigen Stellen in dem globalen Koordinatensystem erzeugt wird, wie 13 zeigt. Die Unter-Bilder werden bei 1434 aggregiert unter Bildung eines Unter-Bildes f₁ – ebenfalls an der richtigen Stelle im globalen Koordinatensystem.
Das durch 1416 dargestellte Unter-Sinogramm wird bei 1436, 1438, 1440 und 1442 ebenso unterteilt. Diese Unter-Sinogramme werden bei 1444, 1446, 1448 bzw. 1450 rückprojiziert und bei 1452 aggregiert unter Bildung des Unter-Bildes f₄. Die gestrichelten Kästen 1454 und 1456 stehen für die entsprechenden Ebene-2-Zerlegungen von f₂ bzw. f₃. In 14 sind alle Zerlegungen annähernde Zerlegungen einer Fächerstrahl-Projektion für L = (N/M) = 2.
15 stellt eine gemischte (eine genaue, eine annähernde) Zwei-Ebenen-Zerlegung einer Fächerstrahl-Projektion für L = (N/M) = 2 dar. In der Figur wird das Sinogramm 22, das P Projektionen umfasst, in die genauen Unter-Sinogramme 1510, 1512, 1514 und 1516 zerlegt, die ebenfalls jeweils P Projektionen umfassen. Die als g₁, g₂, g₃ bzw. g₄ bezeichneten genauen Unter-Sinogramme werden (mit annähernden Parametern) unter Verwendung des in 11 dargestellten annähernden Zerlegungsoperators in annähernde Sub-Sinogramme zerlegt, die jeweils P/2 Projektionen umfassen. Zum Beispiel wird g₁ bei 1518, 1520, 1522 und 1524 in die annähernden Unter-Sinogramme g'_1,1, g'_1,2, g'_1,3, g'_1,4 zerlegt. Diese Unter-Sinogramme werden bei 1526, 1528, 1530 und 1532 in das globale Koordinatensystem zurückprojiziert, so dass die Unter-Bilder f_1,1, f_1,2, f_1,3, f_1,4 an den richtigen Stellen im globalen Koordinatensystem erzeugt werden. Diese Unter-Bilder werden bei 1534 unter Bildung eines Unter-Bildes f₁ aggregiert.
Das in 15 als g₄ bezeichnete Unter-Sinogramm, das bei 1516 erzeugt wird, wird bei 1536, 1538, 1540 und 1542 ebenfalls weiter zerlegt unter Bildung der Unter-Sinogramme g'_4,1, g'_4,2, g'_4,3, g'_4,4. Rückprojektionen bei 1544, 1546, 1548 und 1550 erzeugen die Unter-Bilder f_4,1, f_4,2, f_4,3 und f_4,4, die bei 1552 aggregiert werden unter Bildung eines Unter-Bildes f₄. Die gestrichelten Kästen 1554 und 1556 stehen für die entsprechenden gemischten Zwei-Ebenen-Zerlegungen von g₂ bzw. g₃. Diese Verfahren erzeugen die Unter-Bilder f₂ und f₃, und die Unter-Bilder f₁, f₂, f₃ und f₄ werden bei 1558 aggregiert, wodurch das Bild 26 erzeugt wird, das in 15 mit f bezeichnet ist.
Die Rekursion kann fortgesetzt werden, bis die Unter-Bilder eine gewünschte minimale Größe M_min haben. Dann werden die Rückprojektionen B_Pmin,Mmin in das globale Koordinatensystem unter Verwendung von Gleichung (5) durchgeführt. Das optimale M_min hängt gewöhnlich von der Implementierung ab. Es kann bis zu M_min = 1 gewählt werden, so dass die kleinsten Unter-Bilder eine Größe von 1 Pixel oder Voxel haben. Unter der Annahme, dass N eine zweite Potenz ist und N/M = 2 auf allen Ebenen der Zerlegung, gibt es log₂N solche Ebenen. Nimmt man weiter an: L = (N/M) = 2 und P = kN, wobei k eine kleine ganze Zahl ist, beinhaltet die letzte Ebene N² Ein-Pixel-Rückprojektionen B_P/N,1 mit den Gesamtkosten cPN. Unter der Annahme, dass Interpolatoren mit fester Länge für die Dezimierungs- und Verschiebeschritte bei der Implementierung von O[L,M,δ →] verwendet werden, lässt sich zeigen, dass die Rechenkosten für jede Stufe c₁NP betragen, wobei c₁ eine Konstante ist. Die Gesamtkosten für die log₂N Ebenen betragen daher O(PN logN). Weil P = O(N), werden die Gesamtkosten des rekursiven Zerlegungsalgorithmus zu O(N²logN).
In der Praxis werden aufgrund der Verwendung diskreter Detektoren Projektionsdaten in t mit dem Intervall Δt abgetastet. Deshalb werden die Projektionenverschiebeoperatoren M ^_[δ] und M ^+[δ] in diskreter Form implementiert, erfordern aber Verschiebungen um nicht-ganzzahlige Vielfache von Δt. Weil die Projektionen vor der Rückprojektion gefiltert worden sind, sind sie bandbeschränkt und können unter Verwendung einer Kombination von Interpolation und erneuter Abtastung verschoben werden, wie z.B. im US-Patent 6 287 257 beschrieben. Die Verschiebung um nicht-ganzzahlige Vielfache von Δt können aufgeteilt werden in eine Verschiebung um ganzzahlige Vielfache von Δt in einer Kaskade mit einer Verschiebung um einen Bruchteil von Δt. Die ganzzahligen Verschiebungen entsprechen einer einfachen Neuindexierung, wohingegen die Bruchteil-Verschiebungen mittels Interpolation in der t-Variablen durchgeführt werden können. Diese Interpolation kann zudem durch radiales Filtern ausgeführt werden, was der Konvolution in der radialen Dimension mit dem Filter h_D(p,t) in dem allgemeingültigen Dezimierungs-Filterschritt in Gleichung (14) entspricht.
Es ist entscheidend, dass bei dem hier vorgeschlagenen spezifischen Projektionenverschiebeschema das Abtastintervall der Projektionen gleichmäßig und konstant bleibt, so dass Interpolation und erneute Abtastung als diskret-zu-diskret-Kartierung durchgeführt und als einfacher Digitalfilter effizient implementiert werden können. Man kann den Fehler bei einer solchen Bruchteil-von-t-Verschiebung mit der Länge des Interpolators exponentiell abklingen lassen, so dass kurze Interpolatoren – sogar einfache lineare Interpolation – in der Regel genügen. Die Genauigkeit der t-Interpolation nimmt zudem durch Überabtastung (Oversampling) der Projektionen in t zu, was mittels digitaler Interpolation (höhere Abtastung (Upsampling) und anschließendes Filtern) der Projektionen effizient durchgeführt werden kann, wenn gewünscht.
Die Genauigkeit der Rückprojektion, die unter Verwendung der annähernden Zerlegung erhalten wird, hängt von den verschiedenen Parametern in dem Verfahren ab, einschließlich P, N, M, L, der Auflösung des Bildes und der Bandbreite der Projektionen, dem maximalen Fächerwinkel von Source-Strahlen, die das Objekt schneiden, der Detektorreaktion und den bei dem Algorithmus verwendeten Interpolatoren. Diese Parameter lassen sich in den Divergentstrahl-Algorithmen annähernd optimieren. Ein wichtiger Parameter, der die Genauigkeit beeinflusst, ist das Verhältnis (P/L)/M zwischen der Anzahl der Projektionen, die in den Rückprojektionen verwendet wird, und der Größe der Unter-Bilder. Durch Erhöhen dieses Verhältnisses, auch als Winkel-Überabtastung(-Oversampling) bezeichnet, steigt die Genauigkeit der annähernden Zerlegung.
Es gibt verschiedenerlei Weisen zur Erhöhung des Winkel-Oversampling, u.a. die folgenden: (i) digitale Winkel-Interpolation der Projektionen g zur Erhöhung von P vor Beginn der Zerlegung; (ii) Wählen von L < (N/M) auf verschiedenen Stufen der annähernden Zerlegung; und (iii) Verwenden genauer Zerlegungsschritte, wie z.B. in 15 dargestellt. Zum Verständnis des Verfahrens (iii) erinnere man sich daran, dass die genaue Zerlegung in Gleichung (12) die Unter-Bilder verkleinert, aber nicht die Anzahl der Projektionen in einem Unter-Sinogramm. Dadurch führt sie zu einer Zunahme des Winkel-Oversampling um den Faktor (N/M) und zu größerer Genauigkeit.
Zwar können Besonderheiten der Implementierung einem der Verfahren (ii) oder (iii) einen kleinen Rechenvorteil gegenüber dem anderen geben, aber das Winkel-Oversampling-Verfahren (iii) entspricht dem Verfahren (ii) aufgrund der Identität von K ^_M/2[δ ⇀_ij]K ^_M[δ ⇀_i] = K ^_M/2[δ ⇀_ij], die immer dann gilt, wenn δ ⇀_ij der Mittelpunkt eines Unter-Bildes von dem bei δ ⇀_i zentrierten Unter-Bild ist. Deshalb folgt daraus, dass die Zweistufen-Zerlegung mit dyadischer Unterteilung, einer genauen und einer annähernden, äquivalent zu einer einzelnen annähernden Zerlegung in 16 Unter-Bilder mit jeweils einer Größe von N/4 mit der Dezimierung L = 2 anstelle des maximal möglichen L = 4 ist. Somit liefern beide Verfahren den gleichen Faktor 2 für das Winkel-Oversampling. Wählt man dagegen L = 3, kann die annähernde Ein-Ebenen-Zerlegung in 16 Unter-Bilder einen Faktor von 4/3 für das Winkel-Oversampling liefern. Dadurch erhält man zusätzliche Flexibilität bei der Wahl des Arbeitspunktes. Im Allgemeinen führen Verfahren zur Erhöhung des Winkel-Oversampling zu einer Zunahme der Rechenkosten für den Algorithmus um einen konstanten Faktor. Damit kann die Balance zwischen Berechnung und Genauigkeit gesteuert werden.
Eine bevorzugte Ausführungsform des Algorithmus beinhaltet eine Reihe genauer Unterteilungen und eine Reihe annähernder Unterteilungen oder eine annähernde Unterteilung mit Dezimierungsfaktoren, die kleiner als ihre Unterteilungsfaktoren sind. Die Anzahl der genauen Unterteilungen und/oder die Dezimierungsfaktoren steuern den Over sampling-Faktor. Sie werden derart gewählt, dass die gewünschte Balance zwischen Berechnung und Genauigkeit erhalten wird.
Die einfachste Ausführungsform verwendet eine Unterteilung mit einem konstanten Faktor N/M = 2 auf jeder Ebene, wobei L = 2 in den annähernden Unterteilungsebenen. Es gibt natürlich andere Wahlmöglichkeiten.
Wie dargelegt, erfordert der Algorithmus, dass P durch das Produkt der Dezimierungsfaktoren L auf allen Ebenen teilbar ist. Kann dies durch die Auswahl der Dezimierungsfaktoren nicht leicht oder effizient erreicht werden, kann P mittels Winkelinterpolation und erneuter Abtastung der Projektionen modifiziert werden.
Zur weiteren Verdeutlichung der rekursiven Struktur der bevorzugten Ausführungsform des Algorithmus ist in 16 ein Beispiel für einen als rekursive Funktion geschriebenen Pseudo-Code für den Algorithmus dargestellt, wobei M/N = L = 2. Der Pseudo-Code versteht sich zwar von selbst, wird aber kurz erläutert.
Die Funktion FAN FHBP ist eine rekursive Funktion, die in dem ersten Top-Level-Invocation-Call an den vollständigen Sinogramm G arbeitet. Bei anschließenden Aufrufen von sich selbst zerlegt diese Funktion das Bild in Unter-Bilder und führt Sinogramm-Zerlegungen durch, bei eine minimale Bildgröße N_min erreicht ist. Q dieser Zerlegungen sind genau, der Rest annähernd. Ist die Bildgröße N_min erreicht, erfolgt eine exakte Rückprojektion in das globale Koordinatensystem an dem Unter-Sinogramm, das dem Unter-Bild entspricht. Bei größeren Unter-Bildern als N_min werden die Projektionen abgeschnitten, so dass sie dem Unter-Bild entsprechen. Dies bildet die genaue Sinogramm-Unterteilung. Sind die genauen Zerlegungen beendet (wenn Q < 1), werden die weiteren Schritte für die annähernden Zerlegungen durchgeführt. In 16 werden die Projektionen verschoben, gewichtet-dezimiert und erneut verschoben. Die Unter-Bilder werden bei FAN_FHBP = f = Σ4i=1 fi aggregiert.
Auch die Unter-Bild-Mittelpunkte werden rekursiv berechnet, wobei der Mittelpunkt δ ⇀_i des Unter-Bildes i von einem bei δ ⇀ zentrierten größeren Unter-Bild ausgedrückt wird als δ ⇀_i = δ ⇀ + ξ ⇀_i. Dabei ist ξ ⇀_i der annähernde Positionsvektor des Unter-Bildes f_i relativ zu δ ⇀. Der Mittelpunkt des Bildes mit der vollen Größe liegt per definitionem am Ursprung des globalen Koordinatensystems. Deshalb werden alle folgenden Unter-Bild-Mittelpunkte im globalen Koordinatensystem ausgedrückt.
Bei einer bevorzugten Ausführungsform des schnelle Rückprojektionsalgorithmus werden die rekursiven Zerlegungen durch die folgende Beobachtung vereinfacht. Die Schritte Verschieben der Projektionen in t, Abschneiden der Projektionen und Gewich tung der Projektionen kommutieren (mit geeigneten Koordinatentransformationen). Daher können die Schritte M ^_[δ ⇀_i]W_2,i einer Stufe mit W_1,ijM ^₊[δ ⇀_ij] der folgenden Stufe zu einem Verschiebe- und Gewichtungsschritt kombiniert werden. Ebenso können die Verschiebe- und Gewichtungsschritte M ^_[δ ⇀_i]W_2,i kurz vor der endgültigen Unter-Bild-Rückprojektion B_Mmin,Pmin mit Gewichtungs- und Interpolationsschritten bei der Rückprojektion selbst kombiniert werden. Bei einer rekursiven Implementierung erfordert daher O[L,M,δ ⇀_j] nur einen Projektionsverschiebeschritt und (höchstens) eine Gewichtung pro Stufe.
Für den Fachmann auf dem Gebiet der Gestaltung und Implementierung von Signalverarbeitungsalgorithmen ist ersichtlich, dass die verschiedenen Schritte, welche die erfindungsgemäßen Verfahren definieren, in verschiedenen anderen Reihenfolgen angeordnet werden können, als sie in den hier gezeigten Figuren beschrieben sind, oder parallel ausgeführt werden können, ohne dass sich die Funktion oder die Prinzipien des Verfahrens ändern. Einige dieser Umordnungen können unter dem Gesichtspunkt der Implementierung, des Betriebs oder der Kosten vorteilhaft sein. Zum Beispiel kann der Schritt der Dezimierung der Anzahl an Projektionen für ein Unter-Bild beginnen, sobald zwei oder mehr Projektionen für dieses Unter-Bild abgeschnitten oder verschoben worden sind, ohne dass man auf das Abschneiden und Verschieben aller P Projektionen wartet. Diese Umordnungen können dazu verwendet werden, den Rekonstruktionsalgorithmus weiterzuleiten und mit der Datenerfassung zu überlappen.
Die erfindungsgemäßen Verfahren können zur stetigen Aktualisierung der tomographischen Bildgebung angewendet werden, weil Projektionen in zeitlicher Aufeinanderfolge erhalten werden. Dies kann erzielt werden durch Verwendung des vorgeschlagenen Verfahrens zur Rückprojektion eines Sinogramms, das aus der Differenz zwischen den neuesten erhaltenen Projektionen und einem älteren Satz an Projektionen besteht. Das Rückprojektionsergebnis wird dann zu dem gegenwärtigen Bild hinzuaddiert, wodurch das aktualisierte Bild erhalten wird. Solche Aktualisierungen können sich in aufeinander folgender Weise fortsetzen.
Schnelle Reprojektion
Rückprojektion und Reprojektion sind Tandemschritte oder Adjunkte voneinander. Daher sind die schnellen nativen Fächerstrahl-Algorithmen für ihre Implementierung eng miteinander verbunden. Daher wird der schnelle native Fächerstrahl-Algorithmus etwas weniger ausführlich beschrieben, wobei selbstverständlich sein sollte, dass alle Erläute rungen zum schnellen Rückprojektionsalgorithmus eine Entsprechung beim Reprojektionsalgorithmus haben.
Siehe 8: Man betrachte das M×M-Unter-Bild f' = K_M[δ ⇀]f von f. Seine Fächerstrahl-Projektionen bei den Source-Winkeln θ_q = qΔθ,q = 0 ... Q – 1 sind g(q,t) = PM,Q[δ →]f'(q,t) = (PKM[δ →]f)(q,t). (19)
Wir bezeichnen mit P_M,Q[δ →] den zugehörigen Reprojektionsoperator in das globale Koordinatensystem für das M×M-Unter-Bild mit dem Mittelpunkt an der Stelle δ ⇀. Genauer gesagt, sind (Pf) = P_N,P[0 ⇀]f die Fächerstrahl-Projektionen des gesamten Bildes. Die Rechenkosten von P_M,Q betragen cM²Q bei einer Konstanten c.
Man betrachte als nächstes die in 9 dargestellte Bildunterteilung der Gleichung (11). Aufgrund der Linearität des Projektionsoperators kann die folgende genaue Zerlegung der Reprojektion in die Summe von J Reprojektionen kleinerer M×M-Bilder für die gleiche Anzahl von P Source-Winkeln wie folgt erhalten werden:
Diese genaue Zerlegung für die Reprojektion ist in 17 für N/M = 2 dargestellt. Siehe 17: Das elektronische Bilde 30, in der Figur auch als f bezeichnet, wird mit den Abschneideoperatoren 1710, 1712, 1714 und 1716 in die Unter-Bilder f₁, f₂, f₃ bzw. f₄ unterteilt. Die Unter-Bilder werden durch 1718, 1720, 1722 und 1724 in das globale Koordinatensystem reprojiziert, wodurch die Unter-Sinogramme g₁, g₂, g₃ und g₄ erhalten werden. Die Unter-Sinogramme werden bei 1726 unter Bildung des Sinogramms 34 aggregiert, das in 17 auch als g bezeichnet ist.
Aus den gleichen Gründen wie bei der genauen Zerlegung der Gleichung (12) für die Rückprojektion liefert diese Zerlegung für die Reprojektion keine Beschleunigung verglichen mit der Ein-Schritt-Rückprojektion (Pf) = P_N,P[0 ⇀]f. Sie besitzt aber insbesondere in Kombination mit der annähernden Zerlegung weitere Anwendungsmöglichkeiten.
Der schnelle Reprojektionsalgorithmus verwendet die Zerlegungsidee mit der folgenden, bereits erwähnten zusätzlichen Eigenschaft: Bei einer festen Bildauflösung (Bandbreite) ist die Anzahl der Fächerstrahl-Projektionen, die das Bild charakterisieren, oder ersatzweise die Bandbreite der Projektionen in die θ-Variable proportional zur Bildgröße. Charakterisieren also P Projektionen ein N×N-Bild f, dann charakterisieren P' = P/(N/M) Projektionen ein Unter-Bild f', und die restlichen Projektionen können aus diesem kleineren Satz erhalten werden.
Daher basiert der schnelle Reprojektionsalgorithmus auf der Berechnung eines Satzes von P' Projektionen eines M×M-Unter-Bildes f' und der Verwendung eines "Projektionenvermehrungsverfahrens" zur Gewinnung des größeren Satzes an P Projektionen. Dieses Vermehrungsverfahren erfolgt durch einen Operator ν[L,M,δ →], der seinerseits aus mehreren Operatoren besteht. Die Projektionenverschiebeoperatoren M^_[δ] und M ^+[δ] wurden bereits in Gleichung (13) definiert. Der L-fache Winkel-Interpolationsoperator I_↑L erhöht durch L-fach erhöhte Winkelabtastung und anschließendes Winkelfiltern die Anzahl der Ansichten in einem Satz von Q Projektionen auf QL Projektionen. Die erhöhte Winkelabtastung führt L-1 nullwertige Projektionen zwischen jede zweite Projektion ein. Bezeichnet man mit g(p,t) die Projektionen mit erhöhter Abtastung, erfolgt das Filtern wie in Gleichung (14) und der zugehörigen Erläuterung beschrieben. Der so erhaltene Winkel-Interpolationsoperator I_↑L ist das Adjunkt des Winkel-Dezimierungsoperators D_↓L, der im schnellen Rückprojektionsalgorithmus erscheint. Der Projektionenvermehrungsoperator ν[L,M,δ ⇀], eine Zusammensetzung aus Verschiebe- und Winkelinterpolationsschritten, ist dann definiert durch
wobei P' = QL. Der Operator ν[L,M,δ →] enthält somit eine multiplikative Gewichtung mit den Gewichtungsfunktionen W_k(θ,t),k = 1,2 vor und nach der Interpolation mit ähnlichen Funktions- und Auswahlkriterien wie bei O[L,M,δ →]. Siehe 18: Beispielsweise können auf eine erste Verschiebefunktion 1810 eine erste Gewichtungsfunktion 1812 und eine Interpolationsfunktion 1814, eine zweite Gewichtungsfunktion 1816 und eine zweite Verschiebefunktion 1818 folgen. Falls bevorzugt, können andere Kombinationen von Projektionenerhöhungsoperatoren verwendet werden.
Mit diesen Definitionen ist die annähernde Reprojektion des Unter-Bildes f'
wobei P_M,P/L eine Reprojektion in das globale Koordinatensystem ist, bei der nur P/L Projektionen berechnet werden. Dies führt zu einer annähernden Zerlegung des Reprojektionsschritts für ein unterteiltes Bild analog zu Gleichung (20)
Diese annähernde Zerlegung ist in 19 für L = N/M = 2 dargestellt. In 19 ist der schnelle Reprojektionsprozessor 32 der 1 derart programmiert, dass er das elektronische Bild₂ 30, in der Figur auch als f bezeichnet, akzeptiert und bei 1910, 1912, 1914 und 1916 in die Unter-Bilder f₁, f₂, f₃ bzw. f₄ unterteilt. Reprojektionen in das globale Koordinatensystem werden bei 1918, 1920, 1922 und 1924 durchgeführt, wodurch die Unter-Sinogramme g₁', g₂', g₃' bzw. g₄' mit jeweils P/2 Projektionen erhalten werden. Die Unter-Sinogramme werden bei 1926, 1928, 1930 und 1932 vermehrt, wie zuvor beschrieben, und die erhaltenen Unter-Sinogramme g₁, g₂, g₃ und g₄ mit jeweils P Projektionen werden bei 1934 unter Bildung des Sinogramms 34 aggregiert, das in 19 auch als g bezeichnet ist.
Wie bei der annähernden Zerlegung für die Rückprojektion sind die Rechenkosten L-fach kleiner als bei der Ein-Schritt-Reprojektion P_N,P, wobei die gleiche Erläuterung gilt.
Bei einer bevorzugten Ausführungsform des schnellen Reprojektionsalgorithmus werden die Zerlegungen der Gleichungen (20) und (23) rekursiv angewendet, wobei auf jeder Ebene die Unter-Bilder gemäß Gleichung (11) in kleinere Unter-Bilder unterteilt werden, siehe 13. Die entsprechenden Beispiele für eine annähernde Zwei-Ebenen-Zerlegung und eine Zwei-Ebenen-Zerlegung, die aus einer Kaskade von einer genauen Zerlegung und einer annähernden Zerlegung besteht, sind in den 20 bzw. 21 dargestellt. In 20 beginnt eine annähernde Zwei-Ebenen-Zerlegung einer Fächerstrahl-Reprojektion für L = N/M = 2 bei dem elektronischen Bild₂ 30, in der Figur auch als f bezeichnet. Das Bild f wird bei 2010, 2012, 2014 bzw. 2016 in die Unter-Bilder f₁, f₂, f₃ und f₄ unterteilt, und das Unter-Bild f₁ wird bei 2018, 2020, 2022 und 2024 weiter in die Unter-Bilder f_1,1, f_1,2, f_1,3 und f_1,4 unterteilt. Diese Unter-Bilder werden bei 2026, 2028, 2030 bzw. 2032 in das globale Koordinatensystem reprojiziert, wodurch die Unter-Sinogramme g''_1,1, g''_1,2, g''_1,3 und g''_1,4 mit jeweils P/4 Projektionen erhalten werden. Die Unter-Sinogramme werden bei 2034, 2036, 2038 und 2040 vermehrt, wodurch die Unter-Sinogramme g'_1,1, g'_1,2, g'_1,3 und g'_1,4 mit jeweils P/2 Projektionen erhalten werden, die bei 2042 unter Bildung eines Unter-Sinogramms g'₁ aggregiert werden, das P/2 Projektionen umfasst. Nach einer weiteren Vermehrung bei 2044 wird das Unter-Sinogramm g₁ mit P Projektionen bei 2046 mit den Unter-Sinogrammen g₂, g₃ und g₄ aggregiert.
Das Unter-Bild f₄ wird bei 2048, 2050, 2052 und 2054 weiter zerlegt, und die Unter-Bilder werden bei 2056, 2058, 2060 und 2062 unter Bildung der Unter-Sinogramme g''_4,1, g''_4,2, g''_4,3 bzw. g''_4,4 ebenfalls reprojiziert. Nach Vermehrung bei 2064, 2066, 2068 und 2070 werden die erhaltenen Unter-Sinogramme g'_4,1, g'_4,2, g'_4,3 und g'_4,4 bei 2072 aggregiert, und das erhaltene Unter-Sinogramm wird bei 2074 weiter vermehrt. Das erhaltene Unter-Sinogramm g₄ wird bei 2046 mit den anderen Unter-Sinogrammen unter Bildung des Unter-Sinogramms 34 der 1 aggregiert, das in 20 auch als g bezeichnet ist. Die gestrichelten Kästen 2076 und 2078 stehen für die entsprechenden Zwei-Ebenen-Zerlegungen von f₂ bzw. f₃.
21 verdeutlicht eine gemischte (eine genaue, eine annähernde) Zwei-Ebenen-Zerlegung einer Fächerstrahl-Reprojektion für L = N/M = 2. In der Figur wird das elektronische Bilde 30 der 1, in 21 auch als f bezeichnet, bei 2110, 2112, 2114 bzw. 2116 in die Unter-Bilder f₁, f₂, f₃ und f₄ unterteilt. Das Unter-Bild f₁ wird durch 2118, 2120, 2122 und 2124 in die Unter-Bilder f_1,1, f_1,2, f_1,3 und f_1,4 unterteilt. Die Unter-Bilder werden dann durch 2126, 2128, 2130 und 2132 in das globale Koordinatensystem reprojiziert, wodurch die Unter-Sinogramme g'_1,1, g'_1,2, g'_1,3 und g'_1,4 mit jeweils P/2 Projektionen erhalten werden. Die Unter-Sinogramme g_1,1, g_1,2, g_1,3 und g_1,4 mit jeweils P Projektionen werden nach Vermehrung durch 2134, 2136, 2138 und 2140 erzeugt. Aggregation bei 2142 erzeugt ein Unter-Sinogramm g₁ mit P Projektionen, das bei 2144 mit den Unter-Sinogrammen g₂, g₃ und g₄ aggregiert wird unter Bildung des Sinogramms g mit P Projektionen. Ebenso wird das Unter-Bild f₄ wird bei 2146, 2148, 2150 und 2152 weiter zerlegt, wodurch die Unter-Bilder f_4,1, f_4,2, f_4,3 und f_4,4 erhalten werden. Eine Reprojektion durch 2154, 2156, 2158 und 2160 erzeugt die Unter-Sinogramme g'_4,1, g'_4,2, g'_4,3 und g'_4,4, die durch 2162, 2164, 2166 und 2168 vermehrt werden unter Bildung der Unter-Sinogramme g_4,1, g_4,2, g_4,3 und g_4,4. Diese Unter-Sinogramme werden bei 2170 unter Bildung des Unter-Sinogramms g₄ aggregiert. Die gestrichelten Kästen 2172 und 2174 stehen für die entsprechenden Ebene-2-Zerlegungen von g₂ bzw. g₃.
Die Rekursion kann fortgesetzt werden, bis die Unter-Bilder eine gewünschte minimale Größe M_min aufweisen. Dann können die Reprojektionen P_Pmin,Mmin unter Verwendung von Gleichung (7) berechnet werden.
Die folgenden Überlegungen für den schnellen Reprojektionsalgorithmus sind analog zu den im Fall der Rückprojektion erläuterten:

1. Auswahl verschiedener Parameter beim schnellen Reprojektionsalgorithmus, einschließlich P, N, M, L und M_min.
2. Handhabung der Diskretheit in t.
3. Verwendung von Winkel-Oversampling und der Kombination von genauen und annähernden Zerlegungsschritten, wodurch die Balance zwischen Genauigkeit und Geschwindigkeit gesteuert werden kann.
4. Vereinfachung der Abfolgen der Operatoren M ^_[δ ⇀]W_2,i und W_1,ijM ^₊[δ →_ij], indem Verschiebe- und Gewichtungsschritte zwischen aufeinander folgenden Stufen kombiniert werden, so dass ein einziger Verschiebe- und (höchstens) ein einziger Gewichtungsschritt pro Stufe resultiert.
5. Mögliche Umordnungen der Reihenfolge der Schritte.
6. Die -Rechenkomplexität von O(N²logN) erhaltenen hierarchischen Reprojektionsalgorithmus.
Schnelle native Algorithmen für die gleichwinkelige Fächerstrahl-Geometrie Im Fall eines gleichwinkeligen Detektors kann ein ähnliches Schema für die hierarchische Zerlegung des Bildes in kleinere Unter-Bilder angewendet werden. 22 stellt die Projektionsgeometrie für das Unter-Bild f' in diesem Fall dar.

Unter Verwendung der gleichen Definition für K_M wie in Gleichung (8) ist der Unter-Bild-Rückprojektionsoperator B_M,P[δ →] durch Gleichung (9) gegeben, wobei τ →(θ,r →) wie zuvor bei der gleichwinkeligen Geometrie definiert ist (siehe 6). Der Teil von g(p,·), der zur Rückprojektion auf f' beiträgt, ist wiederum K_M[δ →]g, wobei K_M[δ →] jetzt g(p,·) in Fächerstrahl-t auf den winkeligen Träger abschneidet, der durch die Winkel t_A und t_B der Projektion des Trägers des Unter-Bildes f' = K_M[δ →]f festgelegt wird. Die gleiche Erläuterung und die gleichen Argumente wie bei der colinearen, gleichmäßig beabstandeten Fächerstrahl-Geometrie gelten, was zu den Gleichungen (10)-(12) und zu der in 9 dargestellten genauen Zerlegung führt.
Die annähernde Zerlegung verringert unter Verwendung des folgenden Verfahrens wiederum die Anzahl an Projektionen, die zur Rekonstruktion des Unter-Bildes f' verwendet werden. Sei t_OO' der Winkel zwischen den zwei Strahlen, die durch die Mittelpunkte O und O' von f bzw. des Unter-Bildes f' gehen. In direkter Analogie zur t-Verschiebung der abgeschnittenen Projektionen um das Intervall OO'', das der Projektion des Mittelpunkts O' von f' im Fall des colinearen, gleichmäßig beabstandeten Detektors entspricht, werden die abgeschnittenen Projektionen für das Unter-Bild im Fächerwinkel um t_OO' t-verschoben. Gleichmäßig beabstandete Abtastpunkte mit dem gleichen Winkelabstand Δt werden mit dem Mittelpunkt bei t = 0 erhalten. Somit nimmt im gleich winkeligen Fall die Winkelvariable t die Rolle von t an, das bei der gleichmäßig beabstandeten Geometrie auftaucht, und τ →(θ,r →) ist wie in 7 definiert. Ansonsten gilt hier das gleiche rekursive Zerlegungsverfahren wie zuvor beschrieben. Genauer gesagt, werden abgeschnittene Fächerstrahl-Projektionen, die einem Unter-Bild entsprechen, nachdem sie in t verschoben wurden, durch L in Bezug auf den Source-Winkel θ dezimiert. Dann geht man nach dem übrigen Verfahren vor. Mit Ausnahme der bereits aufgeführten Veränderungen sind die Definitionen der verschiedenen Operatoren wenig modifiziert. Die annähernde Zerlegung wird wiederum von 11 beschrieben.
Ähnliche Argumente und Modifikationen gelten für den schnellen nativen Reprojektionsalgorithmus.
Schnelle native Rückprojektions- und Reprojektionsalgorithmen für die Kegelstrahl-Geometrie mit planarem, gleichmäßig beabstandetem Detektor
Das Verfahren für die allgemeine 3D-Divergentstrahl-Geometrie hat die gleichen Komponenten wie der 2D-Fall, mit den folgenden hauptsächlichen Unterschieden:

1. Das Bild wird entlang aller drei Koordinaten in 3D zerlegt.
2. Weil Projektionen zweidimensional sind, sind Abschneiden, Verschieben und Interpolation einzelner Projektionen zweidimensional.
3. Die Projektionen werden an der Source-Position θ entlang des Orbit und nicht im Source-Winkel auf dem Kreis dezimiert (oder interpoliert).
4. Weil der Source-Orbit oft in 3D nicht isotrop ist, kann die bevorzugte Zerlegung nicht-isotrop oder räumlich variabel sein.

Siehe 23: Man betrachte den Rückprojektionsschritt für ein Unter-Bild f'(r →) von f. K_M[δ →] sei ein durch Gleichung (8) definierter Bildabschneideoperator, so dass f' = K_M[δ →]f ein Unter-Bild von f der Größe M×M×M mit dem Mittelpunkt bei O' = δ → ∊ R³ ist. Die Gleichungen (9) und (10) gelten weiterhin, wobei τ durch den 2D-Vektor τ → ersetzt ist. Der Abschneideoperator K_M[δ →] schneidet die Projektion des gesamten Objektes auf den 2D-Träger Ω der Projektion des M×M×M-Trägers des Unter-Bildes f' = K_M[δ →]f ab (siehe 23).
Man betrachte jetzt eine Unterteilung des Bildes f in J = (N/M)³ nichtüberlappende Unter-Bilder mit jeweils einer Größe von M×M×M
wobei die Vektoren δ_i,j = 1, ... J die Mittelpunkte der Unter-Bilder sind, wie in 24 für den Fall M/N = 2 dargestellt.
Gleichung (12) gilt und liefert eine genaue Zerlegung für die Rückprojektion in Rückprojektionen auf die Unter-Bilder. Dies wird durch eine Figur analog zu 9 beschrieben, nur dass sie acht anstelle von vier "Kanälen" aufweisen würde. Die Erläuterung der Rechenkosten ist ähnlich wie beim 2D-Fall, ausgenommen dass die Rechenkosten von B_N,P und B_M,P cN³P bzw. cM³P bei einem Verhältnis J wie im 2D-Fall betragen.
Wie im 2D-Fächerstrahl-Fall verwendet der schnelle 3D-Kegelstrahl-Rückprojektionsalgorithmus die Zerlegungsidee mit der folgenden zusätzlichen Eigenschaft. Bei einer festen Bildauflösung (Bandbreite) ist die Anzahl der Projektionen, die zur Rekonstruktion des Bildes benötigt wird, proportional zur Bildgröße. D.h. werden P Projektionen zur Rekonstruktion eines N×N×N-Bildes benötigt, dann reichen P' = (M/N)P Projektionen zur Rekonstruktion eines M×M×M-Bildes mit der gleichen Auflösung.
Das Unter-Bild f' = K_M[δ →]f wird wiederum aus einer verringerten Anzahl an Projektionen rekonstruiert, wobei der Operator O[L,M,δ →] zur Durchführung der Reduktion verwendet wird. Er ist wie zuvor durch die Zusammensetzung ähnlicher Operatoren definiert, mit den folgenden Änderungen:

1. Der Projektionenabschneideoperator K_M[δ →] schneidet nun die 2D-Region Ω ab, wie vorstehend erläutert.
2. Die Projektionenverschiebeoperatoren M ^±[δ] führen eine 2D-Verschiebung um den Vektor ±r ⇀(θ,δ ⇀) in den t₁,t₂-Koordinaten durch und keine 1D-Verschiebung wie bei der Fächerstrahl-Geometrie. Diese Verschiebung kann trennbar erfolgen, d.h. zunächst eine Verschiebung in t₁ und dann in t₂.
3. Der Operator D_↓L dezimiert an der Source-Position und nicht im Winkel. Besteht der Orbit aus der Vereinigung von zwei oder mehr getrennten Teilen, wie zwei senkrecht aufeinander stehenden Kreisen, erfolgt die Dezimierung in der Source-Position für jeden Teil getrennt.
4. Das Filtern in der Variablen t ⇀ bei der Dezimierung nach Gleichung (14) ist zwei- und nicht eindimensional. Es kann ebenfalls in t₁ und t₂ getrennt durchgeführt werden, aber das t ⇀-Filtern kann möglicherweise nicht von dem Filtern in θ abgetrennt werden.

Die annähernde Zerlegung wird auch hier durch Gleichung (18) gegeben, wobei B_M,P/L hier eine Divergentstrahl-Rückprojektion in das globale Koordinatensystem auf ein M×M×M-Unter-Bild unter Verwendung von P/L Projektionen ist, wie in Gleichung (9) definiert. Für L = N/M = 2 wird diese Zerlegung durch eine Figur analog zu 10 beschrieben, ausgenommen dass sie acht anstelle von vier "Kanälen" aufweisen würde. Die Erläuterung der Rechenkosten ist ähnlich wie beim 2D-Fall, ausgenommen dass die Rechenkosten von B_N,P und B_M,P/L cN³P bzw. cM³P/L sind und J = (N/M)³. Die Gesamtkosten für J Unter-Bild-Rückprojektionen sind erneut L-mal kleiner als für die herkömmliche Gesamtbild-Rückprojektion B_N,P.
Bei einer bevorzugten Ausführungsform des schnellen Rückprojektionsalgorithmus werden die Zerlegungen der Gleichungen (12) und (18) rekursiv angewendet, wobei auf jeder Ebene die Unter-Bilder gemäß Gleichung (24) in kleinere Unter-Bilder unterteilt werden. Die Beschreibung des Verfahrens und verschiedene Merkmale und Verbesserungen, die vorstehend für die Fächerstrahl-Geometrie dargelegt sind, gelten auch für die allgemeine Kegelstrahl-Geometrie mit ähnlichen Änderungen, wie sie für die genaue und die annähernde Zerlegung selbst aufgeführt sind. Genauer gesagt, werden die Gesamtkosten für die Kegelstrahl-Rückprojektion unter Verwendung des rekursiven Zerlegungsalgorithmus zu O(N³logN) im Vergleich zu O(N⁴) für die herkömmliche Kegelstrahl-Rückprojektion.
Bis hier haben wir die einfachste Ausführungsform einer 3D-Bildzerlegung beschrieben, wobei das Bild auf allen drei Achsen gleichmäßig zerlegt wird. In einigen Fällen der Kegelstrahl-Geometrie ist möglicherweise der Source-Orbit nicht isotrop in 3D, so dass die bevorzugte Zerlegung ebenfalls nicht isotrop ist, und sie kann zudem räumlich variabel sein. Wichtige Beispiele sind der einzelne kreisförmige und der helixförmige Kegelstrahl-Orbit (4, 6 und 3). Diese Fälle zeigen eine zylindrische Symmetrie (oder Fast-Symmetrie) des Source-Orbit um die z-Achse, so dass die auf der z-Achse gewählte Bildunterteilung anders als bei der x- und der y-Achse sein kann. Sogar wenn der gesamte Orbit keine Symmetrieachse hat, kann er manchmal in die Vereinigung von Teilen zerlegt werden, die eine Symmetrieachse haben, wie bei einem Source-Orbit aus zwei senkrecht aufeinander stehenden Kreisen. Der Rückprojektionsschritt kann dann durchgeführt werden als Summe von zwei Rückprojektionen, jeweils mit ihrem eigenen, getrennten, axial symmetrischen Orbit.
In diesen Fällen kann es genügen, wenn Teile des Bildes weniger häufig unterteilt werden oder in Unter-Bilder, die auf der z-Achse größer sind. Bei einem extremen Beispiel kann sich die Unterteilung überhaupt auf die (x,y)-Richtungen beschränken, siehe 25. Eine allgemeinere nicht-isotrope und räumlich ungleichmäßige Unterteilung ist in 26 dargestellt. Dies kann zu Recheneinsparungen führen, indem die Anzahl an erforderlichen Verschiebungen, Interpolationen oder Filterschritten an Projektionen entlang der t₂-Achse verringert wird.
Der Grund, warum eine solche verringerte Unterteilung funktioniert, ist für den Fall eines einzelnen kreisförmigen Source-Orbit aus den Gleichungen (2) und (3) ersichtlich. Die Gleichungen beschreiben die θ-Abhängigkeit der Position der Projektion von einem Punkt auf einem Detektor. Man betrachte einen Punkt r → mit kleinem z und/oder kleiner
D.h. einen Punkt nahe der Ebene des Source-Orbit und/oder nahe der z-Achse. Man beachte, dass für einen solchen Punkt τ₂(θ,r →) viel langsamer mit θ variiert als τ₁(θ,r ⇀). Damit kann man zeigen, dass ein Unter-Bild, dessen Mittelpunkt δ ⇀ derart lokalisiert ist, aus einer verringerten Anzahl an Ansichten genau rekonstruiert werden kann, sogar wenn seine z-Dimension nicht im gleichen Ausmaß verringert wird wie seine (x,y)-Dimensionen. Infolgedessen genügt eine feinere Unterteilung in der (x,y)-Ebene, damit eine genaue Rekonstruktion aus einer geringeren Anzahl an Projektionen gewährleistet ist. Ein ähnliches Argument gilt im Fall eines helixförmigen Orbit, wenn eine globale Verschiebung um hθ(2π) in t₂ zunächst auf die Projektionen angewendet wird, wodurch die Bewegung des Mittelpunkts der Koordinaten auf dem Detektor entlang der z-Achse kompensiert wird, wobei θ variabel ist (siehe 3).
Mit dieser nicht-isotropen und möglicherweise ungleichmäßigen Unterteilung wird die rekursive Zerlegung entsprechend modifiziert. Wird die Unterteilung der 25 rekursiv angewendet, hat die entsprechende Zerlegung vier Kanäle auf jeder Rekursionsebene, wobei nur in der t₁-Variablen verschoben wird. Sie wird durch die 10-15 beschrieben. Wird dagegen die Unterteilung der 25 bei einer annähernden Zwei-Ebenen-Zerlegung eingesetzt, hat die erste Zerlegungsebene acht Kanäle (entsprechend den acht Oktanten bei der groben Unterteilung des Würfels). Dann hat aber jede der Zerlegungen auf der nächsten Ebene nur sieben Kanäle, die den sieben Bildern entsprechen, in die jeder grobe Bild-Oktant unterteilt wird.
Die Gestaltung einer schnellen Reprojektion für die allgemeine Divergentstrahl-Geometrie, einschließlich der oben erläuterten Spezialfälle der Divergentstrahl-Geometrie, ist das Tandem zu den entsprechenden schnellen Rückprojektionsalgorithmen. Sie hat die gleiche Beziehung zur schnellen Rückprojektion bei der Divergentstrahl-Geometrie wie die Fächerstrahl-Reprojektion zur schnellen Fächerstrahl-Rückprojektion.
Weil die letzte Beziehung bereits in gewisser Ausführlichkeit erläutert wurde, wird sie hier nicht wiederholt.
Implementierungen und Experimente
Die neuen Algorithmen wurden in den Programmiersprachen MATLAB^TM und C implementiert und in mehreren numerischen Experimenten erfolgreich getestet, die die Rekonstruktion des Shepp-Logan-Kopfphantoms (eines Standard-Testbildes, das bei numerischen Untersuchungen tomographischer Algorithmen eingesetzt wird) der Größe N×N = 512×512 im 2D-Fall und der Größe N×N×N = 128×128×128 im 3D-Fall beinhalteten. Wir verwendeten das so genannten unmodifizierte Phantom mit einem realistischen (und hohen) Kontrast zwischen Hirnschädel und Hirn. Dies stellt einen herausfordernden Test für die Algorithmen dar, weil sogar kleine Artefakte, die durch die Rekonstruktion des Hirnschädels erzeugt werden, vor dem Hirn-Hintergrund hervorstechen.
Zunächst werden die Ergebnisse der 2D-Fächerstrahl-Rekonstruktion beschrieben. Bei gleichmäßig beabstandeten und gleichwinkeligen Fächerstrahl-Bildgebungsgeometrien stellten wird den Source-Bildmittelpunkt-Abstand auf das 1,25-fache der Bildgröße ein, d.h. D = 1,25N, den Gesamt-Fächerwinkel auf 2α₀ = 1,17 Radian, wobei sich 1025 Detektoren auf dem Fächer befanden. Wir erzeugten Voll-Scan-Fächerstrahl-Projektionsdaten mit P = 1536 gleichmäßig beabstandeten Source-Winkeln θ_p ∊ [0,2π] für das 10-Ellipsen-Phantom unter Verwendung analytischer Ausdrücke für die Geradenintegrale durch eine Ellipse.
Bei unserer Simulation zerlegten wir das Bild auf jeder Ebene in 4 Unter-Bilder mit jeweils der halben Größe. Für die annähernden Zerlegungen verwendeten wir einen Dezimierungsfaktor von L = 2. Das Bild wurde iterativ in Ein-Pixel-Unter-Bilder zerlegt, d.h. M_min = 1. Zur Steigerung des Winkel-Oversampling waren die ersten Q Ebenen der Zerlegung genaue Zerlegungen ohne Winkel-Dezimierung. Wir verwendeten Q = 2 und Q = 3, wobei Q als "Holdoff-Parameter" bezeichnet wird. Die Parameter ("Footprints") für den Projektionenabschneideoperator K_M[δ →] wurden direkt unter Verwendung von einfacher herkömmlicher Annäherung auf Kosten eines gewissen Rechenanstiegs berechnet.
Die 27 und 28 vergleichen den genauen Rückprojektionsalgorithmus und unsere schnellen Fächerstrahl-Rückprojektionsalgorithmen für die beiden Fächerstrahl-Bildgebungsgeometrien. Als Bezug und Vergleichspunkt zeigen die 27(a) und 28(a) das originale Phantom. Die 27(b) und 28(b) sind die rekonstruierten Bilder unter Verwendung von genauer Fächerstrahl-Rückprojektion für gleichmäßig beabstandete und gleichwinkelige Detektoren. Die 27(c), 28(c), 27(d) und 28(d) sind die Rekonstruktionen unter Verwendung unserer schnellen Algorithmen mit dem Holdoff-Parameter Q = 2 bzw. Q = 3. Die Beschleunigungen durch die schnellen Algorithmen betragen 30X und 60X für Holdoff Q = 3 bzw. Q = 2. Die Rekonstruktionen unter Verwendung der genauen und schnellen Algorithmen zeigen sogar für Q = 2 wenig wahrnehmbaren Qualitätsunterschied, wenn überhaupt. Diagramme durch Reihen und Spalten der rekonstruierten Bilder der 28 sind in den 29(a) und 29(b) dargestellt. 29(a) vergleicht Schnitte durch die 410. Reihe. 29(b) vergleicht Schnitte durch die 356. Reihe der 28. In diesen Diagrammen steht die durchgezogene Linie für das ursprüngliche Phantom, die gestrichelte Linie für die genaue Fächerstrahl-Rückprojektion, die Strichpunktlinie für den schnellen Algorithmus mit Holdoff 3 und die gepunktete Linie für den schnellen Algorithmus mit Holdoff 2.
Als nächstes vergleichen wir die Punktverbreiterungsfunktionen (PSF) der genauen und schellen Fächerstrahl-Rückprojektionsalgorithmen für die Geometrie mit gleichmäßig beabstandeten Detektoren. Wir führten die genauen und schellen Algorithmen an einem Bild der Größe 256×256 mit vier Punktzielen an den Positionen (10,10), (128,128), (49,37) und (201,93) durch. Wir verwendeten ähnliche Parameter wie zuvor: Source-Mittelpunkt-Abstand D = 1,25N, Anzahl der Source-Winkel P = 2N und Gesamt-Fächerwinkel von 1,20 Radian mit 737 gleichmäßig beabstandeten Detektoren. Die Punktreaktionen waren an allen vier Positionen sehr ähnlich, was auf eine im Wesentlichen verschiebungsunveränderliche Reaktion hindeutet. Die Punktverbreiterungsfunktionen (PSF) für das Punktziel bei (128,128), die mit den genauen und schnellen Algorithmen (mit Q = 3) erhalten wurden, sind in den 30(a) und 30(b) verglichen. 30(a) stellt eine Mittelpunkt-Pixel-PSF für den genauen Algorithmus dar, 30(b) eine Mittelpunkt-Pixel-PSF für den schnellen Algorithmus mit Q gleich 2. Schnitte durch die beiden Achsen der PSF sind in den 31(a) und 31(b) dargestellt. 31(a) stellt einen Schnitt durch die 128. Reihe für den genauen Algorithmus (durchgezogene Linie) und den schnellen Algorithmus mit Q = 3 (gestrichelte Linie) dar. 31(b) stellt einen Schnitt durch die 128. Spalte für den genauen Algorithmus (durchgezogene Linie) und den schnellen Algorithmus mit Q = 3 (gestrichelte Linie) dar. Die PSF der beiden Algorithmen sind fast identisch, was die Genauigkeit der schnellen Algorithmen bestätigt.
Als nächstes erläutern wir die Ergebnisse von Experimenten mit einer Kegelstrahl-Rekonstruktion mit einzelnem kreisförmigem Source-Orbit in der (x,y-Ebene). Der Abstand von der Source zum Mittelpunkt des Volumens betrug das 1,25-fache der Länge der Lateralen. Kegelstrahl-Projektionen wurden bei 512 Source-Winkeln aufgenommen, die in [0,2π] gleichmäßig beabstandet waren. Der Fächerwinkel betrug (im Azimuth) 1,20 Radian, der Kegelwinkel 1,5 Radian. Wir berücksichtigten die planare, gleichmäßig beabstandete und die zylindrische Detektorgeometrie, die in 4 bzw. 6 dargestellt sind. In jedem Fall bestand die Detektoroberfläche aus 375×375 Detektoren, die entlang der vertikalen (t₂-) Achse gleichmäßig beabstandet waren und sich in der t₁-Koordinate in gleichmäßig beabstandeten oder gleichwinkeligen Intervallen befanden.
Der herkömmliche Rekonstruktionsalgorithmus für den einzelnen kreisförmigen Kegelstrahl-Orbit mit planarem Detektor ist der herkömmliche Feldkamp-Davis-Kress(FDK-)Algorithmus, der umfasst: Gewichten der Projektionen, Filtern entlang der t₁-Koordinate und anschließende Durchführung einer gewichteten Kegelstrahl-Rückprojektion. Gewichtung und Filtern entsprechen denjenigen bei der Fächerstrahl-Geometrie mit planarem Detektor, wobei der Unterschied in der Verwendung von Kegelstrahl-anstelle von Fächerstrahl-Geometrie besteht. Dies ist ein annähernder Algorithmus, der an einigen inhärenten Artefakten bei großen Kegelwinkeln leidet. Trotzdem ist er ein Standard-Bezugspunkt für diese Kegelstrahl-Geometrie. Für die zylindrische Detektorgeometrie verwendeten wird den analogen Algorithmus, der erhalten wird, indem man den Fächerstrahl-Algorithmus für die gleichwinkelige Detektorgeometrie durch Einbeziehen einer gewichteten Kegelstrahl-Rückprojektion auf den zylindrischen Detektor auf die Kegelstrahl-Geometrie ausweitet.
Zu Demonstrationszwecken verwendeten die entsprechenden schnellen Algorithmen Rückprojektion mit den zuvor beschriebenen rekursiven Zerlegungen und dem in 25 gezeigten einfachen Unterteilungsschema. Bei unseren Simulationen beendeten wir das Zerlegungsverfahren, wenn die Unter-Bild-Größe auf 4×4×128 verringert war, weil der Overhead bei unserer besonderen Implementierung des Algorithmus der Beschleunigung der Rückprojektionsoperatoren bei diesen kleinen Volumina entgegenwirkte. Die schnellen Algorithmen verwendeten das gleiche Filtern und Gewichten wie der herkömmliche FDK-Algorithmus.
Die Ergebnisse der herkömmlichen FDK-Rekonstruktion und der schnellen Rekonstruktion unter Verwendung des erfindungsgemäßen Verfahrens sind in 32 und 33 für die planare bzw. die zylindrische Detektorgeometrie verglichen. 32 vergleicht 3D-Kegelstrahl-Rekonstruktionen für einen einzelnen kreisförmigen Orbit und eine Geometrie mit planaren, gleichmäßig beabstandeten Detektoren. Die obere Reihe ist ein xy-Schnitt bei z = –0,25, die mittlere Reihe ein yz-Schnitt bei x = –0,13 und die untere Reihe ein xz-Schnitt bei y = –0,04. Die linke Spalte ist das Rekonstruktionsergebnis bei Verwendung des herkömmlichen FDK-Algorithmus, die mittlere und die rechte Spalte ver wenden den erfindungsgemäßen schnellen Algorithmus mit Holdoff 2 bzw. 1. 33 zeigt einen Vergleich der 3D-Kegelstrahl-Rekonstruktionen für einen einzelnen Source-Orbit und eine zylindrische Detektorgeometrie. Die obere Reihe ist ein xy-Schnitt bei z = –0,25, die mittlere Reihe ein yz-Schnitt bei x = –0,13 und die untere Reihe ein xz-Schnitt bei y = –0,04. Die linke Spalte ist das Rekonstruktionsergebnis bei Verwendung des herkömmlichen FDK-Algorithmus, die mittlere und die rechte Spalte verwenden den erfindungsgemäßen schnellen Algorithmus mit Holdoff 2 bzw. 1. 34 zeigt Schnitte durch Bilder in den entsprechenden Reihen der 32. Das obere Diagramm zeigt Reihe 102, das mittlere Diagramm Reihe 80 und das untere Diagramm Spalte 81. Die durchgezogenen Linien stammen von der herkömmlichen FDK-Rekonstruktion. Der schnelle Algorithmus mit Holdoff 2 ist als Strichpunktlinien dargestellt, der schnelle Algorithmus mit Holdoff 1 als gepunktete Linien. 35 veranschaulicht Schnitte durch Bilder in den entsprechenden Reihen der 33. Das obere Diagramm zeigt Reihe 102, das mittlere Diagramm Reihe 80 und das untere Diagramm Spalte 81. Die durchgezogenen Linien stammen von der herkömmlichen FDK-Rekonstruktion, die Strichpunktlinien zeigen Ergebnisse unter Verwendung des schnellen Algorithmus mit Holdoff 2 und die gepunkteten Linien Ergebnisse unter Verwendung des schnellen Algorithmus mit Holdoff 1. Die 34 und 35 zeigen Diagramme der rekonstruierten Dichte entlang bestimmter Geraden durch die Bilder der 32 und 33 und liefern zusätzliche Einzelheiten. Es sind die Ergebnisse für Holdoff = 1 und Holdoff = 2 angegeben.
Die schnelle Version des FDK-Algorithmus unter Verwendung der Erfindung mit dem Holdoff-Faktor Q = 2 ist etwa dreimal schneller als der herkömmliche Feldkamp-Algorithmus, wobei die Bildqualität wenig schlechter wird, wenn überhaupt. Bei 7-facher Beschleunigung ist jedoch die Bildqualität sogar mit dem Holdoff-Faktor Q = 1 für einige Anwendungen annehmbar. Wir erwarten, dass sich die Beschleunigung auf mehrere Zehnfache für Bildvolumina von 512×512×512 und eine volle 3D-Unterteilung erhöht, wie in den 24 oder 26 und der entsprechenden Zerlegung. (Tatsächlich werden bei Kegelstrahl-Geometrien mit fast-zylindrischer Symmetrie, wie einem kreisförmigen oder helixförmigen Orbit, Beschleunigungen auf mehrere Zehnfache für 512×512×H-Bilder für jede beliebige z-Größe von H erhalten.)
Experimentelle Schlussfolgerung
Bei unseren Experimenten in 2D lieferten die neuen Fächerstrahl-Rückprojektionsalgorithmen praktische Beschleunigungen von 30X-60X für 512×512-Bilder ohne sichtbaren Verlust bei der Genauigkeit. Bei 3D mit einzelnem kreisförmigen Orbit und Kegelstrahl-Geometrie mit planarem Detektor unter Verwendung der einfachsten Form der Zerlegung betrug die Beschleunigung 3X-7X für 128×128×128. Beschleunigungsfaktoren im Bereich des Zehnfachen sollten bei voller Implementierung der vorgeschlagenen erfindungsgemäßen Zerlegung und für 512×512×512-Bilder möglich sein.
Die Experimente zeigten die Effizienz des erfindungsgemäßen Verfahren bei unterschiedlichen Detektorgeometrien. Diese Ergebnisse für die Geometrie mit einzelnem kreisförmigen Orbit wird sich auf die allgemeinere Divergentstrahl-Geometrie, wie den helixförmigen Kegelstrahl, ausweiten. Aufgrund der ähnlichen Prinzipien und der ähnlichen Struktur von Reprojektions- und Rückprojektionsalgorithmen sollten ähnliche Ergebnisse für die erfindungsgemäßen schnellen Divergentstrahl-Reprojektionsalgorithmen gelten. Reprojektions- und Rückprojektionsalgorithmen sind einzeln anwendbar; kombiniert können sie für verschiedene iterative Ansätze zur tomographischen Rekonstruktion oder zur Artefaktkorrektur bei der Tomographie eingesetzt werden, wie zum Beispiel in den US-Patenten 6 263 096 und 6 351 548 beschrieben. Eine Kombination aus Rückprojektion und Reprojektion kann auch zur Konstruktion schneller Reprojektionsalgorithmen für das Problem der langen Objekte bei der helixförmigen Kegelstrahl-Geometrie verwendet werden.
Die Prinzipien der Erfindung sind zwar vorstehend in Verbindung mit speziellen Apparaturen und Anwendungen beschrieben, aber es sollte selbstverständlich sein, dass diese Beschreibung nur als Beispiel und nicht als Beschränkung für den Umfang der Erfindung gegeben wurde.

Claims

Verfahren zum Erzeugen eines elektronischen Bildes aus einem Divergentstrahl-Sinogramm (22), das einer Rückprojektion unterworfen werden kann, umfassend die Schritte: Unterteilen (1010, 1012, 1014, 1016) des Sinogramms (22) in eine Mehrzahl Unter-Sinogramme (g₁, g₂, g₃, g₄ in 10); Durchführen einer gewichteten Rückprojektion der Unter-Sinogramme in das globale Koordinatensystem (1018, 1020, 1022, 1024) zum Erzeugen einer Mehrzahl entsprechender Unter-Bilder (f₁, f₂, f₃, f₄ in 10) an richtigen Stellen im globalen Koordinatensystem; und Aggregieren (1026) der Unter-Bilder (f₁, f₂, f₃, f₄ in 10) zum Erzeugen des elektronischen Bildes (f, 26 in 10).
Verfahren nach Anspruch 1, wobei die Sinogramm-Unterteilung annähernd ist (1210, 1212, 1214, 1216).
Verfahren nach Anspruch 1, wobei das Sinogramm in eine Mehrzahl Unter-Sinogramme rekursiv unterteilt wird, bis jedes Unter-Sinogramm ein Unter-Bild mit gewünschter Größe darstellt, wobei die Unterteilungsschritte eine gewünschte Anzahl genauer Unterteilungen (1510, 1512, 1514, 1516) und eine gewünschte Anzahl annähernder Unterteilungen (1518, 1520, 1522, 1524) umfassen.
Verfahren nach Anspruch 3, wobei die genauen Unterteilungsschritte Approximation ohne Verschiebung beinhalten.
Verfahren nach Anspruch 3, wobei die Unter-Sinogramme Unter-Bildern entsprechen, die eine Größe von nur einem Pixel oder Voxel haben.
Verfahren nach Anspruch 1, wobei der Aggregationsschritt rekursiv durchgeführt wird.
Verfahren nach Anspruch 1, wobei das elektronische Bild eine tomographisches Bild ist.
Verfahren nach Anspruch 1, das zudem eine Vorverarbeitung (16) umfasst, wobei Überabtastung in den Datenkoordinaten zur Verbesserung der Genauigkeit des elektronischen Bildes genutzt wird.
Verfahren nach Anspruch 1, wobei das Sinogramm vorverarbeitete Divergentstrahl-Projektionen enthält.
Verfahren nach Anspruch 1, wobei die Sinogramm-Unterteilungsschritte annähernd sind, Approximierung (1110), Verschiebung (1112, 1120) und Dezimierung (1116) des Sinogramms umfassen, so dass jedes erzeugte Unter-Sinogramm einem Unter-Bild in dem globalen Koordinatensystem entspricht.
Verfahren nach Anspruch 10, wobei die annähernden Unterteilungsschritte zudem Gewichtung (1114, 1118) enthalten.
Verfahren nach Anspruch 10, wobei die Unterteilungen rekursiv durchgeführt werden und wobei Verschiebevorgänge in aufeinander folgenden Stufen der Rekursion zur Verringerung der Rechnerkosten kombiniert werden.
Verfahren nach Anspruch 1, wobei Unterteilungen nicht-isotrop und ungleichmäßig im Raum durchgeführt werden, so dass die Unter-Bilder auf einer gegebenen Rekursionsebene nicht die gleiche Form oder Größe haben müssen.
Verfahren nach Anspruch 13, wobei die Teilung anhand der Symmetrieeigenschaften des Source Orbit gewählt werden, so dass eine gewünschte Balance zwischen Rechnerkosten und Bildgenauigkeit erhalten wird.
Verfahren nach Anspruch 1, das auf mehrere Sinogramme einzeln angewendet wird, die aus mehreren Source Orbits erhalten wurden, so dass Bilder erzeugt werden, die zu dem endgültigen Bild kombiniert werden können.
Verfahren nach Anspruch 1, wobei der genaue Unterteilungsschritt derart eingerichtet wird, dass der Vorgang der Approximierung von Projektionen in einem Si nogramm oder Unter-Sinogramm für ein Unter-Bild beginnt, nachdem ein oder mehr Projektionen für dieses Unter-Bild verfügbar werden.
Verfahren nach Anspruch 1, wobei der Unterteilungsschritt derart eingerichtet wird, dass die Unterteilung eines Sinogramms für ein Unter-Bild beginnt, nachdem ein oder mehr Projektionen für dieses Unter-Bild verfügbar werden.
Verfahren nach Anspruch 1, das zeitlich nacheinander auf mehrere Divergentstrahl-Sinogramme angewendet wird, sobald diese verfügbar sind.
Verfahren nach Anspruch 1, wobei das Divergentstrahl-Sinogramm in einer Konvergentstrahl-Geometrie aufgenommen wird.
Vorrichtung zum Erzeugen eines elektronischen Bildes (36) von einem Objekt, umfassend: Einrichtungen (12) zum Abtasten des Objekts, wodurch Daten erzeugt werden, die ein Bild des Objekts darstellen; Einrichtungen (16, 20) zum Manipulieren der Daten, so dass ein verarbeitetes Sinogramm (22) erzeugt wird; Einrichtungen (24) zum Unterteilen des Sinogramms (22) in eine Mehrzahl Unter-Sinogramme (1010, 1012, 1014, 1016); Einrichtungen (24) zum Rückprojizieren jedes Unter-Sinogramms in das globale Koordinatensystem (1018, 1020, 1022, 1024) zum Erzeugen einer Mehrzahl entsprechender Unter-Bilder (f₁, f₂, f₃, f₄ in 10) an richtigen Stellen im globalen Koordinatensystem; Einrichtungen (24) zum Aggregieren der Unter-Bilder zum Erzeugen des elektronischen Bildes; und Einrichtungen (38) zum Speichern und/oder Anzeigen und/oder Analysieren des elektronischen Bildes.
Vorrichtung nach Anspruch 20, wobei die Sinogramm-Unterteilung annähernd ist (1210, 1212, 1214, 1216).
Vorrichtung nach Anspruch 20, wobei die Sinogramme in eine Mehrzahl Unter-Sinogramme rekursiv unterteilt werden, bis jedes Unter-Sinogramm einem Unter- Bild mit gewünschter Größe entspricht, und die Einrichtung zum Unterteilen eine gewünschte Anzahl genauer Unterteilungen (1510, 1512, 1514, 1516) und eine gewünschte Anzahl annähernder Unterteilungen (1518, 1520, 1522, 1524) durchführt.
Vorrichtung nach Anspruch 22, wobei die genauen Unterteilungen Approximation ohne Verschiebung beinhalten.
Vorrichtung nach Anspruch 22, wobei die Unter-Sinogramme Bildern entsprechen, die eine Größe von nur einem Pixel oder Voxel haben.
Vorrichtung nach Anspruch 20, wobei die Einrichtung zum Aggregieren rekursiv arbeitet.
Vorrichtung nach Anspruch 20, wobei das elektronische Bild eine tomographisches Bild ist.
Vorrichtung nach Anspruch 20, wobei die Einrichtung zum Manipulieren eine Überabtastung in den Datenkoordinaten (16) zur Verbesserung der Genauigkeit des elektronischen Bildes durchführt.
Vorrichtung nach Anspruch 20, wobei die Einrichtung zum Unterteilen derart betrieben werden kann, dass sie eine annähernde Unterteilung inklusive Approximierung (1110), Verschiebung (1112, 1120) und Dezimierung (1116) des Sinogramms durchführt, so dass jedes erzeugte Unter-Sinogramm einem Unter-Bild in dem globalen Koordinatensystem entspricht.
Vorrichtung nach Anspruch 28, wobei die annähernde Unterteilung zudem Gewichtung (1114, 1118) enthält.
Verfahren zum Erzeugen eines Divergentstrahl-Sinogramms (34) aus einem elektronischen Bild, umfassend die Schritte: Unterteilen (1710, 1712, 1714, 1716) des Bildes in eine Mehrzahl Unter-Bilder (f₁, f₂, f₃, f₄ in 17 und in 9); Berechnen (1718, 1720, 1722, 1724) von Unter-Sinogrammen jedes Unter-Bildes in dem globalen Koordinatensystem (g₁, g₂, g₃, g₄); und Aggregieren der Unter-Sinogramme (1726) zum Erzeugen des Sinogramms (34).
Verfahren nach Anspruch 30, wobei die Berechnung von Unter-Sinogrammen annähernd ist (1926, 1928, 1930, 1932).
Verfahren nach Anspruch 30, wobei die Unterteilung rekursiv durchgeführt wird (z.B. wie in 13), bis jedes Unter-Bild eine gewünschte Größe hat, und die Berechnung von Unter-Sinogrammen in einer gewünschten Anzahl an Ebenen der Rekursion annähernd ist und in den verbleibenden Ebenen der Rekursion genau ist.
Verfahren nach Anspruch 30, wobei die Aggregationen rekursiv durchgeführt werden.
Verfahren nach Anspruch 30, wobei die Aggregationen durch Addieren von Unter-Sinogrammen ohne Verschiebung durchgeführt werden.
Verfahren nach Anspruch 30, wobei die Berechnung von Unter-Sinogrammen annähernd ist und durch Verschieben (1810, 1818), Interpolieren (1814) und Addieren von Unter-Sinogrammen durchgeführt wird und jedes verschobene und interpolierte Unter-Sinogramm einem Unter-Bild im globalen Koordinatensystem entspricht.
Verfahren nach Anspruch 35, wobei die Schritte der annähernden Berechnung von Unter-Sinogrammen auch Gewichten (1812, 1816) enthalten.
Verfahren nach Anspruch 32, wobei die genaue Berechnung von Unter-Sinogrammen ohne Verschieben durchgeführt wird.
Vorrichtung zum Erzeugen eines Bildes (38) von einem Objekt, umfassend: einen Scanner (12), der Daten von einem Objekt erzeugt; einen Prozessor (16) zum Erzeugen mindestens einer Divergentstrahl-Projektion von dem Bild; Einrichtungen zum Rekonstruieren (24) oder, bei einer iterativen Rekonstruktionsvorrichtung, zum Rückprojizieren eines Bildes aus mindestens einer Projektion; Einrichtungen (28) zum Ermitteln von Fehlern in dem durch die Rekonstruktionseinrichtung erzeugten Bild; Einrichtungen (32) zum Reprojizieren eines Bildes nach der Fehlerkorrektur, Durchführen zusätzlicher Korrekturen am Sinogramm und Zuführen (20) des korrigierten Sinogramms der Rekonstruktions- oder Rückprojektionseinrichtung (24); und Einrichtungen (38) zum Anzeigen des Bildes, das von der Rekonstruktionseinrichtung nach Korrektur der Fehler erzeugt wird, wobei die Reprojektionseinrichtung das Bild in eine Mehrzahl Unter-Bilder (f₁, f₂, f₃, f₄ in 17) unterteilt, Unter-Sinogramme der Unter-Bilder in dem globalen Koordinatensystem (g₁, g₂, g₃, g₄ in 17) berechnet und die Unter-Sinogramme (1726) zum Erzeugen des Sinogramms (g in 17) aggregiert.
Vorrichtung nach Anspruch 38, wobei die Berechnung von Unter-Sinogrammen annähernd ist.
Vorrichtung nach Anspruch 38, wobei die Unterteilung rekursiv durchgeführt wird, bis jedes Unter-Bild eine gewünschte Größe hat, und die Berechnung von Unter-Sinogrammen in einer gewünschten Anzahl an Ebenen der Rekursion annähernd ist und in den verbleibenden Ebenen der Rekursion genau ist.
Vorrichtung nach Anspruch 38, wobei die Aggregationen rekursiv durchgeführt werden.
Vorrichtung nach Anspruch 38, wobei die Aggregationen durch Addieren von Unter-Sinogrammen ohne Verschiebung durchgeführt werden.
Vorrichtung nach Anspruch 38, wobei die Berechnung von Unter-Sinogrammen annähernd ist und durch Verschieben, Interpolieren und Addieren von Unter-Sinogrammen durchgeführt wird und jedes verschobene und interpolierte Unter-Sinogramm einem Unter-Bild im globalen Koordinatensystem entspricht.
Vorrichtung nach Anspruch 43, wobei die annähernde Berechnung von Unter-Sinogrammen auch Gewichten enthält.
Vorrichtung nach Anspruch 40, wobei die genaue Berechnung von Unter-Sinogrammen ohne Verschieben durchgeführt wird.