DE102007020879A1

DE102007020879A1 - Verfahren und Vorrichtung für die äußerst schnelle Symmetrie- und SIMD- gestützte Projektion/Rückprojektion für die 3D-PET-Bildrekonstruktion

Info

Publication number: DE102007020879A1
Application number: DE102007020879A
Authority: DE
Inventors: Zang Hee Cho; Cheol Ok Lee; In Ki Hong; Young Don Son; Young Bo Seongnam Kim; Sung Tak Yongin Chung; Hang Keun Kim
Original assignee: Industry Academic Cooperation Foundation of Gachon University of Medicine and Science
Current assignee: Industry Academic Cooperation Foundation of Gachon University of Medicine and Science
Priority date: 2006-05-10
Filing date: 2007-05-04
Publication date: 2009-04-02
Also published as: US7456406B2; US20070278410A1; USRE42034E1

Abstract

Es werden ein Verfahren und eine Vorrichtung zum Rekonstruieren eines 3-D-Bildes geschaffen. Das Verfahren kann die folgenden Schritte enthalten: Erfassen mehrerer von einem Objekt ausgesendeter Antwortlinien (LORs); Umwandeln der mehreren LORs in erste Sinogrammdaten; Rückprojizieren der ersten Sinogrammdaten mit mehreren Projektionswinkeln, um Bilddaten für das Objekt zu erzeugen; und Projizieren der erzeugten Bilddaten mit den mehreren Projektionswinkeln, um die Bilddaten in zweite Sinogrammdaten umzuwandeln. Das Rückprojizieren kann das Füllen von Pixeln der Bildebene für jeden der mehreren Projektionswinkel mit den ersten Sinogrammdaten und das Drehen einer Koordinatenachse der Bildebene mit einem entsprechenden Projektionswinkel, um die Bilddaten zu erzeugen, enthalten. Das Projizieren kann das Drehen der Bilddaten mit einem entsprechenden Projektionswinkel in einer Gegenrichtung vor dem Projizieren der Bilddaten mit den mehreren Projektionswinkeln enthalten. Das Projizieren und das Rückprojizieren können Symmetrieeigenschaften im Koordinatenraum verwenden.

Description

Die vorliegende Anmeldung beansprucht die Priorität der Koreanischen Patentanmeldung Nr. 10-2006-0042155 , eingereicht am 10. Mai 2006, der Koreanischen Patentanmeldung Nr. 10-2007-0027305 , eingereicht am 20. März 2007, und der Koreanischen Patentanmeldung Nr. 10-2007-0040515 , eingereicht am 25. April 2007, deren gesamter Gegenstand hier durch Literaturhinweis eingefügt ist.
HINTERGRUND
1. Gebiet
Die vorliegende Erfindung bezieht sich auf ein Verfahren und auf eine Vorrichtung für die äußerst schnelle Symmetrie- und SIMD-gestützte Projektion/Rückprojektion für die 3D-PET-Bildrekonstruktion. Genauer bezieht sich die vorliegende Erfindung auf ein Verfahren und auf eine Vorrichtung, die auf den Symmetrieeigenschaften der Projektions- und der Rückprojektionsprozesse, insbesondere in dem 3D-OSEM-Algorithmus, der mehrere Projektionen und Rückprojektionen erfordert, beruhen.
2. Hintergrund
In den vergangenen Jahren hat es insbesondere auf den Gebieten der Hardware, der Software und der Computerimplementierung der Bildrekonstruktion einen bemerkenswerten Fortschritt der PET-Entwicklung gegeben. Jüngste Entwicklungen bei PET-Scannern (z. B. der von CTI (jetzt Siemens) entwickelte HRRT (High Resolution Research Tomograph – hochaufgelöste Forschungstomograph)) ermöglichen eine stark verbesserte Auflösung und Empfindlichkeit. In solchen PET-Scannern enthält die Menge gesammelter Koinzidenzliniendaten mehr als 4,5·10⁹ Koinzidenzantwortlinien, die durch nicht weniger als 120.000 Kerndetektoren erzeugt werden. Eine so große Datenmenge und die Rekonstruktion dieses Datensatzes stellen in dem HRRT ein echtes Problem dar. Das heißt, sie stellen schwere Probleme bei der Erzielung von Weiterentwicklungen und weiteren Anwendungen hochaufgelöster PET-Scanner dar. Somit erfordert das Erhalten eines Satzes rekonstruierter Bilder in diesen Scannertypen häufig viele Stunden Bildrekonstruktion. Zum Beispiel beträgt die Bildrekonstruktionszeit in dem HRRT mit voller Datenerhebung in normalen Gehirnabtastungen (unter Verwendung von SPAN 3) fast achtzig Minuten. Da die Bildrekonstruktionszeit viele Tage (nicht weniger als 43 Stunden für mehr als 32 Rahmen dynamische Bildkonstruktion) dauert, macht es dies praktisch unmöglich, auf der Grundlage der dynamischen Abbildung irgendeine Listenbetriebsart zu versuchen.
Im Allgemeinen können tomographische Bilder durch zwei Zugänge rekonstruiert werden, von denen einer ein analytisches Verfahren und der andere ein iterativer Zugang ist. Mitte der 1970-er Jahre wurden PET-Scanner verschiedener Typen entwickelt, wobei die Anwendung verschiedener tomographischer Bildrekonstruktionstechniken natürlich in das Gebiet eingeführt wurde. Im Fall eines analytischen Zugangs wie etwa der Rückprojektion und Filterung oder der gefilterten Rückprojektion (FB) wird häufig ein Artefakt erzeugt, der als ein Streifenartefakt bekannt ist. Dies trifft insbesondere dann zu, wenn die Detektoranordnungen wie etwa im Fall des HRRT (z. B. des hochaufgelösten Forschungstomographen von Siemens), bei dem die Detektoren in einem Satz von Blöcken in achteckiger Form angeordnet sind, nicht gleichförmig sind. Wegen der Zwischenräume zwischen den Blöcken verursachen diese Detektoranordnungstypen häufig fehlende Daten und führen im Fall der FB-Technik zu einem schweren Streifenartefakt. Somit wird nach alternativen Zugängen wie etwa nach einem EM-Algorithmus (Erwartungsmaximum) gesucht. Im Allgemeinen erfordert der EM-Zugang in dem Rekonstruktionsprozess mehrere Schritte, von denen die zwei Hauptschritte Folgende sind: Projektion (Vorwärtsprojektion) zum Erzeugen von Projektionsdaten von dem Bild oder Objekt und Rückprojektion in den Bildbereich für die Fertigbildrekonstruktion. In den EM-Algorithmen werden diese zwei Prozesse wiederholt, bis zufriedenstellende Bilder erhalten werden. Offensichtlich sind diese wiederholten Projektions- und Rückprojektionsprozesse zeitaufwändig und haben im Vergleich zum Algorithmus der reinen gefilterten Rückprojektion (FB-Algorithmus) den Hauptnachteil des EM-Zugangs. Außerdem nimmt im Fall der 3D-Bildrekonstruktion die Rechenlast wegen der astronomischen Zunahmen der Koinzidenzlinien oder der Antwortlinien (LOR) überproportional zu. Dies ist ein Haupthindernis im Tagesbetrieb hochaufgelöster PET-Scanner. Somit gibt es, insbesondere bei High-End-PET-Scannern wie etwa dem HRRT, einen starken Bedarf an der Verbesserung der Rechengeschwindigkeit oder der Rekonstruktionszeit in EM-Zugängen.
Projektionsverfahren nutzen üblicherweise eine Systemmatrix, die durch den Geometriefaktor des Scanners bestimmt ist. Während sich die Auflösung des PET-Bildes verbessert und die Anzahl der Schnitte erhöht, erhöht sich die Größe der Matrix proportional zu den Zunahmen der LORs ebenfalls drastisch, was nicht nur zur Notwendigkeit eines großen Speichers, sondern auch der Gesamtrechenzeit führt. Zum Beispiel erfordert die gegenwärtige HRRT außer der Erzeugung von Sinogrammen und geeigneten Daten-Streaming-Prozessen wie etwa Dämpfungs-, Zufalls- und Streukorrekturen einer Menge von Vorläufern zu dem Rekonstruktionsprozess nahezu achtzig Minuten Rekonstruktionszeit. Zur Behebung der Rechenlast der Bildkonstruktion sind eine Anzahl alternativer Vorschläge wie etwa die lineare Integration sowie die Verwendung mehrerer CPUs oder ein Cluster-Computersystem-Zugang vorgeschlagen worden. Da eine solche Clusterberechnung eine große Datenübertragungszeit erfordert, sind die meisten der Techniken allerdings nicht praktisch nutzbar, obgleich die Gesamtberechnung schneller als bei einer Einzeleinheit ist.
Das Erhalten oder Erzeugen von Projektionsdaten kann in zwei Kategorien, d. h. in ein strahlgesteuertes Verfahren und in ein pixelgesteuertes Verfahren, unterteilt werden. Das strahlgesteuerte Verfahren berechnet die lineare Integration direkt entlang des Strahlwegs, der die Mitten der zwei gegenüberliegenden Detektorzellen verbindet, während das pixelgesteuerte Verfahren die lineare Integration entlang des Strahlwegs berechnet, der für die gesamten Projektionswinkel um ein Bildpixel zentriert ist. In der Projektion wird häufig das strahlgesteuerte Verfahren verwendet, während in der Rückprojektion das pixelgesteuerte Verfahren verwendet wird.
In frühen Rekonstruktionstechniken wurde die Projektion dadurch erhalten, dass der durch die Pixelbereiche gehende Strahl unter der Annahme gewichtet wurde, dass der Strahlweg ein Streifen mit einer endlichen Breite ist.
Allerdings umfassen sie eine große Menge an Berechnung sowie die Speicherung einer großen Matrix- oder Datenanzahl. Shepp und Logan schlugen gleichzeitig einen einfachen und rechentechnisch effizienten Algorithmus vor, der die Berechnung der Länge des Strahlwegs, der sich mit jedem Pixel schneidet, (anstelle der Flächen) erfordert.
Um die Berechnung zu beschleunigen, gibt es eine Anzahl von Versuchen, die Rekonstruktionszeit zu verringern. Außerdem ist ein inkrementeller Algorithmus entwickelt worden, in dem die Symmetrieeigenschaft zwischen den Nachbarpixeln betrachtet wird, um die Position des Schnitts eines Strahls zu berechnen. Diese Idee wurde auf die 3D-Rekonstruktion in zylindrischer Geometrie unter Verwendung schiefer Strahlen erweitert. In 3D-Form wurde bei einem Mehrringsystem wie etwa dem HRRT sichtbar, dass echte 3D-Zugänge erforderlich sind, um die schiefen Strahlen vollständig zu nutzen und dadurch die Statistik des Bildes zu verbessern.
KURZBESCHREIBUNG DER ZEICHNUNG
Anordnungen und Ausführungsformen können ausführlich anhand der folgenden Zeichnung beschrieben werden, in der sich gleiche Bezugszeichen auf gleiche Elemente beziehen und in der:
1A ein 3D-Objekt und seine Projektion auf eine 2D-Projektionsebene veranschaulicht.
1B eine y'-z-Ebenen-Ansicht von 1A veranschaulicht.
1C ein Beispiel der Querebene bei θ = 0° und der Kurvenintegrale entlang der auf das x_r projizierten y'-Linie veranschaulicht.
2A eine Drehung des Projektionsstrahls (Strahlwegs) oder des Rahmens auf die Bildebene veranschaulicht, die auf der festen Referenzkoordinate (x, y)-Koordinate ist.
2B den Fall veranschaulicht, dass der Projektionsstrahlrahmen (Strahlwegrahmen) mit der festen (x', y')-Koordinate zusammenfällt.
3 unter Verwendung der y'-z-Ebenen-Ansicht die Beziehungen zwischen z, y, x, θ, x_r, y_r _'
und
veranschaulicht.
4A ein Beispiel der in dem vorgeschlagenen SSP-Verfahren verwendeten "Spiegelsymmetrie" veranschaulicht.
4B ein Beispiel der in dem vorgeschlagenen SSP-Verfahren verwendeten "⌀-Symmetrie" veranschaulicht.
5A die y'-Symmetrie in dem vorgeschlagenen SSP-Verfahren veranschaulicht.
5B die θ-Symmetrie in dem vorgeschlagenen SSP-Verfahren veranschaulicht.
6 das Konzept der ausgeglichenen Aufgabenverteilung auf der Grundlage der Summe der Strahlweglänge veranschaulicht.
7 ein Ablaufplan der Projektion gemäß einer Ausführungsform der vorliegenden Erfindung ist.
8 ein Ablaufplan der Rückprojektion gemäß einer Ausführungsform der vorliegenden Erfindung ist.
9A einen Vergleich der Projektionsdaten zwischen dem vorhandenen Verfahren und dem vorgeschlagenen SSP-Verfahren veranschaulicht.
9B Schnittansichten eines Sinogramms in einer spezifischen Ansicht veranschaulicht.
10A einen Vergleich einfacher Rückprojektionsbilder zwischen dem vorhandenen Verfahren und dem vorgeschlagenen SSP-Verfahren veranschaulicht.
10B Schnittansichten (Profile) auf einer x-Achse (y = 154, z = 103) veranschaulicht.
11A einen Satz mit dem vorhandenen Verfahren und mit dem vorgeschlagenen SSP-Verfahren rekonstruierter Bilder und die Differenzen veranschaulicht.
11B Schnittansichten (Profile) auf einer x-Achse (y = 154, z = 103) veranschaulicht.
AUSFÜHRLICHE BESCHREIBUNG
Anhand der beigefügten Zeichnung kann eine ausführliche Beschreibung gegeben werden. Der Durchschnittsfachmann auf dem Gebiet kann erkennen, dass die folgende Beschreibung nur veranschaulichend und in keiner Weise einschränkend ist. Weitere Ausführungsformen der vorliegenden Erfindung fallen dem Fachmann auf dem Gebiet unter Nutzung dieser Offenbarung leicht ein.
Übersicht über die Projektion, über die Rückprojektion und über Symmetrieeigenschaften
A. Übersicht über die Projektion und über die Rückprojektion in dem ausgerichteten (Referenz-)Rahmen mit gedrehter Projektionsebene
1A veranschaulicht ein 3D-Objekt und seine Projektion auf eine 2D-Projektionsebene. 1B ist eine y'-z-Ebenen-Ansicht von 1A. Es wird eine Beziehung zwischen den Projektionsebenen und dem Weg des Projektionsstrahls, d. h.
gegeben. Das Kurvenintegral wird entlang
ausgeführt. Der Winkel +θ gibt den schiefen Winkel der Bildebenen zu einer Oberseite an. x_r und y_r sind die Koordinaten der 2D-Projektionsebene des 3D-Objekts. 1C ist ein Beispiel der Querebene bei θ = 0° und der Kurvenintegrale entlang der auf x_r projizierten y'-Linie. ϕ gibt in dieser Figur die Drehung des in Bezug auf die Koordinatenachse (x, y) eingestellten Projektionsstrahls an. Die Strichlinien mit Pfeil geben die Projektionsstrahlen an. Jeder Projektionsstrahl ist durch vier Variable, d. h. (x_y, y_r, ϕ, θ), bestimmt. Bei einem gegebenen θ und ϕ besteht die in 1A gezeigte Projektionsebene aus 2D-Projektionsdaten (oder einem Nagelbett) in der Koordinate (x_y, y_r).
Wie in 1 gezeigt ist, hat der Projektionsstrahl in der tomographischen 3D-Bildverarbeitung 4 Dimensionen, d. h. x_r, y_r, ϕ und θ. Die Projektion ist ein Prozess des Umsetzens einer 3D-Objektfunktion oder von 3D-Bildinformationen in 3D-Koordinaten in 2D-Projektionskoordinaten und in die wie in 1 gezeigte Projektionsebene.
Bei einem festen Betrachtungswinkel ϕ und bei einem schiefen Winkel θ ist die Projektionsebene wie in 1A gezeigt bestimmt. In dieser Projektionsebene kann die Definition der Projektionsoperation in Form eines linearen Integrals ausgedrückt werden, das wie folgt gegeben ist:
mit
dem Strahlweg, xr= x' = xcosϕ + ysinϕ, yr = z – (–xsinϕ + ycosϕ)tanθ.
Gleichung (1) beschreibt, dass die Projektion P_ϕ,θ(x_r, y_r) eine Summe der Pixelwerte in einer Bildfunktion I(x, y, z) entlang des Wegs des Projektionsstrahls
in dem Bildbereich ist. δ(.) repräsentiert hier eine Abtastfunktion.
2 veranschaulicht die Konzepte des Drehprojektionsstrahlrahmens (Strahlwegrahmens) mit einer festen Bildebene in (x, y) und des vorgeschlagenen festen (ausgerichteten) Projektionsrahmens mit einer rotierenden Bildebene. 2A veranschaulicht die Drehung des Projektionsstrahls (Projektionsstrahlwegs) oder des Rahmens auf die Bildebene, die auf der festen (x, y)-Referenzkoordinate ist. Dies ist das herkömmliche Schema, das auf die meisten Bildrekonstruktionen angewendet wird. 2B veranschaulicht einen Fall, in dem der Projektionsstrahlrahmen (Projektionsstrahlwegrahmen) mit der festen (x', y')-Koordinate zusammenfällt. Das Letztere bedeutet, dass anstelle des Referenzrahmens, des Projektionsstrahlrahmens, nun die Bildebene gedreht wird. Dieses letztere Schema ist die Grundlage des vorgeschlagenen SSP-Verfahrens (Symmetrie- und SIMD-gestützten Projektions/Rückprojektions-Verfahrens). Dieses Schema ermöglicht, dass die Symmetrieeigenschaft der Bildebene genutzt wird, wodurch die gesamte Projektions- und Rückprojektionszeit verringert wird. Im Fall der Rückprojektion repräsentiert die Bildfunktion I(x', y', z) die Zwischenphase eines zu rekonstruierenden Bildes, wobei sie temporäre Bilddaten sind.
Es ist bekannt, dass eine Projektionsebene um einen Winkel ϕ gegenüber Referenzachsen (x, y) gedreht (oder ausgerichtet) werden kann. Ferner fällt in 2 x' mit x_r zusammen. Die gut bekannte Beziehung zwischen der gedrehten Koordinate (x', y', z) und der Bildkoordinate (x, y, z) in der Zylinderkoordinate ist durch Folgendes gegeben:
mit
R_ϕ der Drehmatrix.
Daraufhin kann Gleichung (1) als die gewichtete Summe entlang des Strahlwegs vereinfacht und wie folgt geschrieben werden:
mit I(x', y', z) = RϕI(x, y, z) x' = xr
In (3) ist angenommen, dass y' die Integrationsvariable ist.
Als ein Spezialfall, wenn der schiefe Winkel θ gleich null ist, sind die Projektionsstrahlen parallel zur y'-Achse (siehe 1). Da die volle Nutzung schiefer Strahlen wegen der erhöhten Statistik die Bildqualität verbessert, werden in der Mehrschicht- oder 3D-PET schiefe Strahlen (θ ≠ 0) gesammelt und für die echte 3D-Rekonstruktion verwendet. Die Erweiterung von (3) auf einen 3D-Fall kann durch θ ≠ 0 repräsentiert werden. Es wird angemerkt, dass die Koordinaten in der Querebene jetzt unabhängig von θ sind und dass der auf die Querebene projizierte Strahl parallel zu der y'-Achse ist. Durch die in 3 gezeigte trigonometrische Beziehung kann z wie folgt geschrieben werden: z = yr + y'tanθ (4)
3 veranschaulicht unter Verwendung der y'-z-Ebenen-Ansicht die Beziehungen zwischen z, y', x', θ, x_r, y_r,
und
. Es wird angemerkt, dass diese Beziehungen in allen y'-z-Ebenen mit irgendeinem x' gültig sind. Außerdem wird y' nun mit einem diskreten Wert n bezeichnet.
Durch Einsetzen von (4) in (3) kann die folgende diskrete Form der Projektionsdaten erhalten werden:
mit

∊ GANZE ZAHL,
∊ REELL, und 0 ≤
< 1 ist.
einem ganzzahligen Wert von

einem Interpolationskoeffizient oder Rest von

n einem diskreten Wert von y'
einer Integrationslänge von y' bei gegebenen y_r, ϕ und θ.
Es wird angemerkt, dass immer eine Interpolation erforderlich ist, da die Projektionsstrahlen und die Koordinaten (x, y, z) oder (x', y', z) nicht immer zusammenfallen. Im Fall der SSP-Implementierung ist, wie in (5) gezeigt ist, eine lineare 1D-Interpolation entlang der z-Achse notwendig. Es wird angemerkt, dass
und
die Interpolationskoeffizienten entlang der z-Achse sind.
Andererseits ist der Rückprojektionsprozess der Transponierte der Projektion. Somit kann auf die Rückprojektion oder Bilddatenerzeugung die gleiche Idee wie in dem Projektionsprozess (oder Projektionsdaten-Erzeugungsprozess) angewendet werden. Im Fall einer iterativen Rekonstruktion wie etwa des OS-EM-Algorithmus sind mehrere Projektionen und Rückprojektionen erforderlich. Das heißt, der Algorithmus erfordert viele wiederholte Bildrekonstruktionen sowie Projektionsdatenerzeugungen.
Nun kann die Rückprojektion von Projektionsdaten P_ϕθ(x_r, y_r) wie im Folgenden gezeigt ausgeführt werden:
Unter Verwendung der in (2) definierten Matrix kann (6) auch wie folgt geschrieben werden:
In dem ausgerichteten Rahmen ist x_r = x' definiert.
Zur Vereinfachung der Schreibweise wird I_ϕ als die Summe der Projektionsebenen um θ für ein festes ϕ definiert, wie sie im Folgenden als ein Zwischenphasenbild vor dem rekonstruierten Fertigbild gegeben ist. Die Rechnung kann dann z. B. unter Verwendung der Symmetrieeigenschaft in θ werter vereinfacht werden:
Unter nochmaliger Verwendung der in (2) gegebenen Drehmatrix kann (8) wie folgt geschrieben werden: Iϕ(x', y', z) = RϕIϕ(x, y, z) oder Iϕ(x, y, z) = R–ϕIϕ(x', y', z) (9)
Unter Verwendung von (9) kann (7) ferner wie folgt geschrieben werden:
Gleichung (10) ist äquivalent der Drehung der Menge teilweise rekonstruierter Bildebenen um einen gegebenen Winkel ϕ und daraufhin der Integration dieser teilweise rekonstruierten Bilder über alle ϕ, um das rekonstruierte Fertigbild I(x, y, z) zu bilden.
B. Symmetrieeigenschaften
Wie oben diskutiert wurde, ist der zentrale Teil des vorgeschlagenen SSP-Verfahrens die volle Nutzung der Symmetrieeigenschaften. Im Folgenden repräsentieren die Symmetriepunkte die ähnlichen Koeffizienten (den gleichen Wert, aber mit anderer Polarität und/oder mit anderen Koordinaten), die in der Rechnung ohne zusätzliche Berechnungen sowie unter Verwendung gleicher Anweisungen gemeinsam genutzt werden. Das heißt, die Symmetriepunkte haben den gleichen Interpolationskoeffizienten oder einen komplementären Wert. Der Anmelder hat festgestellt, dass es sechzehn Symmetriepunkte, d. h. eine "Spiegelsymmetrie", eine "ϕ-Symmetrie", eine "y'-Symmetrie" und eine "θ-Symmetrie", gibt, die praktisch nutzbar sind.
(a) Spiegelsymmetrie (+x' und –x', siehe 4A)
Einer der interessanten Aspekte der Bildrekonstruktion ist die Verwendung der Symmetrieeigenschaften jedes Bildpunkts. Das heißt, wenn ein Punkt mit einem geeigneten Interpolationskoeffizienten berechnet worden ist, kann in dem gegenüberliegenden Punkt bzw. in den gegenüberliegenden Punkten unter Verwendung verschiedener Typen von Symmetrieeigenschaften der gleiche Koeffizientenwert zugewiesen werden. Eine der einfachsten Formen ist die Verwendung der y'-Achsen-Symmetrie, wie sie in 4A gezeigt ist. In dieser Berechnung hat ein Punkt (x'₀, y'₀) bei +x_r nur mit einer Änderung der Polarität bei –x_r den gleichen Wert (x'₁, y'₁). Eine Veranschaulichung dieser "Spiegelsymmetrie" ist in 4A gezeigt. In dieser Veranschaulichung ist ein Beispiel für Symmetriepunkte wegen der "Spiegelsymmetrie", d. h., (x'₀, y'₀) → (x'₁, y'₁) = (–x'₀, y'₀), gegeben. Diese Symmetrieeigenschaft verringert die Rechenzeit um die Hälfte.
(b) ϕ-Symmetrie (ϕ und ϕ + 90°, siehe 4B)
Ähnlich dem Fall der "Spiegelsymmetrie" kann nachgewiesen werden, dass es eine Symmetrie zwischen den Symmetriepunkten bei ϕ und ϕ + 90° gibt. In diesem Fall ist (x'₀, y'₀) → (x'₁, y'₁) = (–y'₀, x'₀), d. h., es werden sowohl die Polarität als auch die Koordinaten geändert. 4B veranschaulicht diesen Fall einer ϕ-Symmetrie mit einem Zahlenbeispiel. Wenn, wie in der Figur gezeigt ist, ein Punkt bei ϕ berechnet wird und für diesen Punkt ein Interpolationskoeffizientenwert erhalten wird, z. B. bei (x'₀, y'₀), kann bei Koordinaten- und Polarisationsänderungen wie (–y'₀, x'₀) nach einer 90°-Drehung, d. h. bei ϕ + 90°, bei (x'₁, y'₁) der gleiche Koeffizientenwert verwendet werden. Diese Symmetrieeigenschaft ermöglicht wiederum, dass die Rechenzeit um die Hälfte verringert wird.
4 veranschaulicht die in dem vorgeschlagenen SSP-Verfahren verwendeten Beispiele der "Spiegelsymmetrie" und der "ϕ-Symmetrie". 4A veranschaulicht ein Beispiel der "Spiegelsymmetrie" gegen die y'-Achse. Das heißt, wenn ein Koordinatenpunkt (x'₀, y'₀) bekannt ist, kann ein symmetrischer Punkt (x'₁, y'₁) ohne zusätzliche Berechnung über die y'-Achse erhalten oder zugewiesen werden. Dies beschleunigt die Berechnung um einen Faktor zwei oder die Rechenzeit um die Hälfte. 4B veranschaulicht die Symmetrieeigenschaft der Projektionsdaten gegen ϕ, d. h., falls ein Interpolationskoeffizient eines Koordinatenpunkts (x'₀, y'₀) bei ϕ bekannt ist, kann bei ϕ + 90° ohne zusätzliche Berechnung ein symmetrischer Punkt (x'₁, y'₁) erhalten werden. Mit anderen Worten, die Projektionsdatenberechnung bei ϕ + 90° nutzt bei einfachen Koordinatenänderungen, d. h. (x', y') bei ϕ in (–y', x') bei ϕ + 90°, die gleichen Koordinatenwerte bei ϕ. Diese Symmetrieeigenschaft verringert die Rechenzeit um die Hälfte.
(c) y'-Symmetrie (+n und –n, siehe 5A)
Anders als die Spiegelsymmetrie und die ϕ-Symmetrie sind die y'-Symmetrie und die θ-Symmetrie in der y'-z-Ebene. Es wird angemerkt, dass die y'-Symmetrie eng mit der θ-Symmetrie verwandt ist. Somit ist es leicht, die zwei wie folgt zusammen zu veranschaulichen. Die y'-Symmetrie-Eigenschaft bewirkt, dass die Interpolationskoeffizienten in komplementäre Werte geändert werden. Tatsächlich hat die y'-Symmetrie, wie im Folgenden diskutiert wird, die gleiche Symmetrieeigenschaft wie die θ-Symmetrie. Zunächst ist die Eigenschaft der y'-Symmetrie in 5A veranschaulicht. Wie in 5A gezeigt ist, kann die Variable z nun wie folgt geschrieben werden (siehe auch 3):
wobei
und
die Interpolationskoeffizienten repräsentieren.
Unter Verwendung dieser Symmetrie können zwei Interpolationsoperationen gleichzeitig ausgeführt werden, wodurch die Anzahl der Schleifen in y'(= n) um die Hälfte verringert werden kann. Somit kann (5) wie folgt umgeschrieben werden:
Um die Anzahl der zur Ausführung der Berechnung erforderlichen Anweisungen zu verringern, kann Gleichung (13) in eine faktorisierte Form umgewandelt werden. Darüber hinaus kann der Interpolationskoeffizient
auf einfache Funktionen nur von n und θ reduziert werden und aus (11) abgeleitet werden. Da y_r immer eine ganze Zahl ist und
der Bedingung 0 =
< 1 genügt, ist
eine Funktion nur von n und θ. Somit kann die Schreibweise
zu r_n,θ vereinfacht werden (siehe und vergleiche mit den 3 und 5A). Ferner ist
gleich
. Ähnlich ist
gleich r_n,θ. Unter Verwendung der vereinfachten Formen kann (13) wie im Folgenden gezeigt in eine faktorisierte Form umgewandelt werden:
(d) θ-Symmetrie (+θ und –θ, siehe 5B)
Nun soll veranschaulicht werden, wie die θ-Symmetrie die Berechnung der Interpolationskoeffizienten mit der oben diskutierten y'-Symmetrie gemeinsam nutzt. Wie in 5B gezeigt ist, ändern nur die Änderungen in den Interpolationskoeffizienten ihre eigenen komplementären Werte, wenn sich die Vorzeichen von θ ändern.
Unter Verwendung von (4), (5), (11) und (12) kann gezeigt werden, dass z und r wie folgt gegeben sind:
Somit ist
Unter Verwendung der in (15) gegebenen Relation kann (14) wie folgt geschrieben werden:
Gleichung (14) ist für P_ϕ,θ(x_r, y_r), während (16) für das P_ϕ,–θ(x_r, y_r) ist. Für die Berechnung von (16) können die in (15) gegebenen Relationen, d. h. r_n,θ, r_n,–θ,
und
noch einmal verwendet werden, wobei ferner gleichzeitig P_ϕ,θ(x_r, y_r) und P_ϕ,–θ(x_r, y_r) berechnet werden können.
Wie oben gezeigt wurde, haben die y' und θ-Symmetrien die üblichen Eigenschaften gemeinsam. Somit wird die y'-Symmetrie verwendet, um die Anzahl der Schleifenzählwerte y' (= n) zu verringern, während die θ-Symmetrie verwendet wird, um die Anzahl der Schleifenzählwerte θ zu verringern.
Zusammengefasst werden für das SSP-Verfahren insgesamt sechzehn Symmetriepunkte, d. h. 2 (Spiegelsymmetrie), 2 (ϕ + 90°-Symmetrie), 2 (y'-Symmetrie) und 2 (θ-Symmetrie), verwendet. Somit führt die Multiplikation von diesen auf 16. Wie genau beschrieben wird, beschleunigen die SIMD zusätzlich zu diesen symmetriegestützten Verringerungen der Rechenzeit (16-fach) die Berechnung um einen weiteren Faktor 4 weiter. Somit kann mit dem neuen SSP-Verfahren ein Gesamtgeschwindigkeitsgewinn von 64 (16 × 4) erzielt werden.
5 veranschaulicht die y'- und/oder die θ-Symmetrie in dem vorgeschlagenen SSP-Verfahren. 5A zeigt die y'-Symmetrie, während 5B die θ-Symmetrie zeigt. Die linke Seite jeder Figur zeigt die z-y'-Ebenen-Schnittansicht, während der rechte Teil eine erweiterte Ansicht des linken zeigt. Jeder Quadrant veranschaulicht die Einzelheiten des Interpolationskoeffizienten r's. Diese y'- oder θ-Symmetrie-Eigenschaft gibt an, dass dann, wenn einer der Interpolationskoeffizienten für einen Punkt berechnet wird, die verbleibenden Interpolationskoeffizienten für die anderen drei Punkte ohne zusätzliche Berechnung erhalten werden können. Somit kann die Menge der Berechnungen unter Verwendung der y'- und der θ-Symmetrie auf ein Viertel (1/4) verringert werden.
Parallelisierung der Berechnung unter Verwendung von SIMD für das SSP-Verfahren
A. SIMD (Einzelanweisungs-Mehrfachdaten) und Symmetrieeigenschaften
Für praktische Zwecke ist es wichtig, die Anzahl der Anweisungen pro Schleife zu verringern, während die Anzahl der Schleifen innerhalb des Projektors/Rückprojektors verringert wird, um einen schnellen Algorithmus zu realisieren. Die Anzahl der Anweisungen pro Schleife wurde unter Verwendung der SIMD-Technik verringert. Die Verwendung der SIMD verringert die Anzahl der Anweisungen pro Schleife und verringert dadurch die Gesamtrechenzeit um einen Faktor vier. Um die Effizienz zu verbessern, sollte der Operand der SIMD sorgfältig bestimmt werden. Wie oben gezeigt wurde, wurde die Projektion/Rückprojektion von sechzehn Punkten gleichzeitig ausgeführt. Für diese Operation wurden sechzehn Symmetriepunkte in zwei Gruppen, eine für +θ und die andere für –θ, aufgeteilt. Da ein einzelner SIMD-Anweisungs-Datensatz aus vier Datensätzen besteht, war jede Gruppe mit dem gleichen (+θ oder –θ mit einem der zwei SIMD-Anweisungs-Datensätze gekoppelt.
Diese zwei Gruppen mit dem gleichen θ beruhen auf (14) und (16). TABELLE I zeigt diese sechzehn Paare der Symmetriebeziehung zwischen gedrehten Bilddatensätzen und Projektionsdatensätzen. Die obere Gruppe (a-1 und a-2) nutzt das gleiche +θ gemeinsam, während die untere Gruppe –θ gemeinsam nutzt. Jede Gruppe ist in zwei Untergruppen unterteilt. Jede obere Untergruppe nutzt die gleichen linearen Interpolationskoeffizienten gemeinsam, während jede untere Gruppe die komplementären Werte zu jenen der jeweiligen oberen Gruppe gemeinsam nutzt. Somit kann jede einzelne Untergruppe als ein einzelner SIMD-Datensatz (oder als ein einzelnes SIMD-Pack) behandelt werden.
und
nutzen den gleichen Wert gemeinsam.
Ähnlich nutzen
und
den gleichen Wert gemeinsam.
TABELLE I DIE SYMMETRIEBEZIEHUNG IN SIMD
Wie oben angemerkt wurde, sind sechzehn symmetrische Paare in zwei Gruppen unterteilt, die jeweils das gleiche θ gemeinsam nutzen. Außerdem wurden diese zwei Gruppen noch einmal in zwei Untergruppen unterteilt. Die Elemente der Untergruppen a-1 und b-2 nutzen den gleichen Interpolationskoeffizienten (z. B. r_y',θ) gemeinsam, während die Untergruppen a-2 und b-1 komplementäre Werte (z. B. 1 – r_y',θ) als Interpolationskoeffizienten haben. Wie im Folgenden gezeigt ist, repräsentiert jede Untergruppe eine gepackte SIMD-Dateneinheit:
wobei
SI einen SIMD-Bilddatensatz bedeutet.
SP einen SIMD-Projektionsdatensatz bedeutet.
Gleichung (17) zeigt die Anordnung des SIMD-Datensatzes für die obere Gruppe aus TABELLE I. Dieser Prozess wird als Packen bezeichnet, wobei ein einzelner SIMD-Anweisungs-Datensatz aus vier Punkten besteht. Dieser Prozess kann auf die unteren Gruppen mit –θ angewendet werden.
Wie oben angemerkt wurde, beruhen die obere Gruppe und die untere Gruppe aus TABELLE I auf (14) und (16). Die Gleichungen (14) und (16) können nun in wie in (17) gegebene SIMD-Versionen umgewandelt werden, wie sie im Folgenden gezeigt sind:
wobei
gemäß der Untergruppe "a-1" aus TABELLE I sind.
gemäß der Untergruppe "a-2" aus TABELLE I sind.
gemäß der Untergruppe "b-1" aus TABELLE I sind.
gemäß der Untergruppe "b-2" aus TABELLE I sind.
Dieses SIMD-Packungskonzept kann gleichfalls auf eine Rückwärtsprojektion angewendet werden.
Sicherheitshalber sei angemerkt, dass es innerhalb der obigen symmetrischen Paare mehrere Spezialfälle gibt. Diese Spezialfälle geben an, dass einige Symmetriepaare den gleichen Bilddrehpunkt gemeinsam nutzen, was zu Fehlern führt. TABELLE II zeigt diese Spezialfälle, die besondere Beachtung und Korrektur erfordern.
TABELLE II DIE SPEZIALFÄLLE OHNE SYMMETRIEGEGENSTÜCK
B. Optimierung der Aufgabenverteilung in einer Mehrprozessorumgebung
Die in letzter Zeit verfügbaren PCs stellen eine Mehrprozessorumgebung bereit. Um diese Parallelrechenfähigkeit zu nutzen, kann eine Anweisung verwendet werden, die "Thread-Programmiertechnik" genannt wird. Das vorgeschlagene SSP-Verfahren kann auf Mehrprozessorsysteme erweitert werden. In diesem Fall ist eine Aufgabenverteilung mit einer gleichen Lastmenge für die optimale Leistung wichtig.
Da das Verfahren gemäß der vorliegenden Erfindung auf der Annahme kreisförmiger oder zylindrischer FOV beruht, ist die Aufgabenlast so beschaffen, dass sie entlang der x'- oder x_r-Achse gleichverteilt ist. Wie in 6 veranschaulicht ist, liegt das daran, dass es eine ungleichmäßige Verteilung der Rechenlast gibt, falls jeder CPU ein gleicher Bereich von x' zugewiesen wird.
6 veranschaulicht das Konzept der ausgeglichenen Aufgabenverteilung auf der Grundlage der Summen der Strahlweglänge, d. h. der Gleichverteilung der Rechenlast für jede Gruppe S1, S2, S3 und S4. In diesem Fall sind die Aufgabenlasten in vier Gruppen mit gleicher Fläche, d. h. S1, S2, S3 und S4, geteilt. Die Umverteilung legt nahe, dass jede Gruppe fast die gleiche Fläche, d. h. die gleiche Menge an Rechenlast, hat.
Beschreibung des Gesamtverfahrens
Die Schemata des neuen SSP-Projektors und -Rückprojektors sind in 7 bzw. 8 gegeben. Für die SSP wurde unter Nutzung der Drehoperation sowie der Symmetrieeigenschaften zum Verringern der Anzahl der Schleifenzählwerte der ausgerichtete Rahmen verwendet. Außerdem wurde die Gesamtzahl der Anweisungen durch Kombinieren der Symmetrieeigenschaften und der SIMD-Technik verringert.
7 ist ein Ablaufplan, der die Projektion zeigt. Der innerste Block 706 innerhalb der Schleife enthält eine z-Schleife, die Interpolationsoperation, die Integration entlang des Strahlwegs 708 und die Abtastung entlang der x'-Achse 710. Nach der Projektion für den gesamten θ-Bereich sollten die Zwischenergebnisse der gedrehten Projektion entpackt werden, damit sie an die ursprüngliche Position zurückgespeichert werden. Es wird angemerkt, dass die Größen der (I-Pack-) und SIMD-Datensätze (I-Pack + P-Pack) kleiner als die Cache-Größen von L1 bzw. L2 sein müssen, um das Cache-Treffer-Verhältnis und die Datenzugriffsgeschwindigkeit zu maximieren.
8 ist ein Ablaufplan, der die Rückprojektion zeigt. Zunächst werden 8 Projektionsdaten mit symmetrischen Beziehungen zur Verwendung in die SIMD gepackt. Diese gepackten Projektionsdaten, P-Pack, werden entlang der y_r-Richtung linear interpoliert, um für alle θ-Richtungen die gepackten rückprojizierten Bilddaten, I-Pack, für das Halbfertigbild I_ϕ(x', y', z) zu berechnen. Als ein Zwischenschritt werden diese gepackten rückprojizierten Bilddaten nach jedem Schritt der θ-Schleife entpackt und auf jede gültige Position in den Halbfertigbilddaten I_ϕ(x', y', z) ausgerichtet. Diese Prozeduren werden über die Schleifen von y_r und x_r wiederholt, bis das Halbfertigbild I_ϕ(x', y', z) vollständig rekonstruiert worden ist. Das Halbfertigbild I_ϕ(x', y', z) wird für jeden Betrachtungswinkel ϕ rekonstruiert. Jedes Halbfertigbild I_ϕ(x', y', z) wird um seinen eigenen Betrachtungswinkel ϕ gedreht und daraufhin summiert, um das Fertigbild I(x, y, z), das Bild in den ursprünglichen kartesischen Koordinaten, zu rekonstruieren. Ähnlich wie bei der Projektionsoperation sollten die Größen der SIMD-Datensätze (I-Pack + P-Pack) und (I-Pack) für die Maximierung der Berechnung kleiner als die Cache-Größen von L2 bzw. L1 sein.
Wie in 7 und 8 gezeigt ist, wurde die Drehoperation verwendet, um den ausgerichteten Rahmen zu nutzen. Vor dem Projektionsschritt wird die Projektionsebene auf den ausgerichteten Rahmen (x', y') ausgerichtet und jede ankommende Bildebene daraufhin gemäß ihrem jeweiligen Betrachtungswinkel gedreht. Nach dem Rückprojektionsschritt wird das ursprüngliche Bild durch Zurückdrehen auf den ursprünglichen Betrachtungswinkel von dem ausgerichteten Rahmen wiedergewonnen, sodass es ermöglicht wird, rekonstruierte Zwischenbilddaten zu bilden.
Es wurde festgestellt, dass die effiziente Verwendung des Cache ein wichtiger Faktor in der Gesamtberechnung ist, um die Verarbeitungszeit zu verringern. Speicherzuweisung und Schleifenreihenfolgen sollten für das Speicherzugriffsmuster optimiert werden. Da für die Variable der äußersten Schleife x_r(= x') eingestellt wird, kann der Speicherschritt beim Lesen und Schreiben auf einen bestimmten Bereich innerhalb der L2-Cache-Größe beschränkt werden. Da sich die meisten der rechenaufwändigen Operationen in dem SSP-Verfahren auf die Interpolation von z oder y_r beziehen, sollte die Variable der innersten Schleife darüber hinaus z oder y_r sein. Somit ist es erwünscht, den Anfangsspeicher für die Rückprojektion über die Reihenfolge [x']:[y']:[z] zuzuordnen, während die Reihenfolge für die Projektion [ϕ]:[x_r]:[θ]:[y_r] sein sollte. Somit bestimmen diese Speicherzuordnungen die Schleifenreihenfolgen. Die Schleifenreihenfolgen des Projektors/Rückprojektors sind ebenfalls in 7 und 8 gezeigt. Im Ergebnis der Verwendung dieser Symmetrien kann die Gesamtzahl der Schleifenzählwerte in jedem Schritt um die Hälfte verringert werden.
Für die SIMD wird vor jedem Projektions/Rückprojektions-Schritt eine Speicherpackung (Schritt) ausgeführt. Nach dem Projektionsschritt wird der gepackte SIMD-Datensatz in dem Projektionsbereich entpackt und zu der ursprünglichen Position zurückgekehrt. Ähnlich wird nach jedem Rückprojektionsschritt der gepackte SIMD-Datensatz in dem Bildbereich entpackt.
Experimentelle Ergebnisse
A. Verfahren
Zur Bewertung der Wirksamkeit des SSP-Verfahrens wurden die Daten oder das Signogramm auf einen HRRT (von CPS/Siemens, Knoxville, Tennessee, USA, entwickelter hochaufgelöster Forschungstomograph) angewendet. Der HRRT besitzt 119.808 Detektoren und ist für die Gehirnabtastung mit äußerst hoher Auflösung ausgelegt. Da der HRRT unter den derzeit verfügbaren PET-Scannern für Menschen (bei der kleinsten Detektorgröße) die höchste Anzahl von Detektoren hat, besitzt er die höchste Rechenkomplexität für die Bildrekonstruktion. Die Anzahl der Berechnungen hat sich (proportional zur Anzahl der LORs) wesentlich erhöht und ist im Vergleich zu anderen vorhandenen PET-Scannern viel größer. Die durch das HRRT-Paket gelieferte Rekonstruktions-Software unterstützt iterative Algorithmen wie etwa OP-OSEM3D und ist in dem HRRT-System der bevorzugte Rekonstruktionsalgorithmus. Einer der Vorteile des OP-OSEM3D ist u. a. das Verhindern des Auftretens einer negativen Vorspannung, wenn die Zufallsrate hoch ist. Dieser Algorithmus wurde für den Leistungstest des neuen SSP-Verfahrens der Erfinder gewählt, da er die größte Speichergröße erfordert und der rechenaufwändigste Algorithmus ist. Sofern nichts anderes angegeben ist, wird somit als das vorhandene Verfahren auf die HRRT-OP-OSEM3D-Paket-Software Bezug genommen.
Wie zu sehen ist, wird das SSP-Verfahren durch kommerziell verfügbare Computersysteme wie etwa durch einen üblichen PC unterstützt. Für den Vergleich des SSP-Verfahrens gegenüber dem vorhandenen Verfahren wurden drei weitere Plattformen gewählt. Die erste Plattform ist ein PC1, der einen Intel Dual Core, 3,0 GHz, 4 GB RAM, 1 MB. L2-Cache pro CPU, besitzt. Die zweite Plattform, PC2, ist ein High-End-PC mit zwei Dual Core AMD 2,4 GHz, 8 GB RAM, 1 MB L2-Cache pro CPU. Die dritte Plattform ist die gegenwärtige HRRT-Computerplattform, die aus einem Acht-Knoten-Clustersystem besteht (dieses wird hier als HRRT-CPS-System bezeichnet), wobei jeder Knoten aus zwei Xeon 3,06 GHz-Prozessoren mit 512 kB-L2-Cache und 2 GB RAM besteht.
Das Konzept der "Spanne" bedeutet in 3D die Betriebsart der axialen Kompression. Der Test wurde in Spanne 3 und in Spanne 9 ausgeführt, die im HRRT üblicherweise verwendet werden. Die Parameter für OP-OSEM3D sind wie folgt: 256 × 256 × 207 Bildpixel; sechs Iterationen; und 16 Teilmengen. Außerdem wird die Thread-Programmierung verwendet.
B. Ergebnisse
Bevor die Rechengeschwindigkeit diskutiert wird, scheint es, dass die Genauigkeit der rekonstruierten Bilder (im Vergleich zu dem ursprünglichen von CPS/Siemens gelieferten HRRT-Paket) hervorgehoben werden muss. Zunächst wurden die Ergebnisse des Projektors und des Rückprojektors der SSP mit dem vorhandenen HRRT-Paket verglichen. Wie in den 9 und 10 gezeigt ist, wurde nur eine kleine Differenz ermittelt. Wie zuvor erwähnt wurde, ist die SSP-Technik oder das SSP-Verfahren grundsätzlich eine trilineare Interpolation, während das vorhandene Verfahren ein bilineares Verfahren ist. 9B zeigt die Profile für Projektionsdaten bei 45° und 90°. Derselbe Vergleich wurde an den Rückprojektionsdaten sowie an den rekonstruierten Fertigbildern ausgeführt, wobei die Ergebnisse in den 10 und 11 gezeigt sind. Wie darin gezeigt ist, sind die projizierten Daten und die Bildqualität der SSP äquivalent dem vorhandenen Verfahren, jedoch bei einem Vorteil von nahezu zwei Größenordnungen in Bezug auf die Rechengeschwindigkeit.
Außerdem wurde der Leistungsaspekt des SSP-Verfahrens (d. h. die Ausführungszeit auf den PC1 und 2) verglichen. Die Rechenzeit wurde für die zwei vorhandenen Kompressionsbetriebsarten, d. h. Spanne 3 und Spanne 9, gemessen. Die Ergebnisse geben an, dass das SSP-Verfahren die Leistung im Vergleich zu dem vorhandenen Verfahren (ursprünglichen HRRT-Paket) erwartungsgemäß tatsächlich fast um das 80-fache erhöhte, 2(ϕ + 90°-Symmetrie) × 2 (x_r- oder Spiegelsymmetrie) × 2(y'-Symmetrie) × 2(θ-Symmetrie) × 4(SIMD) + weitere Optimierungen. Es wird angemerkt, dass diese Leistungsverbesserung mehr oder weniger unabhängig von der Systemarchitektur erhalten werden kann. Wie das Verfahren nahe legt, ist die relative Verbesserung darüber hinaus, wie es aus den TABELLEN III und IV ersichtlich ist, umso höher, je schiefere Strahlen oder Winkel verwendet werden.
9A veranschaulicht einen Vergleich der Projektionsdaten zwischen dem vorhandenen Verfahren und dem vorgeschlagenen SSP-Verfahren. (Oben: vorhanden, Mitte: SSP, unten: die Differenz zwischen der SSP und dem vorhandenen Verfahren). 9B veranschaulicht die Schnittansichten eines Sinogramms in einer spezifischen Sicht. (Oben: Vergleich bei ϕ = 45°, Mitte: Vergleich bei ϕ = 90°, unten: die Differenzen zwischen der SSP und dem vorhandenen Verfahren bei ϕ = 45° bzw. ϕ = 90°).
10A veranschaulicht einen Vergleich einfacher Rückprojektionsbilder zwischen dem vorhandenen Verfahren und dem vorgeschlagenen SSP-Verfahren. (Oben: vorhandenes Verfahren. Mitte: neue SSP-Technik. Unten: die Differenzen.) 10B veranschaulicht die Schnittansichten (Profile) auf einer x-Achse (y = 154, z = 103). (Oben: Vergleich zwischen der vorhandenen und der neuen SSP. Unten: Die Differenzen zwischen den beiden.)

11 veranschaulicht einen Vergleich rekonstruierter Bilder mit sechs Iterationen und ihre Differenzen. 11A veranschaulicht einen Satz rekonstruierter Bilder mit dem vorhandenen Verfahren (oben), mit dem neuen SSP-Verfahren (Mitte) und die Differenzen (unten). 11B veranschaulicht die Schnittansichten (Profile) auf einer x-Achse (y = 154, z = 103). TABELLE III DIE VERGLEICHE DER LAUFZEIT DER VERSCHIEDENEN OPERATIONEN AUF DEM PC1

	vorhanden	neu	Verhältnis
Projektion (Spanne 9)	1539(s)	24(s)	64,1
Rückprojektion (Spanne 9)	1527(s)	25(s)	61,1
Projektion (Spanne (3)	4276(s)	51(s)	83,8
Rückprojektion (Spanne 3)	4255(s)	53(s)	80,3
Projektion + Rückprojektion (Spanne 9)	3066(s)	49(s)	62,6
Projektion + Rückprojektion (Spanne 3)	8531(s)	104(s)	82,0

PC1: Intel 3,0 GHz, 1 Dual-Core-CPU, 4 GB RAM, 1 MB L2-Cache pro CPU.

	vorhanden	neu	Verhältnis
Projektion (Spanne 9)	556(s)	9(s)	61,8
Rückprojektion (Spanne 9)	679(s)	11(s)	61,7
Projektion (Spanne (3)	1600(s)	22(s)	72,2
Rückprojektion (Spanne 3)	1767(s)	22(s)	80,3
Projektion + Rückprojektion (Spanne 9)	1235(s)	20(s)	61,7
Projektion + Rückprojektion (Spanne 3)	3608(s)	44(s)	82

PC2: AMD 2,4 GHz, 2 Dual-Core-CPU, 8 GB RAM, 1 MB L2-Cache pro CPU.

Die Verwendung schieferer Strahlen oder Winkel (z. B. Spanne 3) führt zu besserer axialer Auflösung und Statistik. Da die Rekonstruktion unter Verwendung des EM-Algorithmus aus einer großen Anzahl von Projektionen und Rückprojektionen besteht, ist diese Geschwindigkeitsverbesserung besonders wichtig. Schließlich wurde eine Vergleichsuntersuchung des HRRT-CPS-Systems mit dem ursprünglichen HRRT-CPS-Algorithmus gegenüber der 2-Dual-Core-CPU-Konfiguration (ohne Cluster), PC2, der Erfinder mit dem SSP-Verfahren verglichen und ein Geschwindigkeitsgewinn von einem Faktor 9 bis 11 erhalten. Der Grund für die Unterschiede zwischen der Berechnung (64-fach) und der tatsächlichen systemgestützten Berechnung (9–11-fach) liegt wahrscheinlich an dem ursprünglichen HRRT-CPS-System auf der Grundlage von 16 Xeon-CPUs. Das HRRT-CPS-System mit 16 Xeon-CPUs ist viel leistungsfähiger als der PC2 mit nur 2-Dual-Core-CPUs. Dieser Hardware-Unterschied könnte das Ergebnis dieses verringerten Geschwindigkeitsgewinns gewesen sein. Außerdem kann das SSP-Verfahren leicht auf die Cluster-Version erweitert werden.
Das vorgeschlagene SSP-Verfahren verbessert die Gesamtrechenzeit, insbesondere dann, wenn eine dynamische Funktionsuntersuchung erwünscht ist. Gegenwärtig erfordert die Bildrekonstruktionszeit für einen normalen HRRT-PET-Scan mit einem kommerziell verfügbaren PC (vom PC2-Typ) im Vergleich zu achtzig Minuten mit dem gegenwärtigen HRRT-Rekonstruktionspaket (OP-OSEM3D) unter Nutzung eines 8-Knoten-Clustersystems für eine Spanne 3 nach Abschluss der Sinogrammbildung und der Vorkorrekturprozesse etwa sieben bis acht Minuten.
Schlussfolgerung
Die vorliegende Erfindung schafft ein schnelles Verfahren zum Berechnen des Strahlwegs, das unter Verwendung der Symmetrieeigenschaften der Projektion und der Bilddaten zu den Projektionsdaten und zum rekonstruierten Bild in 3D führt. Unter Nutzung dieser einfachen geometrischen Symmetrieeigenschaften und mit Hilfe der SIMD wurde ein neues, äußerst schnelles Rekonstruktionsverfahren entwickelt, das im Vergleich zu dem vorhandenen Rekonstruktionspaket (insbesondere zu dem des OP-OSEM3D, dem neuesten und insbesondere für den HRRT ohne Gefährdung der Bildqualität am umfassendsten verwendeten PET-Verfahren) einen Rechengeschwindigkeitsvorteil von nahezu zwei Größenordnungen aufweist. Die in dem SSP-Verfahren eingeführten Schlüsselkonzepte können wie folgt zusammengefasst werden. Zunächst werden die Interpolationsoperationen unter Verwendung einer drehungsgestützten Projektion oder des Konzepts des ausgerichteten Rahmens von drei Dimensionen auf eine Dimension verringert. Zweitens wird die drehungsgestützte Projektion mit den symmetrischen Eigenschaften des (θ, ϕ, y', x_r) kombiniert und mit der SIMD-Technik bei Verwendung eines optimierten L1- und L2-Cache gekoppelt. Die sechzehn Symmetriepunkte, die die gleichen Interpolationskoeffizienten (oder ihre komplementären Werte) gemeinsam nutzen, wurden gruppiert, was somit ermöglicht, dass sie in einer Projektions/Rückprojektions-Operation gleichzeitig verarbeitet werden, um dadurch eine wesentlich verbesserte Gesamtrekonstruktionszeit zu erzielen. Um diese Symmetriepunkte in der Projektion/Rückprojektion gleichzeitig zu verarbeiten, sind außerdem SIMD-Operationen integriert worden. Insbesondere ermöglichte das SIMD-Schema, dass die Erfinder gleichzeitig auf vier Daten zugreifen. Außerdem wurden die für eine L2-Cache-Größe geeignete Datengröße pro Schleife sowie die Datenstrukturen optimiert, um den Speicherschritt zu minimieren.
Zusammengefasst wurden Projektions- und Rückprojektionsoperationen in dem ausgerichteten Rahmen unter Verwendung des Symmetriekonzepts und der SIMD mit Cache-Optimierung ausgeführt und wurde die Anzahl der Anweisungen in der Rekonstruktion eines Bildes aus den durch die HRRT-PET gelieferten Sinogrammdaten verringert. Mit anderen Worten, das vorgeschlagene SSP-Verfahren enthält maximal die Symmetrieeigenschaften der Projektion, um dadurch die Rechenzeit um nicht weniger als das 16-fache zu verringern. Zusammen mit der Verwendung des SIMD-Operators und der Cache-Optimierung wurde ein Gesamtrechengeschwindigkeitsgewinn um einen Faktor von nahezu 80 erzielt. Eine solche Verbesserung der Rechenzeit öffnet die Tür für die dynamische Funktionsuntersuchung der hochaufgelösten PET wie etwa der HRRT, die an einer inakzeptabel langen Rekonstruktionszeit litt und ein ernstes unerwünschtes Hindernis bei der Molekülabbildung ist.
Diese Verbesserung schafft für PET-Nutzer, insbesondere für HRRT-Nutzer in der Gemeinschaft der Molekülabbildung, zahlreiche neue Gelegenheiten. Es ist klar, dass das vorgeschlagene Verfahren den PET-Forschern dabei hilft, eine mögliche dynamische Untersuchung zu betrachten, die durch die durch den Rekonstruktionsprozess auferlegte inakzeptable Bildrekonstruktionszeit beschränkt war. Dieses neu entwickelte SSP-Verfahren kann ebenfalls auf Vorkorrekturprozesse wie etwa auf die Dämpfungs- und Streukorrektur angewendet werden. Schließlich ist es mit dem neuen SSP-Verfahren möglich, eine echte interaktive Abtastung zu erzielen. Außerdem kann das neue SSP-Verfahren für dynamische Untersuchungen auf große Klassen vorhandener PET-Scanner verschiedener Typen sowie auf die Clustersysteme angewendet werden.
Irgendeine Bezugnahme in dieser Beschreibung auf "eine Ausführungsform" auf eine "Ausführungsform", auf eine "beispielhafte Ausführungsform" usw. bedeutet, dass ein besonderes Merkmal, eine besondere Struktur oder eine besondere Eigenschaft, das/die in Verbindung mit der Ausführungsform beschrieben wird, wenigstens in einer Ausführungsform der Erfindung enthalten ist. Das Auftreten dieser Begriffe an verschiedenen Stellen in der Beschreibung bezieht sich nicht notwendig auf dieselbe Ausführungsform. Wenn in Verbindung mit irgendeiner Ausführungsform ein besonderes Merkmal, eine besondere Struktur oder eine besondere Eigenschaft beschrieben wird, wird ferner behauptet, dass es im Bereich des Fachmanns auf dem Gebiet liegt, dieses Merkmal, diese Struktur oder diese Eigenschaft in Verbindung mit anderen der Ausführungsformen zu bewirken.
Obgleich Ausführungsformen in Bezug auf eine Anzahl veranschaulichender Ausführungsformen davon beschrieben worden sind, können vom Fachmann auf dem Gebiet selbstverständlich zahlreiche weitere Änderungen und Ausführungsformen konstruiert werden, die im Erfindungsgedanken und im Umfang der Prinzipien dieser Offenbarung liegen. Insbesondere sind im Umfang der Offenbarung, der Zeichnung und der beigefügten Ansprüche verschiedene Änderungen und Abwandlungen an den Einzelteilen und/oder Anordnungen der jeweiligen Kombinationsanordnung möglich. Außer den Änderungen und Abwandlungen an den Einzelteilen und/oder an den Anordnungen sind für den Fachmann auf dem Gebiet außerdem alternative Verwendungen sichtbar.
ZITATE ENTHALTEN IN DER BESCHREIBUNG
Diese Liste der vom Anmelder aufgeführten Dokumente wurde automatisiert erzeugt und ist ausschließlich zur besseren Information des Lesers aufgenommen. Die Liste ist nicht Bestandteil der deutschen Patent- bzw. Gebrauchsmusteranmeldung. Das DPMA übernimmt keinerlei Haftung für etwaige Fehler oder Auslassungen.
Zitierte Patentliteratur

- KR 10-2006-0042155 [0001]
- KR 10-2007-0027305 [0001]
- KR 10-2007-0040515 [0001]

Claims

Verfahren zum Rekonstruieren eines 3D-Bildes, wobei das Verfahren umfasst: Erfassen mehrerer von einem Objekt ausgesendeter Antwortlinien (LORs); Umwandeln der mehreren LORs in erste Sinogrammdaten; Rückprojizieren der ersten Sinogrammdaten mit mehreren Projektionswinkeln, um Bilddaten für das Objekt zu erzeugen; und Projizieren der erzeugten Bilddaten mit den mehreren Projektionswinkeln, um die Bilddaten in zweite Sinogrammdaten umzuwandeln, wobei das Rückprojizieren das Füllen von Pixeln der Bildebene für jeden der mehreren Projektionswinkel mit den ersten Sinogrammdaten und das Drehen einer Koordinatenachse der Bildebene mit einem entsprechenden Projektionswinkel, um die Bilddaten zu erzeugen, enthält, wobei das Projizieren das Drehen der Bilddaten mit einem entsprechenden Projektionswinkel in einer Gegenrichtung vor dem Projizieren der Bilddaten mit den mehreren Projektionswinkeln enthält, und wobei das Projizieren und das Rückprojizieren Symmetrieeigenschaften im Koordinatenraum verwenden.
Verfahren nach Anspruch 1, bei dem die Symmetrieeigenschaften wenigstens eine von Symmetrieeigenschaften in den Betrachtungswinkeln (ϕ-Symmetrie), Symmetrieeigenschaften in den schiefen Winkeln (θ-Symmetrie) und Symmetrieeigenschaften in den Interpolationskoeffizienten entlang der Projektionslinien (y'-Symmetrie) enthalten.
Verfahren nach Anspruch 1, bei dem das Projizieren und das Rückprojizieren das Ausführen mehrerer Projektionsoperationen mit den Symmetrieeigenschaften gleichzeitig unter Verwendung von Einzelanweisungs-Mehrfachdaten (SIMD) enthalten.
Verfahren nach Anspruch 3, bei dem die Verwendung der SIMD das Anordnen eines Datensatzes mit dem gleichen schiefen Winkel in der Amplitude enthält.
Verfahren nach Anspruch 1, bei dem das Projizieren und das Rückprojizieren das Teilen der Bilddaten in mehrere gleiche Teilmengen und das Ausführen der Projektionsoperation für jede der mehreren Teilmengen bei jedem von mehreren Prozessoren enthalten.
Verfahren nach Anspruch 5, bei dem das Teilen der Bilddaten in mehrere gleiche Teilmengen das Verteilen von Aufgaben in der Weise, dass die Summen der Strahlweglängen gleich sind, enthält.
Verfahren nach Anspruch 2, bei dem das Projizieren und das Rückprojizieren das Ausführen mehrerer Projektionsoperationen mit den Symmetrieeigenschaften gleichzeitig unter Verwendung von Einzelanweisungs-Mehrfachdaten (SIMD) enthalten.
Verfahren nach Anspruch 7, bei dem die Verwendung der SIMD das Anordnen eines Datensatzes mit dem gleichen schiefen Winkel in der Amplitude enthält.
Verfahren nach Anspruch 2, bei dem das Projizieren und das Rückprojizieren das Teilen der Bilddaten in mehrere gleiche Teilmengen und das Ausführen der Projektionsoperation für jede der mehreren Teilmengen bei jedem von mehreren Prozessoren enthalten.
Verfahren nach Anspruch 9, bei dem das Teilen der Bilddaten in mehrere gleiche Teilmengen das Verteilen von Aufgaben in der Weise, dass die Summen der Strahlweglängen gleich sind, enthält.
Vorrichtung zum Rekonstruieren eines 3D-Bildes, wobei die Vorrichtung umfasst: Mittel zum Erfassen mehrerer von einem Objekt ausgesendeter Antwortlinien (LORs); Mittel zum Umwandeln der mehreren LORs in erste Sinogrammdaten; Mittel zum Rückprojizieren der ersten Sinogrammdaten mit mehreren Projektionswinkeln, um Bilddaten für das Objekt zu erzeugen; und Mittel zum Projizieren der erzeugten Bilddaten mit den mehreren Projektionswinkeln, um die Bilddaten in zweite Sinogrammdaten umzuwandeln, wobei das Rückprojizieren das Füllen von Pixeln der Bildebene für jeden der mehreren Projektionswinkel mit den ersten Sinogrammdaten und das Drehen einer Koordinatenachse der Bildebene mit einem entsprechenden Projektionswinkel, um die Bilddaten zu erzeugen, enthält, wobei das Projizieren das Drehen der Bilddaten mit einem entsprechenden Projektionswinkel in einer Gegenrichtung vor dem Projizieren der Bilddaten mit den mehreren Projektionswinkeln enthält, und wobei das Projizieren und das Rückprojizieren Symmetrieeigenschaften im Koordinatenraum verwenden.
Vorrichtung nach Anspruch 11, bei der die Symmetrieeigenschaften wenigstens eine von Symmetrieeigenschaften in den Betrachtungswinkeln (ϕ-Symmetrie), Symmetrieeigenschaften in den schiefen Winkeln (θ-Symmetrie) und Symmetrieeigenschaften in den Interpolationskoeffizienten entlang der Projektionslinien (y'-Symmetrie) enthalten.
Vorrichtung nach Anspruch 11, bei der das Projizieren und das Rückprojizieren das Ausführen mehrerer Projektionsoperationen mit den Symmetrieeigenschaften gleichzeitig unter Verwendung von Einzelanweisungs-Mehrfachdaten (SIMD) enthalten.
Vorrichtung nach Anspruch 13, bei der die Verwendung der SIMD das Anordnen eines Datensatzes mit dem gleichen schiefen Winkel in der Amplitude enthält.
Vorrichtung nach Anspruch 11, bei der das Projizieren und das Rückprojizieren das Teilen der Bilddaten in mehrere gleiche Teilmengen und das Ausführen der Projektionsoperation für jede der mehreren Teilmengen bei jedem von mehreren Prozessoren enthalten.
Vorrichtung nach Anspruch 15, bei der das Teilen der Bilddaten in mehrere gleiche Teilmengen das Verteilen von Aufgaben in der Weise, dass die Summen der Strahlweglängen gleich sind, enthält.
Vorrichtung nach Anspruch 12, bei der das Projizieren und das Rückprojizieren das Ausführen mehrerer Projektionsoperationen mit den Symmetrieeigenschaften gleichzeitig unter Verwendung von Einzelanweisungs-Mehrfachdaten (SIMD) enthalten.
Vorrichtung nach Anspruch 17, bei der die Verwendung der SIMD das Anordnen eines Datensatzes mit dem gleichen schiefen Winkel in der Amplitude enthält.
Vorrichtung nach Anspruch 12, bei der das Projizieren und das Rückprojizieren das Teilen der Bilddaten in mehrere gleiche Teilmengen und das Ausführen der Projektionsoperation für jede der mehreren Teilmengen bei jedem von mehreren Prozessoren enthalten.
Vorrichtung nach Anspruch 19, bei der das Teilen der Bilddaten in mehrere gleiche Teilmengen das Verteilen von Aufgaben in der Weise, dass die Summen der Strahlweglängen gleich sind, enthält.