COMPUTERTOMOGRAPHISCHE BILDREKONSTRUKTION

Die Erfindung bezieht sich auf das Gebiet der Computertomographie, insbesondere auf ein Verfahren, ein Datenverarbeitungssystem und ein Computerprogramm zur Bildrekonstruktion gemäss dem Oberbegriff der entsprechenden unabhängigen Patentansprüche.

STAND DER TECHNIK

Ein derartiges Verfahren ist beispielsweise aus WO 02/067201 bekannt. Zur Bildrekonstruktion wird ein statistischer Optimierungsalgorithmus verwendet, welcher das Spektrum der verwendeten polyenergetischen Strahlungsquelle explizite berücksichtigt. Es werden zwei statistische Verfahren unter Verwendung eines Poisson-Modells vorgestellt. In beiden Fällen wird die Optimierung mittels geordneter Untermengen ("ordered subsets") beschleunigt, wobei Summen über den Winkel-Index eines Sinogramms ersetzt werden durch eine Reihe von Summen über Wmkel-Untermengen des Sinogramms. Vor der iterativen Optimierung wird zuerst ein einzelnes Bild rekonstruiert. Man wird eine FBP (filtered back projection) bevorzugen, weil diese Art von Algorithmen die Rekonstruktionszeit verringert [7]. Dann wird eine Segmentierung dieses Bildes durchgeführt, d.h. eine Zuordnung jedes Bildpunktes zu einer bestimmten Materialkategorie vorgenommen. Diese Zuordnung bleibt über mehrere Iterationsschritte hinweg fest. Die Zuordnung kann angepasst werden, falls die Optimierung nicht wie gewünscht konvergiert. Die genannte WO 02/067201 wird hiermit in ihrer Gesamtheit durch Referenz

auf genommen, insbesondere der Grundlagenteil mit Bezug auf die Figuren 1, 2, 3, 8 und 9, sowie die Beschreibung der statistischen Bildrekonstruktion.

Die vorgestellten Verfahren sind auf medizinische Anwendungen ausgelegt. Deshalb wird davon ausgegangen, dass das untersuchte Objekt aus bis zu 90% Wasser besteht, und dass die Gewebedichte zwischen 0.9 und 1.9 variiert. Dies erlaubt vereinfachende Annahmen, welche aber nicht auf andere Anwendungen übertragen werden können.

Ferner verlangt das Verfahren, um den Rechenaufwand in Grenzen zu halten, dass die gesamten Bilddaten der Projektionsmatrix dauernd im Arbeitsspeicher eines Rechners zur Verfügung stehen.

Es zeigt sich ferner, dass dieses und ähnliche Verfahren, insbesondere bei der Verwendung im industriellen Bereich, Bildartefakte erzeugt, welche tatsächliche Bildmerkmale überdecken.

DARSTELLUNG DER ERFINDUNG

Die Terminologie ("Bildpunkt", "Pixel") und die Beispiele sind der Erklärung halber auf zweidimensionale Bilder ausgerichtet. Die Erfindung lässt sich jedoch vom Prinzip her ohne weiteres auf dreidimensionale Darstellungen ("Volumenelemente", "Voxel") anwenden.

Der Begriff "Bildelement" steht deshalb, je nach Zusammenhang, für Pixel oder für Voxel. Analog bezieht sich der Begriff "Schnittbild" wahlweise auf eine zweidimensionale oder eine räumliche Repräsentation von Dichtewerten.

Es ist Aufgabe der Erfindung, ein Verfahren, ein Datenverarbeitungssystem und ein Computerprogranini zur Bildrekonstruktion der eingangs genannten Art zu schaffen, welche die oben genannten Nachteile behebt.

Diese Aufgabe lösen ein Verfahren, ein Datenverarbeitungssystem und ein Computerprogramm zur Bildrekonstruktion mit den Merkmalen der entsprechenden unabhängigen Patentansprüche.

Das Verfahren zur Rekonstruktion eines Röntgenstrahltomogrammes zum Erzeugen eines Schnittbildes ausgehend von Roh-Bilddaten, welche durch eine polyenergetische Strahlungsquelle erzeugt worden sind, weist also die folgenden Schritte auf:

• Bereitstellen eines Anteilmodells, welches den Anteil von mindestens zwei Materialien eines Bildelementes in Funktion der Dichte des Bildelelementes bestimmt;

• Rekonstruktion eines ersten rekonstruierten Schnittbildes anhand der Roh- Bilddaten;

• Iteratives Berechnen von weiteren rekonstruierten Schnittbildern unter Verwendung eines statistischen Optimierungsalgorithmus mit Berücksichtigung des polyenergetischen Spektrums der Strahlungsquelle und unter Verwendung des Anteilmodells.

Die Roh-Bilddaten sind typischerweise als Sinogramme darstellbar. Zur Rekonstruktion eines ersten rekonstruierten Schnittbildes als Initialbild für die Iteration wird ein bekanntes Verfahren wie beispielsweise die gefilterte Rückprojektion (filtered back projection, FBP) verwendet. Wenn die Iteration ein vorgegebenes Gütekriterium oder eine maximale Anzahl von Iterationsschritten erreicht, wird das letzte rekonstruierte Schnittbild als endgültiges Bild oder Ergebnisbild ausgegeben.

Dank der Verwendung des Anteilmodells in der iterativen Optimierung wird insbesondere vor den iterativen Berechnungen keine Segmentierung durchgeführt. Es

besteht also keine festgelegte Zuordnung von Bildelementen zu einzelnen Materialien, die über mehrere Iterationen hinweg beibehalten wird, wie beispielsweise in der eingangs zitierten WO 02/067201. Stattdessen kann ein Bildelement variierende Anteile verschiedener Materialen aufweisen, und werden diese Anteile bei jeder Iteration anhand der Dichte des Bildelementes neu bestimmt.

Der Schritt "Bereitstellen eines Anteilmodells", weist vorzugsweise die folgenden Schritte auf, welche in der Regel automatisch d.h. ohne Benutzerinteraktion durchgeführt werden: • Rekonstruktion eines ersten Schnittbildes anhand der Roh-Bilddaten;

• Berechnen eines Dichtebildes aus dem Abschwächungsbild entsprechend der Roh-Bilddaten;

• Erstellen eines Dichtehistogramms des ersten Schnittbildes;

• Identifizieren von Extrema des Histogramms; und • Bestimmen des Anteilmodells als Funktion der Extrema.

Das Dichtehistogramm kann auch direkt aus dem Histogramm der Werte des Abschwächungsbildes bestimmt werden, ohne dass explizite ein Dichtebild berechnet wird. Das beschriebene Verfahren eignet sich insbesondere für technische Anwendungen, oder medizinische CT, beispielsweise bei der Untersuchung von nichtüberlappenden Objekten welche im wesentlichen homogen mit diskreten Dichtewerten sind. Das Verfahren eignet sich grundsätzlich aber auch für den medizinischen Bereich und für Objekte mit überlappenden respektive gemischten Anteilen. Beispielsweise kann das Verfahren im medizinischen Bereich vorteilhaft eingesetzt werden, um Gewebe mit metallischen oder anderen Fremdkörpern zu untersuchen.

In einer bevorzugten Ausführungsform der Erfindung werden die Werte einer

Systemmatrix zur Repräsentation des Einflusses eines Bildelementes (pixel oder voxel) j auf eine Messung i nicht in ihrer Gesamtheit gespeichert. Insbesondere

werden diese Werte nicht im Arbeitsspeicher eines Computers gespeichert. Stattdessen werden jeweils für die i Messungen korrespondierende Werte γ _t als Summe der Werte von /z,y über alle j Bildelemente berechnet und gespeichert. Diese Werte Y _J werden beim iterativen Berechnen von weiteren rekonstruierten Schnittbildern verwendet. Dabei wird jeweils der Wert von H^y ₁ anstelle von h* verwendet.

Die Anwendung ist damit auch mit grossen Datensätzen möglich, wie sie bei der Verwendung von Fächerstrahlen, Konusstrahlen, Helix-Scans (fan-beam, cone-beam, helical CT) auftreten. Bei solchen Datenmengen kann nicht, wie für herkömmliche Verfahren notwendig, eine Projektionsmatrix im Arbeitsspeicher gehalten werden.

In einer weiteren bevorzugten Ausführungsform der Erfindung werden für die Bildrekonstruktion a priori Informationen über die Zusammensetzung und die interne Struktur der untersuchten Objekte verwendet, die von anderen Messmethoden oder Informationsquellen stammen. Beispielsweise kann mit optischen Messungen der Oberfläche oder aus CAD-Modellen, oder aus Ultraschalluntersuchungen ein Startbüd anstelle des ersten segmentierten Bildes oder in Kombination mit diesem verwendet werden. Es zeigt sich, dass dadurch auch die Konvergenz des Iterationsverfahrens beschleunigt werden kann.

Eine Vorrichtung zur Rekonstruktion eines Röntgenstrahltomogrammes weist Mittel zur Ausführung des beschriebenen Verfahrens auf. Beispielsweise ist die Vorrichtung ein Datenverarbeitungssystem und weisen die Mittel Speichermittel mit darin gespeicherten Computerprogrammcodemitteln auf, welche ein Computerprogramm beschreiben, und Datenverarbeitungsmittel zur Ausführung des Computerprogramms, wobei die Ausfuhrung des Computerprogramms zur Durchführung des Verfahrens gemäss der Erfindung führt. Das Datenverarbeitungssystem kann in einem Tomographiegerät eingebaut sein, oder es kann ein separater Computer mit EuWAusgabegeräten, einem Speicher und einem Prozessor sein. Da

das Verfahren parallelisierbar ist, kann auch ein Datenverarbeitungssystem mit mehrem, parallel arbeitenden Prozessoren verwendet werden.

Ein Computerprogramm zur Rekonstruktion eines Röntgenstrahltomogrammes ist in einen internen Speicher einer digitalen Datenverarbeitungseinheit ladbar und weist

Computerprogrammcodemittel auf, welche, wenn sie in einer digitalen

Datenverarbeitungseinheit ausgeführt werden, diese zur Ausführung des erfindungsgemässen Verfahrens bringen. In einer bevorzugten Ausführungsform der

Erfindung weist ein Computerprogrammprodukt einen Datenträger, respektive ein computerlesbares Medium auf, auf welchem die Computerprogrammcodemittel gespeichert sind.

Weitere bevorzugte Ausführungsformen gehen aus den abhängigen Patentansprüchen hervor. Dabei sind Merkmale der Verfahrensansprüche sinngemäss mit den Vorrichtungsansprüchen kombinierbar und umgekehrt.

KURZE BESCHREIBUNG DER ZEICHNUNGEN hu folgenden wird der Erfmdungsgegenstand anhand von bevorzugten Ausführungsbeispielen, welche in den beiliegenden Zeichnungen dargestellt sind, näher erläutert. Es zeigen jeweils schematisch:

Figur 1 eine Anordnung beim Erzeugen eines Röntgenbildes; Figur 2 eine Anordnung beim Erzeugen eines Röntgentomogramms;

Figur 3 eine Abschwächung bei einer polychromatischen Strahlungsquelle;

Figur 4 einen kontinuierlichen übergang bei einer Segmentierung eines Bildes in verschiedene Materialien;

Figur 5 ein Schnittbild, wobei die verschiedenen Graustufen unterschiedliche Dichten der Bildelemente darstellen;

Figur 6 ein Histogramm dieses Schnittbildes über die Dichte;

Figur 7 aus dem Histogramm hergeleitete Anteilfunktionen für verschiedene

Materialien; und Figur 8 eine eindimensionale Funktion mit verschiedenen Surrogatfunktionen:

Die in den Zeichnungen verwendeten Bezugszeichen und deren Bedeutung sind in der Bezugszeichenliste zusammengefasst aufgelistet. Grundsätzlich sind in den Figuren gleiche Teile mit gleichen Bezugszeichen versehen.

WEGE ZUR AUSFüHRUNG DER ERFINDUNG

Die Terminologie ("Bildpunkt", "Pixel") und die Beispiele sind der Erklärung halber auf zweidimensionale Bilder ausgerichtet. Die Erfindung lässt sich jedoch vom Prinzip her ohne weiteres auf dreidimensionale Darstellungen ("Volumenelemente", "Voxel") anwenden.

Der Begriff "Bildelement" steht deshalb, je nach Zusammenhang, für Pixel oder für Voxel. Analog bezieht sich der Begriff "Schnittbild" wahlweise auf eine zweidimensionale oder eine räumliche Repräsentation von Dichtewerten.

Figur 1 zeigt eine Anordnung beim Erzeugen eines Röntgenbildes. Eine Strahlungsquelle 1 sendet, durch einen Pfeil angedeutet, Röntgenstrahlen durch ein zu untersuchendes Objekt 2 gegen einen Bildaufnehmer oder Sensor 3. Die innere Struktur des Objektes 2 wird auf dem resultierenden Bild sichtbar. Figur 2 zeigt eine Anordnung beim Erzeugen eines Röntgentomogramms. Dabei sind die Strahlungsquelle 1 und der Bildaufnehmer 3 bezüglich des Objektes 2 drehbar angeordnet. So wird eine Mehrzahl von Bildern unter verschiedenen Drehwinkeln aufgenommen. Weitere Anordnungen sind beispielsweise in der WO 02/067201 in den Figuren 1 bis 3 und der dazugehörigen Beschreibung gezeigt.

Es wird im Folgenden ein statistisches Bildrekonstruktionsverfahren für Röntgen- Computertomographie beschrieben Darin wird ein polyenergetisches Quellenstrahlenspektrum modelliert und in einer physikalischen Nachbildung der Messungen berücksichtigt. Bestimmte Teile der damit durchgeführten Optimierung sind beschrieben in Erdogan[2] und Elbakri [3].

Es wird für jedes Bildelement für einen gegebenen Dichtewert eine Zusammensetzung aus einer Menge von Referenzmaterialien berechnet, das heisst ein relativer Anteil jedes dieser Referenzmaterialien am Bildelement. Deshalb ist keine vorgängige benutzergeführte oder automatische Vorsegmentierung des Bildes erforderlich. Bei dieser Vorsegmentierung wird jedes Bildelement einem bestimmten Material oder einer Materialklasse wie beispielsweise Knochen, Luft, Gewebe, Metall, Legierung, ... zugeordnet. Der Zusammenhang zwischen der Dichte und der Zusammensetzung wird durch eine Berechnungsvorschrift dargestellt, die im Folgenden auch Anteilsmodell genannt wird. Das Anteilsmodell erlaubt, dass Bildelemente in übergangsbereichen mit einer Mischung von Referenzmaterialien modelliert werden. Bildelemente in homogenen Bildbereichen weisen einen überwiegenden Anteil eines bestimmten Materials auf und können als aus diesem Material bestehend kategorisiert werden.

Für jedes Büdelement wird der Abschwächungskoeffizient modelliert als Produkt der zu ermittelnden Dichte mit der gewichteten Summe von energieabhängigen und massenbezogenen Abschwächungskoeffizienten. Diese Summe wird über die im Bildelement vorhandenen Materialien gebildet. Gewichtλmgsfunktionen werden durch eine Analyse des Dichtehistogramms eines ersten rekonstruierten Bildes sowie aus vorbekannten Materialeigenschaften gebildet.

Falls keine a priori information über die Materialien vorliegt, kann das Verfahren auch mit einem angepassten Verdrängungsmodell ("adapted displacement model") arbeiten. Dieses "adapted displacement model" wird aus Messungen mit einer Vielzahl von üblichen Materialien ermittelt.

Für das polyenergetische Modell wird eine "penalized-likelihood" Funktion verwendet, und der Iterationsalgorithmus zur Bestimmung der Dichte der Büdelemente wird entsprechend dem Anteilmodell modifiziert. Im Folgenden wird die Erklärung der Einfachheit halber anhand von zweidimensionalen Beispielen vorgenommen. Die Erweiterung auf drei Dimensionen ist ohne weiteres möglich.

Polyenergetische Röntgenstrahl-Computertomographie

Die Messungen werden als unabhängig Poissonverteilte Zufallsvariablen modelliert (Erdogan [2]) welche durch Hintergrundrauschen gestört sind. Das Hintergrundrauschen wird vor allem durch Streuung erzeugt. Für die Messungen wird das folgende statistische Modell angenommen:

M _j ~)Poisson {i(£).exp - l..N

wobei μ(x, y, E) die unbekannte ortsabhängige und energieabhängige Abschwächung durch das Objekt am Ort x,y und für die Energie E ist. Das Integral im Exponenten wird über die Linie Li berechnet, und I(E) repräsentiert die kombinierte Energieabhängigkeit der Strahlungsquelle und der Detektorempfmdlichkeit. Die angenommene mittlere Streuung wird durch s& dargestellt. N ist die Anzahl einzelner Messstrahlen respektive Messpunkte. Die Funktion I(E) kann beispielsweise aus Kalibriermessungen mit verschiedenen Materialien bestimmt werden, oder durch eine Monte-Carlo Simulation.

Figur 3 zeigt eine Abschwächung bei einer polychromatischen Strahlungsquelle. Die Wellenlänge ist entlang der Abszisse aufgetragen, die Intensität entlang der Ordinate. Die oberste Kurve 5 zeigt das Spektrum der Strahlungsquelle, d.h. die Energie in

Abhängigkeit der Wellenlänge. Beispielhaft ist eine vorgefilterte 450 keV-Quelle gezeigt. Die mittlere Kurve 6 zeigt das Spektrum nach 1 cm Aluminum, und die untere Kurve 7 nach 5 cm Aluminium.

Ein Objektmodell wird bei der polyenergetischen Modellierung wie folgt gebildet: Der Objektraum wird mit quadratischen Pixeln quantisiert (in einem dreidimensionalen Modell werden analog Voxel verwendet). Die Abschwächung wird damit berechnet als

P μ(x,y,E)= ∑b _j {x,y)μ _j {E)

wobei b _| (x,y) eine Basisfunktion für quadratische Pixel ist und μ/E) der unbekannte respektive gesuchte energieabhängige lineare Abschwächungskoeffizient des Materials im Pixel j ist. Die Basisfunktion gibt für Werte der Koordinaten x,y innerhalb des Pixels j den Wert 1 und alle anderen Werte von x,y den Wert Null. Für ein bestimmtes Material wird der Abschwächungskoeffizient bestimmt als Produkt des massebezogenen Abschwächungskoeffϊzienten mit der Dichte des Materials. Für das j-te pixel ist der Abschwächungskoeffizient gleich der gewichteten Summe der massebezogenen Abschwächungskoeffϊzienten der Referenzmaterialien, aus denen sich das Material im j-ten Pixel zusammensetzt, mal die Dichte des Pixels/ Voxels. Der Ausdruck für μj, den unbekannten energieabhängigen Abschwächungs- koeffizienten des Materials im j-ten Pixel ist:

M _j (E) = P _j M _n (E) = P _j ∑ω ^k (p _j )μ^(E) k=l wobei ω ^k der dimensionslose Anteil des k-ten Referenzmaterials an der Abschwächung im Pixel j ist, und μ _m ^k jeweils der bekannte energieabhängige massebezogene Abschwächungskoeffizient des £-ten Referenzmaterials ist. Die Gewichtungsfunktionen ω ^k können als Verallgemeinerung einer diskreten Segmentierungsoperation durch einen kontinuierliche Darstellung betrachtet werden.

AIs Illustration dieser Verallgemeinerung zeigt die Figur 4 einen kontinuierlichen übergang bei einer Segmentierung eines Bildes in verschiedene Materialien. Entlang der Abszisse ist eine beliebige Ortskoordinate eines Bildes aufgetragen. Entlang der Ordinate ist der Anteil eines bestimmten Materials am massenbezogenen Abschwächungskoeffizienten für jeweils einen Bildpunkt im Objekt aufgetragen, wobei der Anteil eines weiteren Materials diesen auf 100% ergänzt. Bei einer diskreten Segmentierung ist jeder Punkt entlang dieser Koordinate entweder ganz dem einen Material Mat#l oder ganz dem anderen Material Mat#2 zugeordnet. Bei der kontinuierlichen Darstellung jedoch werden Punkte im Bereich des übergangs zwischen den Materialien jeweils als aus beiden Materialien bestehend betrachtet.

Die ldassische Behandlung der Strahlaufhärtung geht davon aus, dass jedes Pixel einem bestimmten Material fest zugeordnet ist, das heisst dass eine Segmentierung a priori vorliegt. Dieser Ansatz ist erfolgreich, wenn zu Beginn ein Schnittbild guter Qualität vorliegt. üblicherweise wird dabei das Verfahren der gefilterten Rückprojektion (Filtered Backprojection, FBP) verwendet, da damit Bildanteile mit hoher (räumlicher) Frequenz gut rekonstruiert werden. Es können aber auch andere Verfahren verwendet werden, beispielsweise ein statistischer Algorithmus unter Annahme von monochromatischer Strahlung. Das derart rekonstruierte Schnittbild kann auch als erstes Schnittbild oder Initialbild für die iterative Optimierung des Bildes verwendet werden. Für verrauschte Bilder wird dieser Ansatz aber sehr ineffizient, weil bei der Segmentierung eine grosse Anzahl von Büdelementen falsch klassifiziert wird. Aus diesem Grunde wird hier die Segmentierung und eine systematische Kalibrierung ganz vermieden, wobei aber Information über das Strahlspektrum benötigt wird. Der Ansatz gleicht teilweise dem segmentationsfreien Verfahren von Elbakri[3], mit dem Unterschied, dass hier die Modellierung von einer Materialzusammensetzung durch eine endliche Anzahl von nichtüberlappenden Objekten ausgeht, und die Objekte im wesentlichen homogen sind und scharfe Kanten aufweisen.

Das Anteüsmodell wird wie folgt gebildet: Es wird von einem Bild ausgegangen, welches mit einer FBP oder einem statistischen Algorithmus unter Annahme einer monochromatischen Strahlungsquelle erzeugt wird. Typischerweise werden darin Strahlaufhärtungseffekte in der Bildmitte kleinere Abschwächungswerte erzeugen, und es werden streifenförmige oder strahlenförmige Artefakte ("streak artefacts") auftreten. Die Abschwächungskoeffizienten eines Materials im j-ten Pixel entsprechen einem Mittelwert über das Spektrum der Röntgenröhre, also

wobei p _j die Dichte im j-ten Pixel ist. Anhand dieser Beziehung kann aus dem gemessenen Bild der Abschwächungswerte μ ein Dichtebild berechnet werden.

Aus praktischen Gründen kann auch eine andere Variante zum Erstellen eines Dichtebildes aus dem Abschwächungsbild entsprechend der Roh-Bilddaten durchgeführt werden: Dazu wird zuerst ein Histogramm des Abschwächungsbildes erstellt. Jede Spitze des Histogramms von μ entspricht einem a priori bekannten

Material mit einer bekannten Dichte p. Somit ist eine Anzahl von Wertepaaren (μ,p) bekannt, und es kann ein Polynom, beispielsweise 2. oder 3. Ordnung bestimmt werden, das aus der Abschwächung μ die Dichte bestimmt. Mit diesem Polynom, oder allgemein mit einer Zuordnungsfunktion, wird aus dem gemessenen Bild der

Abschwächungswerte μ ein Dichtebild berechnet.

Im Folgenden wird davon ausgegangen, dass wie in industriellen Anwendungen typisch, ein Objekt nicht aus einer Vielzahl unterschiedlicher Gewebe besteht, sondern aus einer begrenzten Anzahl von Teilen aus unterschiedlichen Materialien. In einer Variante der Erfindung sind die Teile zudem mchtüberlappend. Dabei weisen die Materialien in der Regel extrem unterschiedliche Abschwächungskoeffizienten auf, die im Voraus bekannt sind. Auf dieser Basis wird

das Anteilsmodell aus dem Histogramm in der folgenden Weise automatisch gebildet:

Die Figur 5 zeigt ein Schnittbild, wobei die verschiedenen Graustufen unterschiedliche Dichten der Bildelemente darstellen. Figur 6 zeigt das daraus resultierende Dichtehistogramm. Jede der Spitzen im Histogramm entspricht einem bestimmten Material. Jede Spitze wird im Bereich um die Spitze herum durch ein

Gauss'sches Modell angenähert, wozu beispielsweise ein nichtlineares

Optimierungsverfahren wie der Levenberg-Marquardt-Algorithmus verwendet wird. Die Dichte entsprechend der Mitte der derart eingepassten Gauss'schen Verteilung

(pf ,σ% ) wird als gemessene Dichte pf zum Material k bezeichnet. Hohlräume entsprechen dem Material "Luft". Diesen Bildpunkten wird für die weitere Berechnung der Abschwächungskoeffizient von Luft und die (sehr kleine) Dichte zugeordnet.

Anschliessend wird berücksichtigt, dass die Dichtewerte in jedem Bildelelement aufgrund der Strahlaufhärtung kleiner sind als die tatsächliche Dichte pu des k-ten Referenzmaterials. Gemäss dem Anteilsmodell wird deshalb für Dichte werte p zwischen pi _c ^g und p _k der Materialanteil ω ^k des k-ten Materials auf 100% gesetzt. Zwischen diesem Bereich und dem entsprechenden Bereich des nächsten Materials k+1 mit der nächsthöheren Dichte geht der Anteil des k-ten Materials von 1 auf 0 und geht der Anteil des nächsten Materials von 0 auf 1. Die Summe der relativen Anteile aller Materialien ist immer 1, also

∑® ^k (pj)= l, 0 ≤ ω ^k ( _Pj )≤ l Jt=I In diesem übergangsbereich wird der Anteil der Materialien jeweils durch ein Polynom dritten Grades dargestellt, so dass erstens die obigen Bedingungen erfüllt sind und die Funktion beim übergang zu den Bereichen mit den Werten 0 respektive

1 stetig und zweimal differenzierbar ist. Das Anteilsmodell für das k-te Material kann im Bereich zwischen pf und p[ _+i somit ausgedrückt werden als

Im Bereich unter p _k ^S geht analog der Anteil von 1 auf 0 zurück und nimmt der Anteil des k-l-ten Materials von 0 auf 1 zu. Für das Beispiel der Figuren 5 und 6 sind die resultierenden Anteilfunktionen für verschiedene Materialien in der Figur 7 dargestellt. Die Dichte p ist entlang der Abszisse aufgetragen, die Anteile ω der verschiedenen Materialien Mat. #0 bis Mat. #3 entlang der Ordinate. Dabei bezeichnet Mat. #0 die Hohlräume respektive Luft.

In einer anderen Variante des Verfahrens wird der Anteil der Materialien jeweils durch zwei Polynome vierten Grades dargestellt, wobei die Bereiche zwischen dem Wert Null und 0.5 sowie zwischen 0.5 und 1 auf die beiden Polynome aufgeteilt sind. Dabei werden die Polynome so gewählt, dass beim Wert 0.5 ein stetiger übergang zwischen den Polynomen resultiert (die erste und zweite Ableitung sind im übergangspunkt gleich). Dies ist von Vorteil, wenn die Abstände zwischen den Dichtewerten für die Werte Null, 0.5 und 1 sich stark unterscheiden.

In der späteren iterativen Optimierung wird die Materialzusammensetzung jedes Pixels durch das Anteilsmodell entsprechend seiner Dichte bei jedem Iterationsschritt angepasst.

Falls die Materialien nicht a priori genau bekannt sind, kann auch von einem

Standardsatz von Materialien ausgegangen werden, der der Situation angepasst ist. Beispielsweise kann Eisen als Repräsentant für mehrere andere Metalle verwendet werden, oder es können Gewebe als eine Mischung von Wasser und Knochen

betrachtet werden („displacement model"). Anschliessend können verschiedenen Mischungsverhältnissen weitere Gewebe wie Herz. Lunge, Blut etc. zugeordnet werden.

Im Folgenden wird für die Erläuterungen von einem Modell mit wenigen (2-3) Materialien mit signifikant unterschiedlichen Dichten ausgegangen, so dass die Anteilsfunktion wie gezeigt automatisch aus dem Histogramm und seinen Spitzen bestimmbar ist.

Das Linienintegral für die i-te Messung kann als Skalarprodukt ausgedrückt werden:

wobei H = {hj } im Folgenden als Systemmatrix bezeichnet wird. Die einzelnen Elemente hy geben den Einfluss von Pixel j auf Messung i wieder. P ist die Anzahl der Pixel.

Die Matrix H wird üblicherweise auch als Projektionsmatrix oder Projektionsoperator bezeichnet und durch einen Distance-Driven [8] Algorithmus bestimmt. Die Distance-Driven Methode ist gut geeignet für parallelisierbare Rechnungen, und benötigt (zum Gegensatz zu beispielsweise einer Siddon-Methode) nur eine einzige Speicherung des Bildes und des Sinogramms. Die Matrix ist schwach besetzt und kann für eine kleine Anzahl von Messungen im Arbeitspeicher eines Rechners gehalten werden. Für grossere Mengen von Messwerten wie bei konischen Röntgenstrahlen, Helix-Scannern oder Fan-beam Geräten sind jedoch, auch bei effizienter Speicherung, mehrere Gigabyte bis tausende von Terabyte erforderlich.

Deshalb wird in der Röntgenstrahltomographie das Linienintegral jeweils immer neu "on the fly" berechnet, ohne dass die Matrix H gespeichert wird. Gemäss dem vorliegenden Verfahren wird, wie weiter unten vorgestellt, eine Projektion der Matrix gespeichert, was entsprechende Auswirkungen auf die iterative Optimierung hat.

Das Linienintegral kann umgeformt und in der folgenden Form ausgedrückt werden:

mit: ϊ„

Der Vektor s, wird oft als effektive Dicke bezeichnet.

Damit kann der Ausdruck für die log-likelihood-Funktion, unter Verwendung des idealen Erwartungswertes der Messung M; geschrieben werden als:

mit :

E h _t iή ^ Mi logüή-t, h _t {t)= ^--l,

Eine Kostenfunktion für das polyenergetische Modell wird aus den folgenden Gründen eingeführt: Bei einer wohlkonditionierten Matrix H und unverrauschten Messungen würde die Minimierung der negativen Likelihood-Funktion L ein ideales Ergebnis liefern. In der Realität sind diese Umstände nicht gegeben und ist eine Regularisierung des Problems nötig. Dies wird durch Hinzufügen einer

Kostenfunktion (penalty function) zu der Likelihood-Funktion erreicht: Es wird eine paarweise Regularisierung mit einem Regularisierungsterm oder Kostenterm der folgenden Form betrachtet:

wobei φ eine Potentialfunktion ist und Nj eine bestimmte Nachbarschaft des j-ten Pixels ist. Beispielsweise wird dafür ein konvexes kantenerhaltendes Hyperflächen- potential verwendet: φ(u) = Vl + u ²

Diese Norm hat sich in der Regularisierung von Tomogrammen bewährt (Charbonnier [4]). Die Kombination der Likelihood-Funktion mit der Kostenfunktion ergibt die folgende penalized-likelihood (PL) Kostenfunktion:

φ{p)= -L{p)+ßR{p) wobei ß ein Skalar ist welcher die relative Gewichtung zwischen der Datenüber- einstimmung und der Kostenfunlction kontrolliert.

Das Ziel der Bildrekonstruktion wird damit als eine Minimierung des Funktionais φ unter bestimmten Randbedingungen formuliert. Diese Randbedingungen sind insbesondere, dass negative Dichten nicht zulässig sind und dass vom Benutzer eine maximale Dichte vorgegeben werden kann. Somit ist das Ziel der Minimierung der Satz oder das Bild (2D oder 3D) von Dichtewerten p*, für welche p* = min φ(/>) O≤P≤Pmax Dieses multidimensionale nichtlineare Minimierungsproblem wird mit einem iterativen Algorithmus gelöst: Der Einfachheit wird im Folgenden auf den Likelihood-Term L(p) eingegangen; analoge Resultate ergeben sich für den Kostenterm R(p). Die Mininiierung fusst auf der Theorie des Optimierungstransfers ("optimisation transfer"), welche in [2] detailliert beschrieben ist.

Einfach gesagt erlaubt der Ansatz des Optimierungstransfers den Ersatz einer komplexen Likelihood-Funktion durch einfachere Ersatzfunktionen ("Surrogatfunktionen, Surrogate cost functions") welche einfacher zu minimieren sind. Durch den Optimierungstransfer kann dafür gesorgt werden, dass die Kostenfunktion mit jedem Iterationschritt monoton fällt. In Figur 8 ist dieses Prinzip für den eindimensionalen Fall illustriert: Anstelle dass φ(μ) minimiert wird, wird in der n- ten Iteration die Ersatzfunktion φ(μ; μn) minimiert. Hier hat die Ersatzfunktion φ2 eine schwächere Krümmung und ist breiter als φl, somit führt φ2 zu einer grosseren Schrittweite und zu einer schnelleren Konvergenz zum lokalen Minimum.

Das Prinzip des Optimierungstransfers wird zweimal angewandt: erstens unter Verwendung der Eigenschaft der multiplikativen Konvexität; und zweitens unter Verwendung von parabolischen Ersatzfunktionen (parabola Surrogates, siehe [2]). Dies liefert eine separierbare und einfache Ersatzfunktion, welche einfacher als die negative Log-Likelihood-Funktion zu minimieren ist. Diese paraboloide Ersatzfunktion ist [3]:

wobei n die n-te Iteration bezeichnet, und die folgenden Definitionen verwendet werden:

ist eine Ersatzfunktion von g" wobei wiederum

gf [μ _m -feWy, _{+ sc} λ _E

Die Ersatzfunktion ist separierbar, weil eine Summe über alle Pixel P gebildet wird. Dank der Separierbarkeit kann die Berechnung der Ersatzfunktion problemlos 5 parallelisiert werden.

Zur Herleitung des Optimierungsverfahrens wird die erste Ableitung bestimmt und im Sinne der Newtonschen Methode gleich Null gesetzt. Dies ergibt die folgende Vorschrift zur Aktualisierung:

worin Q eine Ersatzfunktion ("Surrogate function") für den Regularisierungsterm R(ρ) ist.

Im Folgenden wird die Bestimmung der Berechnung des Gradienten und der 5 Hesseschen Matrix von P beschrieben (Q ist unabhängig von der Energie und bereitet keine Schwierigkeiten). Der Gradient der Ersatzfunktion P ist für die n-te Iteration ohne weiteres bestimmbar als

_/ = ₁ I Mfö }+ scj) σPj

Für die Hessesche Matrix von P werden die zweiten Ableitungen von q[ benötigt. An der Stelle p" ergibt sich:

Darin ist C, ⁿ die Krümmung der paraboloiden Ersatzfunktion. Die Krümmung wird gemäss der Theorie des „optimisation transfer" so gross gewählt, dass die Ersatzfunktion immer grösser als die zu approximierende Funktion ist. Dieser Wert wird auch als vorausberechnene Krümmung ("precomputed curvature") bezeichnet,

10 da er von den Iterationen unabhängig ist.

Zur Beschleunigung der Optimierung kann, unter Verzicht auf absolute Monotonie des Verfahrens, das Verfahren der geordneten Untermengen (ordered subsets) angewandt werden, siehe zum Beispiel [5]. Darin werden Untermengen der Daten 15 umgeordnet, was die Konvergenz drastisch verbessert aber zum Teil die Monotonizität beeinträchtigt. Da die statistischen Verfahren in ihrer Grundform zu langsam sind, ist eine solche Art der Beschleunigung für eine praktische Anwendung vertretbar.

Damit das Verfahren praktisch umgesetzt werden kann, werden die folgenden drei Vereinfachungen oder Annäherungen gemacht:

1. Die Krümmung Cf (E) ₁ - wird durch eine im voraus berechnete Krümmung M ₁ ersetzt, welche unabhängig von der Iteration und von der Energie ist. Dies wird in Erdogan [2] gezeigt:

Dabei ist Mj die wirkliche, i-te Messung, wohingegen Mj(s) eine Abschätzung

E 2. Es wird ausgenutzt, dass die Variation der massebezogenen Abschwächung über das Spektrum statt dessen für die mittlere (mean) Energie ausgewertet werden kann. Dies ist insbesondere der Fall, wenn 450kV-Röhren mit externen Filtern verwendet werden, beispielsweise einem 4.5 mm Messing-Filter so dass die mittlere Energie des Strahls bei 230 keV liegt. In diesem Bereich andern sich die massebezogenen Abschwächungskoeffizienten nur wenig. Also :

3. Die Berechnung von /^. ist ohne Speichern der Projektionsmatrix H nicht mit akzeptablem Aufwand möglich. Das Speichern von H wiederum ist aus den oben schon vorgestellten Gründen in vielen Fällen nicht möglich. Deshalb werden die Elemente von H über alle Pixel j summiert und diese Summen in einem Summenvektor γ. , typischerweise einem Spaltenvektor, gespeichert. Damit wird als weitere Näherung die folgende Ungleichung eingesetzt: ^h U ^{≤ h} ü-∑ ^h ij = ^h y-ri

Bemerkung: Der Term hj _j taucht in den vorherigen und nachfolgenden Gleichungen auf, aber immer in der Form ^h ₁₁ X ₁ , d.h. im Zusammenhang mit einer Projektion oder Rückprojektion. Es werden nie die einzelnen Werte von hjj gespeichert, sondern nur die (rück)projizierten Werte.

Ferner wird b " (E) über das Spektrum normalisiert, .

Damit ergibt sich eine obere Grenze für die zweite Ableitung von P als:

Diese Funktion ist somit immer noch. Kandidat für die Krümmung der Ersatzfunktion. Prinzipiell ist die Konvergenz für grossere Werte der Krümmung langsamer.

Mit dieser Näherung besteht keine Garantie für Monotonizität. Es hat sich aber gezeigt, dass mit einem ersten rekonstruierten Schnittbild guter Qualität wie aus einer FPB-Rekonstruktion die Kostenfunktion zuverlässig abnimmt.

Der Rechenaufwand beträgt pro Iteration K Vorwärtsprojektionen und K Rückprojektionen, wobei K die Anzahl unterschiedlicher Referenzmaterialien ist. Dies ist also äquivalent zu rund 2K FPB-Durchläufen. Dank der dritten Vereinfachung kann die Krümmung ohne Speichern der Systemmatrix H berechnet werden.

Aus dem Spektrum der Röntgenstrahlen oder aus äquivalenten Messungen kann M ₁ und der Gradient V ^Mj bestimmt werden, beispielsweise in Form einer konischen

Flächenapproximation, (Cardinal and Fenster [I]) oder durch Interpolation von tabellierten Werten von M _j und VJJ-M/ . Die zweite Variante ist bei Verwendung von mehr als zwei Materialien besonders vorteilhaft.

Alternativ kann ein Array mit vorausberechneten Werten von M, und dem Gradienten V^-M/ als Basis für eine Interpolation für weitere Werte von M,- verwendet werden.

Das Verfahren kann zusammengefasst mit dem folgenden Pseudocode dargestellt werden. Bei einzelnen Schritten ist der Aufwand als Anzahl von Projektionen oder Rückprojektionen angegeben. Hier ist auch die Verwendung von geordneten Untermengen mit dargestellt:

• Initialisieren der Dichte der Bildelemente mit p°=p ^FBP aus einer gefilterten Rückproj ektion.

• Berechnen des Dichtehistogramms und Berechnen der Extrema durch Gauss'sche Approximation

• Berechnen der Anteilmodelle jω [p )j _κ

Berechnen der Polynome — ( p ω [p J • Berechnen der Parameter der konischen Approximation von

• Vorausberechnen des zweiten Terms D ₂ der Krümmung

or jede Iteration n=l ..νiter

Berechnen des Gradienten und der Krümmung der paraboloiden Ersatzfunktion

5des Regularisierungsterms: -^- . -^4- dp _t dp-

For jede Untermenge S=l..NbSub

Bereclinen von sf = ^' ∑Oi _j P _j ∞lψ _j ] s" = H'"— ^ ¹ " J ( ^κ Projektionen)

Berechnen von Mj iyf L V^M/ Ly" 1 mit konischen Approximationen

Berechnen von Zj (K Sinogramme)

N

Berechnen von Z : BPZ ^h = υ _J h _i UjZf (K Rückprojektionen) i=l

For jedes Pixel j = 1 ...P

Berechnen von Nj = ∑ _k J'j BPZj , D _χ = J] _k ή BPZ^, Aufdatieren der Dichte gemäss der Iterationsvorschrift:

N _j + + D ₂ ) _f

End of loop - über Pixel End of loop - über Untermengen End of loop - über Iterationen

In einer Variante ohne geordnete Untermengen wird anstelle der Schlaufe über die Untermengen eine Schlaufe über die gesamte Anzahl von Messungen, das heisst über das gesamte Sinogramm ausgeführt. Die neuen Werte von p" ^+ι sind nur von den

Werten von p" abhängig, das heisst, dass das Aufdatierungsverfahren simultan ist.

Referenzen :

[1] H. N. Cardinal and A. Fenster. An accurate method for direct dual-energy calibration and decomposition. Med. Phys.. 17(3):327-41, May 1990.

[2] H Erdogan, J A Fessler. Monotonie algorithms for transmission tomography. IEEE Trans. Med. Imag., 18(9):801-14, Sep. 1999.

[3] I A Elbakri, J A Fessler. Statistical image reconstruction for polyenergetic X-ray computed tomography. IEEE Trans. Med. Imag., 21(2):89-99, Feb. 2002.

[4] P. Charbonnier et al. Deterministic edge-preserving regularization in computed tomography. IEEE Trans. Image. Processing , vol 6, no 2, pp 298-311, 1997

[5] H. M. Hudson and R. S. Larkin. Accelerated image reconstruction using ordered subsets of projection data. IEEE Transactions on Medical Imaging, 13(4):601- 609, December 1994.

[6] A. R. De Pierro, A modified expectation maximization algorithm for penalized likelihood estimation in emission tomography," IEEE Tr. Med. Im., vol. 14, no. 1, pp.l32-137 ₅ March l995.

[7] W. Zbijewski and FJ Beekman "Suppression of intensity transition artefacts in Statistical x-ray CT reconstruction through Radon Inversion initialization" Med.Phys 31(1), January 2004.

[8] Bruno De Man and Samit Basu. Distance-driven projection and backprojection in three dimensions. Phys. Med. Biol. 49 No 11 (7 June 2004) 2463-2475.