Studie en evaluatie van fasecontrast tomografie aan een hoge resolutie X-stralenbuis

Faculteit Wetenschappen Vakgroep Fysica en Sterrenkunde Voorzitter: Prof. Dr. D. Ryckbosch

Studie en evaluatie van fasecontrast tomografie aan een hoge resolutie X-stralenbuis door Wouter Devulder

Promotor: Prof. Dr. L. Van Hoorebeke Scriptiebegeleiders: ir. M. Boone, Dr. ir. M. Dierick

Masterproef ingediend tot het behalen van de graad van Master in de Fysica en Sterrenkunde

Academiejaar 2010–2011

Faculteit Wetenschappen Vakgroep Fysica en Sterrenkunde Voorzitter: Prof. Dr. D. Ryckbosch

Studie en evaluatie van fasecontrast tomografie aan een hoge resolutie X-stralenbuis door Wouter Devulder

Promotor: Prof. Dr. L. Van Hoorebeke Scriptiebegeleiders: ir. M. Boone, Dr. ir. M. Dierick

Masterproef ingediend tot het behalen van de graad van Master in de Fysica en Sterrenkunde

Academiejaar 2010–2011

Voorwoord Toen ik vijf jaar geleden besloten had fysica te studeren, wist ik eigenlijk totaal niet wat me te wachten stond. Het enige dat ik wist — of toch hoopte — was dat het enorm boeiend zou zijn. Nu, vijf jaar later, schrijf ik het laatste deel van mijn thesis en kan ik enkel bevestigen dat mijn vermoeden correct was. Zoals het in een voorwoord hoort, wens ik hier ook enkele mensen te bedanken. In de eerste plaats professor Van Hoorebeke, om mij de kans te geven mijn thesis in zijn onderzoeksgroep te maken. Mijn thesisbegeleiders, Matthieu Boone en Manuel Dierick, om me het hele jaar door met raad en daad bij te staan. Pieter Vanderniepen voor de hulp bij de scans, Yoni De Witte voor de uitleg bij zijn doctoraatsthesis en uiteraard ook de rest van de groep voor de aangename werksfeer. Deze vijf jaar zouden niet geweest zijn wat ze waren zonder mijn medestudenten, bedankt voor de leuke momenten en de aangename groepssfeer! Om dit voorwoord af te sluiten een welgemeend woord van dank naar mijn ouders toe. Zij hebben mij de kans gegeven om dit alles te verwezenlijken, en hebben me steeds gesteund in wat ik deed. Mijn zus wil ik uiteraard ook bedanken voor het nalezen van mijn thesis!

Wouter Devulder, 26 mei 2011

i

Toelating tot bruikleen

“De auteur geeft de toelating deze scriptie voor consultatie beschikbaar te stellen en delen van de scriptie te kopiëren voor persoonlijk gebruik. Elk ander gebruik valt onder de beperkingen van het auteursrecht, in het bijzonder met betrekking tot de verplichting de bron uitdrukkelijk te vermelden bij het aanhalen van resultaten uit deze scriptie.”

Wouter Devulder, mei 2011

ii

Studie en evaluatie van fasecontrast tomografie aan een hoge resolutie X-stralenbuis door Wouter Devulder Masterproef ingediend tot het behalen van de graad van Master in de Fysica en Sterrenkunde Academiejaar 2010–2011 Promotor: Prof. Dr. L. Van Hoorebeke Scriptiebegeleiders: ir. M. Boone, Dr. ir. M. Dierick Faculteit Wetenschappen Universiteit Gent Vakgroep Fysica en Sterrenkunde Voorzitter: Prof. Dr. D. Ryckbosch

Samenvatting In deze masterproef bestuderen we het gebruik van fasecontrast in X-stralen computertomografie aan een hoge resolutie X-stralenbuis. Wanneer een X-stralenbundel zich in een object voortplant wordt hij geattenueerd, maar hij ondergaat ook een faseverschuiving. Dit resulteert in een afbuiging van de bundel. Aan een opstelling met een X-stralenbuis resulteert deze afbuiging in een piek en een dal in de intensiteit. Indien hiervoor niet gecorrigeerd wordt, veroorzaakt deze verandering in het intensiteitspatroon artefacten in de reconstructie die gebaseerd is op de attenuatie. In deze thesis proberen we in een eerste deel dit fenomeen zo goed mogelijk te simuleren. Indien dit fasesignaal met een goede nauwkeurigheid kan gesimuleerd worden, biedt dit mogelijkheden om het te implementeren in een iteratief reconstructie algoritme. Het fasesignaal kan in bepaalde reconstructie algoritmes (bv. het Modified Bronnikov Algoritme (MBA) en de methode van Paganin (MP)) gebruikt worden om op basis hiervan een reconstructie van het object te maken. De hierboven vermelde algoritmes lijken sterk op elkaar, hoewel ze van verschillende veronderstellingen uitgaan. In het tweede deel van deze masterproef worden beide methodes met elkaar vergeleken en wordt nagegaan of de ene methode eventueel betere resultaten oplevert dan de andere. De methodes worden geëvalueerd met gesimuleerde en echte samples.

Trefwoorden Computertomografie, fasecontrast, Modified Bronnikov Algoritme, methode van Paganin

Inhoudsopgave Voorwoord

i

Toelating tot bruikleen

ii

Overzicht

iii

Inhoudsopgave

iv

Lijst van gebruikte afkortingen

vii

1 Inleiding 2 X-stralen tomografie 2.1 X-stralen . . . . . . . . . . . . . 2.1.1 Productie van X-stralen . 2.1.2 Interactie met materie . . 2.2 Computertomografie . . . . . . . 2.2.1 CT-scan . . . . . . . . . . 2.2.2 Analytische reconstructie 2.2.3 Praktisch . . . . . . . . . 2.2.4 Algebra¨ısche reconstructie

1

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

3 Fasecontrast 3.1 Inleiding . . . . . . . . . . . . . . . . . . . 3.2 Stralenoptica . . . . . . . . . . . . . . . . 3.2.1 Complexe brekingsindex . . . . . . 3.2.2 Refractie . . . . . . . . . . . . . . 3.2.3 Verband tussen faseverschuiving en 3.2.4 Afleiding van de TIEW . . . . . . 3.3 Visualisatie van de fase . . . . . . . . . . 3.3.1 Interferometrie . . . . . . . . . . . 3.3.2 Diffraction enhanced imaging . . . 3.3.3 Grating interferometrie . . . . . . 3.3.4 Vrije propagatie . . . . . . . . . . 3.4 Golfoptica . . . . . . . . . . . . . . . . . . iv

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

4 4 5 5 6 6 8 10 10

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . de gemeten intensiteit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

14 14 15 15 17 18 20 21 21 22 22 22 24

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

Inhoudsopgave 3.4.1 3.4.2

v Lineaire shift-invariante systemen . . . . . . . . . . . . . . . . . . . . Vrije propagatie als lineair shift-invariant systeem . . . . . . . . . .

4 Experimentele methode 4.1 Opstelling UGCT . . . 4.2 Gebruikte software . . 4.3 Sample . . . . . . . . . 4.4 Secondary spot . . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

5 Simulatie van fasecontrast 5.1 Inleiding . . . . . . . . . . . . . . . . . . 5.2 Simulatie homogene cilinder . . . . . . . 5.2.1 Vlakke golf versus sferische golf . 5.2.2 Vergelijking met het experiment 5.2.3 Polychromatische bundel . . . . 5.2.4 Variërende SOD . . . . . . . . . 5.3 Bespreking . . . . . . . . . . . . . . . . 6 Fasereconstructie 6.1 Bronnikov algoritme . . . . . . . . . . . 6.1.1 Theorie . . . . . . . . . . . . . . 6.1.2 Modified Bronnikov algoritme . . 6.1.3 Implementatie . . . . . . . . . . 6.2 Methode van Paganin . . . . . . . . . . 6.2.1 Theorie . . . . . . . . . . . . . . 6.2.2 Implementatie . . . . . . . . . . 6.2.3 Vergelijking tussen MBA en MP

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

24 25

. . . .

30 30 31 31 32

. . . . . . .

34 34 35 35 35 36 40 43

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

45 45 45 47 48 48 48 49 50

7 Vergelijkende studie tussen MBA en MP 7.1 Invloed van de attenuatie . . . . . . . . . . . . . . 7.1.1 Simulaties . . . . . . . . . . . . . . . . . . . 7.1.2 Hars . . . . . . . . . . . . . . . . . . . . . . 7.2 Invloed van de homogeniteit . . . . . . . . . . . . . 7.2.1 Invloed van de homogeniteit op MBA . . . 7.2.2 Invloed van de homogeniteit op MP . . . . 7.2.3 Verschillende verhouding bij MBA . . . . . 7.2.4 Inhomogeen sample: hars met grafietstaafje 7.2.5 Bespreking . . . . . . . . . . . . . . . . . . 7.3 Ruis . . . . . . . . . . . . . . . . . . . . . . . . . . 7.3.1 Simulatie van de ruis . . . . . . . . . . . . . 7.3.2 Resultaat . . . . . . . . . . . . . . . . . . . 7.4 Toepassingen . . . . . . . . . . . . . . . . . . . . . 7.4.1 Pilletje . . . . . . . . . . . . . . . . . . . . . 7.4.2 Diamant . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

51 51 51 52 54 54 55 58 61 62 65 65 66 68 68 68

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

Inhoudsopgave 7.4.3 7.4.4

vi Houtschimmel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Vliegenpoot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

69 70

8 Besluit

73

A Fantoombeschrijving

75

Bibliografie

77

Lijst van figuren

82

Lijst van tabellen

84

Lijst van gebruikte afkortingen CT

Computertomografie

FBP

Filtered Backprojection Algorithm

FDK

Cone beam reconstructie algoritme beschreven door Feldkamp, David en Kress

FT

Fouriergetransformeerde (Engels: Fourier Transform)

MBA

Modified Bronnikov Algoritme

MP

Methode van Paganin

NRMSE

Normalized Root Mean Squared Error

ODD

Object-Detector Distance

SART

Simultaneous Algebraic Reconstruction Technique

SDD

Source-Detector Distance

SOD

Source-Object Distance

TIE

Transport of Intensity Equation

TIEW

Transport of Intensity Equation for Weak attenuation

XRT

X-stralenbuis (Engels: X-ray tube)

vii

Hoofdstuk 1

Inleiding Door de jaren heen zijn tal van technieken ontwikkeld voor het visualiseren van objecten. Verschillende probes zoals elektromagnetische golven, neutronen, elektronen, geluidsgolven. . . worden gebruikt. Bij de elektromagnetische golven vinden we o.a. X-stralen terug. Deze zijn sterk penetrerend en ondergaan een attenuatie wanneer ze door een voorwerp gaan. Reeds van bij hun ontdekking werden ze gebruikt om op basis van hun attenuatie projectiebeelden te maken van objecten. In een radiografie bekomt men een projectiebeeld van een object: de X-stralen worden geattenueerd door het sample en op de detector wordt een intensiteitspatroon gemeten. De dimensie die hierbij verloren gaat, wordt in X-stralen computertomografie (CT) gereconstrueerd door het object te roteren en een reeks radiografieën te nemen. Tomografie slaat op de reconstructie van doorsnedes uit projectiedata. De wiskundige basis voor de reconstructie van een doorsnede uit projectiedata dateert van de paper van Radon in 1917 [1]. De eerste praktische toepassing kwam in 1971 met de medische CT-scanner, ontworpen door Godfrey Hounsfield. In 1979 ontving hij hiervoor de Nobelprijs voor de fysica. [2, 3] Attenuatie van X-stralen was lange tijd het enige contrastmechanisme dat gebruikt werd bij CT. Het nadeel is dat de verschillende structuren die men wil kunnen onderscheiden een verschillende en voldoende hoge attenuatie moeten vertonen. Veel (organische) materialen zijn slechts weinig attenuerend en vertonen onvoldoende contrast om er structuren in waar te nemen. Naast attenuatie ondergaat een X-stralenbundel die door materie gaat ook een faseverschuiving, wat resulteert in een héél kleine afbuiging van de bundel. Op het einde van vorige eeuw verschenen de eerste resultaten waarbij naast de attenuatie, ook de faseverschuiving van de X-stralenbundel gebruikt werd [4, 5, 6, 7]. Materialen die weinig attenueren kunnen bv. wel een significante faseverschuiving veroorzaken. Hierdoor kunnen weinig attenuerende materialen — en structuren binnen het object — gevisualiseerd worden. Verschillende 1

2 technieken zijn ontwikkeld om deze faseverschuiving te detecteren en te gebruiken in de reconstructie1 . In figuur 1.1a wordt een klassieke transmissieradiografie van het hart van een rat getoond. In figuur 1.1b is een projectiebeeld weergegeven gebaseerd op de faseverschuiving. Het is duidelijk dat veel meer details zichtbaar zijn op het fasebeeld. Een bijkomend voordeel van deze techniek is dat bij een hogere fotonenergie — daar waar de attenuatie al zeer laag is — de faseverschuiving toch nog significant kan zijn. Bovendien kan de geabsorbeerde dosis straling verlaagd worden vermits deze techniek niet op absorptie gebaseerd is. [4, 8]

Figuur 1.1: Radiografie van een rattenhart. (a) Radiografie gebaseerd op attenuatie, (b) fasecontrast beeld [9]

Het is duidelijk dat attenuatie en faseverschuiving gelijktijdig optreden. De refractie van de bundel resulteert in een piek in de intensiteit waar de afgebogen bundel terechtkomt op de detector, en een dal waar de bundel — in afwezigheid van refractie — had moeten terechtkomen. Men kan zoals reeds vermeld deze faseverschuiving gebruiken om op basis hiervan een reconstructie van het object te maken. Indien men deze faseverschuiving niet wenst te gebruiken in de reconstructie, is dit een storend effect en veroorzaakt het artefacten. 1

Deze worden besproken in sectie 3.3.

3

Inhoud In deze masterproef zullen we ons focussen op dit fasefenomeen. In hoofdstukken 2 en 3 leggen we eerst algemeen uit wat X-stralen tomografie en fasecontrastbeeldvorming is. De gebruikte opstelling wordt besproken in hoofdstuk 4. In hoofdstuk 5 trachten we het fenomeen te simuleren. Een goede simulatie biedt mogelijkheden om voor het fasesignaal te corrigeren in iteratieve reconstructie algoritmes. Dit zou artefacten t.g.v. de faseverschuiving kunnen elimineren. Vervolgens focussen we op twee algoritmes die het fasesignaal gebruiken om een reconstructie te maken van het refractive index decrement δ 2 . Deze twee algoritmes zullen zeer gelijkaardig blijken. In hoofdstuk 6 leiden we de algoritmes af en in hoofdstuk 7 worden beide methodes vergeleken. In hoofdstuk 8 komen we tot een besluit.

2

Deze parameter wordt besproken in hoofdstuk 3.

Hoofdstuk 2

X-stralen tomografie In dit hoofdstuk leggen we de werking van X-stralen tomografie uit. Vermits deze onderzoekstechniek gebruik maakt van X-stralen, bespreken we eerst hoe een X-stralenbundel met materie interageert. Bij een radiografie verliezen we één van de drie dimensies van het object. In dit hoofdstuk zullen we dan ook uitleggen hoe X-stralen tomografie een oplossing biedt om deze verloren dimensie te reconstrueren. We bespreken in het bijzonder de wiskundige basis van analytische reconstructie, vermits we deze nodig hebben wanneer we in hoofdstuk 6 de fasereconstructie algoritmes zullen bekijken. We leggen ook het principe van iteratieve reconstructie uit.

2.1

X-stralen

Wilhelm Conrad R¨ ontgen (1845-1923) wordt algemeen als de ontdekker van X-stralen beschouwd. Hij is ze systematisch gaan onderzoeken in 1895 en het leverde hem in 1901 de Nobelprijs voor de fysica op. Effecten die aan X-stralen konden toegeschreven worden, werden echter reeds vroeger door andere wetenschappers vastgesteld. X-stralen zijn een vorm van elektromagnetische straling met een golflengte λ kleiner dan ongeveer 10 nm. De golflengte is gerelateerd aan de energie volgens E=

hc λ

(2.1)

met h de constante van Planck en c de lichtsnelheid. X-stralen zijn wegens hun korte golflengte sterk penetrerend. Ze worden typisch gebruikt voor radiografie, kristallografie. . . De naam wordt vaak door elkaar gebruikt met γ-straling. De conventie is echter dat γstraling veroorzaakt wordt door het verval van een instabiele kern of door annihilatie processen tussen materie en anti-materie. De term X-straling wordt gebruikt voor de elektromagnetische straling uitgezonden door geladen deeltjes bij verandering van energieniveau in een atoom (karakteristieke X-stralen) of bij het vertragen in een Coulomb krachtveld

4

2.1. X-stralen

5

(remstraling)1 . [10]

2.1.1

Productie van X-stralen

Synchrotron Een lading die een versnelling ondergaat, zendt elektromagnetische straling uit. Als een lading q met snelheid v in een magnetisch veld B beweegt, zal er een kracht uitgeoefend worden op de lading volgens F = q(v × B). De lading is hierdoor aan een versnelling onderhevig en zendt straling uit. Dit principe wordt gebruikt in een synchrotron en men spreekt dan ook over synchrotronstraling. [11, 12] X-stralenbuis Als een lading q in het veld van een atoomkern beweegt, ondervindt ze een kracht F = qE en is ze aan een versnelling onderhevig. De straling die hierbij uitgezonden wordt noemt men bremsstrahlung, of ook wel remstraling. Dit principe wordt gebruikt bij een Xstralenbuis (Engels: X-ray Tube (XRT)). Een filament wordt verhit en elektronen worden door thermionische emissie uitgezonden. De elektronen worden versneld o.i.v. een elektrisch veld en treffen een trefplaatje. Dit is vervaardigd uit materialen met een hoog atoomnummer Z (bv. koper, wolfraam, molybdeen. . . ), vermits het radiatief stoppend vermogen hiervoor groter is. Het grootste deel van de kinetische energie wordt omgezet in warmte en slechts een kleine fractie in straling. [13, 14]

2.1.2

Interactie met materie

Wanneer een X-stralenbundel door een materiaal gaat, wordt deze geattenueerd2 . De wet van Lambert-Beer geeft de intensiteit van een bundel na doorgang door het object: I = I0 e −

R pad

µ(λ,s) ds

.

(2.2)

Hierin is I0 de intensiteit v´ o´ or de doorgang van de bundel door het object en µ(λ, s) is 3 de attenuatiecoëfficiënt . De integraal wordt over het pad door het object genomen. Dit attenuatiebeeld is wat klassiek gebruikt wordt in radiografie en tomografie. Verschillende fenomenen geven een bijdrage tot de attenuatie [15]: Foto-elektrische absorptie: hierbij wordt het foton geabsorbeerd en wordt een fotoelektron gecreëerd. Het elektron heeft een energie gelijk aan de fotonenergie verminderd met de bindingsenergie die het had in het atoom. De werkzame doorsnede 1

Oudere literatuur maakt soms een onderscheid op basis van golflengte, maar tegenwoordig maakt men een onderscheid op basis van hun oorsprong. 2 Naast attenuatie treedt er ook een faseverschuiving op. Dit wordt later uitvoerig besproken. 3 De attenuatiecoëfficiënt µ(λ, s) is afhankelijk van de golflengte λ van de fotonen en is afhankelijk van de positie s in het object.

2.2. Computertomografie

6

is evenredig met Z n /Eγ3.5 , met n gelegen tussen 4 en 5. Het gevormde gat wordt opgevuld door een elektron uit een hogere schil dat terugvalt. De energie komt vrij onder de vorm van karakteristieke X-stralen of Augerelektronen. Compton verstrooiing: hierbij verstrooit een foton aan een elektron van het absorberend materiaal. Bij de verstrooiing wordt energie van het foton naar het elektron overgedragen. De werkzame doorsnede neemt toe naarmate er meer elektronen zijn om aan te verstrooien (ze neemt dus toe met Z) en neemt algemeen af voor stijgende fotonenergie. Rayleigh verstrooiing: hierbij verstrooit een foton aan de elektronenwolk van het materiaal, het is een coherente verstrooiing. Er is geen overdracht van energie van het foton naar de elektronen, maar het foton verandert wel van richting. Deze vorm van verstrooiing is enkel significant voor lage fotonenergieën (lager dan enkele honderd keV) en neemt toe met het atoomnummer Z. Elektron-positron paarvorming: in de buurt van de kern kunnen fotonen met een energie groter dan tweemaal de rustenergie van een elektron omgezet worden in een elektron en een positron. Dit komt niet voor bij de fotonenergieën die beschouwd worden in deze thesis.

De belangrijkste bijdragen tot de attenuatie van de bundel, door verstrooiing van fotonen of absorptie van fotonen uit de bundel, zijn de fotoëlektrische absorptie en comptonverstrooiing. [2, 15, 16]

2.2

Computertomografie

Bij een radiografie krijgen we een tweedimensionale (2D) projectie van een 3D-object. In X-stralen computertomografie wordt deze derde dimensie terug gereconstrueerd door het object te laten roteren en een reeks projecties te nemen gespreid over verschillende hoeken. In wat volgt bespreken we het proces dat doorlopen wordt om een CT-scan te maken en de wiskundige achtergrond van het reconstructie algoritme.

2.2.1

CT-scan

De opstelling is schematisch geschetst in figuur 2.1. Het object wordt opgesteld tussen de X-stralenbron en de detector. We maken gebruik van een conische bundel. De bronobject afstand (Engels: source-object distance (SOD)) en bron-detector afstand (Engels: source-detector distance (SDD)) bepalen de vergroting M : M=

SDD . SOD

Een CT-scan kan dan als volgt samengevat worden:

(2.3)


7

Figuur 2.1: Schematische voorstelling van een opstelling voor X-stralen CT

Scan proces Voordat de reeks projecties genomen wordt, neemt men eerst een reeks dark field images en flat field images. Deze zijn respectievelijk beelden van een niet-belichte detector en een belichte detector waarbij geen object aanwezig is. Men kan ook altijd gebruik maken van dark field images en flat field images die vroeger genomen zijn (of uit een database), maar om de beste kwaliteit te bekomen neemt men deze best elke keer opnieuw. Het object wordt op de gewenste positie tussen bron en detector gepositioneerd. Het vooropgegeven aantal projecties wordt equidistant over 360◦ genomen.

Verwerking De genormeerde beelden Ni worden uit de projecties Pi bekomen volgens

Ni =

Pi − D . F −D

(2.4)

Hierin stellen D en F respectievelijk de dark field images en flat field images voor. Het reconstructie algoritme heeft — om één doorsnede te reproduceren — de data nodig over de verschillende rotatiehoeken. Daarom wordt vaak vooraf alle data van één slice – zijnde een bepaalde rij in alle projecties — onder elkaar geplaatst. Dit beeld noemt men een sinogram. Voor elke slice wordt dus een sinogram gemaakt en deze sinogrammen dienen dan als input voor het reconstructie algoritme. Als alle doorsnedes gereconstrueerd zijn, wordt het 3D-beeld bekomen door de slices op elkaar te “stapelen”.

In wat volgt bespreken we het algoritme dat gebruikt wordt om de doorsnedes te reconstrueren uit de projectiedata.


2.2.2

8

Analytische reconstructie

We bespreken kort de wiskundige achtergrond van het reconstructiealgoritme4 [2]. Om de gedachten te vestigen wordt het ge¨ıllustreerd voor een parallelle bundel. Voor een conische bundel is de uitwerking uitgebreider maar gebaseerd op hetzelfde principe. Het meest gebruikte algoritme is het filtered backprojection algorithm (FBP). Dit maakt gebruik van het Fourier Slice Theorema dat grafisch ge¨ıllustreerd is in figuur 2.2. Theorema 1. Stel dat Pθ (t) de parallelle projectie is van een objectfunctie f (x, y) genomen onder een rotatiehoek θ. Dan is de 1D-Fouriergetransformeerde van Pθ (t) gelijk aan de waarde van de 2D-Fouriergetransformeerde van f (x, y) langs de rechte BB in de frequentieruimte.

Figuur 2.2: Grafische illustratie van het Fourier Slice Theorema [2]

Voor X-stralen CT is de functie f (x, y) de attenuatiecoëfficiënt µ(x, y). De parallelle projectie Pθ (t) is dan de integraal van µ(x, y) langs het pad van de straal. Deze ge¨ıntegreerde waarde bekomt men door vergelijking 2.2 om te vormen tot Z Iθ (t) . (2.5) µ(λ, s) ds = Pθ (t) = − ln Iθ,0 (t) pad Het Fourier Slice Theorema legt dus het verband tussen de Fouriergetransformeerde (Engels: Fourier transform (FT)) van de projectie van ons object, en de FT van de objectfunctie. Elke projectie onder een hoek θ geeft dan een rechte van de functie f˜(ξ, η), met f˜ de FT van de functie f 5 . Praktisch zijn we beperkt tot een eindig aantal projecties, waardoor we slechts een eindig aantal rechten kennen, zoals weergegeven in figuur 2.3. De waarden ertussen zou men kunnen benaderen door interpolatie. De objectfunctie kan dan teruggevonden worden door de inverse FT te nemen van f˜(ξ, η). 4

Een uitgebreide bespreking en afleiding van de formules is terug te vinden in het boek van Kak en Slaney [2]. Het besproken algoritme is een analytische reconstructietechniek. Er zijn ook iteratieve reconstructiemethodes, waarvan we het principe zullen uitleggen in sectie 2.2.4. 5 We zullen volgende conventie hanteren voor de notatie van een FT: F {f (x, y)} = f˜(ξ, η). De co¨ ordinaten in de frequentieruimte duiden we aan met (ξ, η).


9

Figuur 2.3: Door toepassen van het Fourier Slice Theorema kunnen we rechten van de FT van de objectfunctie berekenen [2]

Het FBP algoritme maakt gebruik van dit theorema en kan zeer efficiënt ge¨ımplementeerd worden. Door over te gaan op poolcoördinaten kan men aantonen dat de objectfunctie gegeven wordt door Z π n o f (x, y) = dθF −1 P˜θ (ξ)|ξ| . (2.6) 0

F −1

Hierin is de inverse FT. De naam filtered backprojection kan als volgt begrepen worden. ˜ De functie Pθ (ξ) is een rechte in de Fourier ruimte. Om alle punten hier te reconstrueren zouden we geen rechte, maar een dubbele wigstructuur willen zoals in figuur 2.4a. Dit wordt benaderd door de filterfunctie |ξ| op P˜θ toe te passen: elke frequentie krijgt een gewicht mee evenredig met de frequentie. Vervolgens wordt een deel van de slice gereconstrueerd door het terug te projecteren (dit is de inverse FT nemen).

Figuur 2.4: Data in het frequentiedomein. (a) ideaal geval, (b) de rechte die bekomen wordt volgens het Fourier Slice Theorema, (c) de benadering door een weging met de frequentie [2]

De reconstructies in deze thesis werden uitgevoerd met Octopus [17]. Dit is een reconstructie softwarepakket ontwikkeld aan UGCT dat reconstructies toelaat voor verschillende bundelgeometrieën. We zullen enkel gebruik maken van de cone-beam geometrie. Het algoritme dat daarvoor gebruikt wordt is het FDK-algoritme, genoemd naar de ontwikkelaars Feldkamp, David en Kress. [2]


2.2.3

10

Praktisch

Cone-beam artefacten Zoals reeds vermeld zal in dit werk de conische bundel geometrie gebruikt worden. Het FDK-algoritme is slechts een benadering en zal — op de centrale slice na in het vlak van de bron — niet volledig correct zijn. Dit kan als volgt begrepen worden. Een slice kan exact gereconstrueerd worden als alle stralen door dit vlak gaan en op een rij pixels van de detector komen (zoals bij een parallelle bundel). Bij een conische bundel maakt een straal — op de centrale slice na — steeds een hoek met de te reconstrueren doorsnedes van het object. Men zou dit artefact kunnen reduceren door de SOD en ODD (object-detector distance) te vergroten zodat de hoek van de conische bundel kleiner wordt. Praktisch is dit niet mogelijk vermits de fotonflux met 1/r2 afneemt en de flux die de detector zou bereiken te klein zou zijn. [2, 16] In dit werk hebben we, als we slices vergelijken, steeds een centrale slice gebruikt om deze artefacten te vermijden. Beam hardening Een X-stralenbundel afkomstig van een XRT is polychromatisch. De attenuatiecoëfficiënt µ is afhankelijk van de energie en neemt af met toenemende fotonenergie. De laag energetische fotonen zullen dus sterker geattenueerd worden dan de hoog energetische fotonen. Het gevolg is dat de uitgaande bundel een ander energiespectrum zal hebben dan de invallende bundel (de laag energetische fotonen zullen door de sterke attenuatie in het slechtste geval helemaal verdwenen zijn uit het spectrum). Dit noemt men beam hardening. Het gevolg is dat de gereconstrueerde attenuatiecoëfficiënt aan de rand groter zal zijn dan in het centrum. Dit noemt men een cupping artefact en is ge¨ıllustreerd in figuur 2.5. [16] Men kan dit artefact reduceren door een filter te plaatsen (bv. een loodplaatje waarvan men de dikte kan kiezen) die de laag energetische fotonen uit het spectrum filtert. Een correctie kan ook toegepast worden in de reconstructiesoftware.

2.2.4

Algebra¨ısche reconstructie

In sectie 2.2.2 hebben we besproken hoe analytische reconstructie werkt. Naast dit type reconstructie bestaan er ook algebra¨ısche reconstructie algoritmes. In wat volgt leggen we het basisprincipe hiervan uit6 . De 3D objectfunctie f (x, y, z) die we willen reconstrueren, kunnen we onderverdelen in N 6

Gebaseerd op het boek van Kak en Slaney [2].


11

Figuur 2.5: Reconstructie van een homogeen fantoom met acht holtes, (a) met beam hardening, (b) gecorrigeerd voor beam hardening [16]

voxels7 . Men kan voor elke detectorpixel i schrijven: N X

wij fj = pi , i = 1, 2 . . . , M

(2.7)

j=1

Hierin is pi de negatieve logaritme van de genormeerde projectiedata op pixel i (zie vergelijking 2.5). De som gaat over de N voxels en M is het totaal aantal pixels (dus als n het aantal projecties voorstelt en Mp het aantal pixels per projectie, dan is M = n · Mp ). De factor wij stelt het gewicht voor waarmee de voxel j bijdraagt tot de intensiteit op pixel i. Als de bundel enkel attenuatie ondergaat, is f (x, y, z) de attenuatiecoëfficiënt en de som in vergelijking 2.7 is de lijnintegraal van µ langs het pad van de straal in het object. We bekomen dus een set van M vergelijkingen die we in een matrixvergelijking kunnen schrijven als w · f = p. (2.8) Doordat deze matrix zeer groot is, moet een iteratieve methode gebruikt worden om dit stelsel op te lossen. Vandaar dat deze methodes ook vaak iteratieve reconstructie algoritmes genoemd worden. De Kaczmarz methode stelt de objectfunctie (f1 , f2 , . . . , fN ) voor als een punt in een N-dimensionale ruimte, dat gegeven wordt door de doorsnede van de hypervlakken bepaald door vergelijking 2.8. We illustreren hoe men de doorsnede van deze 7

Een voxel is een driedimensionale pixel.


12

hypervlakken kan vinden in het geval van twee rechten, gegeven door w11 f1 + w12 f2 = p1

(2.9)

w21 f1 + w22 f2 = p2 .

(2.10)

Men start met een beginwaarde f (0) voor de objectfunctie. Dit punt projecteert men loodrecht op de eerste rechte. Vervolgens wordt het geprojecteerde punt loodrecht op de tweede rechte geprojecteerd. Dit is ook weergegeven in figuur 2.6. Men herhaalt dit tot men zo dicht mogelijk tot het snijpunt komt. In het N -dimensionaal geval wordt de objectfunctie f (0) geprojecteerd op het eerste hypervlak (gegeven door de eerste rij uit de matrixvergelijking 2.8) wat f (1) oplevert. Dit punt wordt geprojecteerd op het tweede hypervlak (gegeven door de tweede rij uit de matrixvergelijking 2.8) wat f (2) geeft. Dit wordt herhaald voor alle M vlakken. Algemeen kan men dit wiskundig uitdrukken als

Figuur 2.6: De Kaczmarz methode om een stelsel van vergelijkingen op te lossen [2]

f (i) = f (i−1) +

pi − f (i−1) · wi wi . wi · wi

(2.11)

met wi = (wi1 , wi2 , . . . , wiN ). Nadat alle M vergelijkingen uit 2.7 doorlopen zijn (en dus M projecties uitgevoerd zijn) kan men opnieuw beginnen tot een zekere convergentie bereikt is. Praktisch kunnen we het algoritme als volgt samenvatten:


13

1. Men berekent de term f (i−1) · wi , wat de lijnintegraal van de objectfunctie langs het pad van de straal is. In het geval dat er enkel attenuatie van de bundel optreedt, is de objectfunctie µ. Essentieel berekent men hier dus het beeld die het huidig virtueel object vormt. Men noemt deze stap dan ook de voorwaardse projectie. 2. Het verschil met de echte waarde pi wordt bepaald en de correctie (de tweede term uit het rechter lid van vergelijking 2.11) op de huidige objectfunctie wordt berekend. 3. De oude objectfunctie wordt aangepast tot de nieuwe f (i) . Dit vormt de basis voor verschillende iteratieve reconstructie algoritmes. Het komt er dus op neer dat een virtueel object steeds ge¨ updatet wordt. Het iteratieve reconstructie algoritme SART (Simultaneous Algebraic Reconstruction Technique) werd aan UGCT ge¨ımplementeerd [16]. Het voordeel van dit type reconstructie t.o.v. analytische reconstructie is dat ook nog een goede reconstructie bekomen wordt bij een beperkt aantal projecties of wanneer projecties bij bepaalde hoeken ontbreken. Bovendien kan men in de voorwaardse projectie ook andere factoren in rekening brengen zoals verstrooiing, beam hardening. . . Het is dus duidelijk dat, indien men het fasefenomeen — dat in hoofdstuk 3 uitgegelegd zal worden — correct kan simuleren en in de voorwaardse projectie in rekening kan brengen, men artefacten t.g.v. het fasesignaal kan reduceren. [2, 16]

Hoofdstuk 3

Fasecontrast 3.1

Inleiding

Zoals reeds vermeld wordt een X-stralenbundel geattenueerd als deze door een object gaat. Uit de projecties die men bekomt kan men de attenuatiecoëfficiënt µ reconstrueren. Sommige materialen echter vertonen een lage attenuatie (bv. organisch weefsel), en het verschil in attenuatie van verschillende materialen is soms zeer klein. Naast attenuatie treedt echter nog een tweede fenomeen op: de X-stralenbundel ondergaat een faseverschuiving wanneer hij door het materiaal gaat. Dit wordt weergegeven in figuur 3.1. Wanneer we deze faseverschuiving kunnen detecteren en gebruiken in de reconstructie, komen we bij fasecontrastbeeldvorming terecht. Het voordeel van fasecontrast t.o.v. attenuatie is dat sommige materialen die weinig attenuatie vertonen, wel een significante faseverschuiving kunnen induceren. Ook zal bij hogere energieën de attenuatie lager zijn (en dus ook de stralingsdosis) terwijl er nog een beduidend fasesignaal is. [18, 19, 20] Het is duidelijk dat de faseverschuiving kan gebruikt worden bij beeldvorming van weinig attenuerende objecten (bv. bij mammografie [21, 22, 23] en andere organische materialen [8]). Vermits de techniek niet op attenuatie gebaseerd is, kan de geabsorbeerde dosis sterk verlaagd worden. Bovendien laat het toe, indien het signaal juist ge¨ınterpreteerd wordt, kleinere details te zien [24]. We kunnen uiteraard niet kiezen of er al dan niet refractie optreedt. Ook wanneer we geen gebruik willen maken van het fasesignaal (en dus enkel een reconstructie op basis van de attenuatie willen maken) kan dit sterk aanwezig zijn. Dit zorgt voor artefacten in de reconstructie die, indien niet correct ge¨ınterpreteerd, tot verkeerde conclusies kunnen leiden. In wat volgt leggen we het fenomeen fasecontrast uit. We zullen dit eerst doen aan de hand van de stralenoptica, die een duidelijk beeld geeft van de fysica. Daarna bespreken we de 14

3.2. Stralenoptica

15

meer geavanceerde golfoptica.

Figuur 3.1: Attenuatie en faseverschuiving van een elektromagnetische golf na doorgang door een materiaal met dikte d [19]

3.2

Stralenoptica

In de stralenoptica worden de golffronten (dit zijn vlakken van constante fase) voorgesteld door lijnen loodrecht op de stralen. De stralen zijn dus parallel met de voortplantingsrichting. Dit is een zeer belangrijk model in de fysica en het is toepasbaar zolang de golflengte kleiner is dan de structuren waarmee de golf interageert [25]. Op basis hiervan bespreken we het fenomeen fasecontrast en leiden we ook een vergelijking af die het verband legt tussen de faseverschuiving veroorzaakt in het object en de gemeten intensiteit op de detector.

3.2.1

Complexe brekingsindex

Een object kan gekarakteriseerd worden door een 3D distributie van zijn complexe brekingsindex, gedefinieerd als n(x, y, z) = 1 − δ(x, y, z) + iβ(x, y, z).

(3.1)

Hierin is δ verantwoordelijk voor de faseshift en men noemt δ het refractive index decrement. Deze waarde is voor X-stralen zeer klein. Het imaginair deel β noemt men de absorption index en veroorzaakt de attenuatie. De golf na doorgang door het materiaal wordt gegeven door U (x, y, z) = Ui (x, y, z)e−B(x,y,z)+iφ(x,y,z) . (3.2) Hierin is

en

2π B(x, y, z) = λ

Z

2π φ(x, y, z) = − λ

Z

β(x, y, z) dl

(3.3)

δ(x, y, z) dl.

(3.4)

pad

pad

Dit is dus een uitbreiding op de wet van Lambert-Beer (zie vergelijking 2.2) waar een faseverschuiving is toegevoegd. De parameters δ en β kunnen ook berekend worden volgens

3.2. Stralenoptica

16

[26] δ=

re λ2 ni f1 2π

(3.5)

en

re λ2 ni f2 . (3.6) 2π Hierin zijn f1 en f2 de atomaire verstrooiingsfactoren voor voorwaardse verstrooiing, re is de klassieke elektronenstraal en ni de atomaire dichtheid van het element. De waarden van f1 en f2 zijn gecatalogeerd (zie bv. [27]) en een voorbeeld is weergegeven in figuur 3.2. Dit geeft ons een gevoel voor de grootteorde van de parameters. De parameter δ is typisch van de orde 10−6 en β van de orde 10−10 voor een fotonenergie rond 20 keV. We merken ook op dat bij hogere energieën f2 sterk afneemt, terwijl f1 nagenoeg constant blijft. De attenuatie daalt dan ook veel sterker dan de faseverschuiving bij hogere energieën. β=

a)

b) 10

10

f2 (e/atoom)

f1 (e/atoom)

1

f1

1

0.1 f2

0.01 0.001 0.0001

0.1 0.01

0.1

1 10 energie (keV)

1e-005 0.01

100

c)

1 10 energie (keV)

100

d) 10

100

1

1

0.1

0.01

0.01 0.001

β (rad-1)

δ (rad-1)

0.1

δ

0.0001 1e-005

0.0001 1e-008

1e-006

1e-010

1e-007

1e-012

1e-008 0.01

0.1

1 10 energie (keV)

100

β

1e-006

1e-014 0.01

0.1

1 10 energie (keV)

100

Figuur 3.2: Data van koolstof [27]. (a) De atomaire verstrooingsfactor f1 , (b) de atomaire verstrooiingsfactor f2 , (c) het refractive index decrement δ en (d) de absorptie-index β

3.2. Stralenoptica

3.2.2

17

Refractie

Wanneer een golf door een materiaal beweegt, dan zal typisch de fasesnelheid van de golf (dit is de snelheid waarmee de golffronten zich voortplanten) in het materiaal verschillend zijn van de fasesnelheid in het vacu¨ um. Inderdaad, de reële brekingsindex nr is de verhouding van de lichtsnelheid c en de fasesnelheid v in het materiaal: nr =

c . v

(3.7)

Zoals reeds gezien is nr = 1 − δ en voor X-stralen een weinig kleiner dan 1. De fasesnelheid in het medium is dus een klein beetje groter dan in het vacu¨ um. Van zodra een gradiënt in de faseverschuiving loodrecht op de voortplantingsrichting aanwezig is1 , zal er refractie optreden. Er ontstaat een piek in de intensiteit waar de afgebogen straal terechtkomt en een dal op de plaats waar de straal had moeten terechtkomen als hij niet was afgebogen. Dit is ge¨ıllustreerd in figuur 3.3. [19]

Figuur 3.3: In een object is een gradiënt aanwezig in de faseverschuiving, loodrecht op de voortplantingsrichting van de golf. Er treedt refractie op [28]

Het is duidelijk dat dit effect sterk aanwezig zal zijn wanneer een overgang optreedt tussen gebieden met een sterk verschillende δ. Dit is bijvoorbeeld het geval aan de rand van een object. Er is een grote fasegradiënt en bijgevolg een sterke afbuiging. Dit veroorzaakt een edge enhancement: een piek — en zoals hiervoor besproken dus ook een dal — in de intensiteit wordt waargenomen2 . We zullen in sectie 3.2.3 zien dat dit fasesignaal evenredig is met de tweede afgeleide van de faseverschuiving. In figuur 3.4 zien we dat de tweede afgeleide van φ aan de rand divergent is: er treedt een versterking van de intensiteit op aan de rand. [28] Wanneer een polychromatische bundel gebruikt wordt, zullen de verschillende frequentiecomponenten anders afgebogen worden. Het gevolg is dat de piek breder wordt. Dit is 1

De faseverschuiving van de golf na doorgang door het materiaal wordt gegeven door vergelijking 3.4. Het verband tussen de gemeten intensiteit en de gradiënt van de faseverschuiving wordt gegeven in sectie 3.2.3. 2 In de golfoptica zullen we zien dat we eigenlijk een diffractiepatroon krijgen aan de rand van het object doordat de afgebogen stralen interfereren met de niet-afgebogen stralen. Door de grootte van de spot aan een XRT (orde µm), neemt men niet meer dit diffractiepatroon waar, maar eerder een edge enhancement zoals hier intu¨ıtief uitgelegd in de stralenoptica.

3.2. Stralenoptica

18

Figuur 3.4: Aan de rand van het object is de tweede afgeleide van de faseverschuiving φ divergent: er treedt een edge enhancement op [28]

geschetst in figuur 3.5. [28]

Figuur 3.5: Fasesignaal van een polychromatische bundel: de verschillende frequentiecomponenten worden anders afgebogen en we krijgen een bredere piek [28]

In figuur 3.6 is het edge enhancement getoond van een epoxyhars staafje (een object dat we nog zullen gebruiken in dit werk). Op de projectie (figuur 3.6a) is een witte lijn te zien aan de rand: dit is de piek in de intensiteit zoals hierboven uitgelegd. Dit wordt ook verduidelijkt met het lijnprofiel van de linker rand (figuur 3.6b). Het dal in de intensiteit is niet meer zichtbaar doordat de attenuatie domineert. De invloed op een gereconstrueerde doorsnede van het staafje is getoond in figuur 3.6c: we zien een abnormale stijging en daling in de dichtheid ten gevolge van het fasesignaal.

3.2.3

Verband tussen faseverschuiving en de gemeten intensiteit

Algemeen kan een verband opgesteld worden tussen de faseverschuiving veroorzaakt in het materiaal en de gemeten intensiteit. Dit wordt gegeven door de Transport of Intensity Equation (TIE)3 [29, 30]: 2π ∂ I(z) = −∇⊥ · [I(z)∇⊥ φ] . (3.8) λ ∂z Hierin is I(z) de intensiteit na propagatie over een afstand z en φ de faseshift veroorzaakt door transmissie door het object (zie vergelijking 3.4). De TIE legt dus een verband tussen 3

De afleiding is terug te vinden in het werk van Teague [29].

3.2. Stralenoptica

19

1.06 hars lijnprofiel 1.04

1.02

intensiteit

1

0.98

0.96

0.94

0.92

0.9

0.88 100

110

120

130

140

(a)

150 pixel

160

170

180

190

200

(b)

doorsnede

0.8

0.6

doorsnede

0.4

0.2

0

-0.2

100

(c)

200

300

400

500 pixel

600

700

800

900

(d)

Figuur 3.6: Illustratie van het edge enhancement bij een epoxyhars staafje. (a) Het edge enhancement is te zien als een witte lijn aan de rand. (b) Het lijnprofiel van de linker rand van het harsstaafje. (c) De invloed op een gereconstrueerde slice. (d) Een horizontale centrale doorsnede van de slice

de verandering in de intensiteit in de z-richting en de variatie van de fase in het vlak loodrecht op de voortplantingsrichting. In het geval dat I(z) traag varieert kunnen we dit vereenvoudigen tot λd 2 ∇ φ , (3.9) Id = I0 1 − 2π ⊥ en dit wordt de Transport of Intensity Equation for Weak attenuation (TIEW) genoemd. In wat volgt geven we een afleiding van de TIEW gebaseerd op de stralenoptica [31].

3.2. Stralenoptica

3.2.4

20

Afleiding van de TIEW

Beschouw een object dat op een afstand s van de bron staat en een detector op een afstand d van het object. Het objectvlak kiezen we als z = 0 en het detectorvlak staat op een afstand z = d. We duiden de overeenkomstige coördinaten respectievelijk aan met de subscripten 0 en d. De afleiding vertrekt van een parallelle bundel (later kan dit herschreven worden als er een vergroting M aanwezig is). Verwaarlozen we om te beginnen de attenuatie, dan wordt de faseverschuiving na doorgang door het object gegeven door vergelijking 3.4. De straal wordt lichtjes afgebogen over een hoek α en komt op de detector op positie xd ≈ x0 + dαx (x0 , y0 )

(3.10)

yd ≈ y0 + dαy (x0 , y0 ).

(3.11)

Hierin zijn αx en αy de projecties van α op de vlakken (x0 , z) en (y0 , z) en worden gegeven door λ ∂φ(x0 , y0 ) 2π ∂x λ ∂φ(x0 , y0 ) αy (x0 , y0 ) ≈ . 2π ∂y

αx (x0 , y0 ) ≈

(3.12) (3.13)

Uitdrukkingen 3.10 en 3.11 geven het verband tussen de coördinaten in beide vlakken. Hieruit kan afgeleid worden dat −1 1 + d ∂αx ∂αx d ∂x0 ∂y0 I(xd , yd ) = I0 (3.14) ∂α ∂α y y d ∂x 1 + d ∂y0 0 −1 dλ 2 ≈ I0 1 + ∇ φ(x0 , y0 ) . (3.15) 2π

De determinant in 3.14 is de jacobiaan en in 3.15 werden de producten van de afgeleiden verwaarloosd. Als (dλ/2π)∇2 φ(x0 , y0 ) << 1 dan kunnen we benaderend schrijven dλ 2 ∇ φ(x0 , y0 ) . (3.16) I(xd , yd ) ≈ I0 1 − 2π Als we de vergroting in rekening brengen, dan worden 3.10 en 3.11 vervangen door xd ≈ M x0 + dαx (x0 , y0 )

(3.17)

yd ≈ M y0 + dαy (x0 , y0 ).

(3.18)

3.3. Visualisatie van de fase Als we hiermee de vorige stappen herhalen bekomen we dλ 2 I(xd , yd ) ≈ I0 1 − ∇ φ(x0 , y0 ) , 2πM

21

(3.19)

h R i wat de TIEW is. Als er nog attenuatie aanwezig is, vervangen we I0 door I0 exp − pad µ dl . Het is interessant op te merken dat de intensiteit dus gegeven wordt door een product van de attenuatie en een faseterm. Deze vergelijking legt dus ook een eenvoudig verband tussen de intensiteit op de detector en de faseverschuiving. Dit zal gebruikt worden in verschillende fasereconstructie algoritmes [32, 33, 34].

3.3

Visualisatie van de fase

In onze bespreking hebben we verondersteld dat de golf na doorgang door het object verder propageert en uiteindelijk op de detector terechtkomt. Het fasefenomeen dat dan optreedt is een edge enhancement. De TIEW legt dan een verband tussen de gemeten intensiteit en de faseverschuiving in het object. Er zijn echter nog andere technieken om de faseverschuiving te achterhalen.

3.3.1

Interferometrie

Voor interferometrie heeft men een coherente bundel nodig en een hoge flux. Synchrotrons voldoen hier aan. De bundel moet in minstens twee coherente bundels gesplitst worden. We leggen het principe uit a.d.h.v. de LLL X-stralen interferometer (L staat voor Laueof transmissiegeometrie). Deze is weergegeven in figuur 3.7. Er worden drie kristallen gebruikt. In het eerste kristal (beam splitter S) wordt de bundel gesplitst, in het tweede kristal (mirror M) [35] worden beide bundels nog eens gesplitst. Het object wordt tussen de kristallen M en A geplaatst. De convergerende bundels worden gerecombineerd door het derde kristal (analyser A). Er ontstaat een interferentiepatroon dat gerelateerd is aan de faseshift van de bundel die door het object gegaan is. [35, 36]

Figuur 3.7: LLL X-stralen interferometer [35]

3.3. Visualisatie van de fase

3.3.2

22

Diffraction enhanced imaging

De opstelling voor deze techniek is geschetst in figuur 3.8. Een bundel valt in op een eerste kristal(len) zodat in het geval van een polychromatische bundel, een monochromatische bundel op het object gestuurd wordt. De SOD is groot zodat een nagenoeg parallelle bundel invalt (deze techniek is dus enkel mogelijk aan een synchrotron). De X-stralen ondergaan een faseverschuiving door transmissie door het object. Dit zorgt voor een afbuiging over een hoek λ ∂φ ∆α = − . (3.20) 2π ∂x Een tweede kristal (analyser crystal ) wordt gebruikt om deze hoekafwijking te bepalen. De intensiteit na doorgang door het analyser kristal is voor elke rotatiehoek gekend. Bij één oriëntatie van het object worden projecties gemaakt bij verschillende kristaloriëntaties, waaruit de faseverschuiving bepaald kan worden. [36, 37]

Figuur 3.8: Opstelling bij diffraction enhanced imaging [37]

3.3.3

Grating interferometrie

Soms duidt men deze techniek ook aan met diffraction enhanced imaging. Hier maakt men gebruik van minstens twee gratings4 G1 en G2 zoals getoond in figuur 3.9. Bij doorgang van de bundel door het object treedt refractie op. De transmissie door de roosters is sterk afhankelijk van hun onderlinge positie, maar ook van de invalshoek. Als na doorgang door het object de invalshoek door refractie verandert, zal een ander intensiteitspatroon gemeten worden. Uit de afbuigingshoek kan de faseshift bepaald worden. Het voordeel van deze techniek is dat geen coherente bron nodig is en dat de techniek dus kan gebruikt worden met een XRT [39]. In dat geval wordt een derde rooster G0 na de bron geplaatst (zie figuur 3.9b). Op die manier creëert men als het ware lijnen van coherente bronnen. [18, 40]

3.3.4

Vrije propagatie

De opstellingen tot nu toe besproken zijn eerder complex en sommige leggen strikte voorwaarden op (coherente straling, grote SOD). Wanneer de bundel door het object gaat treedt refractie op. Als we de uittredende bundel verder later propageren, zal deze op een 4

Een methode met één grating werd voorgesteld door Krejci et al. [38]

3.3. Visualisatie van de fase

23

Figuur 3.9: Opstelling bij grating interferometrie. Na de bron worden twee gratings geplaatst. (a) Coherente bron, (b) incoherente bron: er wordt nog een extra grating G0 geplaatst na de bron [40]

afstand z > 0 interfereren met de niet-afgebogen bundel. Afhankelijk van de afstand z onderscheiden we verschillende regio’s, zoals weergegeven in figuur 3.10. Op een afstand < d2 /λ bevinden we ons in z = 0 hebben we een zuiver absorptiebeeld. Op een afstand z ∼ het Fresnel gebied. In het Fresnel gebied noemen we z < d2 /λ het nabije veld (near field ), of ook wel het edge detection regime. Op deze afstand zijn aan de randen van het object (waar de laplaciaan van de faseshift φ divergent is, zie vergelijking 3.19) een piek en een dal in de intensiteit te zien. Een maximaal contrast wordt bekomen op een afstand z = d2 /2λ [24]. Op een zeer grote afstand (z > d2 /λ) van het object komen we in het Fraunhofer gebied terecht. Het beeld is dan de FT van de golf die door het object is gegaan. [36, 37, 41] Deze techniek kan gebruikt worden bij een XRT. De ruimtelijke coherentie van de bundel moet wel voldoende groot zijn. Een XRT heeft een eindige spotgrootte. Al deze punten van de spot zenden incoherente straling uit: de fase tussen de uitgezonden stralen van de verschillende punten is random. Een maat voor de coherentie is de transverse coherence length dT C [36] λ∆ dT C = . (3.21) 2σ Hierin is λ de gebruikte golflengte, ∆ de SOD en σ de spotgrootte. Bij synchrotrons is ∆ groot, terwijl deze bij een XRT klein is. De transversale coherentielengte is bijgevolg veel groter aan een synchrotron dan aan een XRT5 . Als σ klein genoeg is, is dT C toch nog voldoende voor fasecontrast aan een XRT. Bovendien is de polychromatische bundel geen 5

Aan een synchrotron is dT C typisch van de orde 100 µm, bij een XRT aan UGCT van de orde 1 µm [16, 36]

3.4. Golfoptica

24

probleem bij fasecontrastbeeldvorming door vrije propagatie6 . [16, 28, 36]

Figuur 3.10: De verschillende gebieden voorbij het object bij vrije propagatie. Het object heeft een grootte d. De projecties van een fantoom op verschillende afstanden zijn weergegeven [41]

3.4

Golfoptica

Zoals reeds vermeld zal het fasesignaal zichtbaar worden door vrije propagatie van de bundel na doorgang door het object. De stralenbenadering, die een goed inzicht geeft in de fysica van het fenomeen, werd reeds behandeld. De berekening van de intensiteit kan ook vanuit de meer geavanceerdere golfoptica gebeuren. We gaan iets dieper in op de theorie die de propagatie van een monochromatische vlakke golf beschrijft7 en zullen zo trachten de intensiteit op onze detector te berekenen.

3.4.1

Lineaire shift-invariante systemen

Onderstel een systeem waarvan de inputfunctie gegeven wordt door f (t) en de output door g(t). Een systeem is lineair als de respons van de som a + b gelijk is aan de som van de responsen van a en b. Stel dat de output g(t) op tijdstip t gegeven wordt door een gewogen superpositie van f (τ ) op verschillende tijdstippen τ dan kan dit geschreven worden als Z +∞ g(t) = h(t, τ )f (τ ) dτ, (3.22) −∞

met h(t, τ ) het gewicht waarmee f (τ ) tot g(t) bijdraagt. De functie h(t, τ ) wordt de impulse response van het systeem genoemd. Het systeem is shift invariant op voorwaarde dat, wanneer men de input over ∆t verschuift, de output ook over dit tijdsinterval verschuift. De responsefunctie h(t, τ ) wordt dan een functie h(t−τ ) die enkel afhangt van het tijdsverschil. Vergelijking 3.22 wordt dan Z +∞ g(t) = h(t − τ )f (τ ) dτ. (3.23) −∞ 6 7

Een polychromatische bundel zorgt enkel voor een verbreding van de piek [28]. De afleiding is grotendeels gebaseerd op het werk van Saleh en Teich [42].

3.4. Golfoptica

25

Dit is een convolutie van f (t) met de impulsresponsfunctie h(t). Een convolutie in het tijdsdomein is een gewone vermenigvuldiging in het frequentiedomein. Vergelijking 3.23 kan dus geschreven worden als ˜ f˜(ν), g˜(ν) = h(ν) (3.24) ˜ Waar g˜(ν), h(ν) en f˜(ν) respectievelijk de Fouriergetransformeerden zijn van g(t), h(t) en ˜ f (t). De functie h(ν) wordt de transferfunctie van het systeem genoemd. Het komt er dus op neer dat, wanneer men de FT neemt van de input, we na vermenigvuldiging met de transferfunctie de FT van de output vinden. Toepassing in optica Onderstel een vlakke golf met complexe amplitude U (x, y, z) die door een optisch systeem gaat, zoals weergegeven in figuur 3.11. Dit kan gezien worden als een lineair shift-invariant systeem met input f (x, y) = U (x, y, z = 0) en output g(x, y) = U (x, y, z = d). Het optisch systeem kan een lens zijn, maar ook vrije propagatie. Als we de invallende golf ontbinden in zijn frequentiecomponenten (i.e. de FT nemen) dan kunnen we met formule 3.24 de golf in het vlak z = d berekenen. Inderdaad, de transferfunctie geeft de respons op elke frequentie. We moeten wel de transferfunctie van het systeem kennen. In wat volgt leiden we de transferfunctie voor vrije propagatie af.

Figuur 3.11: Een golf U (x, y, z) gaat door een optisch systeem. De mapping met een lineair shift-invariant systeem gebeurt door f (x, y) = U (x, y, 0) en g(x, y) = U (x, y, d) [42]

3.4.2

Vrije propagatie als lineair shift-invariant systeem

Het optisch systeem uit figuur 3.11 vervangen we door vrije propagatie van de vlakke golf, zoals weergegeven in figuur 3.12. Vermits we een lineair shift-invariant systeem hebben8 kan men de amplitude g(x, y) = U (x, y, z = d) in het vlak z = d berekenen als de complexe amplitude f (x, y) = U (x, y, z = 0) in het vlak z = 0 gekend is. Men kan aantonen dat de 8

Dit wordt uitgelegd in het werk van Saleh en Teich [42].

3.4. Golfoptica

26

Figuur 3.12: Een golf U (x, y, z) plant zich voort tussen de vlakken z = 0 en z = d. Het optisch systeem uit figuur 3.11 is hier vervangen door vrije propagatie [42]

complexe amplitude U (xd , yd ) in het vlak z = d gegeven wordt door [36, 42, 43, 44] Z Z 1 eikr0d cos(n, r0d ) dx0 dy0 . (3.25) U (xd , yd ) = U (x0 , y0 ) iλ r0d Dit is de Fresnel-Kirchoff diffractieformule. Hierin is r0d de vector van een punt in het z = 0 vlak tot een punt in het z = d vlak. De integraal gaat over het z = 0 vlak en de vector n is de normaal op dit vlak. Fresnel benadering Onder de paraxiale benadering zijn de diffractiehoeken klein (< 1 mrad, de straal wordt weinig afgebogen) en is cos(n, r0d ) ≈ 1. In de Fresnel benadering kan r0d ontwikkeld worden (men veronderstelt de transversale dimensies veel kleiner dan de propagatie afstand d [43]) en vergelijking 3.25 wordt Z Z eikd iπ 2 2 U (x0 , y0 ) exp (xd − x0 ) + (yd − y0 ) dx0 dy0 . (3.26) U (xd , yd ) = iλd λd We herkennen hierin de convolutie U (x, y) = U (x, y) ∗ hF (x, y).

(3.27)

Hierin is hF de Fresnel propagator : eikd iπ 2 2 hF (x, y) = exp x +y . iλd λd

(3.28)

Analoog aan vergelijking 3.24 kunnen we dit berekenen in het Fourier domein, waarbij de convolutie een vermenigvuldiging wordt: F{U (xd , yd )} = F{U (x0 , y0 )} · exp(ikd) exp −iπλd(ξ 2 + η 2 ) . Hierin stelt F de FT voor en werd de FT van hF expliciet berekend.

(3.29)

3.4. Golfoptica

27

Fresnel propagator in de simulatie van fasecontrast In vorige sectie hebben we gezien dat, gegeven een complexe golfamplitude in het objectvlak, de golfamplitude in het detectorvlak na vrije propagatie berekend kan worden met de Fresnel propagator. We passen dit nu toe op een golf na doorgang door een materiaal met complexe brekingsindexdistributie n(x, y, z) = 1 − δ(x, y, z) + iβ(x, y, z). De golfamplitude in het vlak net na het sample (dit is het z = 0 vlak) wordt gegeven door Uz=0 (x, y) = T (x, y)Ui (x, y),

(3.30)

waarin Ui de invallende vlakke golf voorstelt en T de transmissie functie. T wordt gegeven door [24] T (x, y) = M (x, y)eiφ(x,y) . (3.31) De functies M (x, y) en φ(x, y) beschrijven respectievelijk de attenuatie en de faseverschuiving en worden gegeven door Z 1 M (x, y) = exp − µ(x, y, z) dz (3.32) 2 pad en

2π φ(x, y) = − λ

Z δ(x, y, z) dz.

(3.33)

pad

De intensiteit Iz=d (x, y) wordt dan gegeven door Iz=d (x, y) = |Uz=d (x, y)|2 ,

(3.34)

en Uz=d kan berekend worden uit Uz=0 met de Fresnel propagator (zie formule 3.27). De amplitude van de vlakke golf op de detector is Ud (x, y) = Ui (x, y)T (x, y) ∗ hF (x, y),

(3.35)

en voor de intensiteit op de detector vinden we Id (x, y) = Ii (x, y)|T (x, y) ∗ hF (x, y)|2 .

(3.36)

Aanpassing Fresnel propagator Bij de berekening van de intensiteit in vergelijkingen 3.35 en 3.36 werd verondersteld dat de golfamplitude constant is over het objectvlak (plane wave benadering). Dit is geldig wanneer de SOD veel groter is dan d (ODD) en als SOD >> D met D de grootte van het object. Bij een opstelling aan een synchrotron is dit meestal voldaan, maar in het geval van een opstelling aan een X-stralenbuis is dit in het algemeen niet het geval. Golosio et al. [44] brengen de korte SOD in rekening. De afleiding start met de veronderstelling van een

3.4. Golfoptica

28

puntbron die sferische golven uitzendt (spherical wave benadering). De amplitude in het objectvlak wordt dan eikr0 Ui (x0 , y0 ) = A . (3.37) r0 Als men dit in rekening brengt kan men aantonen dat de golf in het detectorvlak gegeven wordt door (op een constante fasefactor na) [44] Z Z 1 A ik 1 2 2 U (xd , yd ) = + (x0 − Xd ) + (y0 − Yd ) dx0 dy0 . T (x0 , y0 ) exp iλz1 d 2 z1 d (3.38) Hierin stelt z1 de SOD voor. Dit is een convolutie van T (x, y) met de propagatiekernel 1 A ik 1 2 2 + , (3.39) G(x, y) = exp x +y iλz1 d 2 z1 d geëvalueerd in het punt (Xd , Yd ) = (xd /M, yd /M ) met M de vergroting. Het grootste verschil met de Fresnel propagator 3.26 zit in de vergroting en de factor (1/z1 + 1/d) in de exponentiële. Deze factor wordt pas gelijk aan 1/d voor z1 >> d. De convolutie wordt in het Fourier domein berekend. De FT van de propagatorkernel kan berekend worden: 1 A 2 2 ˜ η) = exp −iπλ (ξ + η ) , (3.40) G(ξ, z1 dγ γ met γ = 1/z1 + 1/d. De genormeerde intensiteit wordt dan: I(x, y) =

|T (x, y) ∗ G(x, y)|2 . |Aeikr /r|2

(3.41)

In realiteit heeft men geen ideale puntbron. Hiervoor dient gecorrigeerd te worden door de convolutie te nemen van de intensiteit met een Gaussische functie g(σ, x, y). Hierin is σ de grootte van de geprojecteerde spot op de detector. We vinden dan finaal voor de intensiteit: I(xd , yd ) = I(xd , yd ) ∗ g(σ, xd , yd ). (3.42) Praktisch Golosio et al. [44] bekomen op basis van de sferische golf benadering resultaten die experimenten aan een XRT zeer goed benaderen. Hoewel de theorie vrij complex lijkt, is het formalisme vrij elegant te implementeren: De waarde van de golfamplitude na transmissie (uitdrukking 3.30) wordt berekend op een discreet grid. Als stapgrootte wordt bij voorkeur de pixelgrootte gekozen, die gerelateerd is aan de stapgrootte in het objectvlak via de vergroting M. Vervolgens berekenen we hiervan de FT De transferfunctie (uitdrukking 3.28 of 3.40) wordt berekend in het Fourierdomein

3.4. Golfoptica

29

Het product van de vorige twee wordt genomen en vervolgens neemt men de inverse FT

Vermits hier goede resultaten mee bekomen werden, is het zeker de moeite om op basis hiervan een simulatie te maken. In hoofdstuk 5 zullen we nagaan of de intensiteit gemeten op de detector overeenkomt met de berekeningen met formules 3.41 en/of 3.36.

Hoofdstuk 4

Experimentele methode 4.1

Opstelling UGCT

De metingen voor deze thesis werden uitgevoerd aan de 900-nm scanner van UGCT1 [45]. De X-stralenbron is een dual-head Feinfocus en bij een vermogen van 1 W kan een resolutie van 900 nm gehaald worden. De maximale X-stralen energie bedraagt 160 keV met een maximaal vermogen van 150 W. De samples worden opgesteld op een Micos rotatiemotor met luchtlager. Objecten met een diameter tot ongeveer 35 cm kunnen ingescand worden. Verschillende detectoren zijn ter beschikking. Hier werd gebruik gemaakt van een 40 × 40 cm Perkin-Elmer flatpanel detector. De opstelling is weergegeven in figuur 4.1 2 . [45]

Figuur 4.1: Experimentele opstelling van de 900 nm CT-scanner aan UGCT

Er werden ook metingen uitgevoerd aan de 400 nm CT-scanner (weergegeven in figuur 4.2). Deze opstelling werd enkel gebruikt wanneer zowel de SOD als SDD gevarieerd moesten worden. De maximale X-stralen energie bedraagt hier 160 keV met een maximaal vermogen 1 2

http://www.ugct.ugent.be/ De afgebeelde detector is de Photonic Science VHR-detector i.p.v. de Perkin-Elmer flatpanel detector.

30

4.2. Gebruikte software

31

van 35 W. Een resolutie van 400 nm kan worden gehaald bij een vermogen van 1 W. [45]

Figuur 4.2: Experimentele opstelling van de 400 nm CT-scanner aan UGCT

4.2

Gebruikte software

Voor het verwerken van de projecties en het reconstrueren van de slices werd gebruik gemaakt van Octopus. Dit is een reconstructiepakket ontworpen door UGCT. [17] De algoritmes en simulaties ontwikkeld tijdens dit thesiswerk werden ge¨ımplementeerd in LabVIEW . Verschillende codefragmenten (zoals het inlezen van tiff -bestanden) waren reeds voor handen, wat het programmeerwerk soms een stuk reduceerde. Bovendien geeft LabVIEW een kant en klare grafische gebruikersinterface, wat gebruiksvriendelijk is.

®

®

4.3

Sample

In dit thesiswerk wordt voor experimentele metingen voornamelijk gebruik gemaakt van een epoxyhars (Epofix van de firma Struers). Dit wordt veel door de geologen van UGCT gebruikt om zeer poreuze gesteenten samen te houden tijdens het verzagen. Het voordeel van dit object is dat het zeer homogeen3 blijkt te zijn, wat niet evident is op de schaal die we hier beschouwen. Het hars werd in een rietje gestold om zo tot een cilinder epoxyhars te komen. Een tweede voordeel van dit object is dat het zeer glad is. Hierdoor is de tweede afgeleide van de faseverschuiving aan de rand divergent (zie vergelijking 3.9) en wordt een duidelijk fasesignaal gegenereerd, wat uiteraard gewenst is in dit werk. 3

Een homogeen object is nodig om een fantoom zo goed mogelijk te benaderen.

4.4. Secondary spot

4.4

32

Secondary spot

Bij de metingen werd aanvankelijk veel hinder ondervonden door het secondary spot effect, dat zeer goed zichtbaar was op een zeer homogeen sample. Dit houdt in dat in de XRT ook straling veroorzaakt wordt door elektronen die — eventueel na verstrooiing aan het trefplaatje — de wand van de buis raken. Dit veroorzaakt een tweede spot die veel groter is dan de eigenlijke spot op het trefplaatje. De invloed van de spotgrootte op de projectie is ge¨ıllustreerd in figuur 4.3: de spot wordt ook vergroot en veroorzaakt een waziger beeld. De tweede, grotere spot creëert een bijkomende projectie die veel waziger is en over een groter oppervlak uitgesmeerd is. Dit is ge¨ıllustreerd in figuur 4.4. Hiervoor werd gecorrigeerd door een collimator (vervaardigd uit een dun loodplaatje) te plaatsen net na de buis. De invloed van de secondary spot op een reconstructie wordt weergegeven in figuur 4.5. Dit is de lijndoorsnede van een gereconstrueerde slice van het homogeen staafje epoxyhars. De attenuatiecoëfficiënt van dit homogeen sample is overal hetzelfde. Een lijndoorsnede van een gereconstrueerde slice moet dus een horizontale lijn zijn. Er kan eventueel wel een lichte daling zijn naar het centrum toe door het cupping artefact (zie sectie 2.2.3). In figuur 4.5 is in het centrum een stijging in de gereconstrueerde attenuatiecoëfficiënt te zien. Dit is het gevolg van de storende tweede projectie die gecreëerd wordt. Deze zorgt ervoor dat na normering een lagere intensiteit teruggevonden wordt. Hierdoor zal de gereconstrueerde attenuatiecoëfficiënt groter zijn. Door het plaatsen van een collimator wordt het secondary spot effect zeer goed gereduceerd. Bij zeer hoge energieën echter blijkt het dunne loodplaatje toch nog secundaire straling door te laten. Dit kan verholpen worden door eventueel in de toekomst een dikkere collimator te gebruiken.

Figuur 4.3: Invloed van de spotgrootte op een projectie [16]

4.4. Secondary spot

33

Figuur 4.4: De secondary spot veroorzaakt een tweede, wazige projectie van het object (een epoxyhars staafje). De grijswaarden werden aangepast zodat het effect duidelijk zichtbaar is

1

doorsnede

0.95

0.9

0.85

MBA 80 kV Paganin 80 kV 0.8 0

100

200

300

400

500

600

700

800

900

pixel

Figuur 4.5: Lijnprofiel van een gereconstrueerde slice. De stijging in het centrum is het gevolg van de tweede, wazige projectie

Hoofdstuk 5

Simulatie van fasecontrast 5.1

Inleiding

Zoals reeds vermeld kunnen we niet kiezen of er al dan niet refractie optreedt. Het fasesignaal kan gebruikt worden om een reconstructie van het refractive index decrement δ te maken, maar als we een reconstructie van de attenuatiecoëfficiënt µ willen maken, is het optreden van refractie hinderlijk. Bij een gewone reconstructie ontstaan er dan artefacten. Naast de analytische reconstructietechniek zoals besproken in sectie 2.2.2, bestaan er ook iteratieve reconstructiealgoritmes (besproken in sectie 2.2.4). Essentieel wordt de projectie van een virtuele representatie van het object steeds opnieuw gesimuleerd en aangepast tot de gesimuleerde projectie overeenkomt met de werkelijke projectie. Hoe beter de simulatie, hoe beter de reconstructie. Als men dus de fase effecten ook in rekening kan brengen in de simulatie, kunnen artefacten t.g.v. de fase gereduceerd worden. De simulatie zou dan kunnen ge¨ımplementeerd worden in het iteratief reconstructie algoritme SART [46]. Simulaties op basis van de TIE werden binnen UGCT al uitgevoerd maar overschatten het signaal [16, 47]. Nochtans hebben Peterzol et al. [31] aangetoond dat de TIEW en de Fresnel benadering vrij goede overeenkomst vertonen. Deze experimenten werden echter uitgevoerd aan een synchrotron. El-Ghazaly et al. [19] vonden een goede overeenkomst tussen simulaties op basis van Fresnel diffractie en experiment. Ook hier werd de opstelling gekenmerkt door een grote SOD (X-stralen geproduceerd aan een microtron). Dit is te begrijpen vermits bij grote SOD de coherentielengte groot is (zie vergelijking 3.21). In sectie 3.4 werd uitgelegd hoe het fasesignaal kan berekend worden op basis van de golfoptica. In sectie 3.4.2 werd een aanpassing voorgesteld op de klassieke Fresnel propagator om de korte SOD in rekening te brengen. Vermits Golosio et al. [44] hier goede resultaten mee bekomen en het een elegant formalisme is om te implementeren — dus ook eventueel mogelijkheden biedt in een iteratief reconstructiealgorime — gingen we deze methode na. 34

5.2. Simulatie homogene cilinder

5.2

Simulatie homogene cilinder

5.2.1

Vlakke golf versus sferische golf

35

Eerst werd het lijnprofiel1 van een homogene cilinder berekend met vergelijkingen 3.36 en 3.41. De parameters zijn weergegeven in tabel 5.1. De waarden van de scanneropstelling zijn deze van de opstelling aan UGCT. De waarden voor de energie (polychromatische bron), spotgrootte en δ van het object zijn minder evident. Voor de energie is het gewogen gemiddelde de beste keuze. De spotgrootte is van de orde (enkele) micrometer. Voor een zuiver materiaal (bv. koolstof) kan δ berekend worden uit de atomaire verstrooiingsfactor. Voor samenstellingen is dit niet mogelijk. Het vinden van een geschikt zuiver object is niet evident. De parameterwaarden in de simulatie zijn zo gekozen dat het berekend fasesignaal overeenkomt met het experimenteel gemeten fasesignaal van het epoxystaafje. We merken verder op dat dit realistische parameterwaarden zijn. In figuur 5.1 wordt het resultaat getoond aan de linker rand van de cilinder, voor de vlakke golf2 en sferische golf benadering (beide benaderingen werden besproken in sectie 3.4.2). Er is duidelijk een verschil in hoogte en frequentie van het patroon. In figuur 5.2 wordt een volledig lijnprofiel getoond en is de eindige spotgrootte in rekening gebracht: het diffractiepatroon aan de rand verdwijnt en er blijft enkel een piek en een dal over3 . We vinden dus het resultaat terug zoals voorspeld in de stralenoptica. We merken meteen dat er een groot verschil is tussen de twee benaderingen. Later zullen we zien dat de sferische benadering de experimentele waarden best benadert. De berekeningen werden uitgevoerd met een veel kleinere stapgrootte dan de pixelgrootte4 . Daarom werd het resultaat uitgemiddeld zodat we een waarde krijgen per pixel. Het resultaat is getoond in figuur 5.3. Het uitmiddelen heeft een invloed op de piekhoogte. Deze invloed wordt meer significant als de piek zeer smal wordt.

5.2.2

Vergelijking met het experiment

Zoals reeds vermeld is het moeilijk een sluitend bewijs te geven dat de berekeningen het fasesignaal helemaal correct teruggeven. De vele parameters (polychromatische bundel, spotgrootte, δ, invloed van de detector) maken dit zeer moeilijk. Het kan echter nuttig zijn 1

De simulatie werd beperkt tot 1 dimensie: de intensiteit op één rij pixels van de detector werd berekend. Met een lijnprofiel bedoelen we dan ook de intensiteit op een rij pixels van de detector. 2 De vergroting werd in vergelijkingen 3.26 tot 3.36 niet in rekening gebracht. Hier brengen we de vergroting wel in rekening vermits we met een conische bundel geometrie te maken hebben. 3 Het dal in figuur 5.2b is niet zo goed zichtbaar door de attenuatie. 4 Uit het Nyquist sampling theorema volgt dat, om een analoog signaal accuraat te beschrijven op basis van discrete bemonstering, dit signaal bemonsterd moet worden met een frequentie die minstens tweemaal groter is dan de hoogste frequentie die in dat signaal aanwezig is. Gezien de hoge frequenties in het signaal in figuur 5.1b is een stapgrootte kleiner dan de pixelgrootte nodig om het signaal juist te evalueren.


36

Tabel 5.1: Parameters voor de simulatie van de projectie van een homogene cilinder

breedte detector (pixels) pixelgrootte (mm) straal cilinder (mm) SOD (mm) SDD (mm) µ (mm−1 ) δ λ (nm) spotgrootte (mm)

2048 0.2 2.2 25 1405 0.055 1.70 · 10−7 0.035429 0.002

eens te kijken of we met realistische parameters de piek kunnen bekomen. De berekeningen getoond in vorige sectie waren uitgevoerd met parameters die goede overeenkomsten vertonen met het epoxyhars staafje. De experimentele en berekende lijnprofielen van zo een hars staafje zijn weergegeven in figuur 5.4. De parameters zijn uiteraard zo gekozen dat de lijnprofielen samenvallen5 . Maar het toont dat, met realistische parameters, een goede overeenkomst kan gevonden worden. Met de klassieke Fresnel propagator is het resultaat een grootteorde verschillend. Een detail van het lijnprofiel (linker rand) is weergegeven in figuur 5.5. Hierin is ook een uitgemiddeld en niet-uitgemiddeld lijnprofiel getekend. Het is duidelijk dat het uitmiddelen een invloed heeft op het lijnprofiel.

5.2.3

Polychromatische bundel

In realiteit heeft men te maken met een polychromatische X-stralenbundel. Zoals reeds intu¨ıtief uitgelegd werd in de stralenoptica (zie figuur 3.5) heeft dit een invloed op de breedte van de piek. Dit zou misschien kunnen verklaren waarom de echte piek in figuur 5.5 iets breder is dan de berekende waarde. We gingen na of we dit in onze berekeningen — gebaseerd op de golfoptica — ook waarnemen. Hiervoor werd een spectrum gesimuleerd zoals het waargenomen wordt door de Varian detector bij een XRT buisspanning van 60 kV. Het spectrum is weergegeven in figuur 5.6. Dit spectrum werd opgedeeld per 5 keV en er werd een simulatie gemaakt voor de centrale energiewaarde van elke bin. Elke energiebin krijgt dan een contributie in het lijnprofiel volgens het getoonde spectrum. Het refractive index decrement δ is afhankelijk van de energie, en wordt gegeven door formule 3.5. Hierin kan f1 als constant beschouwd worden6 en dan is δ evenredig met 1/E 2 (zie formule 3.5). De waarden voor δ en de gewichtsfactoren bij de verschillende energieën voor deze simulatie zijn weergegeven in tabel 5.2. De parameters van de opstelling zijn weergegeven in tabel 5.3. 5 6

Het hars staafje wijkt wel lichtjes af van een perfecte cilinder, zoals te zien is op het lijnprofiel. In het energiegebied 1 keV tot 100 keV is dit een juiste veronderstelling, zie hiervoor figuur 3.2a.


37

a) vlakke golf

intensiteit

2.5 2 1.5 1 0.5 -695

-690

-685

-680

-675

-670

2.2

b)

sferische golf

2 1.8 intensiteit

1.6 1.4 1.2 1 0.8 0.6 0.4 -695

-690

-685

-680

-675

-670

pixel

Figuur 5.1: Simulatie van het fasesignaal van een homogene cilinder. De linker rand van de projectie is getoond. Er is nog niet gecorrigeerd voor de eindige spotgrootte. (a) Vlakke golf benadering, (b) sferische golf benadering. De rand bevindt zich rond pixel 681

De berekeningen werden uitgevoerd op een grid met een stapgrootte die veel kleiner is dan de pixelgrootte, daarom moeten we de resultaten finaal uitmiddelen om tot een waarde per pixel te komen. De niet-uitgemiddelde resultaten zijn weergegeven in figuur 5.7, de uitgemiddelde resultaten in figuur 5.8. In figuur 5.7a zijn enkele lijnprofielen getekend die zullen bijdragen tot het totale lijnprofiel dat geschetst is in figuur 5.7b. De verschillende bijdragen zoals getoond in figuur 5.7a zijn nog niet geschaald, vermits sommige bijdragen dan zeer klein worden en moeilijk te plotten zijn. Daarom werd de ongeschaalde versie geplot. In figuur 5.7b is het profiel van de polychromatische bundel geplot. Daarna hebben we ook het lijnprofiel geplot zoals we bij onze vergelijking met de experimenten steeds doen: voor de energie nemen we het gewogen gemiddelde van het spectrum (32 keV) en we passen δ aan zodat de tophoogte overeenkomt (δ = 1.41 · 10−7 ). Het resultaat met de polychromatische bundel blijkt iets breder te zijn, maar het verschil is zeer klein (zeker op het gebind profiel weergegeven in figuur 5.8b). Dit verklaart dus niet waarom de echte piek iets breder is. In sectie 3.2.2 werd uitgelegd dat bij een polychromatische bundel, de piek minder hoog maar breder is dan bij een monochromatische bundel. De oppervlakte onder de piek moet wel dezelfde zijn. Dit werd geverifieerd door een monochromatische bundel te simuleren


38

a) 1.6 vlakke golf

1.5

intensiteit

1.4 1.3 1.2 1.1 1 0.9 0.8 0.7 -800

-600

-400

-200

0

200

400

600

800

200

400

600

800

1.05

b)

sferische golf 1

intensiteit

0.95 0.9 0.85 0.8 0.75 0.7 -800

-600

-400

-200

0 pixel

Figuur 5.2: Simulatie van een lijnprofiel van een homogene cilinder. Er is gecorrigeerd voor de eindige spotgrootte door convolutie met een Gauss. (a) Vlakke golf benadering, (b) sferische golf benadering

a) 1.6 vlakke golf

1.5

intensiteit

1.4 1.3 1.2 1.1 1 0.9 0.8 0.7 -800 b)

-600

-400

-200

0

200

400

600

800

200

400

600

800

1.05 sferische golf 1

intensiteit

0.95 0.9 0.85 0.8 0.75 0.7 -800

-600

-400

-200

0 pixel

Figuur 5.3: Simulatie van een lijnprofiel van een homogene cilinder na uitmiddelen. (a) Vlakke golf benadering, (b) sferische golf benadering


39

1.05

1

intensiteit

0.95

0.9

0.85

0.8

0.75 gebind exp data hoge resolutie 0.7 -800

-600

-400

-200

0 pixel

200

400

600

800

Figuur 5.4: Lijnprofiel van een epoxyhars staafje. De experimentele en berekende waarden (uitgemiddeld en in hoge resolutie) zijn weergegeven

gebind exp data hoge resolutie

1.04

intensiteit

1.02

1

0.98

0.96

-700

-695

-690

-685

-680 pixel

-675

-670

-665

-660

Figuur 5.5: Detail van de linker rand van een epoxyhars staafje. De experimentele en berekende waarden (uitgemiddeld en in hoge resolutie) zijn weergegeven


40

die een zelfde oppervlak onder de piek heeft als de polychromatische bundel. Hiervoor namen we E = 22 keV en δ = 1.55 · 10−7 7 . Dit is ook weergegeven in figuur 5.7b. We zien kwalitatief het juiste gedrag, maar het blijkt ook weer een klein effect te zijn. Tabel 5.2: Gebruikte energieën en δ’s in de polychromatische bundel

E(keV)

λ (nm)

relatief gewicht

δ

2.55E+00 7.55E+00 1.25E+01 1.75E+01 2.25E+01 2.75E+01 3.25E+01 3.75E+01 4.25E+01 4.75E+01 5.25E+01 5.75E+01

0.48627451 0.164238411 0.099598394 0.071060172 0.055233853 0.045173042 0.038212635 0.033110814 0.029210836 0.026132771 0.023641563 0.021583986

1.28E-09 4.15E-02 7.73E-02 1.03E-01 1.46E-01 1.45E-01 1.26E-01 1.19E-01 1.07E-01 8.74E-02 6.16E-02 2.81E-02

1.16E-05 1.32E-06 4.86E-07 2.47E-07 1.50E-07 1.00E-07 7.16E-08 5.37E-08 4.18E-08 3.35E-08 2.74E-08 2.28E-08

Tabel 5.3: Parameters voor de simulatie van een homogene cilinder bij een polychromatische bundel

breedte detector (pixels) pixelgrootte (mm) straal cilinder (mm) SOD (mm) SDD (mm) µ (mm−1 ) spotgrootte (mm)

5.2.4

2048 0.2 2.2 25 1405 0 0.002

Vari¨ erende SOD

We hebben gezien dat het mogelijk is een realistisch fasesignaal te fitten aan experimentele gegevens. Vervolgens werd de SOD aangepast en de invloed op de tophoogte nagegaan. De scan werd bij twee verschillende buisspanningen uitgevoerd (50 kV en 80 kV). De maximale tophoogte werd bepaald met een speciaal ontworpen LabVIEW programma. De gemiddelde tophoogte voor 475 detectorrijen van een projectie werd bepaald, evenals de standaardfout hierop. Vervolgens werd het fasesignaal ook berekend. De parameters voor

®

7

We nemen deze waarden om niet te veel af te wijken van de waarden gekozen in tabel 5.2


41

Varian detector 5e-005

detector respons

4e-005

3e-005

2e-005

1e-005

0 0

10

20

30

40

50

60

Energie (keV)

Figuur 5.6: Simulatie van een X-stralen spectrum zoals waargenomen door de Varian-detector bij een buisspanning van 60 kV

a) 1.12 12.5 keV 17.5 keV 27.5 keV 37.5 keV 47.5 keV

1.1

intensiteit

1.08 1.06 1.04 1.02 1 0.98 0.96

b) 1.05

gewogen gemiddelde mono-energetisch mono-energetisch zelfde piekoppervlak

1.04

intensiteit

1.03 1.02 1.01 1 0.99 0.98 -684

-683.5

-683

-682.5

-682

-681.5

-681 pixel

-680.5

-680

-679.5

-679

-678.5

-678

Figuur 5.7: Detail van de linker rand van het lijnprofiel van een cilinder bij een polychromatische bundel. (a) De bijdragen van enkele energieën (niet gewogen), (b) de totale piek veroorzaakt door de verschillende energieën


42

a) 1.1 12.5 keV 17.5 keV 27.5 keV 37.5 keV 47.5 keV

1.08

intensiteit

1.06 1.04 1.02 1 0.98 0.96

b) 1.05

gewogen gemiddelde mono-energetisch mono-energetisch zelfde piekoppervlak

1.04

intensiteit

1.03 1.02 1.01 1 0.99 0.98 -685

-684.5

-684

-683.5

-683

-682.5

-682 pixel

-681.5

-681

-680.5

-680

-679.5

-679

Figuur 5.8: Detail van de linker rand van het lijnprofiel van een cilinder bij een polychromatische bundel na uitmiddelen. (a) De bijdragen van enkele energieën (niet gewogen), (b) de totale piek veroorzaakt door de verschillende energiebijdragen

de simulaties worden weergegeven in tabel 5.4. De parameters voor de detector en SDD zijn uiteraard deze van de experimentele opstelling. De waarden voor δ en de grootte van de spot werden gekozen opdat het fasesignaal bij een SOD van 25 mm de juiste tophoogte zou geven. Voor de fotonenergie bij de verschillende buisspanningen werd een gewogen gemiddelde genomen. De resultaten zijn weergegeven in figuur 5.9. De fout op het gemeten signaal is de standaardfout over 475 rijen. Bij de berekening van het fasesignaal werd met een stapgrootte gerekend die 100 keer kleiner is dan de pixelgrootte. Dit betekent dat finaal het resultaat moet uitgemiddeld worden8 , zodat er per pixel één waarde voor de intensiteit is. De manier waarop men de bins kiest ter hoogte van de piek blijkt een grote invloed te hebben op de tophoogte. Daarom werd het gemiddelde genomen van 100 verschillende bins9 . We moeten opmerken dat de daling van de piek in de berekeningen het gevolg is van het uitmiddelen. De niet-uitgemiddelde resultaten blijven stijgen, maar de piek wordt zeer dun, zodat piek en dal op dezelfde pixel komen te liggen. Hierdoor wordt globaal toch een daling waargenomen. 8

Een bin waarover uitgemiddeld wordt om één waarde per pixel te bekomen bevat dan 100 waarden De positie waarop we starten voor de eerste bin te maken nemen we achtereenvolgens van 1 tot 100 zodat een piekhoogte berekend wordt voor elke mogelijke manier van het vormen van bins. 9

5.3. Bespreking

43

Tabel 5.4: Parameters voor de simulatie van een homogene cilinder bij verschillende SOD

breedte detector (pixels) pixelgrootte (mm) straal cilinder (mm) SDD (mm) µ (mm−1 ) δ λ (nm) spotgrootte (mm)

50 kV

80 kV

2048 0.2 2.2 1405 0.055 1.70 · 10−7 0.035 0.002

2048 0.2 2.2 1405 0.042 1.33 · 10−7 0.031 0.0025

a) 1.05 exp data simulatie

1.045

tophoogte

1.04 1.035 1.03 1.025 1.02 1.015 1.01 1.005 50 b)

100

150

200

250

300

1.04 exp data simulatie

1.035

tophoogte

1.03 1.025 1.02 1.015 1.01 1.005 1 50

100

150

200

250

300

SOD

Figuur 5.9: De gemeten en berekende tophoogtes bij verschillende SOD. (a) Buisspanning van 50 kV, (b) 80 kV

5.3

Bespreking

In figuur 5.9 vallen meteen de grote fouten op. Bij de berekende waarden is de piekhoogte sterk afhankelijk van de manier waarop men de bins kiest. Er is wel een zelfde trend

5.3. Bespreking

44

zichtbaar. De grote spreiding op de uitgemiddelde waarden is te verklaren doordat de piek zeer smal wordt en de bin waarover uitgemiddeld wordt te breed is10 . Een betere maat dan de tophoogte zou hier bv. de oppervlakte onder de piek zijn, maar deze was experimenteel moeilijk te bepalen. Er werden ook metingen van het fasesignaal uitgevoerd aan de 400 nm CT-scanner. Hier kan ook de SDD aangepast worden. Het fasesignaal was hier wel significant kleiner en er werd veel hinder ondervonden van de secondary spot (hier is geen collimator aanwezig). Trends waren veel minder aanwezig en de overeenkomst was veel minder goed. Zoals we in sectie 5.1 uitgelegd hebben, is het de bedoeling om in een iteratief reconstructiealgoritme rekening te houden met dit fasesignaal. Eén van de voordelen van iteratieve reconstructie t.o.v. analytische reconstructie is dat men meer fysica (i.e. verstrooiing van de bundel, . . . ) kan introduceren in de reconstructie. Als men de faseverschuiving hierin wil implementeren, is een zo goed mogelijke simulatie van het fenomeen nodig. De parameter δ zou dan ook een plaatsafhankelijke variable zijn die aangepast wordt zodat de simulatie het beeld beter benadert. Vandaar dat hier ook kwalitatief een overeenkomst gezocht werd. Uit onze resultaten is gebleken dat we kwalitatief het juiste resultaat terugvinden (i.e. de piek en het dal in de intensiteit aan de rand van een cilinder met realistische parameters). Bij het variëren van de SOD blijkt een zelfde trend zichtbaar te zijn tussen experiment en simulatie, maar de overeenkomst blijkt moeilijker te bereiken aan de 400 nm scanner. Een bijkomend probleem is dat een zeer kleine stapgrootte nodig is om een goede simulatie te bekomen, wat uiteraard de implementatie niet meer praktisch maakt. Hoewel het berekende fasesignaal helemaal niet verkeerd lijkt en de wiskundige constructie met een convolutie aantrekkelijk is, blijkt deze methode niet zo bruikbaar. De vele (onbekende) parameters en de grote afhankelijkheid van het uitmiddelen maken dat er grote fouten zijn en de trends soms niet zo goed overeenkomen.

10

Dit volgt uit het Nyquist sampling theorema: als we een te grote stapgrootte nemen in de evaluatie van de functie, bemonsteren we het signaal met een te lage frequentie en dan zal de vorm (en dus ook de tophoogte) veranderen.

Hoofdstuk 6

Fasereconstructie In de vorige hoofdstukken hebben we de theorie rond fasecontrast gezien. De berekeningen in het vorig hoofdstuk illustreerden kwalitatief het edge-enhancement effect dat door de faseshift veroorzaakt wordt. In dit hoofdstuk bespreken we twee methodes die dit effect gebruiken om het refractive index decrement δ te reconstrueren. Deze zullen zeer gelijkaardig blijken in vorm. De methodes worden vergeleken in het volgende hoofdstuk.

6.1 6.1.1

Bronnikov algoritme Theorie

Het Bronnikov algoritme1 geeft een elegante reconstructiemethode voor een puur fase object die zeer gelijkaardig is aan het FBP-algoritme. Beschouw het coördinatenstelsel zoals weergegeven in figuur 6.1, waarbij de coördinaten in het object gegeven worden door (x1 , x2 , x3 ) en de co¨ ordinaten in het beeldvlak door (x, y). Het object bevindt zich op z = 0 en de detector op z = d. Een monochromatische vlakke golf zal na transmissie door het object beschreven worden door Uθ (x, y) = Tθ (x, y)Ui ,

(6.1)

waarbij Ui en Uθ respectievelijk de invallende en uitgaande vlakke golf voorstellen. De index θ slaat op de rotatiehoek van het object. De functie Tθ is de transmissiefunctie en wordt gegeven door 1 (6.2) Tθ = e− 2 µθ (x,y)+iφθ (x,y) . De eerste term in de exponentiële stelt de attenuatie van de golf voor en de tweede term de faseshift. De argumenten van de exponentiëlen, µθ en φθ , zijn de lijnintegralen langs het pad van de straal van respectievelijk µ en f 2 (met f = nr − 1). Ze worden gegeven 1 2

De volledige afleiding en wiskundige uitwerkingen zijn terug te vinden in [32] en [33]. We noteren het refractive index decrement δ als f om verwarring met de Dirac delta functie te vermijden.

45

6.1. Bronnikov algoritme

46

Figuur 6.1: Co¨ ordinatenstelsel gebruikt bij de afleiding van het Bronnikov algoritme [32]

door

Z µθ (x, y) = R2

en

2π φθ (x, y) = λ

µ(x1 , x2 , y) · δ(x − x1 cos θ − x2 sin θ) dx1 dx2

Z R2

f (x1 , x2 , y) · δ(x − x1 cos θ − x2 sin θ) dx1 dx2 .

(6.3)

(6.4)

De intensiteit op een afstand d van het object wordt, bij een weinig variërende µθ , gegeven door [32] λd 2 Iθ,z=d (x, y) = Iθ,z=0 1 − ∇ φθ (x, y) , (6.5) 2π wat de TIEW is. Deze legt een verband tussen de gemeten intensiteit op de detector en de laplaciaan van de fasefunctie. Het algoritme gebruikt vergelijkingen 6.4 en 6.5 om de fasedistributie f (x1 , x2 , x3 ) terug te vinden. Het maakt gebruik van het volgende theorema. Theorema 2. Stel gθ (x, y) =

Iθ,z=d (x, y) − 1, Iθ,z=0

(6.6)

met Iθ,z=d (x, y) en Iθ,z=0 respectievelijk de intensiteit op de detector en de intensiteit in het objectvlak na doorgang door het materiaal; dan geldt ∂2 ˆ 1 f (s, θ, ω) = − gˆθ (s, ω). 2 ∂s d Hierin worden gˆ(s, ω) en fˆ(s, θ, ω) gegeven door3 Z gˆ(s, ω) = g(x, y)δ(s − x sin ω − y cos ω) dx dy

(6.7)

(6.8)

R2

en fˆ(s, θ, ω) =

Z R3

3

f (x1 , x2 , x3 )δ(s − (x1 cos θ + x2 sin θ) · sin ω − x3 cos ω) dx1 dx2 dx3 . (6.9)

gˆ(s, ω) en fˆ(s, θ, ω) worden de Radon getransformeerden van g(x, y) en f (x1 , x2 , x3 ) genoemd.

6.1. Bronnikov algoritme

47

De functie gˆ is de integratie van g over een rechte s = x sin ω +y cos ω terwijl fˆ de integratie is van f over het vlak met normaal n =(cos θ sin ω, sin θ sin ω, cos ω). Het bewijs is terug te vinden in het artikel van Bronnikov [32]. De fasedistributie f (x1 , x2 , x3 ) wordt bekomen door de inverse Radon getransformeerde te nemen van fˆ(s, θ, ω) en wordt gegeven door Z π Z π ∂2 1 dθ 2 fˆ(s, θ, ω). f (x1 , x2 , x3 ) = − 2 sin(ω) dω (6.10) 4π 0 ∂s 0 Uit het theorema volgt dat de tweede afgeleide van fˆ vervangen kan worden door −ˆ gθ (s, ω)/d. Men kan aantonen dat dit verder vereenvoudigd kan worden tot [32] Z Z Z π 1 |x3 − y| f (x1 , x2 , x3 ) = 2 dθ dx dygθ (x, y) . 4π d 0 (x1 cos θ + x2 sin θ − x)2 + (x3 − y)2 (6.11) Met |y| (6.12) q(x, y) = 2 x + y2 kan dit geschreven worden als Z π 1 f (x1 , x2 , x3 ) = 2 q ∗ gθ dθ. (6.13) 4π d 0 Hierin is q ∗ gθ de convolutie van gθ (x, y) met q(x, y). Deze kan best in het Fourier domein berekend worden. De FT van q(x, y) wordt gegeven door q˜(ξ, η) =

ξ2

|ξ| . + η2

(6.14)

De convolutie q ∗ gθ wordt dan F −1 [F{q} · F{g}] met F de FT.

6.1.2

Modified Bronnikov algoritme

Het Bronnikov algoritme veronderstelt een puur faseobject, een object dat dus geen attenuatie vertoont. Voor samples die wel attenueren, zijn eigenlijk twee metingen nodig: een meting op een afstand d van het object en een meting op d = 0. Op een afstand d = 0 is er enkel attenuatie en nog geen fasesignaal. Vergelijking 6.5 toont dat we het zuiver fasesignaal kunnen overhouden door de projecties op een afstand d te delen door de projecties op d = 04 . Dit wordt ook gedaan in uitdrukking 6.6. Groso et al. [48] stelden een oplossing voor opdat slechts één reeks metingen nodig zou zijn op één bepaalde afstand voor objecten die toch een lichte attenuatie vertonen. Er wordt een correctieparameter α toegevoegd aan de filterfunctie5 6.14 [41, 48]: q˜(ξ, η) = 4

ξ2

|ξ| . + η2 + α

(6.15)

De vergroting moet dan wel voor beide projecties dezelfde zijn, zodat we — op een fasesignaal na — twee identieke beelden hebben. 5 De parameter α is in de implementatie van het algoritme afhankelijk van de dimensie van het beeld.

6.2. Methode van Paganin

48

Het algoritme met deze filterfunctie duiden we aan met het Modified Bronnikov Algoritme (MBA). De invloed van deze parameter kan als volgt begrepen worden. De functies 6.14 en 6.15 zijn low-pass filters: deze onderdrukken de hoge frequentiecomponenten en laten de lage frequentiecomponenten door. De parameter α in 6.15 zorgt er voor dat de lage frequenties (zijnde de traag variërende attenuatie bijdrage) meer onderdrukt worden dan in 6.14.

6.1.3

Implementatie

®

Het algoritme werd ge¨ımplementeerd in LabVIEW . De plaats van het algoritme in de reconstructie is als volgt. De reconstructieformule 6.13 is zeer gelijkaardig aan het FBP algoritme (zie vergelijking 2.6). In de filter q˜(ξ, η) herkennen we de ramp filter |ξ|. We splitsen de filter in twee delen. Eerst wordt op de functie gθ de filter 1/(ξ 2 + η 2 + α) toegepast. Deze inverteert de laplaciaan om de faseverschuiving te krijgen [16]. De waarde die men hieruit bekomt is zeer klein, daarom wordt die vermenigvuldigd met een bepaalde factor A. Daarvan wordt de exponentiële genomen zodat we eAq∗gθ

(6.16)

bekomen. Dat wordt dan ingevoerd in het klassieke FBP algoritme, waar de negatieve natuurlijke logaritme genomen wordt en waar de ramp filter wordt toegepast (zie vergelijkingen 2.5 en 2.6). Door de extra factor in de exponentiële be¨ınvloeden we de absolute waarde van de gereconstrueerde grootheid van ons sample. Daarom zullen we ons in de bespreking van de reconstructies concentreren op de relatieve waarden van de verschillende componenten.

6.2 6.2.1

Methode van Paganin Theorie

De methode van Paganin (MP) [34] vertrekt van de veronderstelling dat het sample homogeen 6 is. Op die manier kan men simultaan µ en δ van een sample bepalen. Het algoritme wordt eerst afgeleid voor een parallelle bundel. Het algoritme gaat uit van de TIE: 2π ∂ I(r⊥ , z) (6.17) ∇⊥ · [I(r⊥ , z)∇⊥ φ(r⊥ , z)] = − λ ∂z die het verband legt tussen de variatie in de fase φ in het vlak loodrecht op de voortplantingsrichting en de verandering in de intensiteit I na propagatie. Stel T (r⊥ ) de geprojecteerde dikte van het sample op het objectvlak. De intensiteit van de bundel na doorgang door het materiaal wordt gegeven door I(r⊥ , z = 0) = I0 e−µT (r⊥ ) . 6

Met homogeen bedoelen we gemaakt uit één materiaal en dus met een uniforme µ en δ.

(6.18)


49

Bij dunne objecten is de faseverschuiving na transmissie evenredig met de geprojecteerde dikte: 2π φ(r⊥ , z = 0) = − δT (r⊥ ). (6.19) λ Als we deze uitdrukkingen substitueren in de TIE bekomen we −

δ ∂ I0 ∇2⊥ e−µT (r⊥ ) = I(r⊥ , z = 0). µ ∂z

(6.20)

Het rechter lid van deze vergelijking kunnen we benaderen door I(r⊥ , z = d) − I0 e−µT (r⊥ ) ∂ I(r⊥ , z = 0) ≈ ∂z d met d de afstand tussen het sample en de detector. Vergelijking 6.20 wordt nu dδ 2 − ∇⊥ + 1 I0 e−µT (r⊥ ) = I(r⊥ , z = d). µ

(6.21)

(6.22)

De intensiteiten in het objectvlak en het beeldvlak worden als Fourierintegralen geschreven: Z Z I0 −µT (r⊥ ) F{e−µT (r⊥ ) ]}eik⊥ · r⊥ dk⊥ (6.23) I0 e = 2π Z Z 1 (6.24) I(r⊥ , z = d) = F{I(r⊥ , z = d)}eik⊥ · r⊥ dk⊥ . 2π Als we dit substitueren in vergelijking 6.22, rekening houden met de vergroting [49] en oplossen naar T (r⊥ ) bekomen we: 1 F{M 2 I(M r⊥ , z = d)}/I0 −1 T (r⊥ ) = − ln F µ . (6.25) µ dδ|k⊥ |2 /M + µ Uit de geprojecteerde dikte kan met vergelijkingen 6.18 en 6.19 de intensiteit en de faseverschuiving teruggevonden worden.

6.2.2

Implementatie

®

Het algoritme werd in LabVIEW ge¨ımplementeerd. Het algoritme levert de geprojecteerde dikte. Relatief gezien is het reconstrueren van µ equivalent met het reconstrueren van δ (want zowel de attenuatie als de faseverschuiving zijn na transmissie door het object evenredig met de geprojecteerde dikte T , dus de functie Pθ (t) uit vergelijking 2.5 is op een factor na dezelfde voor beide grootheden). Voor een efficiënte implementatie wordt de geprojecteerde dikte met −µ vermenigvuldigd en neemt men daarvan de exponentiële. Het komt er dus op neer dat we gewoon het argument van de logaritme moeten uitrekenen. Dit dient als invoer voor het klassiek FBP-algoritme waarmee de verdere reconstructie uitgevoerd wordt. Het is dus duidelijk dat we µ reconstrueren, maar die is evenredig met δ. De verhouding µ/δ is een input parameter die in de reconstructie manueel aangepast wordt.


6.2.3

50

Vergelijking tussen MBA en MP

Beide algoritmes worden op een totaal verschillende manier afgeleid, maar de vorm van beide is zeer gelijkaardig. Ze kunnen als volgt samengevat worden. 1. Men neemt de FT van een projectie. Bij MBA moeten we theoretisch wel de FT nemen van de functie gθ (waar de attenuatie weggedeeld is, zie formule 6.6), maar om de attenuatie te onderdrukken wordt de parameter α ingevoerd. 2. Er wordt een filter toegepast van de vorm f (ξ, η) =

a . ξ2 + η2 + b

(6.26)

Bij MP wordt daarna nog de natuurlijke logaritme genomen. 3. Hetgeen we bekomen uit de algoritmes dient als invoer voor het FBP-algoritme waarmee de reconstructie verder uitgevoerd wordt. In het volgende hoofdstuk vergelijken we beide methodes. We onderzoeken of de voorwaarden die de methodes opleggen (een niet-attenuerend object bij MBA en een homogeen sample bij MP) effectief noodzakelijk zijn en of de ene methode dan ook beter presteert dan de andere.

Hoofdstuk 7

Vergelijkende studie tussen MBA en MP Zoals reeds vermeld lijken beide algoritmes sterk op elkaar hoewel ze van verschillende veronderstellingen vertrekken. De methode van Paganin veronderstelt een homogeen sample. MBA daarentegen veronderstelt een niet-attenuerend object. Voor de weinige resterende attenuatie wordt de correctieparameter α ge¨ıntroduceerd. We gaan na of deze voorwaarden ook noodzakelijk zijn voor een goede werking van de desbetreffende methode. Hiertoe maken we gebruik van een echt object en gesimuleerde objecten. Het gebruikte fantoom wordt beschreven in appendix A. Het voordeel van gesimuleerde objecten is dat er veel meer vrijheid is om parameters te variëren, wat in de praktijk niet haalbaar is.

7.1

Invloed van de attenuatie

We onderzoeken eerst de invloed van de attenuatie. Hiervoor nemen we een homogeen sample om aan de voorwaarde van MP te voldoen. De attenuatiecoëfficiënt µ van het object wordt gevarieerd en de reconstructies met beide algoritmes worden vergeleken.

7.1.1

Simulaties

Er werden simulaties uitgevoerd om het effect van toenemende attenuatie op de reconstructie van een homogeen object te verifiëren. Het voordeel hiervan is dat we de gereconstrueerde doorsnedes kunnen vergelijken met de theoretische doorsnede. Hiervoor werd een homogene bol gesimuleerd. De parameters van de simulatie worden weergegeven in tabel 7.1. De fase werd gesimuleerd op basis van [32] λd ∂µθ ∂φθ ∂µθ ∂φθ λd 2 ∇ φθ + + , (7.1) Iθ,z=d ≈ Iθ,z=0 1 − 2πM 2πM ∂x ∂y ∂y ∂x wat een uitbreiding op de TIEW is. Vervolgens werden op de gesimuleerde beelden beide algoritmes toegepast. De gefilterde projecties werden daarna gereconstrueerd. De gerecon51

7.1. Invloed van de attenuatie

52

strueerde beelden werden vergeleken met de theoretische doorsnedes. De fout op de reconstructie werd gekwantificeerd d.m.v. de Normalized Root Mean Squared Error (NRMSE) [50]: q Pn

i=1 (x1,i −x2,i )

N RM SE =

n

xmax − xmin

2

.

(7.2)

Hierin wordt de som genomen over de n pixels, x1,i stelt de theoretische waarde voor en x2,i de gereconstrueerde waarde. De NRMSE wordt weergegeven voor beide filters als functie Tabel 7.1: Parameters voor de simulatie van een homogene bol

SOD SDD µ δ

10 mm 250 mm 0.05 mm−1 0.25 · 10−6

van de attenuatie in figuur 7.1. Hierin zien we dat de NRMSE toeneemt voor MBA als de attenuatie toeneemt. MP daarentegen ondervindt nagenoeg geen hinder van de toenemende attenuatie. In figuur 7.2 tonen we de genormeerde lijnprofielen1 van de gereconstrueerde centrale slice bij stijgende attenuatie (4.7 %, 17.7 % en 28.8 %). Deze figuur toont mooi hoe het cupping artefact toeneemt bij MBA met toenemende attenuatie, terwijl MP hier geen hinder van lijkt te ondervinden. We merken verder op dat MP een zekere attenuatie nodig heeft [34] om tot een goede reconstructie te komen: als de attenuatie te laag wordt neemt de fout toe. Dat volgt direct uit formule 6.25 die de geprojecteerde dikte geeft: bij µ = 0 krijgen we een singulariteit.

7.1.2

Hars

Experiment In praktijk is het zeer moeilijk een object te vinden dat homogeen is op de schaal die we hier beschouwen. Als homogeen object werd een epoxyhars gebruikt dat gestold werd in een cilindrische vorm. De attenuatiecoëfficiënt werd gevarieerd door verschillende X-stralen energieën te gebruiken. Het sample werd gescand bij buisspanningen van 50 kV, 80 kV en 110 kV. De attenuatie bedroeg respectievelijk 22 %, 17.5 % en 11 %. Aanvankelijk werd veel hinder ondervonden door het secondary spot effect (zie sectie 4.4), dat zeer goed zichtbaar was op het homogeen sample. Met de collimator werd dit effect zeer goed gereduceerd. Bij hogere energieën lijkt de collimator een kleine fractie van de secundaire straling door te laten, maar dit kan nog gecorrigeerd worden. 1

Met een lijnprofiel van een slice bedoelen we de waarden langs een lijn van een gereconstrueerde slice.

7.1. Invloed van de attenuatie

53

0.12

0.1

NRMSE

0.08

0.06

0.04

0.02 MBA MP 0 0

5

10

15

20

25

30

35

attenuatie (%)

Figuur 7.1: NRMSE van de centrale slice van een gesimuleerde homogene bol als functie van de attenuatie

Resultaat Op de projecties worden MBA en MP toegepast. De genormeerde lijnprofielen van een slice voor verschillende energieën worden weergegeven in figuur 7.3. Op de projecties bij 110 kV werd wel nog gecorrigeerd voor het resterend secondary spot effect. Het sample is homogeen, heeft een uniforme densiteit. Dit werd nagegaan door het sample in twee te zagen en te reconstrueren. We zien dat het cupping artefact bij elke energie groter is bij MBA dan bij MP. Het cupping artefact is te wijten aan de fasefilter. Als de attenuatie afneemt zien we dat de cupping bij MBA minder erg wordt, zoals we verwachten. Op de doorsnede gereconstrueerd met MP is ook een cupping artefact zichtbaar, maar dit is vermoedelijk door beam hardening (zoals uitgelegd in sectie 2.2.3), vermits op een gewone reconstructie ook een cupping artefact aanwezig is (zie doorsnede figuur 3.6d). We kunnen besluiten dat we kwalitatief hetzelfde resultaat terugvinden zoals in de simulaties.

7.2. Invloed van de homogeniteit

doorsnede

a)

54

1 0.9 0.8 MBA 4.7% attenuatie MP 4.7% attenuatie 0.7

doorsnede

b)

1 0.9 0.8 MBA 17.7% attenuatie MP 17.7% attenuatie 0.7

c) doorsnede

1 0.9 0.8

MBA 28.8% attenuatie MP 28.8% attenuatie

0.7 0

100

200

300

400

500

pixel

Figuur 7.2: Genormeerd lijnprofiel van de centrale slice van een gesimuleerde homogene bol bij verschillende attenuaties: (a) 4.7 %, (b) 17.7 % en (c) 28.8 %

7.2 7.2.1

Invloed van de homogeniteit Invloed van de homogeniteit op MBA

Vermits MBA een niet-attenuerend object veronderstelt, simuleren we een niet-attenuerende bol met 4 inclusies die een andere δ hebben. De doorsnede van de centrale slice wordt getoond in figuur 7.4. De bulk heeft δ0 = 0.25 · 10−6 . De meer dense inclusies hebben respectievelijk een δ van 4/3δ0 en 5/3δ0 . De minder dense inclusies hebben 2/3δ0 en 1/3δ0 als δ. De parameters zijn weergegeven in tabel 7.2. De reconstructie werd uitgevoerd voor 1 tot 4 inclusies. Dit stelt het minder homogeen worden voor. De reconstructies zijn weergegeven in figuur 7.5. We zien dat de reconstructie van δ uitgevoerd wordt en dat de inhomogeniteit dus geen probleem vormt voor MBA.


doorsnede

a)

55

1 0.9 0.8 0.7

MBA 50 kV MP 50 kV

0.6

doorsnede

b)

1 0.9 0.8 0.7

MBA 80 kV MP 80 kV

0.6

doorsnede

c)

1 0.9 0.8 0.7

MBA 110 kV MP 110 kV

0.6 0

100

200

300

400

500

600

700

800

900

pixel

Figuur 7.3: Lijnprofiel van een slice van een epoxyhars staafje bij verschillende buisspanningen: (a) 50 kV, (b) 80 kV en (c) 110 kV Tabel 7.2: Parameters van het gebruikte fantoom bij de simulatie van een inhomogeen, nietattenuerend faseobject

7.2.2

inclusie

δ

µ (mm−1 )

0 (bulk) 1 2 3 4

0.25 · 10−6 4/3δ0 5/3δ0 2/3δ0 1/3δ0

0 0 0 0 0

Invloed van de homogeniteit op MP

De methode van Paganin veronderstelt een homogeen object. Uit formule 6.25, die de geprojecteerde dikte geeft, volgt dat vooral de parameter µ/δ van belang is. Bij een homogeen object is dit uiteraard een constante doorheen het hele object, maar dit wil ook zeggen dat een inhomogeen object met een constante µ/δ correct zou gereconstrueerd


56

Figuur 7.4: De centrale slice van het gebruikte fasefantoom bij de evaluatie van de invloed van de homogeniteit op MBA

Figuur 7.5: Reconstructie van inhomogene faseobjecten met MBA

moeten worden. Hiervoor werd een fantoom analoog als in figuur 7.4 gesimuleerd, maar waarbij de bollen ook een attenuatie hebben zodat µ/δ een constante is (zie tabel 7.3). De attenuatiecoëfficiënt van de bulk µ0 bedraagt 0.15 mm−1 . Er werd vervolgens een simulatie uitgevoerd waarbij µ/δ niet constant was. Hiervoor werd δ0 = 0.25 · 10−6 gekozen en alle inclusies hebben δ = 3 δ0 (de parameters voor het fantoom zijn weergegeven in tabel 7.4). De reconstructies van het fantoom met een constante en verschillende µ/δ worden weergegeven in respectievelijk figuren 7.6a en 7.6b. In figuur 7.6a is te zien dat bij constante µ/δ de fasedistributie δ gereconstrueerd wordt, net zoals we verwachten. Bij verschillende µ/δ wordt de fasedistributie δ niet gereconstrueerd. We krijgen een reconstructie van µ (en niet δ) met faseringen. We kunnen dit als volgt begrijpen. Zolang het object homogeen is of µ en δ evenredig zijn doorheen het object, gaat de veronderstelling op dat voor beiden de ge¨ıntegreerde waarde gegeven wordt door een geprojecteerde dikte


57

Tabel 7.3: Parameters van het fantoom bij de simulatie van een inhomogeen object waarbij µ/δ constant is

inclusie

δ

µ (mm−1 )

0 (bulk) 1 2 3 4

0.25 · 10−6 4/3δ0 5/3δ0 2/3δ0 1/3δ0

0.15 4/3µ0 5/3µ0 2/3µ0 1/3µ0

Tabel 7.4: Parameters van het fantoom bij de simulatie van een inhomogeen object waarbij de verhouding µ/δ niet constant is

(a)

inclusie

δ

µ (mm−1 )

0 (bulk) 1 2 3 4

0.25 · 10−6 3δ0 3δ0 3δ0 3δ0

0.15 4/3µ0 5/3µ0 2/3µ0 1/3µ0

(b)

Figuur 7.6: Invloed van de homogeniteit op de reconstructiemethode van MP. (a) Reconstructie van fantoom met constante µ/δ, (b) reconstructie van een fantoom waarbij µ/δ niet constant is

T maal een constante µ of δ. Van zodra µ en δ niet meer evenredig variëren, kan er geen T gevonden worden die zowel de faseverschuiving als de attenuatie juist beschrijft, zoals weergegeven in vergelijkingen 6.18 en 6.19. Blijkbaar domineert de attenuatie dan in de berekening voor de geprojecteerde dikte en wordt die grootheid gereconstrueerd (zij het wel


58

met faseartefacten). We kunnen besluiten dat MP de grootheden µ en δ juist reconstrueert zolang de verhouding µ/δ constant is.

7.2.3

Verschillende verhouding bij MBA

Vervolgens werd ook de MBA-filter toegepast op een object waarbij de verhouding µ/δ constant is en waarbij µ/δ niet constant is, en µ verschillend van 0. Hiervoor werd hetzelfde fantoom met 4 inclusies als bij MP gebruikt (zie tabellen 7.3 en 7.4). Als de verhouding constant is, lukt de reconstructie zeer goed (zie figuur 7.7a). Wanneer de verhouding niet constant is, verschijnen de faseranden (figuur 7.7b). Het valt opnieuw op te merken dat de fasedistributie niet teruggevonden wordt (δbulk = 0.25 · 10−6 en δinclusies = 3 δbulk ). We zien dus dat ook MBA faalt bij een verschillende µ/δ. Het is duidelijk dat de attenuatie domineert in het signaal en de gereconstrueerde grootheid2 is µ i.p.v. δ.

(a)

(b)

Figuur 7.7: Invloed van de verhouding µ/δ op MBA. (a) Reconstructie van fantoom met constante µ/δ, (b) reconstructie van fantoom met verschillende µ/δ met dezelfde reconstructieparameter α

Absolute waarde gereconstrueerde grootheid Zoals uitgelegd in sectie 6.1.3, wordt bij de implementatie een attenuatieparameter A ge¨ıntroduceerd omdat de bekomen grootheid anders veel te klein is om in de reconstructie te gebruiken (zie vergelijking 6.16). De waarden die we bekomen in de reconstructie zijn dus geschaald. Het is echter nuttig eens te kijken naar de absolute waarde van de gereconstrueerde grootheid en of deze nog gerelateerd is aan δ. Bij de reconstructies in 2

µ.

Merk op dat dit niet de absolute waarde van µ is, maar de gereconstrueerde grootheid is evenredig met


59

figuren 7.7a en 7.7b bekomen we voor de bulk een waarde gelijk aan 1.15. Gegeven dat de attenuatieparameter A = 5000 en rekening houdend met de factor (4π 2 d)−1 in formule 6.13 die niet in de uitdrukking van het klassieke FBP-algoritme aanwezig is (zie formule 2.6), bekomen we δ = 0.24 · 10−6 . Als µ/δ constant is, wordt δ juist gereconstrueerd. Als de verhouding niet constant is, domineert de attenuatie en de gereconstrueerde δ’s verhouden zich volgens µ. Vermits we ons vooral baseren op de relatieve verhoudingen van de gereconstrueerde grootheid van de verschillende componenten, en deze volgens µ blijkt te zijn, zullen we dit aanduiden als een reconstructie van µ. We zien dus hetzelfde gedrag als bij MP, enkel de absolute waarde van de gereconstrueerde grootheid is anders. Homogene bol met fase-inclusie Het is eerder verrassend dat de gereconstrueerde grootheden van de verschillende inclusies zich verhouden volgens µ i.p.v. δ, vermits in de afleiding van het Bronnikov Algoritme expliciet δ gereconstrueerd wordt. We gaan nu na wanneer de attenuatie voldoende laag wordt opdat we toch δ reconstrueren i.p.v. µ. Hiervoor werd een bol gesimuleerd die één fase-inclusie heeft, en homogeen is in zijn attenuatie (parameters weergegeven in tabel 7.5). De attenuatie werd achtereenvolgens 5 %, 2 %, 1 %, 0.5 %, 0.1 % en 0.01 % gekozen. De gereconstrueerde centrale slice evenals het lijnprofiel zijn weergegeven in figuur 7.8. Bij een attenuatie van 5 % verschijnen opnieuw de fase artefacten in de reconstructie. We zien dat de attenuatie domineert en dat de gereconstrueerde grootheid µ is i.p.v. δ. Vanaf een attenuatie van minder dan 1 % lijkt het fasesignaal steeds belangrijker te worden in de reconstructie en bij een attenuatie van 0.1 % verhouden de inclusies zich al vrij goed zoals het refractive index decrement δ. Bij hogere attenuaties echter is de attenuatie dominant en bepaalt die onze reconstructie. Tabel 7.5: Parameters van het fantoom bij de simulatie van een homogeen attenuerende bol met 1 fase-inclusie

inclusie

δ

µ (mm−1 )

0 (bulk) 1

0.25 · 10−6 4/3δ0

µ0 µ0

Kleinere variatie in de parameters De waarde voor δ gekozen in de voorgaande simulaties voor de inclusies en de bulk waren een factor drie verschillend. De simulatie werd nog eens herhaald voor een kleiner verschil tussen de bulk en de inclusies: δinclusie = 1.3δbulk . De verhouding µ/δ is nog steeds niet constant, maar de afwijking t.o.v. de parameters van de bulk zullen kleiner zijn. De resultaten voor MBA en MP zijn respectievelijk weergegeven in figuren 7.9 en 7.10. In


60

Figuur 7.8: Homogene bol met 1 fase-inclusie gereconstrueerd met MBA. De attenuatie bedraagt maximaal (a) 5 %, (b) 2 %, (c) 1 %, (d) 0.5 % , (e) 0.1 % en (f) 0.01 %. Het getoonde lijnprofiel is de verticale door het midden van de slice


61

figuren 7.9a en 7.10a is de reconstructie te zien van het object met δinclusie = 3δbulk . Buiten het feit dat de attenuatie in plaats van de fase teruggevonden wordt, zijn er storende ringen rond de inclusies, wat ook geen interpretatie op basis van de attenuatie toelaat. Hieruit kan geen structuur afgeleid worden van het object. In figuren 7.9b en 7.10b zijn de parameters aangepast om de storende ringen weg te werken. Dit lukt maar gedeeltelijk (zie ook de doorsnede) en het beeld wordt zeer wazig. In figuren 7.9c en 7.10c is de reconstructie voor δinclusie = 1.3δbulk weergegeven. We zien dat het toch nog lukt om een degelijke reconstructie te maken met de bulkparameters. Men moet wel in gedachten houden dat de gereconstrueerde grootheid de attenuatie3 is i.p.v. de fasedistributie zoals de theorie zegt. We merken wel op dat bij MBA (zie figuur 7.9) de verhouding langs de randen lichtjes afwijken van deze verhouding.

7.2.4

Inhomogeen sample: hars met grafietstaafje

Vervolgens werd een opzettelijk inhomogeen sample gemaakt: hetzelfde epoxyhars staafje als hierboven maar deze keer met een grafiet staafje in het centrum. Het object werd weer ingescand bij drie verschillende energieën, en telkens gereconstrueerd met MBA en MP. De gereconstrueerde slices zijn weergegeven in figuur 7.11. De lijnprofielen van de gereconstrueerde slices met MBA en MP zijn weergegeven in figuur 7.12. De doorsnedes werden genormeerd op de attenuatiecoëfficiënt van het centrale koolstofstaafje. Wat opvalt is dat het cupping effect weer groter is bij MBA en afneemt met toenemende energie. Wat nog opvalt is dat bij toenemende energie de attenuatie van het koolstofstaafje relatief meer afneemt dan die van het omgevende hars. Het object is niet homogeen, wat een verschillende µ/δ kan betekenen (hoewel we dit niet met zekerheid kunnen zeggen), bovendien evolueert de attenuatiecoëfficiënt van de beide materialen als functie van de energie anders. In de veronderstelling dat δ weinig verandert over het energie interval dat we hier beschouwen en µ zichtbaar verandert, is het dus zeker niet uitgesloten dat er effectief een verschillende µ/δ is. Toch zien we dat de slice vrij goed gereconstrueerd wordt. Uit onze simulaties is gebleken dat dit mogelijk is zolang µ/δ niet te sterk fluctueert. Opnieuw zien we toch weer dat vooral aan de rand van het object de verhouding bij MBA lichtjes afwijkt. Het dal in de dichtheid rond pixel 200 is een luchtbel. We zien dat dit niet perfect als intensiteit 0 reconstrueerd wordt. Dit is een goede indicator voor de kwaliteit van een reconstructie. Bij MP blijkt dit iets beter te zijn dan bij MBA. 3

Bij MBA vindt men niet de absolute waarde van µ terug, maar de verhoudingen van de gereconstrueerde grootheid van de verschillende componenten verhouden zich volgens µ (zoals hiervoor reeds besproken werd). Bij MP vinden we effectief de waarde van µ terug, wat ons niet hoeft te verwonderen gezien onze implementatie (zoals uitgelegd in paragraaf 6.2.2).


62

doorsnede

a) 1.2 1 0.8 0.6 0.4 0.2 0 0

100

200

300

400

500

300

400

500

300

400

500

pixel

b) doorsnede

1 0.8 0.6 0.4 0.2 0 0

100

200 pixel

c) doorsnede

1 0.8 0.6 0.4 0.2 0 0

100

200 pixel

Figuur 7.9: MBA toegepast op een inhomogeen object. Het getoonde lijnprofiel is een horizontale lijn door het midden van de slice. (a) δinclusie = 3δbulk en reconstructie met α = 6000, (b) δinclusie = 3δbulk en α zo gekozen dat de artefacten verminderen, (c) δinclusie = 1.3δbulk en reconstructie met α = 6000

7.2.5

Bespreking

Uit het voorgaande blijkt dat de basisveronderstellingen bij beide methodes effectief zichtbaar zijn. MBA veronderstelt een niet-attenuerend object en geeft dan ook goede resultaten voor niet-attenuerende samples, de fasedistributie wordt teruggevonden. Uiteraard komt een dergelijk object in de praktijk niet voor. In theorie moet men dan twee reeksen projecties nemen om de attenuatie weg te delen. Om slechts één reeks projecties te moeten uitvoeren, wordt een correctieparameter α ingevoerd om te corrigeren voor de resterende attenuatie. Bij de reconstructies wordt echter wel een cupping artefact waargenomen dat groter wordt naarmate er meer attenuatie is. De methode van Paganin veronderstelt een homogeen sample, dat een zekere attenuatie nodig heeft om tot een goede reconstructie te komen. Deze methode ondervindt dit cupping artefact veel minder, in de simulaties is het zelfs niet aanwezig. MP reconstrueert de fase


63

a) doorsnede

1 0.8 0.6 0.4 0.2 0 0

100

200

300

400

500

300

400

500

300

400

500

pixel

b) doorsnede

1 0.8 0.6 0.4 0.2 0

100

200 pixel

c) doorsnede

1 0.8 0.6 0.4 0.2 0 0

100

200 pixel

Figuur 7.10: MP toegepast op een inhomogeen object. Het getoonde lijnprofiel is een horizontale lijn door het midden van de slice. (a) δinclusie = 3δbulk en reconstructie met de bulkparameters, (b) δinclusie = 3δbulk en de parameters zo gekozen dat de artefacten verminderen, (c) δinclusie = 1.3δbulk en reconstructie met de bulkparameters

(a)

(b)

(c)

Figuur 7.11: Gereconstrueerde slice van een epoxyhars met grafietstaafje. Reconstructie met (a) FBP, (b) MBA en (c) MP


64

a) 1

cross sectie

0.8 0.6

MBA 50 kV MBA 80 kV MBA 110 kV

0.4 0.2 0 b) 1

cross sectie

0.8 0.6

MP 50 kV MP 80 kV MP 110 kV

0.4 0.2 0 0

100

200

300

400 pixel

500

600

700

800

Figuur 7.12: Een grafietstaafje in epoxyhars als inhomogeen object. Het object werd ingescand bij 50 kV, 80 kV en 110 kV. (a) Doorsnede (volgens de getoonde lijn) van een gereconstrueerde slice met MBA, (b) doorsnede van dezelfde slice gereconstrueerd met MP

zeer goed zolang de verhouding µ/δ constant is. Dit volgt direct uit formule 6.25 die de geprojecteerde dikte geeft. We merken op dat het reconstrueren van µ dus equivalent is met het reconstrueren van δ (ze zijn evenredig). Wanneer de verhouding µ/δ niet constant is en sterk fluctueert, ontstaan er artefacten in de reconstructie. We zien dat de fase niet gereconstrueerd wordt, maar dat we een reconstructie bekomen van µ met artefacten. Wanneer µ/δ niet te sterk varieert kan toch nog een degelijke reconstructie van µ gemaakt worden. Hieruit blijkt dus dat de attenuatie domineert en dat we eigenlijk de grootheid µ reconstrueren. Wanneer µ/δ niet constant is blijkt ook MBA hier moeilijkheden mee te hebben. Er kan geen goede correctieparameter α gevonden worden om de fase te reconstrueren. Dit is eerder verrassend vermits dit geen voorwaarde is bij MBA. Een mogelijke verklaring voor dit gedrag is de volgende [16]. Beide methodes maken gebruik van een low-pass filter van

7.3. Ruis

65

de vorm f (ξ, η) =

ξ2

a . + η2 + b

(7.3)

Bij MP weerspiegelt b de verhouding µ/δ, bij MBA de parameter α. De simulaties impliceren dat indien er geen constante µ/δ is, er ook geen correcte parameter α kan gevonden worden. Bovendien vinden we niet het refractive index decrement δ terug, maar de gereconstrueerde grootheid is evenredig met µ en er zijn artefacten. Indien µ/δ niet te sterk varieert, kan ook hier toch nog een degelijke reconstructie gemaakt worden van µ. In de afleiding van het Bronnikov algoritme (zie sectie 6.1) wordt expliciet het refractive index decrement δ berekend. Het lijkt dan ook vreemd dat, wanneer er toch attenuatie is, de attenuatiecoëfficiënt µ gereconstrueerd wordt. We mogen echter niet vergeten dat we de methode eigenlijk toepassen op een verkeerd argument: de methode vraagt de functie gθ (zie formule 6.6), die een puur fasesignaal verminderd met 1 voorstelt. De parameter α die ingevoerd werd, is bedoeld om de attenuatie te onderdrukken. Het gewicht van de lage frequentiecomponenten van het signaal wordt inderdaad lager dan wanneer de parameter α nul is, maar in de filter 7.3 die toegepast wordt krijgen de lage frequentiecomponenten nog steeds het hoogste gewicht (zie figuur 7.13), waardoor het te begrijpen is dat de attenuatie toch nog een belangrijke rol speelt. Slechts bij zeer lage attenuatie wordt δ gereconstrueerd, in de andere gevallen domineert de attenuatie en zolang de verhouding µ/δ niet te sterk fluctueert worden goede reconstructies bekomen van µ. Voor praktische toepassingen waar meestal een zekere attenuatie aanwezig is, lijkt MP dus in het algemeen betere resultaten te geven. De gereconstrueerde grootheid is µ en het cupping artefact is veel minder aanwezig dan bij MBA.

7.3

Ruis

Zoals eerder vastgesteld werd [32, 34, 51, 52], vertonen beide algoritmes een goede stabiliteit onder ruis. Dit werd nagegaan door een bol te simuleren met twee inclusies en twee holtes. De bulk van het materiaal had een attenuatiecoëfficiënt µ0 = 0.3 mm−1 , de inclusies 3µ0 . De maximale attenuatie bedroeg 26 %. De verhouding µ/δ werd constant gehouden; δ0 = 3 · 10−6 voor de bulk en 3δ0 voor de inclusies.

7.3.1

Simulatie van de ruis

De ruis werd gesimuleerd door random getallen te genereren volgens een Poissonverdeling. We lichten dit even toe. Als het verwachte aantal tellen in een pixel van de detector λ is, dan wordt de kans op k tellen gegeven door p(k, λ) =

λk e−λ . k!

(7.4)

7.3. Ruis

66

1/(ξ 2+η2+0.001)

1000 900 800 700 600 500 400 300 200 100 0

1000 900 800 700 600 500 400 300 200 100 0 0.4 0.2 -0.4

0 -0.2

0 ξ

-0.2 0.2

η

-0.4

0.4

Figuur 7.13: De filterfunctie in het Bronnikov algoritme

Als we zonder object op de detector een intensiteit I0 meten en in dat geval λ tellen verwachten, dan is het verwachte aantal tellen bij een gemeten intensiteit Im gelijk aan λ · Im . Het aantal tellen in de detectorpixel wordt dus gevonden door een random getal te genereren volgens een Poissondistributie met verwachtingswaarde λ · Im . Als we een hoge ruis willen simuleren komt dit overeen met een klein aantal tellen in een pixel. Een random getal genereren volgens een Poisson distributie kan gebeuren met de inversie methode. In de implementatie voor het simuleren van de ruis werd echter gebruik gemaakt van een ingebouwde functie in LabVIEW die een getal genereert volgens een Poissondistributie met verwachtingswaarde λ · Im .

®

7.3.2

Resultaat

De simulaties werden uitgevoerd voor verscheidene ruisniveaus. Er werd eerst ruis gesimuleerd op projecties zonder fase en die werden vervolgens gereconstrueerd met het FBPalgoritme. Vervolgens werd ruis gesimuleerd op gemengde projecties, waarbij attenuatie en fase aanwezig was. Hier werd de fase gereconstrueerd met MBA en MP. De fout op de reconstructie werd gekwantificeerd door de NRMSE te berekenen. Het resultaat is weergegeven in figuur 7.14. We zien duidelijk dat MBA en MP een beter resultaat geven bij weinig tellen. Het valt op te merken dat de fout bij MBA steeds een weinig hoger ligt dan bij MP. Dit is te wijten aan het cupping artefact (zoals reeds besproken) door de

7.3. Ruis

67

0.6 FBP MBA MP

NRMSE

0.4

0.2

0 100

1000

10000

100000

tellen

Figuur 7.14: De NRMSE van een gereconstrueerde slice als functie van het gemiddeld aantal tellen op een detectorpixel

resterende attenuatie. De gereconstrueerde slices zijn weergegeven in figuur 7.15. Opnieuw

Figuur 7.15: De gereconstrueerde slices bij gemiddeld 316 tellen per pixel: (a) theoretisch fantoom, (b) reconstructie met FBP, (c) reconstructie met MBA, (d) reconstructie met MP

zien we dus vrij gelijkaardig gedrag tussen MBA en MP. Een verklaring voor de stabiliteit bij MBA werd reeds gegeven in [16]: de methode integreert de waargenomen pieken van de fase4 en is zo minder afhankelijk van de ruis. Deze integratie volgt niet direct uit de theorie van MP. Wat we wel kunnen opmerken is dat de filterfunctie 7.3 die beide algoritmes ge4

zie vergelijkingen 6.7 en 6.11.

7.4. Toepassingen

68

bruiken een lowpass filter is, die ruis onderdrukt. Dit verklaart de relatieve ongevoeligheid van beide algoritmes tegen ruis.

7.4

Toepassingen

In een laatste gedeelte passen we beide algoritmes toe op enkele objecten uit de praktijk.

7.4.1

Pilletje

Een ingescand pilletje heeft een maximale attenuatie van ongeveer 30 % en bevat een zeer sterk fasesignaal. De parameters van de scan zijn weergegeven in tabel 7.6. Bij een klassieke reconstructie werkt dit fasesignaal zeer storend, het lijkt alsof er een extra laagje, een coating, op de buitenkant ligt. Dit is weergegeven in figuur 7.16a. Vervolgens werden ook MBA en MP toegepast, waarvan voor elk een gereconstrueerde slice is weergegeven in respectievelijk figuren 7.16b en 7.16c. Het is duidelijk dat hier geen storende rand meer aanwezig is. De structuren zijn nog duidelijk onderscheidbaar bij de reconstructies met MBA en MP, maar de doorsnedes lijken toch iets minder scherp dan in de klassieke reconstructie. Onderaan figuur 7.16 is de doorsnede gegeven van de reconstructies met MBA en MP. Hierin — en ook op de slices — is duidelijk te zien dat bij MBA een sterker cupping artefact aanwezig is dan bij MP: er is een maximale variatie van 5 % in de gereconstrueerde waarde. Tabel 7.6: Parameters voor de scan van het pilletje

# projecties # detector rijen # detector kolommen pixelgrootte (µm) SOD (mm) SDD (mm) Voxelgrootte (µm)

7.4.2

800 644 947 38 10 270 1.4

Diamant

Een diamant werd ingescand en heeft een maximale attenuatie van ongeveer 30 % en een sterk fasesignaal. De parameters van de scan zijn weergegeven in tabel 7.7. Dit sample is zeer homogeen, wat dus voldoet aan de voorwaarde voor MP. De reconstructies met het FBP-algoritme, MBA en MP zijn respectievelijk weergegeven in figuren 7.17a, 7.17b en 7.17c. De doorsnedes van de reconstructies zijn weergegeven onderaan figuur 7.17. Het is opnieuw duidelijk dat MBA een véél sterker cupping artefact vertoont: de gereconstrueerde

7.4. Toepassingen

69

(a)

(b)

(c)

Figuur 7.16: Reconstructie van het pilletje met (a) het FBP algoritme, (b) MBA en (c) MP. Onderaan is een lijndoorsnede (volgens de getoonde lijn) weergegeven voor MBA en MP.

waarde bij MBA ligt maximaal 10 % lager dan bij MP. We merken opnieuw dat MBA en MP de neiging hebben een iets waziger beeld op te leveren. Zo zijn de randen bij MBA en MP uitgesmeerd over 10 pixels op de reconstructie, tegenover 6 bij het klassiek algoritme. Tabel 7.7: Parameters voor de scan van de diamant


7.4.3

1000 1820 1450 127 41.5 857.5 6.15

Houtschimmel

Een scan van een stukje dennenhout met een houtschimmel is weergegeven in figuur 7.18. De parameters van de scan zijn weergegeven in tabel 7.8. Opnieuw werden de klassieke

7.4. Toepassingen

(a)

70

(b)

(c)

Figuur 7.17: Reconstructie van een diamant met (a) het FBP algoritme, (b) MBA en (c) MP. Onderaan is een doorsnede weergegeven volgens de getoonde lijn.

reconstructie (figuur 7.18a), MBA (7.18b) en MP (7.18c) vergeleken. Het is duidelijk dat de schimmel beter zichtbaar is bij MBA en MP dan bij het klassieke reconstructie algoritme. Opnieuw lijken MBA en MP iets waziger, maar het contrast is wel beter. Onderaan figuur 7.18 is het aangeduide stukje doorsnede zoals in 7.18a weergegeven. Dat deze structuren nog kunnen onderscheiden worden, heeft uiteraard te maken met het feit dat we het fasesignaal gebruiken, maar ook met het feit dat de ruis veel meer onderdrukt wordt bij MBA en MP.

7.4.4

Vliegenpoot

Een sample dat al meer gebruikt werd in het faseonderzoek aan UGCT is een vliegenpoot. Op dit object is een duidelijk fasesignaal aanwezig. De instellingen van de scan zijn weergegeven in tabel 7.9. In figuur 7.19 is links een 3D-beeld van een reconstructie met het FBP-algoritme. Het faseartefact zorgt ervoor dat het lijkt alsof de wand van de poot hol

7.4. Toepassingen

(a)

71

(b)

(c)

Figuur 7.18: Reconstructie van een stukje dennenhout met (a) FBP, (b) MBA en (c) MP. Onderaan is voor elke methode het stuk doorsnede weergegeven zoals in (a) aangeduid. Tabel 7.8: Parameters voor de scan van het stukje dennenhout


800 800 970 127 3.29 164.28 0.6

is. De reconstructie met MBA (figuur 7.19 rechts) toont dan dat dit niet zo is. We hebben dan ook MP toegepast op de poot. De reconstructies met MBA en MP zijn weergegeven in figuur 7.20. Hierin zien we dat beide methodes zeer gelijkaardige resultaten geven, die iets minder scherp blijken dan de reconstructie met het FBP-algoritme.

7.4. Toepassingen

72

Tabel 7.9: Parameters voor de scan van de vliegenpoot


600 483 322 0.076 5 245 0.001551

Figuur 7.19: Reconstructie van een vliegenpoot met FBP (links) en MBA (rechts)

(a) MBA

(b) MP

Figuur 7.20: Reconstructie van een vliegenpoot met MBA en MP. De resultaten zijn zéér gelijkaardig

Hoofdstuk 8

Besluit In deze masterproef hebben we fasecontrast aan een hoge resolutie X-stralenbuis bestudeerd. In een eerste deel hebben we geprobeerd het fenomeen te simuleren op basis van de golfoptica. Een goede simulatie van het fenomeen is noodzakelijk om het te kunnen implementeren in een iteratief reconstructie algoritme. De simulatie gaf kwalitatief goede resultaten met realistische parameterwaarden. De stapgrootte in de berekeningen moet echter veel kleiner zijn dan de pixelgrootte, wat enorm veel rekentijd zou vragen en de implementatie in een iteratief reconstructie algoritme ook een stuk moeilijker zou maken. De praktische toepasbaarheid van de simulatie is dus beperkt. In een tweede deel hebben we twee fasereconstructie algoritmes, met name het Modified Bronnikov Algoritme (MBA) en de methode van Paganin (MP), ge¨ımplementeerd en vergeleken. Deze methodes lijken sterk op elkaar, hoewel ze van verschillende veronderstellingen vertrekken en op een verschillende manier afgeleid worden. MBA veronderstelt een niet-attenuerend object. Voor de resterende attenuatie wordt bij MBA gecorrigeerd door een extra parameter α in te voeren in de filter. Bij een puur faseobject wordt het refractive index decrement δ gereconstrueerd. Bij een sample dat nog attenuatie vertoont blijkt de attenuatie nog steeds te domineren en is de gereconstrueerde grootheid evenredig met de attenuatiecoëfficiënt µ. Bovendien wordt een cupping artefact waargenomen dat toeneemt met toenemende attenuatie. Zolang de verhouding µ/δ niet te sterk varieert worden geen noemenswaardige artefacten waargenomen in de reconstructie van µ. MP veronderstelt een homogeen object. Uit de filter blijkt dat ook inhomogene objecten met een constante verhouding µ/δ correct zouden moeten gereconstrueerd worden. Bij dergelijke objecten resulteert het reconstrueren van µ in hetzelfde resultaat als het reconstrueren van δ vermits ze evenredig zijn. We zien dat de reconstructie van δ (en dus ook µ) goed uitgevoerd wordt. Deze methode ondervindt geen hinder van de attenuatie, boven73

74 dien is een zekere attenuatie nodig om tot een goede reconstructie te komen. Wanneer de verhouding µ/δ niet constant is en sterk varieert, worden artefacten waargenomen in de reconstructie. Dan zien we ook dat de grootheid µ gereconstrueerd wordt. Waneer deze verhouding niet te sterk fluctueert doorheen het object kan toch nog een goede reconstructie van µ gemaakt worden. Beide methodes vertonen een betere stabiliteit tegen ruis dan een klassieke reconstructie. De twee methodes werden succesvol toegepast op enkele samples uit de praktijk. Algemeen kunnen we besluiten dat beide methodes vrij gelijkaardig zijn, maar dat MP beter overweg kan met de attenuatie dan MBA. Dit is uiteraard interessant vermits een niet-attenuerend object in de praktijk niet voorkomt. Bij verschillende µ/δ hebben beide methodes wel problemen met de reconstructie, maar wanneer deze verhouding niet te sterk varieert kunnen toch nog goede reconstructies van µ bekomen worden. Bovendien onderdrukken beide algoritmes de ruis zeer goed.

Bijlage A

Fantoombeschrijving In hoofdstuk 7 wordt gebruik gemaakt van gesimuleerde data van een sferisch fantoom. De parameters van het fantoom zijn weergegeven in tabel A.1. Er kunnen tot vier inclusies geplaatst worden. Een homogeen object bekomen we dus door geen inclusies toe te voegen. Een ellipso¨ıde heeft als co¨ ordinaten een middelpunt (x0 , y0 , z0 ) en drie halve assen a,b en c (zie figuur A.1). Hier zijn de inclusies sferisch gekozen (a = b = c). In figuren A.2a en A.2b zijn respectievelijk één projectie en de gereconstrueerde centrale slice van het fantoom getoond.

Figuur A.1: De coördinaten van een ellips [16]

Tabel A.1: Parameters voor het sferisch fantoom met vier inclusies

Inclusie

a

b

c

x0

y0

z0

µ

δ

bulk

0.95

0.95

0.95

0

0

0

µ0

δ0

1

0.2

0.2

0.2

0

0.6

0

µ1

δ1

2

0.2

0.2

0.2

-0.6

0

0

µ2

δ2

3

0.2

0.2

0.2

0

-0.6

0

µ3

δ3

4

0.2

0.2

0.2

0.6

0

0

µ4

δ4

75

76

(a)

(b)

Figuur A.2: Projectie (a) en gereconstrueerde centrale doorsnede (b) van de inhomogene sfeer. De attenuatiecoëfficiënten van de inclusies bedragen µ1 = 4/3µ0 , µ2 = 5/3µ0 , µ3 = 2/3µ0 en µ4 = 1/3µ0

Bibliografie ¨ [1] J. Radon. Uber die Bestimmung von Funktionen durch ihre Integralwerte l¨ angs gewisser Mannigfaltigkeiten. Akad. Wiss., 69:262–277, 1917. [2] A. C. Kak en Malcolm Slaney. Principles of Computerized Tomographic Imaging. IEEE Press, New York, 1988. [3] Wikipedia. Godfrey Hounsfield, 2011. http://en.wikipedia.org/wiki/Godfrey_Newbold_Hounsfield. [4] T.J. Davis, D. Gao, T.E. Gureyev, A.W. Stevenson, en S.W. Wilkins. Phase-contrast Imaging of Weakly Absorbing Materials using Hard X-rays. Nature, 373(6515):595– 598, FEB 16 1995. ISSN 0028-0836. [5] A. Snigirev, I. Snigireva, V. Kohn, S. Kuznetsov, en I. Schelokov. On the possibilities of x-ray phase contrast microimaging by coherent high-energy synchrotron radiation. Review of Scientific Instruments, 66(12):5486–5492, DEC 1995. ISSN 0034-6748. [6] S.W. Wilkins, T.E. Gureyev, D. Gao, A. Pogany, en A.W. Stevenson. Phase-contrast imaging using polychromatic hard X-rays. Nature, 384(6607):335–338, NOV 28 1996. ISSN 0028-0836. [7] P. Cloetens, R. Barrett, J. Baruchel, J.P. Guigay, en M. Schlenker. Phase objects in synchrotron radiation hard x-ray imaging. Journal of Physics D-Applied Physics, 29 (1):133–146, JAN 14 1996. ISSN 0022-3727. [8] Oliver Betz, Ulrike Wegst, Daniela Weide, Michael Heethoff, Lukas Helfen, Wah-Keat Lee, en Peter Cloetens. Imaging applications of synchrotron X-ray phase-contrast microtomography in biological morphology and biomaterials science. 1. General aspects of the technique and its advantages in the analysis of millimetre-sized arthropod structure. Journal of Microscopy-Oxford, 227(1):51–71, JUL 2007. ISSN 0022-2720. [9] F. Pfeiffer, T. Weitkamp, en C. David. X-ray phase contrast imaging using a grating interferometer. Europhysics News, 37(5):13–15, 2006. ISSN 0531-7479.

77

Bibliografie

78

[10] Wikipedia. X-ray, 2011. http://en.wikipedia.org/wiki/X-ray. [11] Wikipedia. Synchrotron radiation, 2011. http://en.wikipedia.org/wiki/Synchrotron_radiation. [12] Wikipedia. Synchrotron, 2011. http://en.wikipedia.org/wiki/Synchrotron. [13] Wikipedia. Bremsstrahlung, 2011. http://en.wikipedia.org/wiki/Brehmstrahlung. [14] Wikipedia. X-ray tube, 2011. http://en.wikipedia.org/wiki/X-ray_tube. [15] Glenn F. Knoll. Radiation Detection and Measurement. Wiley, 2000. [16] Yoni De Witte. Improved and practically feasible reconstruction methods for high resolution X-ray tomography. PhD thesis, Universiteit Gent, 2010. [17] J. Vlassenbroeck, M. Dierick, B. Masschaele, V. Cnudde, L. Van Hoorebeke, en P. Jacobs. Software tools for quantification of X-ray microtomography. Nuclear Instruments & Methods in Physics Research Section A-Accelerators Spectrometers Detectors and Associated Equipement, 580(1):442–445, SEP 21 2007. ISSN 0168-9002. 10th International Symposium on Radiation Physics, Coimbra, PORTUGAL, SEP 17-22, 2006. [18] Christian David Franz Pfeiffer, Timm Weitkamp. X-ray phase contrast imaging using a grating interferometer. Europhysicsnews, 37:13–15, 2006. [19] M. El-Ghazaly, H. Backe, W. Lauth, G. Kube, P. Kunz, A. Sharafutdinov, en T. Weber. X-ray phase contrast imaging at mami. European Physical Journal A, 28:197–208, 2006. [20] M. Langer, P. Cloetens, J. P. Guigay, en F. Peyrin. Quantitative comparison of direct phase retrieval algorithms in in-line phase tomography. Medical Physics, 35(10):4556– 4566, 2008. [21] F. Arfelli, M. Assante, V. Bonvicini, A. Bravin, G. Cantatore, E. Castelli, L. Dalla Palma, M. Di Michiel, R. Longo, A. Olivo, S. Pani, D. Pontoni, P. Poropat, M. Prest, A. Rashevsky, G. Tromba, A. Vacchi, E. Vallazza, en F. Zanconati. Lowdose phase contrast x-ray medical imaging. Physics in Medicine and Biology, 43(10): 2845–2852, 1998.

Bibliografie

79

[22] Liberato De Caro, Cinzia Giannini, Roberto Bellotti, en Sabina Tangaro. A theoretical study on phase-contrast mammography with Thomson-scattering x-ray sources. Medical Physics, 36:4644–4653, 2009. ISSN 0094-2405. [23] Konic Minolta Medical en Inc Graphic. Phase Contrast Technology, 2011. http://www.konicaminolta.com/healthcare/technology/phasecontrast. [24] P. Cloetens, M. PateyronSalome, J. Y. Buffiere, G. Peix, J. Baruchel, F. Peyrin, en M. Schlenker. Observation of microstructure and damage in materials by phase sensitive radiography and tomography. Journal of Applied Physics, 81(9):5878–5886, 1997. [25] M. Alonso en E.J. Finn. Fundamentele natuurkunde, deel 3 Golven. Delta Press, 1997. [26] Center for X-ray Optics en Lawrence Berkeley National Laboratory Advanced Light Source. X-ray Data Booklet, 2009. http://xdb.lbl.gov/. [27] NIST. X-Ray Form Factor, Attenuation and Scattering Tables, 2011. http://physics.nist.gov/PhysRefData/FFast/html/form.html. [28] Nikolay Kardjilov. Further developments and applications of radiography and tomography with thermal and cold neutrons. PhD thesis, Technishce Universität M¨ unchen, 2003. [29] Michael Reed Teague. Deterministic phase retrieval: a green’s function solution. Journal of Optical Society of America, 73(11), 1983. [30] T. E. Gureyev, A. Roberts, en K. A. Nugent. Partially coherent fields, the transportof-intensity equation, and phase uniqueness. Journal of the Optical Society of America a-Optics Image Science and Vision, 12(9):1942–1946, 1995. [31] A. Peterzol, A. Olivo, L. Rigon, S. Pani, en D. Dreossi. The effects of the imaging system on the validity limits of the ray-optical approach to phase contrast imaging. Medical Physics, 32(12):3617–3627, 2005. [32] A. V. Bronnikov. Theory of quantitative phase-contrast computed tomography. Journal of the Optical Society of America a-Optics Image Science and Vision, 19(3):472– 480, 2002. [33] A. V. Bronnikov. Reconstruction formulas in phase-contrast tomography. Optics Communications, 171(4-6):239–244, 1999.

Bibliografie

80

[34] D. Paganin, S. C. Mayo, T. E. Gureyev, P. R. Miller, en S. W. Wilkins. Simultaneous phase and amplitude extraction from a single defocused image of a homogeneous object. Journal of Microscopy-Oxford, 206:33–40, 2002. [35] A. Momose. Phase-sensitive imaging and phase tomography using X-ray interferometers. Optics Express, 11(19):2303–2314, SEP 22 2003. ISSN 1094-4087. [36] Peter Cloetens. Contribution to Phase Contrast Imaging, Reconstruction and Tomography with Hard Synchrotron Radiation. PhD thesis, Vrije Universiteit Brussel, ESRF, 1999. [37] E. Pagot, S. Fiedler, P. Cloetens, A. Bravin, P. Coan, K. Fezzaa, J. Baruchel, en J. Hartwig. Quantitative comparison between two phase contrast techniques: diffraction enhanced imaging and phase propagation imaging. Physics in Medicine and Biology, 50(4):709–724, FEB 21 2005. ISSN 0031-9155. [38] F. Krejci, J. Jakubek, en M. Kroupa. Single grating method for low dose 1-d and 2-d phase contrast x-ray imaging. 12th International workshop on radiation imagin detectors, 2010. [39] Zhili Wang, Kun Gao, Peiping Zhu, Qingxi Yuan, Wanxia Huang, Kai Zhang, Youli Hong, Xin Ge, en Ziyu Wu. Grating-based x-ray phase contrast imaging using polychromatic laboratory sources. Journal of Electron Spectroscopy and Related Phenomena, In Press, Corrected Proof:–, 2010. ISSN 0368-2048. [40] C. David, J. Bruder, T. Rohbeck, C. Gruenzweig, C. Kottler, A. Diaz, O. Bunk, en F. Pfeiffer. Fabrication of diffraction gratings for hard X-ray phase contrast imaging. Microelectronic Engineering, 84(5-8):1172–1177, MAY-AUG 2007. ISSN 0167-9317. 32nd International Conference on Micro- and Nano-Engineering, Barcelona, SPAIN, SEP 17-20, 2006. [41] Andrei V. Bronnikov. Phase-contrast ct: Fundamental theorem and fast image reconstruction algorithms. In Proceedings of the Society of Photo-Optical Instrumentation Engineers (SPIE), 2006. [42] Bahaa E.A. Saleh en Malvin Carl Teich. Fundamentals of Photonics. John Wiley & Sons, Inc., 1991. [43] Joseph W. Goodman. Introduction to Fourier Optics. The McGraw-Hill Comanies, 1968. [44] B. Golosio, P. Delogu, I. Zanette, M. Carpinelli, G. L. Masala, P. Oliva, A. Stefanini, en S. Stumbo. Phase contrast imaging simulation and measurements using polychromatic sources with small source-object distances. Journal of Applied Physics, 104(9), 2008.

Bibliografie

81

[45] M. Dierick, D. Van Loo, B. Masschaele, J. Vlassenbroeck, M.N. Boone, L. Brabant, E. Pauwels, S. Peetermans, P. Vanderniepen, en L. Van Hoorebeke. Scanner design and construction at ugct. In First UGCT seminar, X-ray tomography as a multidisciplinary research tool: Exploring new frontiers, 2010. [46] Y. De Witte, J. Vlassenbroeck, en L. Van Hoorebeke. A multiresolution approach to iterative reconstruction algorithms in x-ray computed tomography. Image Processing, IEEE Transactions on, 19(9):2419 –2427, September 2010. ISSN 1057-7149. [47] Matthieu Boone. Fasecontrast beeldvorming aan een nanofocus x-stralenbuis. Master’s thesis, Universiteit Gent, 2006-2007. [48] A. Groso, R. Abela, en M. Stampanoni. Implementation of a fast method for high resolution phase contrast tomography. Optics Express, 14(18):8103–8110, 2006. [49] A. Pogany, D. Gao, en S.W. Wilkins. Contrast and resolution in imaging with a microfocus x-ray source. Review of Scientific Instruments, 68(7):2774–2782, JUL 1997. ISSN 0034-6748. [50] Wikipedia. Root mean square deviation, 2011. http://en.wikipedia.org/wiki/Root_mean_square_deviation. [51] S. C. Mayo, T. J. Davis, T. E. Gureyev, P. R. Miller, D. Paganin, A. Pogany, A. W. Stevenson, en S. W. Wilkins. X-ray phase-contrast microscopy and microtomography. Optics Express, 11:2289–2302, 2003. [52] B. D. Arhatari, W. P. Gates, N. Eshtiaghi, en A. G. Peele. Phase retrieval tomography in the presence of noise. Journal of Applied Physics, 107(3), FEB 2010. ISSN 00218979.

Lijst van figuren 1.1

Radiografie van het hart van een rat: attenuatie versus fasecontrast . . . . .

2.1 2.2 2.3 2.4 2.5 2.6

Schematische voorstelling van een CT-opstelling . . . Grafische illustratie van het Fourier Slice Theorema Rechten in het Fourier domein van de objectfunctie . Filterfunctie van het FBP algoritme . . . . . . . . . Beam hardening . . . . . . . . . . . . . . . . . . . . De Kaczmarz methode . . . . . . . . . . . . . . . . .

. . . . . .

. . . . . .

7 8 9 9 11 12

3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 3.10 3.11 3.12

Attenuatie en faseshift van een EM golf . . . . . . . . . . . . . . . . . . . Atomaire verstrooiingsfactoren, δ en β van koolstof . . . . . . . . . . . . . Refractie van een EM-golf na transmissie door een object . . . . . . . . . Edge enhancement aan de rand van een object . . . . . . . . . . . . . . . Fasesignaal van een polychromatische bundel . . . . . . . . . . . . . . . . Edge enhancement bij een epoxyhars staafje . . . . . . . . . . . . . . . . . LLL X-stralen interferometer . . . . . . . . . . . . . . . . . . . . . . . . . Opstelling bij diffraction enhanced imaging . . . . . . . . . . . . . . . . . Grating interferometer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . De verschillende gebieden voorbij het object bij vrije propagatie . . . . . . Verband tussen een lineair shift-invariant systeem en een optisch systeem Vrije propagatie als een shift-invariant systeem . . . . . . . . . . . . . . .

. . . . . . . . . . . .

15 16 17 18 18 19 21 22 23 24 25 26

4.1 4.2 4.3 4.4 4.5

Experimentele opstelling: 900 nm scanner . . . . . Experimentele opstelling: 400 nm scanner . . . . . Invloed van de spotgrootte op het beeld . . . . . . Secondary spot bij een epoxyhars staafje . . . . . . Invloed van de secondary spot op een reconstructie

. . . . .

30 31 32 33 33

5.1

Simulatie van een lijnprofiel van een homogene cilinder vóór convolutie met een Gauss . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

37

82

. . . . .

. . . . . .

. . . . .

. . . . . .

. . . . .

. . . . . .

. . . . .

. . . . . .

. . . . .

. . . . . .

. . . . .

. . . . . .

. . . . .

. . . . . .

. . . . .

. . . . . .

. . . . .

. . . . . .

. . . . .

. . . . . .

. . . . .

. . . . . .

. . . . .

. . . . .

2

Lijst van figuren 5.2

83

5.3 5.4 5.5 5.6 5.7 5.8 5.9

Simulatie van een lijnprofiel van een homogene cilinder na convolutie met een Gauss . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Simulatie van een lijnprofiel van een homogene cilinder na uitmiddelen . . Vergelijking van een gemeten lijnprofiel met een berekening . . . . . . . . . Vergelijking van een gemeten lijnprofiel met een berekening, detail linker rand Simulatie van een X-stralen spectrum en de detector respons . . . . . . . . Lijnprofiel van een cilinder bij een polychromatische bundel . . . . . . . . . Lijnprofiel van een cilinder bij een polychromatische bundel na uitmiddelen Tophoogtes bij verschillende SOD . . . . . . . . . . . . . . . . . . . . . . . .

38 38 39 39 41 41 42 43

6.1

Co¨ ordinatenstelsel gebruikt bij de afleiding van het Bronnikov algoritme . .

46

7.1

MBA versus MP: NRMSE van de centrale slice van een gesimuleerde homogene bol als functie van de attenuatie . . . . . . . . . . . . . . . . . . . . . . MBA versus MP: lijnprofiel van een gesimuleerde homogene bol bij verschillende attenuaties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . MBA versus MP: gemeten lijnprofiel van epoxyhars bij verschillende buisspanningen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Fasefantoom inhomogeen faseobject . . . . . . . . . . . . . . . . . . . . . . Invloed van de homogeniteit op MBA . . . . . . . . . . . . . . . . . . . . . Invloed van de homogeniteit op MP . . . . . . . . . . . . . . . . . . . . . . Invloed van de verhouding µ/δ op MBA . . . . . . . . . . . . . . . . . . . . Homogene bol met 1 fase-inclusie gereconstrueerd met MBA . . . . . . . . . Invloed van een kleinere spreiding in de verhouding µ/δ op MBA . . . . . . Invloed van een kleinere spreiding in de verhouding µ/δ op MP . . . . . . . Hars met grafiet staafje . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Lijndoorsnede van een inhomogeen object: grafiet in epoxyhars . . . . . . . De Bronnikov filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . NRMSE van een gereconstrueerde slice als functie van het gemiddeld aantal tellen op een detectorpixel . . . . . . . . . . . . . . . . . . . . . . . . . . . . Effect van MBA en MP onder ruis . . . . . . . . . . . . . . . . . . . . . . . FBP, MBA en MP toegepast op een pilletje . . . . . . . . . . . . . . . . . . MBA en MP toegepast op een diamant . . . . . . . . . . . . . . . . . . . . MBA en MP toegepast op een stukje dennenhout met houtschimmel . . . . Reconstructie van een vliegenpoot met FBP en MBA . . . . . . . . . . . . . MBA en MP toegepast op een vliegenpoot . . . . . . . . . . . . . . . . . . .

7.2 7.3 7.4 7.5 7.6 7.7 7.8 7.9 7.10 7.11 7.12 7.13 7.14 7.15 7.16 7.17 7.18 7.19 7.20

A.1 De co¨ ordinaten van een ellips . . . . . . . . . . . . . . . . . . . . . . . . . . A.2 Projectie en gereconstrueerde centrale doorsnede van het gebruikte fantoom

53 54 55 56 56 57 58 60 62 63 63 64 66 67 67 69 70 71 72 72 75 76

Lijst van tabellen 5.1 5.2 5.3 5.4 7.1 7.2

Parameters voor de simulatie van de projectie van een homogene cilinder . Gebruikte energieën en δ’s in de polychromatische bundel . . . . . . . . . . Parameters voor de simulatie van een homogene cilinder bij een polychromatische bundel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Parameters voor de simulatie van een homogene cilinder bij verschillende SOD

36 40 40 43

Parameters voor de simulatie van een homogene bol . . . . . . . . . . . . . Parameters van het gebruikte fantoom bij de simulatie van een inhomogeen, niet-attenuerend faseobject . . . . . . . . . . . . . . . . . . . . . . . . . . . Parameters van het fantoom bij de simulatie van een inhomogeen object waarbij µ/δ constant is . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Parameters van het fantoom bij de simulatie van een inhomogeen object waarbij de verhouding µ/δ niet constant is . . . . . . . . . . . . . . . . . . . Parameters van het fantoom bij de simulatie van een homogeen attenuerende bol met 1 fase-inclusie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Parameters voor de scan van het pilletje . . . . . . . . . . . . . . . . . . . . Parameters voor de scan van de diamant . . . . . . . . . . . . . . . . . . . . Parameters voor de scan van het stukje dennenhout . . . . . . . . . . . . . Parameters voor de scan van de vliegenpoot . . . . . . . . . . . . . . . . . .

59 68 69 71 72

A.1 Fantoombeschrijving . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

75

7.3 7.4 7.5 7.6 7.7 7.8 7.9

84

52 55 57 57

Studie en evaluatie van fasecontrast tomografie aan een hoge resolutie X-stralenbuis

Recommend Documents