Kurzfassung

Vorgestellt wird ein neuer iterativer Algorithmus zur Strahlaufhärtungskorrektur in der Kegelstrahl-Computertomographie (3D CT). Die Computertomographie ist eine zerstörungsfreie Untersuchungsmethode um Form und Position von Fehlstellen in Gussteilen zu analysieren.

Polychromatische Röntgenquellen verursachen in den rekonstruierten Volumen unter anderem Strahlaufhärtungsartefakte, die die Bildqualität erheblich vermindern. Die Linearisierungsmethode ist ein Korrekturalgorithmus, der unter Verwendung einer nichtlinearen Korrekturkennlinie diese Artefakte reduziert. Die hier vorgestellte Methode extrahiert die Korrekturkennlinie aus dem Untersuchungsobjekt selbst. Ein Referenzobjekt wird nicht benötigt. Die Tauglichkeit der Methode wird an geeigneten Aluminium-Gussteilen demonstriert.

Schlüsselwörter:
Computertomographie, Artefakte, Strahlaufhärtungskorrektur

1. Einführung

Computertomographie ist ein digitales Bildgebungsverfahren, welches auf der Schwächung von Röntgenstrahlung beim Durchgang durch ein Objekt beruht. Das rekonstruierte Volumen des Objektes besteht aus Voxeln, deren Grauwert den linearen Schwächungskoeffizienten m repräsentiert.

Artefakte sind künstliche Strukturen in Bildern, für die es in der Struktur des abgebildeten Objektes keine Entsprechung gibt. Ursache sind entweder technische Fehler des Bilderzeugungssystems oder es werden Aufnahmebedingungen nicht erfüllt oder verletzt.

Die CT-Rekonstruktion setzt monoenergetische Röntgenquellen voraus. Polyenergetische Strahlung wird beim Durchgang durch Materie aufgehärtet, weil niederenergetische Photonen stärker absorbiert werden. Dieser Effekt heißt Strahlaufhärtung (beam hardening BH) und führt unkorrigiert zu störenden Artefakten im rekonstruierten Volumen.

Die Linearisierungstechnik ist eine bewährte Korrekturmethode, die mittels einer Korrekturfunktion polyenergetische auf monoenergetische Projektionsdaten transformiert. Die Korrekturfunktion wird durch Messung eines geeigneten Referenzobjektes ermittelt, was mit zusätzlichen Aufwand verbunden ist.

Ziel dieser Arbeit ist es, diese Kalibriermessung durch eine Methode zu ersetzen, die unabhängig von einem Referenzobjekt funktioniert. Die Iterative Artefakt Reduzierung (IAR) ist ein mehrstufiger Prozess, der aus dem rekonstruierten Volumen die benötigte Korrekturkennlinie berechnet. Der Beitrag ist wie folgt aufgebaut:

Im nächsten Abschnitt wird der physikalische Hintergrund der CT-Bildgebung beschrieben. Abschnitt 3 stellt die Linearisierungstechnik vor, die Abschnitte 4 und 5 stellen die neue Methode zur Berechnung der Korrekturkennlinie vor. Die Ergebnisse werden in Abschnitt 6 präsentiert, in Abschnitt 7 folgt eine Zusammenfassung.

2. Theoretischer Hintergrund

Die Wechselwirkung von elektromagnetischer Strahlung mit Materie wird durch ein exponentielles Schwächungsgesetz beschrieben. Für monoenergetische, kurzwellige Röntgenstrahlung gilt die Transmissionsgleichung von Lampert-Beer:

(1)

und bezeichnen die emittierte bzw. detektierte Intensität des Röntgenstrahls nach Durchstrahlung des Objektes. Der lokale Schwächungskoeffizient ist abhängig von der Energie der einfallenden Strahlung. l bezeichnet die Strecke vom Brennfleck der Quelle zu einem Detektorelement.

Der zugehörige Projektionswert ist eine lineare Funktion der Absorberdicke L:

(2)

Die letzte Gleichung ist nur für homogene Objekte richtig. Eine CT-Rekonstruktion basiert auf folgendem Prinzip: sind für verschiedene Winkel q die Projektionen bekannt, läßt sich die unbekannte Funktion unter Verwendung der inversen Radon-Transformation berechnen.

(3)

Für eine ausführliche Beschreibung der CT-Rekonstruktionsalgorithmen sei z.B. auf [1] hingewiesen. Industrielle Computertomographen sind meistens mit Röntgenröhren ausgestattet. Das erzeugte Energiespektrum S(E) ist polyenergetisch und verursacht in Gleichung (2) einen nichtlinearen Zusammenhang zwischen Projektionswert und Schwächungskoeffizient. Aus dieser Nichtliniarität resultieren Fehler, die als sog. Strahlaufhärtungsartefakte bekannt sind. Unkorrigiert führt eine Rekonstruktion des Volumens aus den gemessenen Projektionen zu Abbildungsfehlern wie z.B. dunklen Streifen entlang langer Absorbtionskanten. Selbst bei homogenen Objekten erscheinen Zentralbereiche weniger dicht als Randzonen. Dieser störende Effekt heißt cupping-Artefakt.

3. Strahlaufhärtungskorrektur für homogene Objekte

Eine weit verbreitete und effektive Korrekturmethode von Strahlaufhärtungsartefakten für homogene Objekte basiert auf der Linearisierungstechnik. Es handelt sich um einen Vertreter der vorverarbeitenden Korrekturtechniken. Trägt man in einem Diagramm monoenergetische und polyenergetische Projektionswerte (y-Achse) gegen die durchstrahlte Weglänge im Objekt (x-Achse) auf, erhält man eine Gerade bzw. eine konkave Kurfe durch den Ursprung. Das Prinzip der Linearisierung besteht darin, einen gemessenen polyenergetischen Projektionswert von der nichtlinearen Kurfe auf eine korrespondierende monoenergetische Strahlsumme auf der linearen Kurfe zu transformieren.

In der Praxis wird die nichtlineare Kurfe experimentell aus der Radiographie eines homogenen Eichkörpers mit bekannten Längen ermittelt. Die Materialeigenschaften von Referenz- und Untersuchungsobjekt sollten möglichst übereinstimmen. In [2] wird für den Fit der Messpunkte ein Polynom N. Grades vorgeschlagen.

Allerdings wird der Zusammenhang der detektierten Schwächungsverhältnisse und der Weglängen eines Röntgenstrahls durch das Objekt physikalisch angemessener durch eine Summe von Exponentialfunktionen beschrieben.

(4)

Die Funktion H_s(.) bezeichnet die Strahlaufhärtungskennlinie bezüglich des Energie-spektrums S(E). Man erhält Datenpunkte der Funktion H_s, indem man für verschiedene bekannte Längen L des Eichkörpers das Schwächungsverhältniss bestimmt. Häufig werden Stufenkeile als Referenzobjekt verwendet. Die gesuchte Korrekturkennlinie H_s^-1 erhält man durch Inversion der approximierten Funktion H_s:

(5)

Wenn die Korrekturkennlinie H_s^-1 ermittelt ist, funktioniert die Strahlaufhärtungskorrektur wie folgt: Für jeden Messwert einer Projektion und für alle Projektionen des Objekts wird das Schwächungsverhältnis durch die Weglänge L aus Gleichung (5) ersetzt. Nach diesem Vorverarbeitungsschritt erfolgt die inverse Radontransformation mit den modifizierten Projektionsdaten, um ein strahlaufhärtungskorrigiertes Volumen zu rekonstruieren.

Der Nachteil der Linearisierungsmethode liegt auf der Hand: für jedes Untersuchungs-objekt wird ein Referenzobjekt aus dem gleichen Material benötigt. Ausserdem ist die Kalibriermessung nur für ein spezielles Energiespektrum S(E) gültig. Der Vorteil unserer Arbeit liegt darin, die benötigte Korrekturkennlinie aus dem Untersuchungsobjekt selbst zu bestimmen. Im nächsten Abschnitt werden ausführlich die Nachverarbeitungsschritte erläutert, die zur Berechnung der Korrekturkennlinie H_s^-1 erforderlich sind.

4. Berechnung der Korrekturkennlinie

Ziel der im folgenden dargestellten 3D Bildverarbeitungsschritte ist die Ermittlung der Korrekturkennlinie aus dem untersuchten, rekonstruierten Volumen. Diese Nachverarbeitungschritte ersetzen die aufwendige Kalibriermessung mit einem Eichkörper. Es handelt sich im Einzelnen um folgende Operationen: Segmentation des Untersuchungsobjektes, Ermittlung der Strahlsummen (Ray Summing), also die Berechnung der durchstrahlten Weglänge L durch das Objekt und anschliessender Fit der erhaltenen Datenpunkte.

A. Segmentation
Bevor Weglängen durch das Objekt berechnet werden können, ist es wichtig, die Geometrie des Objektes genau zu kennen. Dazu müssen alle Voxel des rekonstruierten Volumens möglichst eindeutig in Objekt oder Hintergrund separiert werden. Zu diesem Zweck wenden wir eine 3D Segmentierung auf das Volumen an. Im vorliegenden Fall eines homogenen Objektes, bedeutet Segmentierung die Erzeugung eines Binärvolumens V^bin. Eine schwellwertbasierte Binarisierung funktioniert nur dann zufriedenstellend, wenn die Bildqualität hoch ist. Jedoch verursachen in der Kegelstrahl-CT Strahlaufhärtung, Rauschen und Streustrahlung erhebliche Fehler bei der Kantendetektion. Deshalb ist eine einfache Schwellwertmethode für unseren Anwendungsbereich ungeeignet. Stattdessen erzeugen wir das Binärvolumen in einem zweistufigen Prozess. Nach einer hysteresegesteuerten Vorbinarisierung vollendet ein Konturverfolgungsalgorithmus die Segmentierung. Der Segmentierungsprozess benötigt das rekonstruierte Volumen V^m und ein Konturvolumen V^c, um daraus das Binärvolumen V^bin zu berechnen.

Vorbinarisierung
Betrachten wir ein Voxel im rekonstruierten Volumen und zwei vorgegebene Schwellwerte Th_u, Th_l. Ist der Grauwert des Voxels größer als die obere Schwelle Th_u, wird das Voxel dem Objekt zugerechnet. Liegt der Grauwert unterhalb des unteren Schwellwertes Th_l gehört das Voxel zum Hintergrund. Alle Voxel, deren Grauwerte zwischen diesen beiden Grenzen
Th_l<gv < Th_u liegen, werden weiter untersucht und zwar mit dem

Konturverfolgungsalgorithmus
Die Idee dieses Algorithmuses läßt sich an einem einfachen, zweidimensionalen Beispiel illustrieren: Betrachten wir eine geschlossene Kontur C, zum Beispiel einen Kreis in einer Ebene und einen Testpunkt t. Um zu entscheiden, ob der Testpunkt innerhalb oder ausserhalb der Kontur liegt, verfolgen wir einen Strahl R beginnend am Testpunkt in eine beliebige Richtung ins Unendliche und ermitteln die Anzahl der Schnittpunkte von R mit C. Der Testpunkt t befindet sich innerhalb der Kontur, wenn die Anzahl der Schnittpunkte ungerade ist. Der Algorithmus funktioniert im Prinzip, allerdings existieren mehrere Aus-nahmen, zum Beispiel Tangentenschnittpunkte. Es ist daher erforderlich für einen Testpunkt mehrere Strahlen zu untersuchen und einen Mehrheitsentscheid durchzuführen.

Um diesen Algorithmus auf Voxeldaten anwenden zu können, muß zunächst die Oberfläche des betrachteten Objektes extrahiert werden. Ein dreidimensionaler Kantenoperator liefert auf das rekonstruierte Volumen angewendet ein Konturvolumen V^c. In der dreidimensionalen Bildverarbeitung ist die diskrete Repräsentation eines Strahls eine geordnete Folge von Voxeln n_o...n_n. Beispielsweise kann als Startvoxel n_o eines Strahls einer der 26 Nachbarn des Testvoxels n_t gewählt werden. Im folgenden Schritt, "Ray Summing" genannt, wird die Weglänge durch das nun segmentierte Objekt unter Verwendung des Binärvolumens V^bin "gemessen".

B. Ray Summing
Aufgabe des Ray Summers ist die Ermittlung von "Messpunkten" ( , L) für die gesuchte Funktion H nach Gleichung (4). Die voxelbasierte Ray-Tracing-Technik arbeitet wie folgt: Der Algorithmus verfolgt alle Strahlen von der Punktquelle zu allen Rasterpunkten der Detektorfläche. Dabei berechnet er die Weglänge des Strahls durch das Objekt und registriert die zugehörige Schwächung der Strahlung.

C. Fit
Wie bereits weiter oben aufgeführt werden die ermittelten Datenpunkte ( , L) mit einer Summe von Exponentialfunktionen approximiert. Es handelt sich um ein nichtlineares Ausgleichsproblem, das mit Standardroutinen, vergleiche z.B. das Levenberg-Marquard-Verfahren [3] gelöst werden kann. Die Korrekturfunktion für die Linearisierung erhält man durch Inversion der Approximierten.

5. Iterative Artefakt Reduzierung (IAR)

In den beiden vorangegangenen Abschnitten haben wir die erforderlichen Vor- und Nachverarbeitungsschritte dargelegt, die zur Berechnung und Verwendung einer Korrekturkennlinie für die Strahlaufhärtungskorrektur nötig sind.

Die Iterative Artefakt Reduzierung (IAR) ist ein mehrstufiger Prozess, der diese Vor- und Nachverarbeitungsmethoden kombiniert. Nach jeder Iteration des IAR-Alogrithmuses liegt eine verbesserte Korrekturfunktion H_s^-1 vor. Der Algorithmus funktioniert wie folgt:

1. INITIALISIERUNG oder LINEARISIERUNG
2. REKONSTRUKTION
3. NACHVERARBEITUNG
4. ABBRUCH oder nächste ITERATION

Initialisierung
Die Korrekturfunktion H_s^-1 für die erste Rekonstruktion wird unter der (falschen) Annahme initialisiert, dass die Aufnahme der Projektionsdaten mit monoenergetischer Röntgenstrahlung erfolgt ist:

(6)

Rekonstruktion
Mittels dreidimensionaler inverser Radontransformation rekonstruiert die erste Iteration aus einem gemessenen Satz Rohdaten P₀ die initialen 3D Voxeldaten. Dieses Ausgangsvolumen V₀(r) ist unkorrigiert und enthält somit alle Artefakte.

(7)

Nachfolgende Iterationen rekonstruieren verbesserte Volumen V_N(r), in dem sie die modifizierten Projektionsdaten P_N aus dem Vorverarbeitungsschritt, d.h. der Linearisierung verwenden.

(8)

Nachverarbeitung
Nach jeder Rekonstruktion werden die Nachverarbeitungsschritte Segmentierung, Ray Summing und Approximation durchgeführt, um die Korrekturkennlinie für die nächste Iteration zu berechnen.

Abbruch
Wenn sich die Differenz zweier nachfolgender Korrekturkennlinien nur noch geringfügig unterscheidet terminiert der IAR-Prozess, ansonsten startet die nächste Iteration mit der Linearisierung.

6. Ergebnisse

Um unser IAR-Verfahren zu testen, haben wir einen Testkörper und einen Eichkörper, beide aus der gleichen Aluminium-Legierung konstruiert. Bei dem Testkörper handelt es sich um einen Quader mit den Abmessungen 20 x 10 x 4 cm. Er enthält Schrägen und Bohrungen mit unterschiedlichen Durchmessern. Unser Referenzobjekt ist ein Stufenkeil bestehend aus 30 Stufen von 0.05 cm bis 20 cm Länge und jeweils 0.5 cm Höhe.

Testkörper und Stufenkeil haben wir mit einer Seifert DP 424 Anlage mit folgenden Parametern aufgenommen: 225 kV, 3 mA, 512³ Voxel, Voxelgröße 373 m m, 800 Projektionen, Abstand Quelle-Detektor 3 m.

Zur Überprüfung der IAR-Methode haben wir die "gleiche" Korrekturkennlinie auf drei verschieden Arten berechnet. Wenn unsere Methode funktioniert, sollten die drei erhaltenen Funktionen gut übereinstimmen:

• H^-1_SW,C nach der herkömmlichen Methode, Fit der gemessenen Daten für den Stufenkeil,
• H^-1_SW,IAR nach der IAR-Methode für den Stufenkeil und
• H^-1_TP,IAR nach der IAR- Methode für den Testkörper.

Um die Resultate vergleichen zu können, tragen wir, wie in der Literatur üblich, die Projektionswerte gegen die Durchstrahlungslängen auf. Wie in Abbildung 1 zu sehen ist, unterscheiden sich die drei Funktionen in akzeptablen Grenzen.

Abb 1: Projektionswerte gegen Durchstrahlungslängen: monoenergetische Linie , punktweise Vermessung des Aluminium Stufenkeiles, Stufenkeil: gefittete Daten -In( HSW,C ), Stufenkeil: -In( HSW,IAR ) mit IAR-Methode, Testkörper: -In( HTP,IAR ) mit IAR-Methode.

Die Steigung der monoenergetischen Linie wurde so gewählt, dass sie mit der Steigung der 1. Ableitung von-In( H_SW,C ) für Objektdicke Null übereinstimmt.

a)	b)
Abb 2: Rekonstruierte Schicht des Testkörpers a) vor und b) nach Strahlaufhärtungskorrektur mit IAR-Methode ( H-1TP,IAR ). Beachtenswert sind die dunklen Streifen in a).

Abbildungen 2 mit 4 zeigen eine rekonstruierte Schicht und Profilschnitte unseres Testkörpers vor und nach der Strahlaufhärtungskorrektur. Auftretende Artefakte sind deutlich reduziert worden.

a)	b)
Abb 3: Profilschnitte des Testkörpers bei P1, siehe Abbilung 1. a) cupping Artefakt und b) mit IAR-Korrektur.

a)	b)
Abb 4: Profilschnitte des Testkörpers bei P2, siehe Abbilung 1. a) vor und b) nach IAR- Korrektur. Beachtenswert ist, dass in a) die Grauwerte nicht auf Null absinken..

7. Zusammenfassung

Mit der Iterativen Artefakt Reduzierung (IAR) haben wir ein effektives Verfahren zur Strahlaufhärtungskorrektur vorgestellt. Die Berechnung der Korrekturkennlinie erfolgt aus dem rekonstruierten Volumen selbst. Die IAR-Methode benötigt die Projektionsdaten und funktioniert derzeit nur bei homogenen Objekten. Eine Kenntnis des Energiespektrums oder des untersuchten Materials ist nicht erforderlich. Wir haben unsere Methode an einen Testkörper überprüft und einen Stufenkeil verwendet, um den IAR-Algorithmus mit bisherigen Korrekturmethoden zu vergleichen. Wir haben demonstriert, dass die Bildqualität durch die iterative Korrekturmethode offensichtlich verbessert wird. Basierend auf den Testdaten haben wir gezeigt, dass die IAR ein geeignetes Verfahren ist, um die zeitaufwendige Aufnahme und Auswerung von Eichkörpern zu ersetzen.

Referenzen

1. Avinash C. Kak and Malcolm Slaney, Principles of Computerized Tomographic Imaging, IEEE Press, New York, 1988.
2. Gabor T. Herman, "Correction for Beam Hardening in Computed Tomography", Phys. Med. Biol., Vol. 24, No. 1, 1979, pp 81-106.
3. W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, Numerical Recipes in FORTRAN, Cambridge University Press, Cambridge, pp 678 - 683, 1992