DE102010029281A1

DE102010029281A1 - Verfahren und Bildrekonstruktionseinrichtung zur Rekonstruktion von Bilddaten

Info

Publication number: DE102010029281A1
Application number: DE201010029281
Authority: DE
Inventors: Dr. Popescu Stefan
Original assignee: Siemens AG
Current assignee: Siemens Healthcare GmbH
Priority date: 2010-05-25
Filing date: 2010-05-25
Publication date: 2011-12-01
Also published as: US20110293158A1; US8965091B2

Abstract

Die Erfindung betrifft ein Verfahren und eine Bildrekonstruktionseinrichtung (20) zur Rekonstruktion von Bilddaten (X) auf Basis von Messdaten (Y) eines bildgebenden Systems (1). Dabei werden aus den Messdaten (Y) zunächst rekonstruierte Start-Bilddaten (X0) in einem iterativen Optimierungsverfahren unter Nutzung eines Rauschregularisierungsterms und eines zusätzlichen Sparsity-Regularisierungsterm optimiert.

Description

Die Erfindung betrifft ein Verfahren sowie eine Bildrekonstruktionseinrichtung zur Rekonstruktion von Bilddaten auf Basis von Messdaten eines bildgebenden Systems. Bei dem bildgebenden System kann es sich um ein medizintechnisches System wie beispielsweise ein Magnetresonanz-System, ein Röntgen-Computertomographie-System, PET-System (PET = Positronen-Emissions-Tomographie), ein SPECT-System (SPECT = Single-Photon-Emissions-Computer-Tomographie) oder Ähnliches handeln. Je nach verwendetem bildgebenden System sind die Messdaten dementsprechend Magnetresonanz-Rohdaten, d. h. ortskodierte Magnetresonanzsignale, Röntgenprojektionsdatensätze, PET- oder SPECT-Datensätze etc. Darüber hinaus betrifft die Erfindung ein Verfahren zur Erzeugung von Bilddaten vom Inneren eines Objekts unter Verwendung eines solchen Rekonstruktionsverfahrens.
Zur Rekonstruktion von Bilddaten gibt es heutzutage grundsätzlich zwei Klassen von Verfahren.
Bei der ersten Verfahrensklasse handelt es sich um analytische Methoden wie z. B. bei Computertomographie-Projektionsdaten das sogenannte gefilterte Rückprojektionsverfahren (Filtered Back Projection; FBP) oder bei Magnetresonanz-Rohdaten um eine zweidimensionale inverse Fourier-Transformation. Hierbei handelt es sich jeweils um einschrittige Bildberechnungsprozeduren, die auf idealisierten mathematischen Systemmodellen basieren, welche die Ableitung von analytischen inversen Lösungen der zu berechnenden Gleichung erlauben. Geht man beispielsweise von einem Satz X von unbekannten Objekteigenschaften (z. B. örtliche Gewebeeigenschaften, welche letztlich durch das zu rekonstruierende Bild repräsentiert werden) aus, so werden diese Objekteigenschaften durch die Systemübertragungsfunktion A (die sog. Transferfunktion) des bildgebenden Systems auf einen Satz Y von Messwerten abgebildet: Y = A·X (1)
Die Datensätze X und Y sind dabei in Form von Matrizen geschrieben. In einem solchen idealisierten mathematischen System-Modell ist die Transferfunktion A und auch deren Inverse A^–1 bekannt, so dass sich umgekehrt gemäß X = A^–1·Y (2) aus den erhaltenen Messwerten Y der Satz X der unbekannten Objekteigenschaften, d. h. das zu rekonstruierende Bild, berechnen lässt.
Bei diesen Verfahren handelt es sich in der Regel um sehr zeitsparende Rekonstruktionsmethoden. Das Problem dieser Methoden besteht jedoch darin, dass hierbei die nicht-idealen Eigenschaften von realen bildgebenden Systemen nicht berücksichtigt werden. Dieses nicht-ideale Verhalten zeigt sich beispielsweise bei einem realen Magnetresonanz-Systemen in dem grundsätzlich vorhandenen Gauß'schen Rauschen, in inhomogenen Magnetfeldern, in Defekten hinsichtlich der Aufnahmeempfindlichkeit, Hochfrequenzinhomogenitäten etc. Bei CT-Scannern führen solche realen, bei der Rekonstruktion zu berücksichtigenden nicht-idealen Eigenschaften z. B. zum Quantenrauschen, zur Röntgenstreustrahlung, zu Strahlaufhärtungseffekte etc.
Aus diesem Grunde werden in einer zweiten Gruppe von Rekonstruktionsmethoden relativ zeitintensive rekursive Verfahren verwendet, die das Erfordernis einer vorher zu entwickelnden analytischen Lösung umgehen.
Eine solche iterative Rekonstruktionsmethode beginnt üblicherweise mit einer anfänglichen Schätzung von Start-Bilddaten X₀. Hierbei können beispielsweise mit einer klassischen analytischen Methode gemäß Gleichung (2) ermittelte Bilddaten verwendet werden. Ausgehend von diesen Start-Bilddaten X₀ werden dann iterativ die Bilddaten sukzessive gemäß folgender Gleichung verbessert: X_k+i = X_k – ΔX_k (3)
Das heißt, es werden in jedem Schritt k Bildkorrekturdaten ΔX_k berechnet und diese dazu verwendet, um aus den aktuellen Bilddaten X_k aktualisierte bessere Bilddaten X_k+1 zu erzeugen. Dabei geht es darum, Datenabweichungen zwischen den realen Messwerten Y und den synthetisierten Messwerten Y_k zu minimieren, die jeweils in dem laufenden Iterationsschritt gemäß Y_k = A·X_k (4) berechnet werden.
Um die Bildkorrekturdaten ΔX_k in jedem Iterationsschritt zu berechnen, können verschiedene Algorithmen verwendet werden, wie sie dem Fachmann aus dem Stand der Technik bekannt sind, beispielsweise der Gradient des steilsten Abstiegs (Deepest Descend Gradient), ein konjugierter Gradient (Conjugate Gradient) oder andere ähnliche Methoden. Z. B. können die Bildkorrekturdaten ΔX_k gemäß der Gleichung ΔX_k = α·grad(f_Z) (5) berechnet werden. Der Parameter α bezeichnet einen Relaxationsparameter, der die Geschwindigkeit der Konvergenz steuert. Es kann sich hierbei um einen beliebigen Wert, vorzugsweise im Bereich von 0 bis 1, beispielsweise um einen Wert 0,7, handeln. f_Z ist eine Zielfunktion wie z. B. f_z = ||A·X – Y||₂ (6)
Hierbei wird die Notation ||Z||₂ für die L₂-Norm von Z verwendet.
Das zuvor beschriebene Verfahren kann auch als ein durch Randbedingungen eingeschränktes, nicht lineares mathematisches Optimierungsproblem ausgedrückt werden gemäß: min||A·X – Y₂|| mit X > 0 (7)
Das heißt, es wird die Zielfunktion f_Z gemäß Gleichung (6) unter der Randbedingung minimiert, dass die Werte der einzelnen Bildpixel bzw. -voxel des Bilddatensatzes X alle größer als 0 sind.
Die Gleichungen (6) und (7) gelten für den Fall einer sog. „Least Square-Minimierung”. Grundsätzlich können aber auch andere Zielfunktionen verwendet werden, wie später noch er läutert wird. Die jeweils ideale Zielfunktion hängt u. a. von dem verwendeten bildgebenden System ab.
Der Hauptvorteil dieser iterativen Verfahren besteht darin, dass sie es erlauben, auch Transferfunktionen A für nichtlineare und nicht-ideale Messsysteme zu verwenden.
Um die iterativen Verfahren noch weiter zu verbessern und die realen Messbedingungen noch besser zu berücksichtigen, ist es möglich, innerhalb des Optimierungsverfahrens ein statistisches Modell für das die Messergebnisse beeinflussende Rauschen zu berücksichtigen. Eine Methode hierzu wird beispielsweise in Jeffrey A. Fessler, Statistical Methods for Image Reconstruction, John Hopkins University: „Short" Course, Mai 11, 2007, beschrieben.
Hierzu wird beispielsweise in die Zielfunktion f_Z im oben beschriebenen Optimierungsverfahren ein sogenannter Rauschregularisierungsterm eingesetzt, der mathematisch derart formuliert wird, dass höheres Bildrauschen unterdrückt wird. Sinnvollerweise wird dabei der Rauschregularisierungsterm so aufgebaut, dass er darauf abzielt, innerhalb der Optimierung den Bildkontrast, insbesondere Kanten und Bildschärfe, so gut wie möglich zu erhalten.
Seit einiger Zeit gibt es zudem ein neues Rekonstruktionsverfahren, welches auch als „Compressed Sensing”-Verfahren (CS-Verfahren) bekannt ist. Hierbei handelt es sich um eine Rekonstruktionsmethode, mit der recht gute Rekonstruktionen auch aus unterabgetasteten und/oder komprimierten Messdaten erzeugt werden können. Bei den unterabgetasteten Daten handelt es sich um Messdaten, die unterhalb der Nyquist-Grenze abgetastet wurden. Es hat sich herausgestellt, dass unter bestimmten Bedingungen aus unter der Nyquist-Grenze abgetasteten Messdaten sinnvolle Bilder erzeugt werden können, wenn die Unterabtastung nicht kohärent, sondern zufällig erfolgt. Beispielsweise können bei einem Computertomographen pro Umlauf zufällig gesteuert weniger Projektionen aufgenommen werden. Dies führt zu einer Reduzierung der Röntgenbelastung des Patienten. Bei einem Magnetresonanz-System können beispielsweise inkohärent weniger Linien im Fourier-Raum abgetastet werden, so dass auf diese Weise Messzeit und letztlich auch die Aufwärmung der Gewebe durch Hochfrequenzbelastung reduziert werden können.
Ein geeignetes Rekonstruktionsverfahren wird beispielsweise bei Michael Lustig et al., Compressed Sensing for Rapid MR Imaging, University of Stanford, beschrieben. Auch die Rekonstruktion solcher unterabgetasteten Messdaten erfolgt mit einer iterativen Methode, wobei jedoch die Erzeugung eines Bilds mit einer hohen „Sparsity” (Lückenhaftigkeit) forciert wird, weshalb das Verfahren auch als „Sparsity Enforcing Reconstruction” (SER) bezeichnet wird. Dabei wird versucht, das Bild aus den unterabgetasteten Daten so zu rekonstruieren, dass es nicht zu viele Details hat. In Abhängigkeit von der Untersuchungsfrage macht dies vielfach Sinn. So wird z. B. bei angiographischen Aufnahmen in der Regel nur darauf abgezielt, die mit Kontrastmittel gefüllten Blutgefäße möglichst gut zu sehen, wobei die Details um die Blutgefäße herum eher uninteressant und störend sind. In Bildern mit zu viel Bildrauschen bildet z. B. auch das Rauschen Details, welche jedoch keine Informationen für eine spätere Befundung bieten, sondern eher störend sind. Auch in solchen Fällen ist also ein Bild mit höherer Sparsity erwünscht.
Eine Sparsity Enforcing Reconstruction lässt sich beispielsweise dadurch realisieren, dass in dem oben genannten Optimierungsverfahren, beispielsweise innerhalb der Zielfunktion, ein sogenannter „Sparsity-Regularisierungsterm” genutzt wird. Verschiedene Möglichkeiten zur Bildung von Sparsity-Regularisierungstermen werden später noch genannt.
Ein grundsätzliches Problem des CS besteht aber darin, dass das allgemeine Bildrauschen erhöht wird. Dies liegt beispielsweise bei Magnetresonanz-Systemen daran, dass eine Reduzierung der Messzeit automatisch mit einem erhöhten Bildrauschen verbunden ist. Ebenso wird durch eine Reduktion der Gesamtstrahlungsdosis in Röntgen-Computertomographen und bei ähnlichen tomographischen Methoden wie SPECT oder PET, bedingt durch den bekannten Effekt der genannten „Photonen Starvation”, das Bildrauschen zwangsläufig erhöht. Durch geeignete Wahl eines Sparsity-Regularisierungsterms, beispielsweise durch Verwendung der Variation (Total Variation), kann zwar das Rauschen etwas reduziert werden. Diese Art der Rauschreduzierung ist jedoch ungünstigerweise mit einem Verlust an Bildschärfe verbunden, da es mit einem Sparsity-Regularisierungsterms nicht möglich ist, gleichzeitig die Bild-Sparsity zu maximieren und die Flanken im Bild möglichst gut zu erhalten. Es ist folglich immer zwischen einer Rauschreduktion und einer akzeptablen Bildschärfe abzuwägen.
Es ist daher eine Aufgabe der vorliegenden Erfindung, ein verbessertes Rekonstruktionsverfahren sowie eine Bildrekonstruktionseinrichtung zu schaffen, welche es ermöglichen, bei einer Verwendung von unterabgetasteten oder entsprechend komprimierten Messdaten in den rekonstruierten Bilddaten eine höhere Bildschärfe als mit dem bisherigen SER-Verfahren zu erreichen.
Diese Aufgabe wird durch ein Verfahren gemäß Patentanspruch 1 sowie durch eine Bildrekonstruktionseinrichtung gemäß Patentanspruch 13 gelöst.
Bei dem erfindungsgemäßen Rekonstruktionsverfahren werden zunächst aus den Messdaten Start-Bilddaten rekonstruiert, die dann in einem iterativen Optimierungsverfahren unter Nutzung eines Rauschregularisierungsterms und eines zusätzlichen, d. h. separaten Sparsity-Regularisierungsterms optimiert werden.
Im Rahmen umfangreicher Simulationsverfahren hat der Erfinder überraschenderweise herausgefunden, dass es möglich ist, innerhalb eines iterativen Optimierungsverfahrens nicht nur wie bisher einen Rauschregularisierungsterm oder alternativ einen Sparsity-Regularisierungsterm einzusetzen, um hiermit eine Sparsity Enforced Reconstruction auf unterabgetasteten Daten durchzuführen, sondern dass durch die Verwendung zweier separater Regularisierungsterme nebeneinander ein Synergieeffekt erreicht wird, und damit aus unterabgetasteten Messdaten Bilddaten mit einer hohen Bildschärfe erzeugt werden können, wie sie mit den bisher bekannten Verfahren nicht erreicht wird. Die Auswahl des Rauschregularisierungsterms und des Sparsity-Regularisierungsterms erfolgt dabei in Abhängigkeit von dem verwendeten bildgebenden System, der Art der Messdaten sowie ggf. der Zielvorgabe, wofür die rekonstruierten Bilddaten verwendet werden sollen.
Bei einem erfindungsgemäßen Verfahren zur Erzeugung von Bilddaten vom Inneren eines Objekts mittels eines medizintechnischen Systems werden zunächst mit Hilfe des bildgebenden Systems, beispielsweise einem MRT-Scanner, einem CT-Scanner, einem PET- oder SPECT-Scanner des bildgebenden Systems, Messdaten des Untersuchungsobjekts akquiriert, welche Bildinformationen des betreffenden Objekts repräsentieren. Darauf erfolgt eine Rekonstruktion der Bilddaten aus den erhaltenen Messdaten mit Hilfe des zuvor genannten erfindungsgemäßen Rekonstruktionsverfahrens.
Eine erfindungsgemäße Bildrekonstruktionseinrichtung zur Rekonstruktion von Bilddaten eines Objekts benötigt neben einer Messdatenschnittstelle zur Übernahme von mittels eines bildgebenden Systems gewonnenen Messdaten eine geeignete Rekonstruktionseinheit, welche ausgebildet ist, um Bilddaten basierend auf den Messdaten unter Verwendung eines iterativen Rekonstruktionsverfahrens zu rekonstruieren, wobei aus den Messdaten zunächst rekonstruierte Start-Bilddaten in einem iterativen Optimierungsverfahren unter Nutzung eines Rauschregularisierungsterms und eines zusätzlichen Sparsity-Regularisierungsterms optimiert werden. Darüber hinaus weist die Bildrekonstruktionseinrichtung eine Bilddatenschnittstelle zur Ausgabe der rekonstruierten Bilddaten auf.
Eine solche Rekonstruktionseinheit kann beispielsweise aus einer Start-Rekonstruktionseinheit sowie einer nachgeschalteten Optimierungseinheit bestehen, wobei die Start-Rekonstruktionseinheit beispielsweise aus den Messdaten zunächst mit einem analytischen Verfahren die Start-Bilddaten rekonstruiert und die nachgeschaltete Optimierungseinheit auf Basis dieser Start-Bilddaten die Optimierung durchführt.
Die Bildrekonstruktionseinheit kann im Übrigen selber Teil eines bildgebenden Systems sein, d. h. sie kann beispielsweise in üblicher Weise auf einem Steuer- und Auswerterechner des bildgebenden Systems installiert sein. Grundsätzlich kann eine solche Bildrekonstruktionseinrichtung aber auch in Form von bzw. auf einer anderen Rechnereinheit realisiert sein, die beispielsweise mit einem bildgebenden System über ein Netzwerk zur Datenübernahme verbunden ist oder in sonstiger Weise mit Messdaten versorgt werden kann.
Insbesondere kann die Rekonstruktionseinheit (beispielsweise die Startrekonstruktionseinheit und die Optimierungseinheit) in Form von Software auf einem geeigneten Rechner mit entsprechenden Speichermöglichkeiten realisiert sein. Die Messdaten-Schnittstelle sowie die Bilddaten-Schnittstelle können ebenfalls in Form von reiner Software realisiert sein, sofern nur eine Übernahme der Projektionsmessdaten bzw. eine Ausgabe der Bilddaten von anderen bzw. an andere auf der gleichen Rechnereinheit realisierten weiteren Messdaten-Vorverarbeitungseinheiten bzw. Bilddaten-Weiterverarbeitungseinheiten erforderlich ist. Grundsätzlich können diese Schnittstellen aber auch als kombinierte Hardware-/Software-Schnittstellen realisiert sein, um eine externe Ein- und Ausgabe zu realisieren, beispielsweise mit Hilfe von Softwarekomponenten speziell konfigurierte Hardware-Schnittstellen. Unter einer Ausgabe der Bilddaten ist dabei nicht nur eine externe Ausgabe auf einen Bildschirm, einen Drucker oder dergleichen zu verstehen, sondern jede Ausgabe der Bilddaten durch die Bildrekonstruktionseinrichtung, beispielsweise eine Hinterlegung der Bilddaten für eine spätere Sichtung oder Weiterverarbeitung in einen Speicher.
Eine weitgehend softwaremäßige Realisierung hat den Vorteil, dass auch schon bisher verwendete Bildrekonstruktionseinrichtungen auf einfache Weise durch ein Software-Update nachgerüstet werden können, um auf die erfindungsgemäße Weise zu arbeiten. Insofern wird die Aufgabe auch durch ein Computerprogrammprodukt gelöst, welches direkt in einen Speicher einer programmierbaren Bildrekonstruktionseinrichtung ladbar ist, mit Programmabschnitten, um alle Schritte des erfindungsgemäßen Verfahrens auszuführen, wenn das Programm in der Bildrekonstruktionseinrichtung ausgeführt wird. Erfreulicherweise ist die erfindungsgemäße Rekonstruktionsmethode gegenüber den bekannten iterative Rekonstruktionsverfahren nicht mit einem erhöhten Aufwand verbunden, sondern lässt sich ohne Weiteres auch auf bisher schon genutzten Bildrekonstruktionseinrichtungen implementieren.
Weitere vorteilhafte Ausgestaltungen und Weiterbildungen der Erfindung ergeben sich aus den weiteren abhängigen Ansprüchen sowie der nachfolgenden Beschreibung. Dabei kann die erfindungsgemäße Bildrekonstruktionseinrichtung auch analog zu den abhängigen Verfahrensansprüchen weitergebildet sein.
Wie bei den bisher bekannten Iterationsverfahren werden auch bei der erfindungsgemäßen Methode in den einzelnen Iterationsschritten unter Verwendung von zuvor aktuell erzeugten Bilddaten jeweils aktuelle Bildkorrekturdaten ermittelt, welche dann zur Berechnung von Bilddaten einer nachfolgenden Iterationsstufe genutzt werden. In der ersten Iterationsstufe sind dies beispielsweise die Start-Bilddaten, auf deren Basis dann Bildkorrekturdaten ermittelt werden, die – wie dies oben beispielsweise anhand von Gleichung (3) dargestellt ist – verwendet werden können, um aktualisierte Bilddaten zu erzeugen. In der nächsten Stufe werden dann diese aktualisierten Bilddaten genutzt, um neue Bildkorrekturdaten zu berechnen, die dann wiederum genutzt werden, um die Bilddaten erneut zu aktualisieren. Das Iterationsverfahren wird schließlich abgebrochen, wenn ein vorgegebenes Konvergenzkriterium erreicht ist, beispielsweise wenn die Bildkorrekturdaten kleiner als ein bestimmter ε-Wert sind, was damit gleichbedeutend ist, dass die Zielfunktion innerhalb des Optimierungsverfahrens ausreichend stark konvergiert ist. Ebenso kann sie nach einer maximalen Anzahl von Iterationsschritten abgebrochen werden.
In dem erfindungsgemäßen Iterationsverfahren kann jedoch in jeder Iteration die Erzeugung der aktuellen Bildkorrekturendaten unter Nutzung des Rauschregularisierungsterms und des zusätzlichen Sparsity-Regularisierungsterms erfolgen. Dabei können besonders bevorzugt der Rauschregularisierungsterm und der Sparsity-Regularisierungsterm einfach als zusätzliche Zielfunktionsterme, d. h. als Teile innerhalb einer im iterativen Verfahren zu optimierenden Zielfunktion, eingesetzt werden.
Die Zielfunktion und optional weitere Randbedingungen werden dabei vorzugsweise so gewählt, dass es sich um konvexe Funktionen handelt, die sicher gegen ein eindeutiges Minimum konvergieren und nicht weitere lokale Maxima oder Minima aufweisen. Es existieren spezielle Klassen von mathematischen Optimierungsproblemen, welche sowohl eine konvexe Zielfunktion als auch konvexe Randbedingungen haben, von denen bekannt ist, dass sie immer gegen ein eindeutiges Minimum konvergieren. Es können daher auch in dem erfindungsgemäßen Rekonstruktionsverfahren unter Berücksichtigung diesen Klassen geeignete Zielfunktionen und Randbedingungen gewählt werden.
Vorteilhafterweise werden der Rauschregularisierungsterm und der Sparsity-Regularisierungsterm jeweils mit einem Gewichtungsfaktor bei der Optimierung, beispielsweise innerhalb der Zielfunktion, gewichtet. Durch diese Gewichtungsfaktoren lässt sich einstellen, welchen Einfluss die einzelnen Regularisierungsterme jeweils auf die Optimierung haben. Beispielsweise lässt sich mit dem Gewichtungsfaktor für den Rauschregularisierungsterm das Verhältnis zwischen dem Rauschen und der räumlichen Auflösung variieren, wogegen der Gewichtungsfaktor für den Sparsity-Regularisierungsterm den Grad der möglichen Unterabtastung kontrolliert.
Besonders bevorzugt werden die Gewichtungsfaktoren des Rauschregularisierungsterms und des Sparsity-Regularisierungsterms so gewählt, dass die Summe der beiden Gewichtungsfaktoren zwischen 0 und 1 liegt, beispielsweise um 0,5.
Die Gewichtungsfaktoren können vorzugsweise interaktiv von einem Benutzer ausgewählt werden. Dadurch hat der Benutzer des bildgebenden Systems bzw. der Rekonstruktionseinheit die Möglichkeit, nacheinander mehrfache Rekonstruktionen mit beliebigen verschiedenen Gewichtungsfaktoren durchzuführen.
Die Messdaten werden bevorzugt mit einer bestimmten Unterabtastungsrate unterabgetastet. Alternativ zu einer solchen direkten Unterabtastung ist es auch möglich, unterabgetastete Daten durch eine nachträgliche Komprimierung der Messdaten zu erzeugen, beispielsweise indem innerhalb des Datensatzes inkohärent Messdaten verworfen werden. Das Verfahren eignet sich also auch dann für einen Einsatz, wenn die Messdaten mit einem hohen Bildrauschen verbunden sind und es sinnvoll wäre, diese Messdaten mit einem Sparsity Enforced Reconstruction-Verfahren zu bearbeiten, um so Bilder mit einer hohen Sparsity zu erreichen, in denen unnötige Bilddetails entfernt sind, die ohnehin keine Informationen liefern. Zudem eröffnet dies die Möglichkeit, Messdaten auch in Form von Rohdaten in komprimierter Form mit verringertem Speicherplatzbedarf abzuspeichern.
Unter der Unterabtastungsrate ist dabei der Grad der Unterabtastung im Verhältnis zu einer vollen Abtastung an der Nyquist-Grenze zu verstehen. Das heißt, eine Unterabtastrate von 1 bedeutet, dass eine volle Abtastung an der Nyquist-Grenze erfolgt ist, und eine Unterabtastung von 0,5 bedeutet, dass nur 50% der Daten im Verhältnis zu einer vollen Abtastung akquiriert wurden. Im Bereich der beschleunigten parallelen Bildgebungsverfahren in der MRT wird auch der Begriff Reduktionsfaktor oder Beschleunigungsfaktor R verwendet, wobei R = 2 eine Unterabtastung von 0,5 bedeutet. Das heißt, dass der Grad der Unterabtastung gleich 1/R ist. Besonders bevorzugt wird der Gewichtungsfaktor des Sparsity-Regularisierungsterms in Abhängigkeit von der Unterabtastungsrate gewählt. Dabei ist es bei einer ganz besonders bevorzugten Variante möglich, für die Gewichtungsfaktoren des Rauschregularisierungsterms und des Sparsity-Regularisierungsterms einen Summenwert vorzugeben, der – wie oben erläutert – vorzugsweise zwischen 0 und 1 liegt, und dann den Anteil des Gewichtungsfaktors des Sparsity-Regularisierungsterms am Summenwert in Abhängigkeit von der Unterabtastungsrate zu wählen. Dies ermöglicht eine optimale Wahl der Gewichtungsfaktoren und somit des Einflusses des „normalen” Rauschregularisierungsterms und des Sparsity-Regularisierungsterms.
Wie bereits oben erläutert, werden vorzugsweise in jedem Iterationsschritt zunächst aus den aktuellen Bilddaten unter Nutzung des Rauschregularisierungsterms und des Sparsity-Regularisierungsterms aktuelle Bildkorrekturdaten ermittelt. Bei einer ersten Verfahrensvariante erfolgt dies in einem Schritt und es werden dann, basierend auf den aktuell erzeugten Bilddaten und den aktuellen Bildkorrekturdaten, z. B. durch eine Subtraktion wie in Gleichung (3), neue aktualisierte Bilddaten für den nächsten Iterationsschritt erzeugt. Das heißt, hier findet eine einstufige Bildaktualisierung unter gleichzeitiger Berücksichtigung des Rauschregularisierungsterms und des Sparsity-Regularisierungsterms statt.
Bei einer alternativen Verfahrensvariante werden in einem Iterationsschritt zunächst in einer ersten Stufe unter Verwendung von zuvor aktuell erzeugten Bilddaten unter Nutzung des Rauschregularisierungsterms oder des Sparsity-Regularisierungsterms aktuelle erste Bildkorrekturdaten ermittelt. Basierend auf diesen aktuell erzeugten Bilddaten und den aktuellen ersten Bildkorrekturdaten werden dann neue aktuelle Interims-Bilddaten erzeugt. Erst in einer zweiten Stufe des Iterationsschritts werden dann unter Verwendung der aktuellen Interims-Bildkorrekturdaten unter Nutzung des Sparsity-Regularisierungsterms oder des Rauschregularisierungsterms (je nachdem, welcher dieser Regularisierungsterme in der ersten Stufe schon genutzt wurde) aktuelle zweite Bildkorrekturdaten ermittelt. Basierend auf den aktuellen Interims-Bildkorrekturdaten und den aktuellen zweiten Bildkorrekturdaten werden dann neue aktuelle Bilddaten für den nächsten Iterationsschritt erzeugt. Bei diesem Verfahren erfolgt also eine sukzessive Bildaktualisierung, wobei in der ersten Stufe die Aktualisierung der Bildkorrekturdaten hinsichtlich eines der beiden Regularisierungsterme erfolgt und in der zweiten Stufe dann die Aktualisierung hinsichtlich des anderen Regularisierungsterms.
Bei der Anwendung des erfindungsgemäßen Verfahrens ist es auch möglich, die Messdaten in beliebiger Weise vorzubearbeiten, wie dies auch sonst bei den herkömmlichen Rekonstruktionsverfahren der Fall ist. Beispielsweise können röntgencomputertomographische Messdaten zunächst einer Strahlaufhärtungskorrektur unterzogen werden. Ebenso ist eine Nachbearbeitung der rekonstruierten Bilddaten auf beliebige Weise möglich.
Die Erfindung wird im Folgenden unter Hinweis auf die beigefügten Figuren anhand von Ausführungsbeispielen noch einmal näher erläutert. Es zeigen:
1 eine schematische Darstellung eines Ausführungsbeispiels eines medizintechnischen bildgebenden Systems mit einer Bildrekonstruktionseinrichtung,
2 ein Flussdiagramm für einen möglichen Ablauf eines erfindungsgemäßen iterativen Rekonstruktionsverfahrens gemäß einem ersten Ausführungsbeispiel,
3 ein Flussdiagramm für einen möglichen Ablauf eines erfindungsgemäßen iterativen Rekonstruktionsverfahrens gemäß einem zweiten Ausführungsbeispiel.
In 1 ist zunächst schematisch ein medizintechnisches bildgebendes System mit einer Bildrekonstruktionseinrichtung 20 dargestellt, wobei es sich hierbei um eine erfindungsgemäße Bildrekonstruktionseinrichtung 20 handelt, die – wie später noch erläutert wird – innerhalb eines Optimierungsverfahrens einen Rauschregularisierungsterm und zusätzlich einen Sparsity-Regularisierungsterm verwendet.
Bei dem bildgebenden System 1 handelt es sich beispielsweise um ein Magnetresonanzsystem mit einem üblichen MR-Scanner 2, in welchem in einem Untersuchungstunnel 3 eine Patientenliege 4 untergebracht ist, auf der ein Untersuchungsobjekt O, beispielsweise ein Patient oder Proband, positioniert werden kann. Das MR-System 1 steht hier aber nur als ein Beispiel für andere medizintechnische bildgebende Systeme. Handelt es sich beispielsweise um ein CT-System, so würde anstelle des MR-Scanners 2 ein CT-Scanner verwendet, in welchem an einer um ein Untersuchungsraum umlaufenden Gantry ein Detektorsystem mit einem Detektor und einer dem Detektor gegenüberliegenden Röntgenquelle umläuft. Ebenso wird dort auf einer Patientenliege innerhalb der Gantry ein Untersuchungsobjekt positioniert.
Angesteuert werden der Scanner 2 sowie die Patientenliege 4 durch eine Steuereinrichtung 10. Hierbei kann es sich um einen üblichen zentralen Steuerrechner 10 handeln. Dieser ist mit einer entsprechenden Steuerschnittstelle 14 sowie einer Messsteuereinheit 11 ausgestattet. Diese Messsteuereinheit 11 steuert den Scanner 2 gemäß einem vorgegebenen Messprotokoll P zur Durchführung einer Messung an.
Die Steuerschnittstelle 14 ist hier nur als ein Block dargestellt. Es ist aber klar, dass bei einem medizintechnischen System diese Schnittstelle 14 i. d. R. aus einer Vielzahl von einzelnen Schnittstellen besteht. Die Steuerschnittstelle 14 kann beispielsweise bei einem Magnetresonanzgerät Schnittstellen zur Ansteuerung des Gradientensystems, eine Hochfrequenzsendesystem zur Aussendung der Hochfrequenzpulse etc. umfassen. Bei einem CT-Scanner wäre eine Schnittstelle zur Ansteuerung der Röntgenquelle sowie der Motoren in der Gantry erforderlich. Weiterhin umfasst diese Steuerschnittstelle 14 auch eine Schnittstelle, um Steuersignale an die Patientenliege zu übergeben.
Die vom Scanner 2 akquirierten Messdaten Y werden über eine Bilddaten-Empfangsschnittstelle 15 vom Steuerrechner 10 übernommen. Auch diese Schnittstelle 15 kann aus mehreren Teilschnittstellen bestehen. Die Messdaten Y werden dann an die Bildrekonstruktionseinrichtung 20 übergeben, welche daraus die Bilddaten X rekonstruiert, die dann beispielsweise auf einem Bildschirm eines Terminals 5 vor Ort angezeigt und/oder in einem Speicher 16 hinterlegt und/oder über eine Netzwerkschnittstelle 13 an einen Datenbus 6 übergeben werden, über welchen das bildgebende System 1 mit anderen Einrichtungen innerhalb eines Netzwerks verbunden ist, um beispielsweise Messdaten oder fertig rekonstruierte Bilddaten auf großen Massenspeichern zu hinterlegen oder an Befundungsstationen oder dergleichen zu übersenden.
Die Ansteuerung des gesamten bildgebenden Systems 1 durch einen Bediener erfolgt über das Terminal 5, welches über eine Terminalschnittstelle 12 mit der zentraler Steuereinheit 10 verbunden ist. Mit Hilfe dieses Terminals 5 und einer dort realisierten Benutzerschnittstelle ist es beispielsweise möglich, dass der Bediener ein Messprotokoll P aus dem Speicher 16 auswählt, ggf. anpasst und dafür sorgt, dass die Messsteuereinheit 11 auf Basis dieses Protokolls entsprechende Steuersignale über die Steuerschnittstelle 14 an den Scanner 2 ausgibt, so dass die Messung durchgeführt werden kann.
Die vom Scanner 2 akquirierten Messdaten Y werden an die Bildrekonstruktionseinrichtung 20 über eine Messdaten-Schnittstelle 21 übergeben. Sie gelangen dann an eine Rekonstruktionseinheit 22, welche dafür sorgt, dass aus den Messdaten Y Bilddaten X rekonstruiert werden, die dann über eine Bilddaten-Schnittstelle 23 ausgegeben, z. B. auf dem Terminalmonitor dargestellt, im Speicher 16 hinterlegt und/oder über eine Netzwerk-Schnittstelle 13 und einen Datenbus 6 an andere Komponenten innerhalb eines Netzwerks versandt werden.
Die Rekonstruktionseinheit 22 umfasst hier eine Start-Rekonstruktionseinheit 24, welche zunächst auf Basis der Messdaten Y mit einem analytischen Verfahren Start-Bilddaten X₀ erzeugt. Diese Start-Bilddaten X₀ werden an eine Optimierungseinheit 25 übergeben, in der dann innerhalb eines iterativen Verfahrens die Bilddaten optimiert werden, wie dies nachfolgend noch erläutert wird. Die Optimierungseinheit 25 erhält u. a. verschiedene Rekonstruktionsparameter, z. B. die zu verwendende Zielfunktion, die Regularisierungsterme und/oder Gewichtungsfaktoren, wie sie später noch erläutert werden.
Die zentrale Steuereinheit 10 muss nicht notwendigerweise, wie hier dargestellt, als eine integrale Einheit aufgebaut sein, sondern kann auch aus vielen separaten Einheiten gebildet sein, welche in geeigneter Weise untereinander vernetzt sind. Eine Vielzahl der Komponenten kann auch in Form von Software auf geeigneten Mikrocontrollern realisiert sein. Dies gilt insbesondere für die Rekonstruktionseinheit 22 bzw. deren Untermodule, die Start-Rekonstruktionseinheit 24 und die Optimierungseinheit 25. Sie können ggf. gemeinsam mit weiteren Komponenten auf einem oder mehreren Rechnereinheiten (z. B. Microcontrollern), die zur Steuereinheit 10 gehören, realisiert sein.
Wie eingangs bereits erläutert, kann eine solche Rekonstruktionseinrichtung gemäß der Erfindung auch außerhalb einer solchen zentralen Ansteuereinheit zur Ansteuerung des Scanners existieren. In diesem Fall werden die Messdaten Y beispielsweise direkt über ein Netzwerk an die Bildrekonstruktionseinrichtung übergeben. Insbesondere kann die Bildrekonstruktionseinrichtung Teil des Terminals 5 sein, sofern dies eine geeignete Rechenkapazität aufweist.
Es wird auch darauf hingewiesen, dass die zentrale Steuereinrichtung 10 sowie der Scanner eine Vielzahl von weiteren Komponenten aufweisen kann, die üblicherweise ein Magnetresonanzsystem bzw. ein CT-, PET- oder SPECT-System aufweisen. All diese Komponenten sowie die grundsätzliche Funktionsweise derartiger bildgebender Systeme sind dem Fachmann aber bekannt und brauchen daher im Folgenden nicht näher erläutert zu werden.
Im Folgenden wird zunächst die mathematische Basis für die Durchführung der in der Rekonstruktionseinheit 22 bzw. der Optimierungseinheit 25 ablaufenden iterativen Rekonstruktionen erläutert:
Wie oben erläutert, werden von einer Start-Rekonstruktionseinheit 24 zunächst mit Hilfe eines analytischen Verfahrens, beispielsweise bei CT-Daten mit einer gefilterten Rückprojektion oder, wie in dem vorliegenden Ausführungsbeispiel, bei MR-Messdaten durch eine zweidimensionale inverse Fourier-Transformation, Start-Bilddaten X₀ aus den Messdaten Y rekonstruiert.
Diese gehen dann in das iterative Optimierungsverfahren ein. Bei der Iteration erfolgt in jeder Iterationsstufe k = 0, 1, 2, ... eine Aktualisierung der Bilddaten X_k. Dies kann mit der bereits oben genannten Gleichung (3) beschrieben werden.
Wie oben ebenfalls bereits erwähnt (vgl. Gleichung (5)), sind die Bildkorrekturdaten ΔX_k jeweils eine Funktion der Zielfunktion, beispielsweise der Gradient der Zielfunktion multipliziert mit einem Relaxationsparameter α. In der Optimierung soll diese Zielfunktion f_z minimiert werden, was oben an einem Beispiel durch Gleichung (7) beschrieben ist.
Erfolgt in der erfindungsgemäßen Weise eine Optimierung unter Nutzung eines Rauschregularisierungsterms sowie eines zusätzlichen separaten Sparsity-Regularisierungsterms, so lässt sich dies allgemein in der folgenden statistischen Formulierung schreiben:
Hierbei wird also die Zielfunktion f_z in der Form: f_Z ={L(A·X, Y) + β₁·R(X) + β₂·||ψ(X)||₁} (8) dargestellt. Diese Zielfunktion kann wie in Zusammenhang mit den Gleichungen (6) und (7) oben dargestellt im Rahmen des erfindungsgemäßen Optimierungsverfahrens minimiert werden. Bei dem Term L(A·X, Y) handelt es sich um einen statistisch gewichteten Datenfehlanpassungsterm. Dieser gleicht die Fehlerabweichung entsprechend ihrer Glaubwürdigkeit an die Messdaten an, d. h. es wird hierdurch berücksichtigt, wie verrauscht jeder reale Messdatenwert ist. Der Datenabweichungsterm L(A·X, Y), welcher ein statisches Rauschmodell enthält, kann als nicht negative logarithmische Wahrscheinlichkeit in einer von dem geeigneten mathematischen Rauschmodell abhängigen Form ausgedrückt werden.
Für ein Magnetresonanztomographie-System ist ein geeignetes Rauschmodell ein Gauß'sches Rauschen. In diesem Fall kann der Datenfehlanpassungsterm L wie folgt dargestellt werden: L(A·x, Y) = ||A·X – Y||₂ (9)
Dies entspricht der Zielfunktion gemäß Gleichung (6), nur dass dort weder ein Rauschregularisierungsterm noch ein Sparsity-Regularisierungsterm eingesetzt wurde. Alternativ zur L₂-Norm in Gleichung (9) kann es auch bei Magnetresonanztomographen zu Situationen kommen, beispielsweise bei paralleler Bildgebung, in denen es besser ist, anstelle der L₂-Norm eine gewichtete „Least-Square”-Form für die Zielfunktion zu verwenden, um so die Rauschstatistik besser modellieren zu können, indem die Sensitivitätsprofile der Empfangsspulen berücksichtigt werden.
Handelt es sich bei den Messdaten um CT-, PET- oder SPECT-Messdaten, sollte beispielsweise ein mathematisches Modell gewählt werden, welches das Poissont-Rauschen beschreibt. In diesem Fall lässt sich L wie folgt schreiben: L(A·X, Y) = (Y – A·X)^TW(Y – A·X) (10)
Hierbei ist W der Faltungskern, der die Faltung der Projektionsdaten beschreibt. Man kann W aber auch als Gewichtungskoeffizienten sehen, der hier die Tatsache widerspiegelt, dass das Poisson-Rauschen von der Signalamplitude abhängt.
Der Term R(X) in Gleichung (8) ist der Rauschregularisierungsterm. Hiermit werden unerwünschte Bildcharakteristiken wie z. B. Rauhigkeit oder Rauschen unterdrückt. Für R(X) können verschiedenste Terme verwendet werden, einschließlich quadratischer Funktionen, die das Rauschen reduzieren, aber leider auch die räumliche Auflösung reduzieren. Eine spezielle nicht-quadratische Form ist die häufig verwendete flankenerhaltende Regularisierung gemäß:
Hierbei ist
eine Potenzialfunktion, welche so ausgewählt ist, dass sie weniger schnell ansteigt als eine Parabel. Der Wert δ ist ein frei wählbarer Parameter, der die Halbwertsbreite (FWHM) der Kurve bestimmt. Die Werte x_j sind die einzelnen Pixel oder Voxel der Bilddaten-Matrix X (je nachdem, ob es sich um zweidimensionale oder dreidimensionale Bilddaten handelt) und N ist die Anzahl aller Pixel/Voxel im Bild. Dieses Verfahren wird beispielsweise bei Jeffrey A. Fessler, Statistical Methods for Image Reconstruction, (siehe Zitat oben) beschrieben.
Mit dem Parameter bzw. Gewichtungsfaktor β₁ wird die Regularisierungsstärke hinsichtlich des Rauschregularisierungsterms R(X) und damit auch die Balance zwischen Rauschen und räumlicher Auflösung kontrolliert.
Der Sparsity-Korrekturterm ||Ψ(X)||₁ in Gleichung (8) dient dazu, die Bild-Sparsity zu maximieren, indem die L₁-Norm der „Sparsity-Transformierten” Ψ(X) minimiert wird. Als Sparsity-Transformierte Ψ(x) können verschiedene Funktionen gewählt werden. Besonders bevorzugt handelt es sich hierbei um eine Funktion mit einer Diskreten Cosinus-Transformierten (DCT), mit Wavelets, mit Curvelets, finiten Differenzen, einem Gradientenoperator. Ganz besonders bevorzugt kann als Sparsity-Transformierte Ψ(x) die Variation TV(X) des Bildes X gemäß
genutzt werden. m, n, k sind hierbei die Positionen der Pixel bzw. Voxel in einer Raumrichtung eines dreidimensionalen kartesischen Koordinatensystems. Die Variation gemäß Gleichung (13) entspricht dem vektoriellen Abstand der Bildwerte bzw. der Intensitäten zwischen zwei in jeder Raumrichtung benachbarten Pixeln bzw. Voxeln. Ist beispielsweise ein Bild nur gleichmäßig grau, so ist die Variation TV(X) = 0. Ein anderes Extrembeispiel wäre ein Bild, dessen Voxel schachbrettartig mal weiß, mal schwarz verteilt sind. In diesem Fall wäre die Variation TV(X) maximal.
Durch die Verwendung beider Regularisierungsterme gemäß Gleichung (8) innerhalb der Zielfunktion wird überraschenderweise synergistisch erreicht, dass einerseits eine Rekonstruktion von unterabgetasteten Daten mit einem optimierten Compressed-Sensing-Rekonstruktionsverfahren möglich ist und andererseits dennoch das Bildrauschen unter möglichst guter Bewahrung des Bildkontrastes reduziert werden kann.
Vorzugsweise werden die einzelnen Regularisierungsterme R(X) und ||Ψ(X)||₁ wie in Gleichung (10) dargestellt, als lineare Kombination durch Aufsummieren mit zwei separaten Gewichtungsfaktoren β₁, β₂ gebildet. Diese Art der Kombination garantiert eine konvexe Form der Zielfunktion, sofern die einzelnen Regularisierungsterme R(X) und ||Ψ(X)||₁ und der Datenabweichungsterm L(A·X, Y) sowie die Randbedingungen ebenfalls konvex sind, wodurch auch die Konvergenz des iterativen Optimierungsprozesses zu einem eindeutigen Minimum gewährleistet werden kann (vgl. Stephen Boyd und Lieven Vandenberghe, Book an Convex Opimization, Cambridge University Press, 2004).
Dabei kann der Rauschregularisierungsterm R(X) frei gewählt werden, um die beste Rauschreduzierung mit gleichzeitiger Kantenerhaltung zu erreichen, wogegen die Sparsity-Transformierte Ψ(x) hiervon unabhängig gewählt werden kann, um eine maximale Unterabtastung zu erreichen, die wiederum von vorbekannten Informationen über die notwendige Bilddiversität, beispielsweise auf Basis von Informationen über die abzubildende Region bzw. das Organ abhängig sein kann.
Der Gewichtungsfaktor β₁ kann wie erwähnt die Rauschregularisierungsstärke steuern, wogegen der Gewichtungsfaktor β₂ unabhängig davon den Grad der zulässigen Unterabtastung steuert. Wird β₂ = 0 gewählt, so entspricht die Methode dem bekannten Verfahren der Rauschreduktion mit Kantenerhaltung, für den Fall β₁ = 0 wird eine klassische CS-Rekonstruktion für eine optimale Unterabtastung mit hoher Sparsity in den Bilddaten erreicht. Die beiden Gewichtungsfaktoren β₁, β₂ werden vorzugsweise so gewählt, dass 0 < β₁ + β₂ < 1 gilt, um eine Überregularisierung zu vermeiden, die den iterativen Prozess insgesamt zu stark verlangsamen würde. Zum Beispiel könnte der Summenwert β = β₁ + β₁ ungefähr zwischen 0,5 und 0,7 liegen.
Bei einem besonders bevorzugten Ausführungsbeispiel werden die Gewichtungsfaktoren gemäß β₁ = β·λ und β₂ = ß·(1 – λ) gewählt. Hierbei ist β der o. g. Summenwert und λ ist eine Unterabtastungsrate, so dass λ = 1 eine volle Abtastung an der Nyquist-Grenze bedeutet und λ = 0,5 einen Wert von 50% Unterabtastung relativ zur Nyquist-Grenze repräsentiert. Auf diese Weise kann z. B. der Anteil der Sparsity-Erzwingung innerhalb des Optimierungsverfahrens automatisch von der tatsächlichen Unterabtastungsrate abhängig gemacht werden. Ebenso kann der Benutzer aber auch unabhängig davon, ob tatsächlich unterabgetastete Daten vorliegen, durch geeignete Wahl der Gewichtungsfaktoren β₁, β₂ und/oder der Werte β bzw. λ eine Rekonstruktionsmethode unter Nutzung einer SER einstellen. In diesem Fall würde die Anwendung einer SER Methode eine Reduzierung von Rekonstruktionsartefakten wie Bilddetails, die durch Inkonsistenzen in den Eingangsdaten auftreten (z. B. durch Röntgenstreuung um Metallobjekte in CT-Aufnahmen oder Bewegungsartefakte in MR-Aufnahmen), verstärken.
Der Benutzer des Systems kann dabei bevorzugt interaktiv auf besonders einfache Weise die Gewichtungsfaktoren bestimmen, und nacheinander mehrfache Rekonstruktionen mit verschiedenen Gewichtungsfaktoren durchführen und die Ergebnisse betrachten, indem er beispielsweise mit einem Mauszeiger auf einer virtuellen Skala auf einer graphischen Bedieneroberfläche für λ einen Wert zwischen 0 und 1 einstellt. Der Benutzer kann so stufenlos beliebig zwischen in einem Extremfall vollständig rauschbeschränkt (full noise-constricted) und im anderen Extremfall vollständig sparsity-beschränkt (full sparsity-constricted) rekonstruierten Bildern navigieren, um so spezielle Aspekte im interessierenden Bereich deutlich zu machen.
Bei der Optimierung können im Übrigen neben der Randbedingung, dass keine negativen Bilddaten entstehen können. (d. h. X > 0), auch weitere Randbedingungen eingeführt werden, wie z. B. obere Grenzwerte für die Bildwerte, Strahlaufhärtungs-Randbedingungen, Randbedingungen auf Basis von vorherigen Bildern etc.
2 zeigt ein Flussdiagramm für den Ablauf des Optimierungsverfahrens gemäß einer besonders bevorzugten Methode, wobei in einem Iterationsschritt eine einstufige Aktualisierung des Bilds erfolgt.
Das Verfahren beginnt mit der Übernahme von Messdaten Y, aus denen im Schritt 2.I zunächst Start-Bilddaten, beispielsweise mit einem klassischen analytischen Verfahren (bei MR-Rohdaten z. B. mit einer zweidimensionalen inversen Fourier-Transformation), erzeugt werden. Diese so abgeschätzten Start-Bilddaten X₀ werden dann in dem nachfolgenden iterativen Optimierungsprozess verbessert, wobei i. d. R. mehrfach eine Iterationsstufe A durchlaufen wird.
Hierzu wird zunächst in einem Schritt 2.II die Sparsity-Transformierte Ψ(X_k) für die aktuellen Start-Bilddaten X₀ berechnet. Dies kann beispielsweise mit Hilfe der Gleichung (13) erfolgen, wenn als Sparsity-Transformierte Ψ(x) die Variation TV(X) verwendet wird.
In einem weiteren Schritt 2.III wird dann für die aktuellen Bilddaten, d. h. hier im ersten Iterationsschritt für die Start-Bilddaten X₀, der Rauschregularisierungsterm R(X_k) berechnet, was beispielsweise mit den Gleichungen (11) und (12) erfolgen kann.
Anschließend erfolgt dann im Schritt 2.IV die Berechnung der Bildkorrekturdaten ΔX_k, beispielsweise gemäß der Gleichung (5) in Verbindung mit Gleichung (8). Da die Bildkorrekturdaten ΔX_k somit eine Funktion des Sparsity-Regularisierungsterm ||Ψ(X)||₁ (bzw. der Sparsity-Transformierten Ψ(X_k)), des Rauschregularisierungsterms R(X_k) sowie den zugehörigen Gewichtungsfaktoren β₁ und β₂ ist, ist diese Berechnung in 2 vereinfacht als Funktion f(β₁, β₂, Ψ(X_k), R(X_k)) dargestellt. Dies bedeutet jedoch nicht, dass die Bildkorrekturdaten ΔX_k nicht auch noch von anderen Parametern bzw. Funktionen abhängen, beispielsweise von dem statistisch gewichteten Datenfehlanpassungsterm L(A·X, Y) gemäß Gleichung (9) oder (10) das ist f(Y, X_k, β₁, ß₂, Ψ(X_k), R(X_k).).
Im Schritt 2.V erfolgt dann das Update der Bilddaten, indem von den aktuellen Bilddaten X_k (im ersten Iterationsschritt den Start-Bilddaten X₀), die berechneten Bildkorrekturdaten ΔX_k gemäß Gleichung (3) abgezogen werden.
In einem nachfolgenden Schritt 2.VI wird dann ein Konvergenzkriterium überprüft, beispielsweise, ob die Zielfunktion oder die Bildkorrekturdaten unter einem vorgegebenen ε-Wert liegen. Ist dies nicht der Fall (Abzweig N), so erfolgt ein neuer Iterationsschritt mit den aktualisierten Bilddaten X_k+1 als Eingangsdaten. Andernfalls (Abzweig J) ist im Schritt 2.VII das Rekonstruktionsverfahren beendet und die aktuellen Bilddaten werden als endgültige Bilddaten X akzeptiert.
3 zeigt ein etwas abgewandeltes Verfahren, bei dem die Bildaktualisierung bezüglich des Sparsity-Regularisierungsterms ||Ψ(X)||₁ und bezüglich des Rauschregularisierungsterms R(X_k) sukzessive in zwei Teilstufen B₁, B₂ erfolgt, d. h. hier besteht ein Iterationsschritt B aus zwei Teilstufen B₁ und B₂.
Das Verfahren beginnt im Schritt 3.I (ebenso wie im Schritt 2.I nach dem Verfahren gemäß 2) zunächst damit, dass aus den Messdaten Y Start-Bilddaten X₀ rekonstruiert werden. Das heißt, die Start-Bilddaten X₀ werden mit einem einfachen Rekonstruktionsverfahren abgeschätzt.
Dann wird – anders als bei der Verfahrensvariante nach 2 – noch vor der Iterationsschleife im Schritt 3.II ein erster Rauschregularisierungsterm R(X_Ψ,–1) berechnet, wobei wie im Schritt 3.I dargestellt die Bilddaten X_Ψ,–1 den Start-Bilddaten X₀ gleichgesetzt sind.
Erst dann beginnt die erste Teilstufe B₁ der Iterationsschleife im Schritt 3.III mit der Berechnung der Sparsity-Transformierten Ψ(X_k) für die aktuellen Start-Bilddaten X_k (d. h. zunächst mit den Start-Bilddaten X₀), wobei dieser Schritt wieder analog zum Schritt 2.II im Verfahren gemäß 2 ausgebildet ist.
Anschließend werden im Schritt 3.IV zunächst erste Bildkorrekturdaten ΔX_Ψ,k ermittelt. Dabei wird im Prinzip die gleiche Funktion wie im Schritt 2.IV gemäß der Verfahrensvariante nach 2 genutzt, jedoch mit dem Unterschied, dass für die Berechnung der ersten Bildkorrekturdaten ΔX_Ψ,k lediglich der aktuelle Sparsity-Regularisierungsterm ||Ψ(X)||₁, wie er in Schritt 3.III berechnet wurde, jedoch der Rauschregularisierungsterm R(X_Ψ,k–1) aus dem vorhergehenden Iterationsschritt verwendet wird. Beim ersten Iterationsschritt ist dies der im Schritt 3.II extra vor Beginn der Iteration berechnete Rauschregularisierungsterm R(X_Ψ,–1).
Im Schritt 3.V wird dann analog zum Schritt 2.V eine Aktualisierung der Bilddaten vorgenommen, wobei hier jedoch aktuelle Interims-Bilddaten X_Ψ,k ermittelt werden, die dann innerhalb der zweiten Teilstufe B₂ des Iterationsschritts B noch einmal hinsichtlich des Rauschregularisierungsterms aktualisiert werden. Hierzu wird zunächst im Schritt 3.VI der aktuelle Rauschregularisierungsterm R(X_Ψ,k) auf Basis der aktuellen Interims-Bilddaten X_Ψ,k berechnet.
Danach erfolgt dann im Schritt 3.VII eine Berechnung von zweiten Bildkorrekturdaten ΔX_R,k, wobei auch hier die gleiche Funktion wie im Schritt 2.IV verwendet wird, jedoch nun die bereits in der ersten Stufe B₁ im Schritt 3.III berechnete Sparsity-Transformierte Ψ(X_k), sowie der in der zweiten Stufe B₂ im Schritt 3.VI berechnete aktuelle Rauschregularisierungsterm R(X_Ψ,k) verwendet wird.
Schließlich erfolgt im Schritt 3.VIII ein Update der aktuellen Interims-Bilddaten X_Ψ,k unter Nutzung der im Schritt 3.VII ermittelten zweiten Bildkorrekturdaten ΔX_R,k.
Im Schritt 3.IX wird dann analog zum Schritt 2.VI geprüft, ob die Zielfunktion und/oder die Bildkorrekturdaten (ggf. jeweils auch beide einzelnen Bildkorrekturdaten ΔX_Ψ,k, ΔX_R _, _k) unter einem vorgegebenen ε-Wert liegen. Ist dies nicht der Fall (Abzweig N), wird der nächste Iterationsschritt durchlaufen, wobei die Laufvariable k um 1 inkrementiert wird. Ist das Konvergenzkriterium erfüllt (Abzweig J), so ist das Optimierungsverfahren im Schritt 3.X beendet und die aktuellen Bilddaten werden als endgültigen Bilddaten X angenommen.
Das Verfahren und die Konstruktionseinrichtung wurden in erster Linie anhand einer Rekonstruktion von medizinischen Bilddaten erläutert. Die Erfindung ist jedoch nicht auf eine Anwendung im medizinischen Bereich beschränkt, sondern es können auch grundsätzlich magnetresoananztomographische, computertomographische, PET- oder SPECT-Bilder etc. für andere Zwecke, beispielsweise für Materialprüfungen oder dergleichen, erzeugt und bearbeitet werden.
Es wird abschließend noch einmal darauf hingewiesen, dass es sich bei den zuvor beschriebenen Verfahren und Vorrichtungen lediglich um bevorzugte Ausführungsbeispiele der Erfindung handelt und dass die Erfindung vom Fachmann variiert werden kann, ohne den Bereich der Erfindung zu verlassen, soweit er durch die Ansprüche vorgegeben ist. Es wird der Vollständigkeit halber auch darauf hingewiesen, dass die Verwendung der unbestimmten Artikel „ein” bzw. „eine” nicht ausschließt, dass die betreffenden Merkmale auch mehrfach vorhanden sein können. Ebenso schließt der Begriff „Einheit” bzw. „Modul” nicht aus, dass diese aus mehreren Komponenten besteht, die gegebenenfalls auch räumlich verteilt sein können.
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 Nicht-Patentliteratur

Jeffrey A. Fessler, Statistical Methods for Image Reconstruction, John Hopkins University: „Short” Course, Mai 11, 2007 [0015]
Michael Lustig et al., Compressed Sensing for Rapid MR Imaging, University of Stanford [0018]
Jeffrey A. Fessler, Statistical Methods for Image Reconstruction [0068]
Stephen Boyd und Lieven Vandenberghe, Book an Convex Opimization, Cambridge University Press, 2004 [0072]

Claims

Verfahren zur Rekonstruktion von Bilddaten (X) auf Basis von Messdaten (Y) eines bildgebenden Systems (1), bei dem aus den Messdaten (Y) zunächst rekonstruierte Start-Bilddaten (X₀) in einem iterativen Optimierungsverfahren unter Nutzung eines Rauschregularisierungsterms und eines zusätzlichen Sparsity-Regularisierungsterm optimiert werden.
Verfahren nach Anspruch 1, wobei in einzelnen Iterationsschritten unter Verwendung von zuvor aktuell erzeugten Bilddaten (X_k, X_Ψ,k) unter Nutzung des Rauschregularisierungsterms und des Sparsity-Regularisierungsterms aktuelle Bildkorrekturdaten (ΔX_k, ΔX_Ψ,k, ΔX_R,k) ermittelt werden.
Verfahren nach Anspruch 1 oder 2, wobei der Rauschregularisierungsterm und der Sparsity-Regularisierungsterm Teil einer innerhalb des iterativen Verfahrens zu optimierenden Zielfunktion sind.
Verfahren nach Anspruch 3, wobei die Zielfunktion und optionale weitere Randbedingungen konvexe Funktionen sind.
Verfahren nach einem der Ansprüche 1 bis 4, wobei der Rauschregularisierungsterm und der Sparsity-Regularisierungsterm jeweils mit einem Gewichtungsfaktor (β₁, β₂) bei der Optimierung gewichtet werden.
Verfahren nach einem der Ansprüche 5, wobei die Gewichtungsfaktoren des Rauschregularisierungsterms und des Sparsity-Regularisierungsterms so gewählt werden, dass die Summe der beiden Gewichtungsfaktoren (β₁, β₂) zwischen 0 und 1 liegt.
Verfahren nach einem der Ansprüche 1 bis 6, wobei die Messdaten (Y) mit einer bestimmten Unterabtastungsrate unterabgetastete Messdaten (Y) sind.
Verfahren nach einem der Ansprüche 5 bis 7, wobei die Gewichtungsfaktoren (β₁, β₂) interaktiv von einem Benutzer ausgewählt werden.
Verfahren nach einem der Ansprüche 5 bis 7, wobei der Gewichtungsfaktor (β₂) des Sparsity-Regularisierungsterms in Abhängigkeit von einer Unterabtastungsrate gewählt wird.
Verfahren nach Anspruch 9, wobei für die Gewichtungsfaktoren (β₁, β₂) des Rauschregularisierungsterms und des Sparsity-Regularisierungsterms ein Summenwert vorgegeben wird und der Anteil des Gewichtungsfaktors (β₂) des Sparsity-Regularisierungsterms am Summenwert in Abhängigkeit von der Unterabtastungsrate gewählt wird.
Verfahren nach einem der Ansprüche 1 bis 10, wobei in einem Iterationschritt (A) unter Verwendung von zuvor aktuell erzeugten Bilddaten (X_k) unter Nutzung des Rauschregularisierungsterms und des Sparsity-Regularisierungsterms aktuelle Bildkorrekturdaten (ΔX_k) ermittelt werden und basierend auf den aktuell erzeugten Bilddaten (X_k) und den aktuellen ersten Bildkorrekturdaten (ΔX_k) neue aktuelle Bilddaten (X_k+1) erzeugt werden.
Verfahren nach einem der Ansprüche 2 bis 10, wobei in einem Iterationschritt (B) in einer ersten Stufe (B₁) unter Verwendung von zuvor aktuell erzeugten Bilddaten (X_k) unter Nutzung des Rauschregularisierungsterms oder des Sparsity-Regularisierungsterms aktuelle erste Bildkorrekturdaten (ΔX_Ψ,k) ermittelt werden und basierend auf den aktuell erzeugten Bilddaten (X_k) und den aktuellen ersten Bildkorrekturdaten (ΔX_Ψ,k) neue aktuelle Interims-Bilddaten (X_Ψ,k) erzeugt werden und dann in einer zweiten Stufe (B₂) unter Verwendung der aktuellen Interims-Bildkorrekturdaten (ΔX_Ψ,k) unter Nutzung des Sparsity-Regularisierungsterms oder des Rauschregularisierungsterms aktuelle zweite Bildkorrekturdaten (ΔX_R,k) ermittelt werden und basierend auf den aktuellen Interims-Bildkorrekturdaten (ΔX_Ψ,k) und den aktuellen zweiten Bildkorrekturdaten (ΔX_R,k) neue aktuelle Bilddaten (X_k+1) erzeugt werden.
Verfahren zur Erzeugung von Bilddaten (X) vom Inneren eines Objekts (O) mittels eines medizintechnischen Systems (1) bei dem zunächst mit Hilfe des bildgebenden Systems (1) Messdaten (Y) des Objekts (O) akquiriert werden und nachfolgend mit einem Verfahrens nach einem der Ansprüche 1 bis 12 aus den erhaltenen Messdaten (Y) Bilddaten (X) rekonstruiert werden.
Bildrekonstruktionseinrichtung (20) zur Rekonstruktion von Bilddaten (X) eines Objekts (O) mit – einer Messdaten-Schnittstelle (21) zur Übernahme von mittels eines bildgebenden Systems (1) gewonnenen Messdaten, – mit einer Rekonstruktionseinheit (22), welche ausgebildet ist, um Bilddaten (X) basierend auf den Messdaten (Y) unter Verwendung eines iterativen Rekonstruktionsverfahrens zu rekonstruieren, wobei aus den Messdaten (Y) zunächst rekonstruierte Start-Bilddaten (X₀) in einem iterativen Optimierungsverfahren unter Nutzung eines Rauschregularisierungsterms und eines zusätzlichen Sparsity-Regularisierungsterm optimiert werden, – und einer Bilddaten-Schnittstelle (23) zur Ausgabe der rekonstruierten Bilddaten (X).
Bildgebenden Systems (1) mit einer Messdaten-Akquisitionseinheit (2), zur Akquisition von Messdaten (Y), welche Bildinformationen eines Objekts (O) repräsentieren, und mit einer Bildrekonstruktionseinrichtung (20) nach Anspruch 14.
Computerprogrammprodukt, welches direkt in einen Speicher einer programmierbaren Bildrekonstruktionseinrichtung (20) ladbar ist, mit Programmcodeabschnitten, um alle Schritte eines Verfahrens nach einem der Ansprüche 1 bis 12 auszuführen, wenn das Programm in der Bildrekonstruktionseinrichtung (20) ausgeführt wird.