Segmentatie van transmissie tomografie-data voor kwantitatieve analyse op basis van een discreet grijswaardenmodel

Segmentatie van transmissie tomografie-data voor kwantitatieve analyse op basis van een discreet grijswaardenmodel. Wim van Aarle IWT projectvoorstel 1e termijn

Promotor: Prof. Dr. Jan Sijbers

Inhoudsopgave 1 Probleemstelling

2

2 Doelstelling

3

3 Projectbeschrijving 3.1 Situering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Inleiding en Literatuurstudie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 X-stralen Attenuatie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Reconstructie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Segmentatie van tomografische reconstructies . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Fase 1: Automatische thresholding van cone-beam CT data op basis van gemeten projectie-data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.2 Fase 2: Nieuwe automatische segmentatietechnieken op basis van gemeten projectiedata. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.3 Fase 3: Combinatie van reconstructie en segmentatie d.m.v. discrete tomografie. 3.3.4 Fase 4: Discrete tomografie op sterk absorberende objecten. . . . . . . . . . 3.4 Conclusie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4 4 4 4 5 7

10 11 12 13

4 Planning

14

5 Toepassingsmogelijkheden

15

Referenties

16 1

9

Segmentatie van transmissie tomografie-data voor kwantitatieve analyse

1

Probleemstelling

Transmissietomografie is een techniek waarbij een driedimensionaal voorwerp virtueel gereconstrueerd wordt aan de hand van tweedimensionale projectiebeelden. Computed Tomography (CT), de meest voorkomende vorm van tomografie, gebruikt X-stralen om deze projectiedata te meten. Een wiskundig algoritme reconstrueert nadien het originele object in een driedimensionaal volume, het tomogram. Daarin wordt aan elke voxel een grijswaarde toegekend die afhangt van de gemeten eigenschap van het materiaal in het originele voorwerp (zoals bv. de densiteit). Via tomografie kan de interne structuur van een voorwerp gevisualiseerd worden. Daarnaast is echter ook extractie van kwantitatieve gegevens interessant. Dit omvat afmetingen, volume en morfologische aspecten zoals het aantal kenmerken of de connectiviteit tussen bepaalde delen. De extractie van dergelijke informatie uit een tomogram met een continu grijswaardenspectrum is echter niet triviaal. Daarom is het noodzakelijk dat een accurate en betrouwbare segmentatie kan worden toegepast die een discreet grijswaardenmodel oplegt aan het tomogram. Dit ogenschijnlijk eenvoudig probleem wordt gekenmerkt door een aantal moeilijkheden die stuk voor stuk aangepakt moeten worden: • Vaak wordt voor segmentatie van tomografische reconstructies beroep gedaan op globale thresholding. Deze methode is echter zeer subjectief. De keuze van de drempelwaarden is afhankelijk van de gebruiker. Verschillende scanners geven soms verschillende metingen en reconstructies zijn doorgaans onderhevig aan ruis en artefacten. Uitsluitend op basis van het gereconstrueerde beeld (zoals bij segmentatie op basis van het histogram) is het daarom niet duidelijk hoe de drempelwaarden optimaal bepaald kunnen worden. Daarom is er een objectief kwaliteitscriterium nodig dat efficiënt kan worden geëvalueerd. • Meer geavanceerde segmentatietechnieken (zoals lokale thresholding of edge detection) zijn mogelijk maar kampen met dezelfde problemen [21]. Ook hier wordt momenteel enkel gebruik gemaakt van de beschikbare reconstructie en wordt er geen objectieve kwaliteitsmaatstaf toegekend aan de segmentatie. Bovendien zijn deze technieken afhankelijk van meerdere parameters, wat een snelle optimalisatie zeer moeilijk maakt. • Doorgaans is een groot aantal projecties vereist om tot een kwalitatief goede reconstructie, en dus ook tot een kwalitatief goede segmentatie, te komen. Dit staat in de praktijk synoniem voor een grote stralingsbelasting, wat een veelvuldige longitudinale opvolging van het gescande voorwerp onmogelijk maakt. Het gebruik van belangrijke voorkennis (zoals bv. materiaaleigenschappen) kan het aantal benodigde projecties drastisch terugschroeven zonder in te moeten boeten aan nauwkeurigheid. Deze voorkennis is vaak aanwezig maar wordt momenteel niet of nauwelijks benut. • Reconstructies van objecten met sterk absorberende materialen (zoals bv. protheses of tandvullingen) vertonen vaak ernstige artefacten. Verscheidene Metal Artifact Reduction technieken (MAR) werden reeds voorgesteld om het effect van deze artefacten te verminderen [4,23,24]. Geen van hen zijn er echter in geslaagd om de artefacten volledig te verdringen. Dit maakt een betrouwbare, automatische segmentatie bijzonder lastig, zo niet onmogelijk.

IWT-projectvoorstel

Wim van Aarle

Pagina 2 van 17


2

Doelstelling

In de onderzoekseenheid Visielab van de Universiteit Antwerpen is er recent een nieuwe deelgroep ontstaan met als doel het verbreden van expertise in beeldverwerkingalgoritmen voor transmissie tomografie. Dit project spitst zich daarbij toe op het ontwikkelen van nieuwe beeldsegmentatietechnieken specifiek voor tomografische reconstructies. Hierdoor zal een kwantitatieve analyse mogelijk worden. In eerste instantie zullen nieuwe segmentatiealgoritmes worden ontwikkeld die specifiek gebruik maken van beschikbare projectiedata. In tweede instantie zal worden gewerkt aan de ontwikkeling van algoritmen die reconstructie en segmentatie combineren. Tijdens het ontwikkelen zal er nadruk gelegd worden op de kwaliteit van de segmentaties en op de computationele efficiëntie van de nieuwe algoritmen. • Vooreerst moeten er nieuwe segmentatietechnieken komen voor globale thresholding op tomografische reconstructies. Daarbij zal een objectieve vergelijkingsmaatstaf worden ontwikkeld op basis van de gemeten projectiedata. Door deze te vergelijken met een voorwaartse projectie van een segmentatie kan de kwaliteit van de segmentatie worden berekend. Dit geeft aanleiding tot een optimalisatieprocedure met de meest nauwkeurige segmentatie als doel. Omdat daarvoor een groot aantal segmentaties moeten worden beoordeeld, is het essentieel dat de kwaliteitsmaatstaf zo efficiënt mogelijk te evalueren is. • Nadat een objectief kwaliteitscriterium is toegepast voor globale thresholing, moet onderzocht worden in welke mate het mogelijk is om dezelfde maatstaf ook te gebruiken bij meer geavanceerde segmentatietechnieken. Speciale aandacht zal daarbij nog steeds gaan naar de computationele efficiëntie. Een vergelijkende studie moet aantonen wat de vooren nadelen zijn van de nieuw ontwikkelde segmentatiemethodes in vergelijking met reeds bestaande technieken. • Het aantal projecties dat nodig is om een segmentatie van een reconstructie te maken, moet worden terug teruggebracht. Discrete Tomografie (DT) kan hiervoor een oplossing bieden. In plaats van een klassiek tweestaps proces, waarbij eerst een continue reconstructie berekend wordt en nadien pas wordt gesegmenteerd, wordt bij DT de segmentatie en de reconstructie gecombineerd in één algoritme. We voeren daartoe één extra aanname in: gescande objecten bestaan uit een zeer beperkt aantal materialen of weefsels. Dit moet het mogelijk maken om met substantieel minder projecties een minstens even nauwkeurige segmentatie te bekomen als bij klassieke, continue tomografische technieken. Recent werden enkele nieuwe, veelbelovende DT-algoritmen voorgesteld [5, 7]. Deze kunnen echter nog niet direct toegepast worden op CT-data. We zullen deze daarom eerst doorontwikkelen voor gebruik met experimentele data. De kwaliteit van de segmentaties zal nadien onderzocht worden. • Segmentatieproblemen bij continue reconstructies van voorwerpen die sterk absorberende materialen bevatten, moeten worden voorkomen. De nieuwe discrete reconstructiealgoritmen worden daarom gecombineerd met de klassieke, continue tomografietechnieken. Via voorkennis kan discrete tomografie zorgen voor een nauwkeurige reconstructie en segmentatie van dichte materialen. Tegelijkertijd kan een klassiek algoritme de rest van het volume continu reconstrueren met het verschil van de gemeten projectiedata en de voorwaartse projectie van de discrete reconstructie.

IWT-projectvoorstel

Wim van Aarle

Pagina 3 van 17


3 3.1

Projectbeschrijving Situering

De Universiteit Antwerpen heeft een sterke onderzoeksgeschiedenis in biomedische- en materiaalwetenschappen. Deze maken typisch veelvuldig gebruik van beeldverwerking- en beeldanalysetechnieken. Het in 1991 opgerichte Visielab zorgt voor een multidisciplinaire groepering van dit onderzoeksdomein. Met de recente aanwerving van Joost Batenburg (in 2006 verkozen tot on” derzoeker van het jaar” door de Faculteit Wiskunde en Natuurwetenschappen van de Universiteit leiden) en enkele nieuwe doctorandi is er een nieuwe deelgroep ontstaan met als doel het verbreden van expertise in het construeren en verbeteren van beeldverwerkingalgoritmen voor transmissie tomografie. Het Visielab bevindt zich in een zeer gunstige positie om dit onderzoek te ondersteunen. Zo is er uitstekend contact met het UZA en met tal van partners in de Vlaamse industrie. Via deze samenwerking kunnen nieuwe technieken en ontwikkelingen getoetst worden met data uit medische of industriële toepassingen. Eén partner verdient extra vermelding. SkyScan is een spin-off bedrijf van het Visielab. Het commercialiseert thans X-stralen microtomografen en is daarbij wereldwijd marktleider. De ontwikkeling, constructie en optimalisatie van microtomografie is één van de kernactiviteiten van het Visielab.

3.2

Inleiding en Literatuurstudie

Transmissietomografie is een techniek waarbij een driedimensionaal voorwerp virtueel gereconstrueerd wordt aan de hand van tweedimensionale projectiebeelden. Daarbij wordt het voorwerp onder verschillende hoeken bestraald. Nadien reconstrueert een wiskundig algoritme het originele object in een driedimensionaal volume, het tomogram. Daarin wordt aan elke voxel een grijswaarde toegekend die afhangt van de gemeten eigenschap van het materiaal in het originele voorwerp (zoals bv. de densiteit). In de volgende sectie wordt de acquisitie bij Computed Tomography (CT), de meest voorkomende vorm van transmissietomografie, verduidelijkt [20]. 3.2.1

X-stralen Attenuatie

Een bundel fotonen (X-stralen) van zekere energie E wordt afgevuurd vanuit een stralenbron in de richting van een fotonendetector, door een object (figuur 1a). Als een foton binnen het object daarbij in contact komt met een elektron zal het de detector waarschijnlijk niet bereiken. Men spreekt dan van attenuatie. Bij energieniveaus die gebruikt worden in CT-scanners (typisch tussen 20 en 150 keV) zijn hiervoor twee hoofdoorzaken: bij lage energie de fotoelektrische absorptie en bij gematigde energie het Compton effect. Bij het eerste effect absorbeert een sterk gebonden elektron alle energie van de invallende X-straal. Een deel daarvan wordt gebruikt om zijn binding met de kern te verbreken, de rest wordt kinetische energie. Bij het Compton effect wordt een X-straal gerefracteerd op een vrij of zwak gekoppeld elektron. De X-straal gaat dan, met verlies van energie, verder in een andere richting. De kans dat één van deze effecten optreedt is afhankelijk van de dikte en de aard van het materiaal. Zoals blijkt uit figuur 1a, is elke straal afhankelijk van twee parameters: de projectiehoek θ en de positie van de stralenbron en detector op de rechte t. De verzameling projecties kan dus

IWT-projectvoorstel

Wim van Aarle

Pagina 4 van 17

Segmentatie van transmissie tomografie-data voor kwantitatieve analyse y

s

t

de t e c t or

f ( x , y) x

θ

s t r a l e nbr on θ

( a )

t

( b)

Figuur 1: (a) Opstelling parallel-beam CT-scan (b) Sinogram van het fantoomobject uit (a) geschreven worden als een functie Pθ (t), de zogenaamde sinogramfunctie of de Radon-transformatie (figuur 1b). Eenvoudigheidshalve wordt er daarbij gewerkt in het poolcoördinatenstelsel (t, s). Z +∞ Pθ (t) = f (t, s)ds (1) −∞

In figuur 1 zijn alle X-stralen van een projectiebeeld parallel aan elkaar. We spreken dan over parallel-beam CT. Deze methode is heel eenvoudig, maar omdat het in praktijk zeer lastig is om een brede stralenbundel op te wekken moet de stralenbron en de detector voor elke meting opschuiven in de richting van t. Er kan dus maar één detectormeting tegelijkertijd plaatsvinden. Bij fan-beam CT is de scanner uitgerust met een ééndimensionale array van detectoren. De X-stralen zijn niet meer parallel, maar liggen in de vorm van een waaier. Alle metingen die nodig zijn om één volledige projectie van een slice te maken, kunnen zo simultaan uitgevoerd worden. De acquisitietijd wordt zo sterk gereduceerd. Cone-beam CT gaat nog een stapje verder en breidt de scanner uit met een tweedimensionale array van detectoren. De stralenbundel heeft dan de vorm van een kegel. Zo worden de X-stralen nog efficiënter gebruikt. Vermits er in cone-beam CT gewerkt wordt met tweedimensionale projecties en er per opname dus meerdere sneden van het voorwerp worden bestraald, is het niet mogelijk om de reconstructie snede per snede af te werken. Dit heeft als nadeel dat de ruimtecomplexiteit van het reconstructiealgoritme sterk toeneemt [17, 18]. 3.2.2

Reconstructie

Is de Radon-transformatie gekend, dan kan een reconstructiealgoritme het oorspronkelijk voorwerp volumetrisch nabouwen. Hiervoor zijn verschillende methodes beschikbaar, elk met hun eigen vooren nadelen: Direct Fourier Reconstruction Het Fourier Slice Theorema zegt dat een Fourier getransformeerde Sθ (w) van een projectie Pθ (t) gelijk is aan een doorsnede in diens Fourierruimte onder IWT-projectvoorstel

Wim van Aarle

Pagina 5 van 17


dezelfde hoek θ en door de oorsprong: Z +∞ Z +∞ f (t, s) ds e−j2θwt dt Sθ (w) = −∞ −∞ Z +∞ Z +∞ f (x, y)e−j2θw(x cos θ + y sin θ) dxdy = −∞

(2) (3)

−∞

= F(w cos θ, w sin θ)

(4)

Deze getransformeerde projecties kunnen dus gezien worden als monsters waarmee men de Fourierruimte kan benaderen. Nadien kan men via een inverse Fouriertransformatie tot het originele beeld komen. Daar men in de praktijk beperkt is in het aantal projecties, is er nood aan interpolatie in het frequentiedomein. Omdat elke projectie in de Fourierruimte een lijn door de oorspong is, is de interpolatiefout niet overal gelijk en zijn de lage frequenties van het beeld veel beter bemonsterd dan de hoge frequenties. Dit maakt de reconstructie onscherp. Filtered Backprojection In wezen slechts een wiskundige herformulering van de Direct Fourier Approach, is Filtered Backprojection (FBP) een eenvoudig en efficiënt evalueerbaar algoritme om een beeld te reconstrueren uit diens Radon-ruimte. In [10] tonen Kak en Slaney aan dat f (x, y) in poolcoördinaten (w, θ) geschreven kan worden als: Z π Z +∞ f (x, y) = Sθ (w) |w| ej2π(x cos θ+y sin θ) dwdθ (5) 0

−∞

In deze formule zijn de stappen van het FBP-proces duidelijk zichtbaar: 1. De Radon-ruimte wordt gefilterd met een rampfilter |w|. 2. Het beeld wordt gevormd via een sommatie van alle gefilterde terugprojecties van de Radonruimte (w, θ). Dit zijn beide eenvoudige operaties en dus is FBP een aantrekkelijke en populaire keuze voor tomografische reconstructie. Algebra¨ısche Technieken De Algebraic Reconstruction Technique (ART) werd ontworpen door Gordon, Bender en Herman in [3] en biedt enkele belangrijke voordelen ten opzichte van Filtered Backprojection. Conceptueel is ART veel eenvoudiger, wat uitbreidingen zoals het verwerken van a priori informatie toelaat. Ook zijn er doorgaans minder projecties nodig om een evenwaardige kwaliteit te bekomen, van belang indien de X-stralen schadelijk zijn voor het gescande voorwerp [18]. Het is ook mogelijk om, zonder diepgaande veranderingen in het algoritme aan te brengen, ondersteuning te bieden voor niet-equidistante of niet-equiangulaire projecties. De reconstructie wordt geschreven als een probleem uit de lineaire algebra: Wx = p

(6)

Hierbij is p ∈ RM de gemeten projectiedata (M = aantal detectoren × aantal projecties), x ∈ RN is de reconstructie (N = aantal pixels in het gereconstrueerde beeld) en W ∈ RM ×N een gewichtsmatrix waarbij elke wij weergeeft in welke mate voxel vj in het oorspronkelijke volume een invloed IWT-projectvoorstel

Wim van Aarle

Pagina 6 van 17


( a )

( b)

( c )

Figuur 2: SART-Reconstructies van het fantoomobject uit figuur 1 na 10 iteraties. (a) 4 projecties (b) 20 projecties (c) 90 projecties. heeft op de projectiewaarde pi . Gezien de enorme proporties die de matrix W kan aannemen is het oplossen van een dergelijk stelsel niet triviaal. Technieken die gebruik maken van een matrix inversie zoals LU-factorisatie zijn niet mogelijk. Er wordt daarom gekozen voor benaderende technieken zoals de Kaczmarz methode [9]. Daarbij wordt elke rij wi1 x1 + ... + wiN xN = pi in het stelsel beschouwd als een hypervlak in een N -dimensionale ruimte. De oplossing van het stelsel (8) ligt dan net op het snijpunt van deze M vlakken. Dit snijpunt (of een benadering daarvan) wordt bepaald door, startend vanuit een zekere beginoplossing, iteratief de huidige oplossing loodrecht te projecteren op een willekeurig nieuw hypervlak. Elk stapje van het reconstructieproces kan wiskundig uitgedrukt worden als volgt [18]: P win vnk pi − N 0 wij (7) xj = xj + λ PNn=1 2 n=1 win Bij standaard ART zijn de gewichten wij binair. Zo kunnen ze tijdens de evaluatie worden berekend. Dit is gunstig voor de ruimtecomplexiteit van het algoritme. In de Simultaneous Iterative Reconstruction Technique (SIRT) worden tussentijdse oplossingen pas bijgewerkt nadat elk hypervlak éénmaal werd toegepast. Zo wordt convergentiesnelheid opgeofferd voor een superieure kwaliteit. Simultaneous Algebraic Reconstruction Technique (SART) zoekt een compromis tussen ART en SIRT door de correcties op het beeld projectie per projectie uit te voeren. Verder wordt er gebruik gemaakt van een andere interpolatiekernel, een bilineaire functie, die ruis op het beeld drastisch doet dalen [10]. Figuur 2 toont aan dat de kwaliteit van continue reconstructies daalt indien het aantal projecties afneemt. Dit is een groot nadeel daar het aantal projecties direct proportioneel is met de scantijd en de stralingsbelasting.

3.3

Segmentatie van tomografische reconstructies

Tomografische reconstructies zijn vaak onderhevig aan ruis en onnauwkeurigheden, bijvoorbeeld ten gevolge van een beweging van het object, een slechte signal-to-noise-ratio van de scanner of simpelweg een tekort aan projecties. Dit maakt dat een beeld dat slechts enkele grijswaarden IWT-projectvoorstel

Wim van Aarle

Pagina 7 van 17


( a )

( b)

( c )

Figuur 3: Segmentaties van de reconstructies van figuur 2. (a) 4 projecties (b) 20 projecties (c) 90 projecties. zou mogen bevatten, alsnog uit tientallen grijswaarden kan zijn opgebouwd. Om kwantitatieve informatie uit een dergelijk beeld te halen, is daarom eerst nood aan een nauwkeurige segmentatie. Een veelgebruikte segmentatietechniek is thresholding. De classificatie van pixels is daarbij afhankelijk van hun intensiteit ten opzichte van een verzameling drempelwaarden t1 , t2 , ..., tn−1 . Elke partitie krijgt een eigen, unieke grijswaarde g1 , g2 , ..., gn toegekend:   g1 f (x, y) ≤ t1    g2 t1 < f (x, y) ≤ t2 (8) fS (x, y) = ..  .    g t < f (x, y) n

n−1

De kwaliteit van de segmentatie staat of valt met de keuze van ti en gi . Het manueel aanpassen van de drempelwaarden op basis van visuele inspectie is mogelijk maar moet vermeden worden. Om de keuze te automatiseren wordt vaak gebruik gemaakt van een histogram, een diagram van de grijswaardendistributie [8,19]. Elke grijswaarde wordt dan toegekend aan een piek in het histogram (lokaal maximum) en de drempelwaarden worden geplaatst in de vallei tussen twee dergelijke pieken (lokaal minimum). Deze toekenning kan gebeuren via een regressie van gaussische functies op het histogram. Dit resulteert echter doorgaans in een kwalitatief matige segmentatie: 1. Bij klassieke tomografische beeldvorming bevat het histogram zelden perfecte gaussische distributies. 2. De densiteit van twee verschillende materialen kan zo dicht bij elkaar liggen dat hun grijswaardendistributies met elkaar overlappen. In vier fasen wordt er getracht het segmentatieprobleem te overwinnen. In eerste instantie beperken we ons daarbij tot globale thresholding, maar gaandeweg wordt dit uitgebreid naar andere segmentatieparadigma’s. Uiteindelijk kijken we ook naar een nieuw type reconstructieproces waarbij segmentatie reeds een eigenschap is van het resultaat. Dit alles moet de kwantitatieve analyse van een object mogelijk maken.

IWT-projectvoorstel

Wim van Aarle

Pagina 8 van 17


3.3.1

Fase 1: Automatische thresholding van cone-beam CT data op basis van gemeten projectie-data.

Een segmentatie van een tomografische reconstructie kan, net zoals het oorspronkelijke voorwerp, onderworpen worden aan een voorwaartse projectie. Het verschil met de gemeten projectiedata uit de scanner is dan een goede maatstaf voor de kwaliteit van de segmentatie. Via een optimalisatieproces kunnen de optimale parameters van een segmentatie gevonden worden. Essentieel daarbij is dat daarvoor zo weinig mogelijk evaluaties worden uitgevoerd, het voorwaarts projecteren is namelijk een computationeel zeer onaantrekkelijke operatie. In de sectie over algebra¨ısche reconstructietechnieken werd reeds vermeld dat een tomografische projectie gemodelleerd kan worden als een lineair stelsel Wx = p

(9)

met x het gereconstrueerde beeld, p de projectiedata en W een lineaire projectieoperator. Elke pixel xi wordt slechts geprojecteerd op een klein aantal detectorpixels pi . Bijgevolg is W een ijle matrix. Dit is essentieel voor het voorgestelde algoritme en ligt aan de basis van diens efficiënte evalueerbaarheid. Stel dat het gescande object uit precies l verschillende materialen bestaat. Een beeld moet dan opgedeeld worden in een verzameling partities S = {S1 , ..., Sl } die bepaald worden door de drempelwaarden t1 , ..., tl−1 . Een segmentatiefunctie rS (g) wijst elke pixel xi toe aan een partitie St en aan diens bijhorende grijswaarde gt ∈ R. Voor de optimale segmentatie geldt dat de volgende functie minimaal is: |W rS (g) − p|2 (10) Stel eerst dat de drempelwaarden ti , en dus ook de partities Si , onveranderlijk zijn. Voor elke m×l partitie apart P moet men dan op zoek naar de meeste correcte grijswaarde. Zij A = (ait ) ∈ R met ait = j:s(j)=t wij , het gewicht van de pixels in partitie St op de detectorwaarde pi , dan geldt: [W rS (g)]i =

l X

ait gt = ai g

(11)

t=1

Het totale projectieverschil d ∈ Rm is dan |d|2 = |Ag − p|2 =

X

(aTi g − pi )2

(12)

=

X

(g T ai aTi g − 2pi ai g + p2i )

(13)

= pT Qg + cT g + |p|2

(14)

P P met c = −2pi aTi en Q = ai aTi . De kern van het algoritme rust op het feit dat deze twee termen niet afhankelijk zijn van g en dus bij een optimalisatieproces slechts éénmaal berekend hoeven te worden. Stel nu dat de grijswaarden onveranderlijk zijn en de partitie S in enkele pixels gewijzigd wordt. Het herbereken van c en Q blijkt dan noodzakelijk. Indien er slechts enkele pixels van partitie veranderen, zijn er heel weinig rijen ai die aangepast moeten worden. Dit maakt dat de volgende twee formules ervoor zorgen dat de optimale grijswaarden voor een nieuwe, licht aangepaste segmentatie,

IWT-projectvoorstel

Wim van Aarle

Pagina 9 van 17


zeer efficiënt kunnen herberekend worden: X

c0 = c +

(−2p0i a0i ) − (−2pi ai )

(15)

i:wi,j 6=0

Q0 = Q +

X

(a0i aTi 0 ) − (ai aTi )

(16)

i:wi,j 6=0

In dit onderzoek zullen algoritmen worden ontwikkeld die gebruik maken van dit kwaliteitscriterium. Het bepalen van een optimale segmentatie zal bestaan uit een iteratief optimalisatieproces van de parameters ti en gi . Met kleine stapjes kunnen de drempelwaarden worden verlaagd of verhoogd en kunnen pixels, indien nodig, herverdeeld worden in andere partities. Elke aangepaste partitie krijgt dan vervolgens een nieuwe grijswaarde toegewezen aan de hand van formule (16). Is er slechts één drempelwaarde, dan is het mogelijk om het grijswaardenspectrum stapsgewijs af te gaan en de optimale drempelwaarde te bewaren. Zijn er meerdere drempelwaarden, dan wordt een segmentatieruimte te groot om stapje per stapje te overlopen. Vanuit een startsituatie moet dan door het aanpassen van drempelwaarde per drempelwaarde het verschil d zo klein mogelijk gemaakt worden. Om lokale minima te vermijden, is het belangrijk dat de startsituatie kwalitatief goed is. Geavanceerde optimalisatietechnieken zoals simulated annealing bieden geen soelaas, daar zij het proces trachten te versnellen door, daar waar het kan, grotere stappen in de optimalisatieruimte te nemen. In dit algoritme staan grotere stappen echter gelijk met langzamere herberekening. In de eerste fase van het project wordt de kwaliteitsmaatstaf verder geoptimaliseerd en wordt ervoor gezorgd dat er naast parallel-beam ook ondersteuning is voor fan-beam en cone-beam CT. Planning Automatische globale thresholding op basis van gemeten projectiedata. Subtaak 1: Literatuurstudie. (3 maanden) Subtaak 2: Ontwerp van globale thresholding voor parallel-beam CT met objectieve kwaliteitsmaatstaf. (3 maanden) Subtaak 3: Ontwerp van globale thresholding voor cone-beam CT met objectieve kwaliteitsmaatstaf. (6 maanden) Subtaak 4: Validatie met experimentele data of simulaties. (2 maanden) 3.3.2

Fase 2: Nieuwe automatische segmentatietechnieken op basis van gemeten projectiedata.

Globale thresholding is niet altijd even geschikt om een nauwkeurige segmentatie te bekomen. Zo is het bv. moeilijk om zeer kleine, niet dichte voorwerpen van elkaar te onderscheiden. Het toevoegen van de kwaliteitscontrole uit de eerste fase zal hiervoor geen oplossing bieden want deze voorwerpen hebben slechts een kleine impact op de projectiedata. In deze fase van het project wordt het objectief kwaliteitscriterium daarom uitgebreid naar enkele meer geavanceerde segmentatietechnieken. Enkele voorbeelden: Lokale Adaptieve Thresholding In tegenstelling tot globale thresholding wordt bij lokale thresholding een drempelwaarde ti in elke pixel afzonderlijk gekozen. Bij lokale variantie thresholding hangt deze af van de mediaan en de standaard afwijking van de grijswaarden in een venster rondom de pixel, bij lokale contrast thesholding van het gemiddelde van de grijswaarden [21]. De omvang van het venster rondom de pixel is bepalend voor de kwaliteit van de segmentatie. Via voorwaartse projecties van een tussentijdse segmentatie kunnen de pixels IWT-projectvoorstel

Wim van Aarle

Pagina 10 van 17


die niet tot hun optimale partitie behoren, gelokaliseerd worden. De drempelwaarden in deze pixels kunnen in een volgende iteratie opnieuw bepaald worden. Omdat een pixel slechts invloed heeft op enkele detectorwaarden, kan de volgende voorwaartse projectie zeer efficiënt verlopen. Randdetectie Bij randdetectie gaat men er vanuit dat de precieze vorm van een partitie afhangt van zijn randen [25]. Deze kunnen worden bepaald via randdetectoren zoals de Sobel operator, de Laplaciaanse edge detector of de Marr-Hildreth edge operator. Het resultaat daarvan is een grijswaardenbeeld waarbij hoge intensiteit een sterke rand voorstelt, waarbij lage intensiteit een zwakke rand voorstelt en waarbij geen intensiteit geen rand voorstelt. Via adaptive thresholding kan daaruit een binair beeld worden opgesteld waarbij de verschillende partities duidelijk zichtbaar zijn. Dit is echter niet altijd even nauwkeurig. Bij aanwezigheid van ruis of onnauwkeurigheden in het tomogram kunnen randen verkeerd geplaatst worden. Door de resulterende segmentatie te vergelijken met de projectiedata, is het mogelijk om fouten in de randdetectie iteratief te corrigeren. Nadat is onderzocht hoe de nieuwe, objectieve kwaliteitsmaatstaf past in deze segmentatietechnieken, zal een vergelijkende studie een uitgebreid verslag uitbrengen van de vooren nadelen van de nieuwe segmentatietechnieken in vergelijking met de klassieke, niet-objectieve methoden. Planning Subtaak 1: Subtaak 2:

Nieuwe segmentatietechnieken op basis van gemeten projectiedata. Literatuurstudie. (3 maanden) Ontwerp van nieuwe segmentatietechnieken met objectieve kwaliteitsmaatstaf. (5 maanden) Subtaak 3: Validatie van de nieuwe technieken met experimentele data of simulaties. (2 maanden) Subtaak 4: Een vergelijkende studie die de nieuw ontwikkelde segmentatietechnieken vergelijkt met niet objectieve methoden. (4 maanden) 3.3.3

Fase 3: Combinatie van reconstructie en segmentatie d.m.v. discrete tomografie.

Een heel andere aanpak om een gesegmenteerd beeld te bekomen uit tomografische projectiedata is om de segmentatie reeds toe te passen tijdens de reconstructie. We spreken dan over Discrete Tomografie (DT). DT is fundamenteel verschillend van de klassieke continue tomografische reconstructietechnieken. Naast calculus en lineaire algebra wordt er bij DT ook gebruik gemaakt van combinatoriek, artificiële intelligentie en optimalisatietechnieken [6]. Men gaat er bij DT vanuit dat het te reconstrueren object slechts uit een beperkt aantal materialen of weefsels bestaat. A priori informatie over de materialen van het gescande voorwerp kan gebruikt worden om de kwaliteit van de reconstructie te verbeteren en om het aantal projecties dat daarvoor nodig is, te verlagen. Als er weinig projecties zijn, is het DT-probleem sterk onderbepaald en kunnen er verschillende reconstructies zijn die voldoen aan de projectiedata. Daarom zijn extra regels nodig die een voorkeur geven aan de meest intu¨ıtieve van de verschillende correcte oplossingen [14]. Een heuristische DT-methode is de Discrete Algebraic Reconstruction Technique (DART) [2]. Deze is gebaseerd op het continue SIRT-algoritme. Na een aantal iteraties van SIRT wordt de tussentijdse reconstructie gesegmenteerd volgens thresholdingsregels opgesteld aan de hand van IWT-projectvoorstel

Wim van Aarle

Pagina 11 van 17


a priori informatie. Segmentatiefouten in deze reconstructies komen voornamelijk voor in het grensgebied van twee partities. Pixels die niet in een dergelijke grensregio liggen, zijn waarschijnlijk correct en worden daarom tijdelijk vast gehouden. De volgende SIRT-iteratie werkt dan uitsluitend op pixels die toebehoren aan een grensgebied. Is een pixel in het centrum van een regio toch foutief, dan is de kans erg groot dat deze na enkele iteraties alsnog in een grensgebied zal belanden. Voor elke pixel i die vast wordt gezet, moet formule (8) aangepast worden:   x1    ...   w1,1 . . . w1,i−1 w1,i+1 . . . w1,n    x  .. .. .. ..   i−1  = p − v w (17)   . i i . . .  x  i+1   wm,1 . . . wm,i−1 wm,i+1 . . . wm,n   ...  xn Deze aanpak verlaagt de tijdsduur van de SIRT-iteraties en zorgt ervoor dat de grijswaarden sterk gegroepeerd blijven. Hoewel er in de laatste decennia reeds veel onderzoek werd verricht naar oplossingen voor het DT-probleem [1, 5, 7], staat het alsnog in zijn kinderschoenen. Zo is DT nog nooit toegepast op cone-beam scanner data. Dit is dan ook het hoofddoel van de derde fase van het project. Planning Subtaak 1: Subtaak 2: Subtaak 3: 3.3.4

Combinatie van reconstructie en segmentatie d.m.v. discrete tomografie. Literatuurstudie. (3 maanden) Ontwerp van het DART algoritme voor cone-beam CT data. (4 maanden) Validatie met experimentele data of simulaties. (2 maanden)

Fase 4: Discrete tomografie op sterk absorberende objecten.

Een sterk absorberend object, zoals een metalen voorwerp, heeft door zijn hoger atoomgetal een attenuatiefactor die ettelijke malen groter is dan die van een zacht weefsel. Bij X-stralen tomografie zullen er onvoldoende fotonen de detectoren bereiken en zal de projectiedata “onvolledig” zijn. Als er gebruik gemaakt wordt van FBP of ART zorgt dit voor streaking artefacten die de kwaliteit van de reconstructies aantasten [16]. Verschillende technieken voor metaal artefact reductie (MAR) werden reeds voorgesteld om het effect van deze artefacten te verminderen. Zo zijn er interpolatiemethoden die foutieve projectiedata trachten te corrigeren [11, 15, 22], iteratieve methoden die metaalartefacten trachten te voorkomen tijdens een iteratieve reconstructie [23] en hybride methoden die zowel gebruik maken van projectiedata als van reconstructies [4]. Bij interpolatiemethoden gaat men eerst op zoek naar de projectieregio van het dichte object. Dit bv. door thresholding van projectiedata of van een reconstructie. Nadien wordt de “ontbrekende projectiedata” in deze regio’s artificieel opgevuld. Hiervoor kan men gebruik maken van bv. polynomiale interpolatie [15], lineaire interpolatie [11] of wavelet interpolatie [22,24]. De keuze van de projectieregio van het dichte object is van cruciaal belang voor de kwaliteit van de reductie. Is de regio te groot, dan wordt informatie van het omliggende weefsel onterecht vervangen. Is de regio te klein, dan is de interpolatie gebaseerd op foutieve data. Het toepassen van discrete tomografie om de projectieregio van het dichte object te bepalen is veelbelovend. Dit is het eerste doel in deze fase van het project. IWT-projectvoorstel

Wim van Aarle

Pagina 12 van 17


Een strategie die als tweede bestudeerd zal worden in de laatste fase van het project, is een combinatie van continue tomografie met discrete tomografie. Daarvoor zal het DART-algoritme uit de vorige fase worden uitgebreid. Bij elke iteratie zullen alle pixels die via een tussentijdse segmentatie worden toegewezen aan het sterk absorberende materiaal en die geen grenspixels zijn, in de volgende iteratie worden vastgezet. Via formule (17) kan men dan een klassiek algebra¨ısch reconstructieprobleem formuleren voor het niet dichte deel van het volume. Na elke iteratie worden de randpixels van het dichte materiaal gesegmenteerd en blijven de overige pixels hun grijswaarden behouden. Zo worden de voordelen van continue tomografie en discrete tomografie gecombineerd om elkaars nadelen op te heffen. Sterk absorberende materialen zullen duidelijk te onderscheiden zijn en de streaking artefacten in de reconstructie van het zachte weefsel zullen worden verminderd. Planning Discrete tomografie op sterk absorberende objecten. Subtaak 1: Literatuurstudie. (3 maanden) Subtaak 2: Bepalen van projectieregio’s van dichte objecten om projectiedata te corrigeren. (2 maanden) Subtaak 3: Combinatie van continue en discrete tomografie in een iteratief algoritme (5 maanden) Subtaak 4: Validatie met experimentele data of simulaties. (2 maanden)

3.4

Conclusie

In vier fasen zal dit project een theoretische en praktische basis leggen voor nieuwe technieken die tomografische reconstructies nauwkeurig kunnen segmenteren. In de eerste en de tweede fase zullen segmentaties van tomogrammen automatisch en objectief worden geëvalueerd met behulp van beschikbare projectiedata. De parameters van een dergelijke segmentatie (zoals de drempelwaarden bij thresholding) kan men zo optimaliseren. Dit moet de segmentatie nauwkeuriger en robuuster maken. In de derde fase zal er getracht worden het aantal projecties dat nodig is om een accurate segmentatie te bekomen, terug te dringen. Zo kan de scantijd en de stralingsbelasting geminimaliseerd worden. Hiervoor wordt discrete tomografie voorgesteld. Daarbij wordt voorkennis over de materialen van het voorwerp gebruikt om uit slechts enkele projecties een gesegmenteerde reconstructie op te stellen. De vierde fase van het project heeft als doel om reconstructies van voorwerpen met sterk absorberende materialen, nauwkeurig te segmenteren. Er zal daarvoor gebruik gemaakt worden van klassieke continue reconstructiealgoritmen in combinatie met de, in de derde fase ontworpen, discrete tomografietechnieken. Elke nieuwe techniek zal ondersteuning bieden voor parallel-beam, fan-beam en cone-beam projectiedata. Bovendien zal er nadruk gelegd worden op de snelle evalueerbaarheid en een gunstige ruimtecomplexiteit. Dit alles moet een nauwkeurige en computationeel aantrekkelijke segmentatie van tomografische reconstructies mogelijk maken.

IWT-projectvoorstel

Wim van Aarle

Pagina 13 van 17


4

Planning

Jaar 1: Globale Thresholding: Automatische globale thresholding van CT-beelden, op basis van gemeten projectiedata. • Literatuurstudie. (3 maanden) • Ontwerp van globale thresholding voor parallel-beam CT met objectieve kwaliteitsmaatstaf. (3 maanden) • Ontwerp van globale thresholding voor cone-beam CT met objectieve kwaliteitsmaatstaf. (6 maanden) • Validatie met experimentele scannerdata of met simulaties. (2 maanden)

Jaar 2: Gevorderde Segmentatietechnieken: Automatische segmentatie van CT-beelden met meer geavanceerdere technieken, op basis van gemeten projectiedata. • Literatuurstudie. (3 maanden) • Ontwerp van nieuwe segmentatietechnieken met objectieve kwaliteitsmaatstaf. (5 maanden) • Validatie van de nieuwe technieken met experimentele scannerdata of met simulaties. (2 maanden) • Vergelijkende studie van de nieuwe segmentatietechnieken met niet-objectieve technieken. (3 maanden)

Jaar 3: Discrete Tomografie: Combinatie van reconstructie en segmentatie. • Literatuurstudie. (3 maanden) • Ontwerp van het DART algoritme voor cone-beam CT data. (4 maanden) • Validatie met experimentele scannerdata of met simulaties. (2 maanden)

Jaar 4: DT op sterk absorberende objecten: Technieken om dichte objecten te reconstrueren met DT. • Literatuurstudie. (3 maanden) • Bepalen van projectieregio’s van dichte objecten om projectiedata te corrigeren. (2 maanden) • Combinatie CT en DT in een iteratief algoritme. (5 maanden) • Validatie met experimentele scannerdata of met simulaties. (2 maanden)

Figuur 4: Gantt chart van de planning. IWT-projectvoorstel

Wim van Aarle

Pagina 14 van 17


5

Toepassingsmogelijkheden

De toepassingsmogelijkheden van nauwkeurige en efficiënte segmentatiealgoritmen voor tomografische reconstructies zijn legio. Op biomedisch vlak, vanouds een vakgebied waar tomografie veel gebruikt is, zijn er zeer veel toepassingen. Gesegmenteerde CT-beelden laten toe om een kwantitatieve analyse uit te voeren op scans van een patiënt. Dit kan een diagnose of een karakterisatie van een ziekte of aandoening vereenvoudigen. Ook op andere gebieden, zoals bij materiaalwetenschappen of in industriële toepassingen kan een segmentatie op een tomografisch beeld van belang zijn. Enkele voorbeelden: • Het trabeculair bot bestaat uit een driedimensionaal, roosterachtig netwerk van plaatjes en staafjes. Via een kwantitatieve analyse daarvan kan osteoporose gekarakteriseerd worden. Via een CT-reconstructie kan dit niet-destructief en volledig driedimensionaal gebeuren [2]. Het is wel belangrijk dat er zuinig wordt omgesprongen met het aantal projecties, dit teneinde de stralingsbelasting voor de patiënt te beperken. • Zachte weefsels, zoals het bloedvatennetwerk, kunnen driedimensionaal gereconstrueerd worden met behulp van discrete CT-reconstructies en angiografie [12]. Daarbij wordt het weefsel ingespoten met een contrastvloeistof van hoge atomaire waarde. • Bij radiotherapie tracht men kwaadaardige cellen te vernietigen door ze te bestralen. Sterk absorberende objecten zoals tandvullingen en protheses kunnen de plaatsbepaling van de schadelijke stralenbundels echter ernstig verstoren [13]. Via simulaties aan de hand van CTbeelden kan men broninstellingen kiezen die de bestraling zo goed mogelijk stuurt. Er is dan wel nood aan uiterst nauwkeurige reconstructies zonder metaal artefacten. • In industriële toepassingen kan men een gesegmenteerd CT-beeld gebruiken om de nauwkeurigheid van een industrieel bouwproces te controleren. Zo kan een elektronische eenheid zoals een microprocessor gereconstrueerd worden via micro-CT. De kwaliteit van de gefabriceerde interconnecties kan men zo controleren en vergelijken met een oorspronkelijke blauwprint. • Er worden steeds hogere eisen gesteld aan de karakteriseringstechnieken die nanosystemen op atomaire schaal moeten onderzoeken op structureel, chemisch en elektronisch vlak. Reconstructies van deze nanosystemen met elektrononentomografen zijn vaak van slechte kwaliteit, wat een segmentatie en een analyse onbetrouwbaar maakt. Bovendien zijn sommige materialen niet bestand tegen de stralingsbelasting van transmissietomografie. Door het toepassen van discrete tomografie kan zowel de kwaliteit van reconstructies worden verhoogd als de stralingsbelasting worden verlaagd. • Het zodanig slijpen van ruwe diamantstenen dat de edelsteen zo waardevol mogelijk is, is een niet-triviaal probleem. De ruwe steen wordt daarom best eerst nagebouwd in een computermodel waarmee een optimalisatieproces de optimale slijplijnen kan lokaliseren. Een tomografische micro-CT-scan en reconstructie zijn hiervoor uitermate geschikt. Het is essentieel dat de segmentatie van een dergelijke reconstructie bijzonder nauwkeurig is. Kleine afwijkingen kunnen namelijk grote gevolgen hebben voor de waarde van de geslepen steen. Discrete tomografie is hiervoor uiterst geschikt.

IWT-projectvoorstel

Wim van Aarle

Pagina 15 van 17


Referenties [1] K.J. Batenburg. Network Flow Algorithms for Discrete Tomography. PhD thesis, Universiteit Leiden, 2006. [2] K.J. Batenburg and J. Sijbers. Discrete tomography from micro-ct data: Application to the mouse trabecular bone structure. Proc. of SPIE Medical Imaging 2006, 6142:1325–1335, 2006. [3] R. Gordon, R. Bender, and G.T. Herman. Algebraic reconstruction techniques (art) for threedimensional electron microscopy and x-ray photography. Journal of theoretical Biology, 29:471– 481, 1970. [4] J. Gu, L. Zhang, Z. Chen, Y. Xing, and Z. Huang. A method based on interpolation for metal artifacts reduction in ct images. Journal of X-Ray Science and Technology, 14:11–19, 2006. [5] G.T. Herman and A. Kuba. Discrete Tomography: Foundations, Algorithms, and Applications. Birkhauser, 1999. [6] G.T. Herman and A. Kuba. Discrete tomography in medical imaging. Proceedings of the IEEE, 91:1612–1626, 2003. [7] G.T. Herman and A. Kuba. Advances in Discrete Tomography and Its Applications. Birkhauser, 2007. [8] C. Johnson and C. Hansen. Visualization Handbook. Academic Press, Inc., Orlando, FL, USA, 2004. [9] S. Kaczmarz. Angenherte auflsung von systemen linearer gleichungen. Bull. Internat. Acad. Polon.Sci. Lettres A, pages 335–357, 1937. [10] A.C. Kak and M. Slaney. Principles of Computerized Tomographic Imaging. IEEE Press, 1988. [11] W.A. Kalender, R. Hebel, and J. Ebersberger. Reduction of ct artifacts caused by metallic implants. Radiology, 164:567–577, 1987. [12] A. Kanitsar, D. Fleischmann, R. Wegenkittl, D. Sandner, P. Felkel, and E. Groller. Computed tomography angiography: A case study of peripheral vessel investigation. Visualization Proceedings, pages 477–593, 2001. [13] Y. Kim, W. Tom, M. Bal, T. McNutt, and L. Spies. The impact of dental metal artifacts on head and neck imrt dose distributions. Radiotherapy and Oncology, 79:198–202, 2006. [14] A. Kuba, L. Rusko, L. Rodek, and Z. Kiss. Preliminary studies of discrete tomography in neutron imaging. IEEE Transactions on Nuclear Science, 52:380–385, 2005. [15] R.M. Lewitt and R.H. Bates. Image reconstruction from projections (iii): Projection completion methods. Optik, 50:189–204, 1978. [16] B. De Man, J. Nuyts, P. Dupont, G. Marchal, and P. Suetens. Metal streak artifacts in x-ray computed tomography: A simulationstudy. Nuclear Science Symposium Conference Record., 3:1860–1865, 1998. IWT-projectvoorstel

Wim van Aarle

Pagina 16 van 17


[17] R.M. Manzke. Cardiac Cone Beam CT. PhD thesis, Kings College London, United Kingdom, 2004. [18] K. Mueller. Fast and Accurate Three-Dimensional Reconstruction from Cone-beam Projection Data using Algebraic Methods. PhD thesis, The Ohio State University, 1998. [19] N. Otsu. A threshold selection method from gray level histograms. IEEE Trans. Systems, lan and Cybernetics, 9:62–66, 1979. [20] R.P.C. Schram. X-ray attenuation, application of x-ray imaging for density analysis. NRG report 20002/01.44395/I, 2001. [21] M. Sezgin and B. Sankur. Survey over image thresholding techniques and quantitative performance evaluation. Journal of Electronic Imaging, 12:146–165, 2004. [22] S.Zhao, D.D. Robertson, G. Wang, B. Whiting, and K.T. Bae. X-ray ct metal artifact reduction using wavelets: An application for imaging total hip protheses. IEEE Transactions on Medical Imaging, 19:1238–1247, 2000. [23] G. Wang, D.L. Snyder, J.A. O’Sullivan, and M.W. Vannier. Iterative deblurring for ct metal artifact reduction. IEEE Transactions on Medical Imaging, 15:657–664, 1996. [24] S. Zhao, K.T. Bae, B. Whiting, and G. Wang. A wavelet method for metal artifact reduction with multiple metallic objects in the field of view. Journal of X-Ray Science and Technology, 10:67–76, 2002. [25] D. Ziou and S. Tabbone. Edge detection techniques - an overview. International Journal of Pattern Recognition and Image Analysis, 8:537–559, 1998.

IWT-projectvoorstel

Wim van Aarle

Pagina 17 van 17

Segmentatie van transmissie tomografie-data voor kwantitatieve analyse op basis van een discreet grijswaardenmodel

Recommend Documents