CT-scan met zo min mogelijk projecties

Technische Universiteit Delft Faculteit Elektrotechniek, Wiskunde en Informatica Delft Institute of Applied Mathematics

CT-scan met zo min mogelijk projecties

Verslag ten behoeve van het Delft Institute for Applied Mathematics als onderdeel ter verkrijging

van de graad van

BACHELOR OF SCIENCE in TECHNISCHE WISKUNDE

door

Melissa de Koning Delft, Nederland Mei 2011

c 2011 door Melissa de Koning. Alle rechten voorbehouden. Copyright

2

BSc verslag TECHNISCHE WISKUNDE

“CT-scan met zo min mogelijk projecties ” (Engelse titel: “CT-scan with least amount of projections )”

Melissa de Koning

Technische Universiteit Delft

Begeleider Dr.ir. M.B. van Gijzen

Overige commissieleden Dr. J.L.A. Dubbeldam

Prof.dr.ir. C. Vuik

Dr. J.G. Spandaw

Mei, 2011

Delft

Inhoudsopgave 1 Voorwoord

5

2 Inleiding 2.1 Indeling verslag . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6 7

3 Medische tomografie

8

4 Geschiedenis van computer tomografie 10 4.1 6 generaties van scanners . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 5 Wiskundige probleembeschrijving en onderzoeksvraag 15 5.1 Reeksontwikkeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 5.2 Kleinste kwadraten methode en de pseudo-inverse . . . . . . . . . . . . . . . . . . . . 21 5.3 ART . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 6 Numerieke experimenten 6.1 Testproblemen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 Horizontale projecties . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Horizontale en verticale projecties . . . . . . . . . . . . . . . . . . . . 6.3.1 Random of gericht horizontale en verticale projecties . . . . . . 6.4 Waaierstralen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.5 Waaierstraal met de pseudo-inverse van de projectiematrix A voor een 6.5.1 Waaierstraal vanuit beginpunt x = 0, y = 0.5n . . . . . . . . . 6.5.2 Waaierstraal vanuit beginpunten x ∈ [−3, ..., 0], y = 0.5n . . . . 6.5.3 Conclusie voor de pseudo-inverse met de waaierstraal. . . . . . 6.6 Waaierstraal met ART . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.6.1 Volgorde van de projecties . . . . . . . . . . . . . . . . . . . . . 6.6.2 Nauwkeurigheid met ART . . . . . . . . . . . . . . . . . . . . . 6.6.3 Benodigde precisie van de fout voor een scherpe foto. . . . . . . 6.6.4 Overeenkomsten tussen de pseudo-inverse en ART. . . . . . . . 6.6.5 Conclusie voor waaierstraal met ART . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . foto van . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . 3×3 . . . . . . . . . . . . . . . . . . . . . . . . . . .

24 24 26 28 29 30 31 33 36 39 40 40 42 44 45 47

7 Validatie 7.1 Strategie om de parameters van te voren te kiezen. 7.2 Projectiestrategie voor een foto van 11x11 . . . . . 7.3 Projectiestrategie voor een foto van 100x100 . . . . 7.4 Projectiestrategie voor een foto van 390 × 390 . . .

. . . .

. . . .

48 48 49 50 51

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

8 Conclusie A Waaierstraal met de pseudo-inverse A.1 Waaierstraal vanuit beginpunten x = 0, y ∈ [1, 2, 3] . . . . . A.2 Waaierstraal vanuit beginpunten x ∈ [−3, ..., 0], y ∈ [1, 2, 3] A.3 Waaierstraal vanuit beginpunt x = 0, y = 1 . . . . . . . . . A.4 Waaierstraal vanuit beginpunten x ∈ [−3, ..., 0], y = 1 . . .

53

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

55 55 57 59 61

B Matlab code B.1 Grafieken met de pseudo-inverse: 1 . B.2 Grafieken met de pseudo-inverse: 2 . B.3 Grafieken met ART: 1 . . . . . . . . B.4 Grafieken met ART: 2 . . . . . . . . B.5 Reconstructie met de pseudo-inverse B.6 Reconstructie met ART . . . . . . .

. . . . . .

. . . . . .

4

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

63 63 66 69 72 76 79

1

Voorwoord

Voor het Bacheloreindproject van de studie Technische Wiskunde mogen de studenten zelf een richting kiezen en heb ik voor de Numerieke Wiskunde gekozen. Dit verslag is dan ook gemaakt bij de Numerieke Wiskunde op de afdeling van Toegepaste Wiskundige Analyse op de faculteit Elektrotechniek, Wiskunde en Informatica. De begeleiding van dit project werd gedaan door Martin van Gijzen. Ik wil hem heel erg bedanken voor zijn uitgebreide begeleiding. Hij heeft me heel veel vrijheid gegeven en hij heeft veel geduld gehad om mij op mijn eigen manier de tijd te laten nemen om te kunnen ontdekken, experimenteren en een eigen kant op te gaan met het onderwerp. Hierdoor heb ik heel veel geleerd over het maken van keuzes bij het aanpakken, verdiepen en uitvoeren van een onderwerp voor een project. Ook heeft dit project opeens de gehele bacheloropleiding een plaats en afgerond geheel gegeven en nog meer mijn interesse gewekt in de studie. Ik heb als voorkennis vooral de vakken Numerieke Wiskunde 1 en Lineaire Algebra 1 nodig gehad uit het eerste en tweede jaar van de studie Technische Wiskunde.

5

2

Inleiding

CT is de afkorting van Computerized Tomography en wordt gebruikt in de medische wereld om foto’s te maken van dwarsdoorsneden van het menselijk lichaam. Dit is erg handig omdat het lichaam van binnen kan worden bekeken zonder dat deze hoeft worden open gemaakt. Dit apparaat werkt echter met röntgenstraling wat schadelijk is voor de menselijke gezondheid. Het is dus van extreem belang dat de röntgenstraling minimaal wordt gehouden zonder dat dit de scherpte van de foto (reconstructie) van de dwarsdoorsnede van het lichaam aantast. Als deze foto niet scherp genoeg is kan deze belangrijke details missen waar iemands leven vanaf kan hangen. Bijvoorbeeld bij het vaststellen van de exacte plaats van een tumor om de bestraling op de juiste plek te richten. Het is dus belangrijk dat er zo min mogelijk röntgenstraling wordt toegepast met een zo hoog mogelijke scherpte van de reconstructie van de dwarsdoorsnede en dit is waar de focus van deze scriptie ligt. Bij het onderzoek maken we gebruik van reeksontwikkeling om het probleem te discretiseren. En bij het reconstrueren gebruiken we de pseudo-inverse en het ART algoritme. We gaan er van uit dat de gemeten data geen ruis bevatten.

6

2.1

Indeling verslag

In Hoofdstuk 3 wordt er verteld over het medisch gebruik van de CT-scan. In Hoofdstuk 4 wordt er verteld over de geschiedenis van de CT-scan met de verschillende projectietechnieken die door de jaren heen zijn gebruikt. In Hoofdstuk 5 wordt uitgelegd hoe we de numerieke oplossing als reeks ontwikkelen met behulp van discretisatie dat leidt tot een lineair stelsel. Ook wordt uitgelegd hoe we de reeks kunnen oplossen met de pseudo-inverse die voortkomt uit de kleinste kwadraten methode. En tot slot wordt uitgelegd hoe en waarom we voor de reeks het ART algoritme gebruiken. In Hoofdstuk 6 zijn de numerieke experimenten voor verschillende projectiestrategiën. Hier zullen we het hebben over de waaierstraal in combinatie met de pseudo-inverse en daarna in combinatie met het ART algoritme om tot een projectiestrategie te komen. In Hoofdstuk 7.1 passen we een projectiestrategie toe op afbeeldingen.

7

3

Medische tomografie

(De volgende informatie komt uit [1] en [6].) Tomografie is het proces om met een tomograaf (= apparaat) een tomogram (foto) te maken. Een CT-scan (Computerized Tomography) wordt ook wel een CAT-scan (Computerized Axial Tomography) genoemd. Het woord tomografie komt uit het Grieks van • tomos: een stuk, sectie of laag; • graphein: schrijven, opnemen, tekenen of omschrijven, en axial komt uit het Latijn van • om de as bewegend. De CT-scan heeft van alle toepassingen van beeldreconstructie uit projecties het grootste effect gehad op de wereld en heeft een revolutie teweeggebracht in de radiologie, vooral omdat de rekencapaciteit van computers toenam. Binnen de radiologie wordt naast röntgenstraling (röntgenfoto, CT-scan) ook gebruik gemaakt van beeldvorming door middel van geluidsgolven (echografie) en magnetische velden (Magnetic Resonance Imaging-scan) (zie Figuur 1).

Figuur 1: Echografie (links) en MRI-scan (rechts). Een tomograaf maakt een tomogram van een dwarsdoorsnede van een lichaam door röntgenstralen door het lichaam te sturen. Een röntgenstraal wordt ook wel een projectie genoemd. Het verschil in energie tussen de uitgaande röntgenstraal, uit de röntgenbron, en de gemeten energie, door de röntgendetector in de detectoreenheid, is de absorptie waarmee gerekend wordt. Elk weefsel heeft een andere absorptiecoëfficiënt, wat de mate van absorptie aangeeft. De absorptiecoëfficiënt wordt uitgedrukt in Hounsfieldeenheden (hu). Wit in Figuur 7 correspondeert met een hoge Hounsfieldeenheid (bot) en zwart correspondeert met een lage Hounsfieldeenheid (lucht). In Figuur 2 is te zien dat lucht een CT-getal van −1000 hu heeft, zacht weefsel een CT-getal in de buurt van 0 hu en bot een CT-getal van 1000 hu. Gegeven de verzwakkingscoëfficiënt µ, dan is het CT-getal, gedefinieerd als

CTgetal =

µweef sel − µwater (hu). µwater

In een tomogram worden vooral de grenzen van de verschillende weefsels goed gezien, waardoor bijvoorbeeld gezonde weefsels worden onderscheiden van tumoren. Een tomogram kan hierdoor 8

Figuur 2: Hounsfieldschaal van CT-nummers. worden gebruikt om de bestraling te plannen voor de vernietiging van een tumor. Het doel hierbij is om een zo hoog mogelijke dosis straling te geven aan het beoogde volume zonder dat de risicovolle organen deze schadelijke dosis krijgen. Ook wordt het verschil tussen weefsel en bot goed gezien, waardoor botbreuken kunnen worden vastgesteld. Door de hoge snelheid van de 6e generatie scanners (33 seconden per omwenteling, 120 plaatjes per seconde) wordt de beweging van het hart zeer nauwkeurig weergegeven (3D-beeld). De preciese locatie in het hart van de hartafwijking kan hierdoor worden vastgesteld. Hartafwijkingen kunnen erfelijk zijn en patiënten krijgen preventief medicijnen toegedient. Als er geen hartafwijking bij de patiënt is vastgesteld hoeft deze niet langer medicijnen te slikken. Doordat er naast breedte ook diepte wordt gezien wordt met een druk op de knop bepaald weefsel, zoals bijvoorbeeld bot (met een CT-getal van 1000 hu), onzichtbaar gemaakt, zodat er achter het verwijderde weefsel gekeken kan worden (zie Figuur 3).

Figuur 3: Het CT-getal van omliggend weefsel kan worden verwijderd. Een nadeel van een CT-scan is dat de röntgenstraling, dat wordt gebruikt voor de verzameling van de benodigde data, elektronen losslaat van atomen. Hierdoor kan DNA beschadigd worden en eiwitten kunnen worden afgebroken. Nog een nadeel is dat een CT-scan niet nauwkeurig genoeg details in één bepaald stuk weefsel kan weergeven. Hiervoor wordt de MRI-scan gebruikt.

9

4

Geschiedenis van computer tomografie

(De volgende informatie komt uit [3], [5] en [6].) Op 8 november 1895 ontdekte Conrad Röntgen Röntgenstraling (X-ray) waardoor röntgenfoto’s konden worden gemaakt. Röntgenstraling is ioniserende straling, omdat het genoeg energie heeft om een elektron uit de buitenste schil van een atoom weg te slaan. De straling is het gevolg van radioactiviteit. Een van de eerste foto’s maakte Röntgen van de hand van zijn vrouw (zie Figuur 4).

Figuur 4: Röntgenfoto. In 1930 publiceerde de Italiaanse radioloog Alessandro Vallebona het eerste fotomateriaal van een klinische dwarsdoorsnede van een lichaam met een Axiale Transversale Scanning-methode op een röntgenfilm (zie Figuur 5). De röntgenbuis en röntgenfilm waren aan een as verbonden met een staaf en ze bewogen synchroon, om de as, in tegenovergestelde richting. Het pivot-punt van de staven (de as) was het focuspunt van het te reconstrueren plaatje. Er werd röntgenstraling uitgezonden van de röntgenbuis naar de röntgenfilm en met de absorptie, die door de tussenliggende weefsels optrad, werd een plaatje van het focuspunt (de weefsels) gereconstrueerd. Het proces van deze methode stond bekend als Tomografie.

Figuur 5: Scan van Vallebona. Betere Transversale Axiale Tomografie-methoden werden hierna onafhankelijk ontwikkeld door: Bernhard Ziedses des Plantes (1932), William Watson (1937), Jean Kieffer (1938) en Shinji Takahashi (1947). De Zuid-Afrikaan Allan McLeod Cormack van Tufts Universiteit in Massachusets, had in 1963 een wiskundig model beschreven voor absorptie van ioniserende straling door weefsel. Sir Godfrey Hounsfield werkte voor EMI (Electrical and Music Industries) Central Research Laboratories in Hayes in Engeland. Daar ontwikkelde hij de Transversale Axiale Tomografie methoden tot de Algebra¨ısche Reconstructie Techniek (ART) (zie Paragraaf 5.3) door gebruik te maken van de computer. ART is een heruitvinding van Kaczmarz’s methode gepubliceerd in 1937. Hounsfield beschreef ook de röntgenabsorptiecoëfficiënt voor weefsel en met het model van Cormack ontwikkelde hij in 1967 de CT-machine. De absorptiecoëfficiënt werd van nu af aan gemeten in Hounsfield eenheden (zie Hoofdstuk 3) en de naam Tomography veranderde in Computerized Axial Tomography (CAT). 10

Het originele prototype (1e generatie scanner) van 1968 nam 160 parallelle projecties in 180 hoeken, elk 1 graad verschil. Het duurde 9 dagen om de informatie van de scans te verzamelen en het duurde 2.5 uur om van de verzamelde data, met ART, een foto te reconstrueren. De eerste CT machine werd ge¨ınstalleerd en in gebruik genomen in het Atkinson Morley Hospital in Wimbledon, Engeland op 1 oktober 1971 en het werd de EMI-scanner genoemd (zie figuur 6). Hiermee werd de eerste hersenscan gemaakt (zie Figuur 7: links).

Figuur 6: EMI-scanner Het duurde 4 minuten om de data voor twee naast elkaar liggende plakken te scannen en de rekentijd was 7 minuten per foto van 80 × 80. Het gebruik van ART werd al snel opgevolgd door Filtered Back Projection. FBP werd uitgevonden door Chris Lemay en werd gepatenteerd door Hounsfield, waardoor een foto van 160 × 160 met dezelfde computer werd berekend in 30 seconden. Nu wordt ART weer meer populair door krachtiger computers.

Figuur 7: Hersenscan gemaakt met EMI-scanner (links) en de Siemens Somatom Emotion (rechts). In 1979 deelden Hounsfield en Cormack de Nobelprijs voor medicijnen. In juni 1980 konden, met de 3e generatie scanner, volume scans worden gemaakt van 1200×1200×270 pixels (= vierkantjes waaraan een uniforme grijswaarde is gekoppeld). Er wordt gezegd dat EMI Central Research Laboratories dankzij het succes van The Beatles via EMI-records het onderzoek kon bekostigen en de eerste EMI-scan modellen kon bouwen voor medisch gebruik. 11

In de 80-er jaren werd door Willie Kalender de spiraal-scan CT (4e generatie, zie Paragraaf 4.1) ontwikkeld toen hij voor Siemens Medical Systems werkte, waardoor 3D-beelden sneller werden weergegeven met minder straling. In 1995 was het bekend dat de Mayo Clinic een röntgen-CT (Dynamic Spatial Reconstructor) had gebouwd waarmee in 20 seconden meer dan een miljard (109 ) projecties werden verzameld, waarmee de grijswaarden van meer dan 106 punten (reconstructies voor 10 × 10 × 10 cm3 ) 60 keer per seconde werden gereconstrueerd [6]. Om nog gedetailleerdere en snellere beelden te maken nam het aantal rijen detectoren (= aantal plakken per scan) toe tot 16 detectoren in 2001 en 64 in 2004. Als we naar Figuur 7 kijken zien we de ontwikkeling van de CT-scan terug in het verschil in scherpte van een hersenscan gemaakt met de EMI-scanner (1971, links) en de Siemens Somatom Emotion (2010, rechts).

12

4.1

6 generaties van scanners

Om absorptiecoëfficiënten in de pixels van een tomogram te berekenenen hebben we meerdere projecties nodig. Het systeem om de projecties te nemen is steeds veranderd, waardoor er een onderscheid kan worden gemaakt in 6 generaties van systemen: • 1e generatie. Een enkele röntgenbron en een enkele detector verzamelen alle data voor 1 enkele plak. De bron en de detector bewegen samen en de straal wordt getransleerd over de patiënt om een set parallelle projecties te krijgen in 1 hoek. Het bron-detector paar wordt dan een klein beetje geroteerd en soortgelijke metingen worden opnieuw uitgevoerd tijdens het opnieuw transleren over de patiënt. Dit proces wordt voor elke projectiehoek herhaald tot een maximale hoek van 180 graden. Dit wordt ook wel de translatie/rotatie-scanner genoemd. • 2e generatie. Er wordt een array van detectoren gebruikt die tegenover de röntgenbron liggen, waardoor de röntgenbron een waaier van stralen uit kan zenden. Door hetzelfde translatie/rotatie principe van de 1e generatie krijgen we door de translaties meerdere parallelle projecties. De maximale rotatie is 180 graden. Door de meerdere projecties per translatie is deze translatie/rotatie-scanner efficiënter.

Figuur 8: 1e generatie CT-scan (links) en 2e generatie CT-scan (rechts). • 3e generatie. Er wordt in één keer een waaierstraal projectie van de gehele dwarsdoorsnede van de patiënt gemaakt. Door de grootte van de detectorarray hoeft deze niet meer te transleren, maar alleen te roteren rond het te reconstrueren object. De maximale rotatie is 360 graden. Omdat de röntgenbron en de detector beiden roteren wordt deze scanner ook de rotatie/rotatie-scanner genoemd.

Figuur 9: 3e generatie CT-scan.

13

• 4e generatie. De röntgenstralen worden uitgezonden naar een stationaire, 2-dimensionale, array van detectoren die geheel rondom de patiënt zijn. De detectorring ligt buiten het circulaire pad van de röntgenbron. De maximale rotatie is 360 graden. De stralen worden, in plaats van in een waaiervorm, ook wel in een kegelvorm (spiraal) uitgezonden. Het lichaam wordt met een constante snelheid door het focuspunt bewogen. Deze scanner staat bekend als de rotatie/stationaire-scanner. • 5e generatie. Hetzelfde als de 4e generatie, maar de detectorring ligt binnen het circulaire pad van de röntgenbron en wijkt uit als de röntgenbron langs komt. Deze scanner staat ook bekend als de rotatie/uitwijk-scanner.

Figuur 10: 4e en 5e generatie CT-scan (links, midden) en 3D opname darm (rechts). • 6e generatie. Ook wel de cine CT, milliseconde CT, ultrasnelle CT of elektronenstraal CT genoemd. Deze methode heeft geen mechanische scan beweging meer. Een elektronenstraal maakt dezelfde beweging als de röntgenbron bij de 4e generatie. De elektronen worden versneld richting een anode die om de patiënt heen is en in de anode worden de elektronen tot stilstand gebracht waardoor er röntgenstraling ontstaat. De röntgenstralen maken dezelfde beweging als de röntgenstralen bij de 4e generatie. De tijd om data te verzamelen voor 1 plak is 0.05−0.1 seconde.

Figuur 11: 6e generatie CT-scan.

14

5

Wiskundige probleembeschrijving en onderzoeksvraag

(De volgende informatie komt uit [6], [9], [1].) Er zijn verschillende computer-algoritmen voor de reconstructie van foto’s voorgesteld door vele auteurs. Ze hebben allen als input de beschikbare projectiedata en als output een schatting van de originele structuur (dichtheid door middel van absorptiecoëfficiënten uitgedrukt in grijswaarden). We kijken naar de 3e generatie CT-scanners. Hier worden 2-dimensionale secties uit het 3-dimensionale object genomen en per 2-dimensionale sectie (één plak) worden alle straal-sommen verzameld voor een serie röntgenstralen. Elke plak wordt onafhankelijk van de andere verschillende plakken gereconstrueerd en de verschillende plakken worden op elkaar gestapeld om het 3-dimensionale object te reconstrueren (zie Figuur 5).

Figuur 12: links: 3e generatie scanner (zie hst. 4.1), rechts: projecties door eéń plak We zijn dus ge¨ınteresseerd in het reconstrueren van een klinische dwarsdoorsnede van bijvoorbeeld 3 mm dik. Deze doorsnede kunnen we onderverdelen in kleine vierkante blokjes die allemaal dezelfde afmeting hebben. Deze blokjes worden meestal naar verwezen als volume elementen of afgekort “voxels”. Een CT-getal (zie Hoofdstuk 3) is proportioneel aan de gemiddelde lineaire absorptie in een voxel. Elke voxel heeft dus een absorptiecoëfficiënt met een CT-getal waaraan een uniforme grijswaarde is gekoppeld. Deze kleine vierkantjes worden foto elementen genoemd of afgekort “pixels” (van picture elements). Het probleem is om de grijswaarde per pixel van de geprojecteerde beelden te reconstrueren. Het reconstructieprobleem wordt: gegeven een set projecties van een 2-dimensionaal object, schat de dichtheid per pixel. Het algemene gebied van het bepalen van de distributie van parameters in een structuur die invloed heeft op de metingen wordt binnen de toegepaste wiskunde een “inverse probleem” genoemd [6]. Het reconstrucieprobleem is dus een “inverse probleem”. Eén klinische dwarsdoorsnede noemen we één plak en deze heeft de volgende aannamen: • De plakken zijn oneindig dun; • Alle r¨ ontgenfotonen bewegen in dezelfde rechte lijn (die in de oneindig dunne plak liggen)[6]; • De absorptiecoëfficiënt hangt alleen af van het weefsel waar het doorheen gaat; • Er is geen onderscheid tussen voxels en pixels; • Door de oneindig dunne plak is de grijswaarde in elke punt y 0 (x0 ) proportioneel aan de relatieve absorptiecoëfficiënt in dat punt.

15

Daarom wordt de theorie achter de reconstructiealgoritmen ook wel fotoreconstructie van projecties genoemd. Nu bestaat ons reconstructieprobleem uit het bepalen van de n × n digitale versie van een foto uit de straal-sommen van de grijswaarden van de foto. Er komen verscheidene wiskundige en rekenkundige moeilijkheden naar boven in de reconstructie, doordat de projectiedata ontoereikend zijn en ruis bevatten. De precisie van al de reconstructie algoritmen worden hierdoor be¨ınvloed. In Hoofdstuk 4.1 hebben we gezien dat er bij de verschillende generatie scanners verschillende projectiestrategiën worden toegepast. De projecties worden parallel, in een waaierstraal of in een kegelvorm genomen. Als de projecties zijn genomen moet uit deze projecties de foto worden gereconstrueerd. Onderzoeksvraag: Wat is de beste strategie om de projecties te kiezen? Dat wil zeggen, hoe moeten de projecties worden gekozen om met zo min mogelijk projecties een reconstructie van voldoende kwaliteit te maken?

16

5.1

Reeksontwikkeling

Bij het bepalen van de grijswaarden in pixel xj van een foto gebruiken we reeksontwikkeling [6], waarbij de (x, y)-coördinaten worden uitgedrukt in de j e pixelcoördinaat. De ie projectie-som is dan: 2

projectie-somi =

n X

ai,j xj

j=1

Hierbij is xj de te bepalen grijswaarde van de j e pixel uit de gemeten projectie-sommen (1 ≤ i ≤ n2 ). En ai,j is de lengte door een pixel j in de ie meting, die 1 is als deze door de pixel xj gaat en anders 0 is. Bekijk het volgende voorbeeld van een foto van 2 × 2: 2 3 F= . 4 5 Als we aannemen dat de beschikbare projecties twee rij-sommen en twee kolom-sommen zijn dan krijgen we 5 (= 2 + 3), 9 (= 4 + 5), 6 (= 2 + 4) en 8 (= 3 + 5). Deze “input data” noemen we de “projectie-som-vector” b. Met b (= gemeten) reconstrueren we a b . F= c d De 4 onbekenden in F vinden we door het systeem a + b = 5, c + d = 9, a + c = 6, b + d = 8, in matrixvorm te schrijven: 

1  0   1 0

1 0 0 1

0 1 1 0

 a 0  b 1   0  c 1 d

 5   9  =    6  8 



(1)

en op te lossen. De 4 × 4 matrix is de projectiematrix en als we deze A noemen krijgen we in het algemeen de vergelijking Ax = b, waarbij alle pixels van de fotomatrix F per rij onder elkaar gezet zijn in wat de kolomvector x oplevert. [10] Bekijken we nu als voorbeeld een foto van 3 × 3:   2 3 4 F= 5 6 7  8 9 10 17

en nemen we aan dat de beschikbare projecties 3 rij-sommen en 3 kolom-sommen zijn dan reconstrueren we met b opnieuw matrix F, ditmaal dus:   a b c F =  d e f . g h i We lossen weer het stelsel Ax = b op om de 9 onbekenden in F te vinden. In matrixvorm hebben we ditmaal:         

1 0 0 1 0 0

1 0 0 0 1 0

1 0 0 0 0 1

0 1 0 1 0 0

0 1 0 0 1 0

0 1 0 0 0 0

0 0 1 1 0 0

0 0 1 0 1 0

0 0 1 0 0 1

            

a b c d e f g h i

           =         

9 18 27 15 18 21

    .   

(2)

Als we nu een foto van 4 × 4 bekijken en we nemen aan dat de beschikbare projecties 4 rij-sommen en 4 kolom-sommen zijn dan krijgen we een systeem met 8 beschikbare projecties en 4 · 4 = 16 onbekenden in F. De projectiematrix A heeft in dit geval 8 rijen en 16 kolommen. De projectiematrix A kan in geen van de 3 gevallen volgens de normale manier ge¨ınverteerd worden om een unieke oplossing te krijgen, want deze is 1. Onderbepaald (rang < n2 ), de oplossing is niet uniek. De 3 systemen hebben een unieke oplossing als ze consistent zijn en de rang van A is n2 . Als we de projectiematrix A in echelonvorm zetten zien we dat in systeem (1) rang 3 < 4. Omdat het systeem consistent is heeft deze wel een oplossing, maar geen unieke oplossing. Zoals we in Tabel 1 kunnen zien neemt bij een grotere fotomatrix F het aantal pixels sneller toe dan de rang. De rang van deze systemen zullen dus nooit in de buurt komen van rang n2 om een unieke oplossing te krijgen voor F. Foto’s met n ≥ 2 pixels hebben geen unieke oplossing, maar oneindig veel oplossingen. 2. Niet vierkant (voor systeem (2): determinant bestaat niet). beschikbare projecties: n rij-sommen en n systeem Foto pixels projecties rang (1) 2×2 4 4 3 (2) 3×3 9 6 5 4×4 16 8 7 5×5 25 10 9 6×6 36 12 11

kolom-sommen A onderbepaald en wel vierkant onderbepaald en niet vierkant onderbepaald en niet vierkant onderbepaald en niet vierkant onderbepaald en niet vierkant

Tabel 1: Eigenschappen systeem Ax=b voor verschillende afmetingen foto’s.

18

Het aantal beschikbare projecties van de foto’s breiden we uit door het toevoegen van diagonaalsommen, die precies door het middelpunt van de pixels gaan. De onbekenden in F voor een foto 2 × 2 en 3 × 3 vinden we door de systemen in matrixvorm (3) en (4) te schrijven en op te lossen: 

Ax = b :



Ax = b :

                   

1 0 0 1 0 0 1 0 0 0 0 0

1 0 0 0 1 0 0 0 1 0 1 0

      

1 0 0 0 0 1 0 1 0 0 0 0

1 0 1 0 1 0

0 1 0 1 0 0 0 0 0 1 1 0

1 0 0 1 0 1

0 1 0 0 1 0 1 1 0 0 0 0

0 1 1 0 0 1

0 1 0 0 0 1 0 0 1 0 0 1

0 1 0 1 1 0

0 0 1 1 0 0 0 1 0 0 0 0



 

      

0 0 1 0 1 0 0 0 0 1 0 1



 a   b  =   c   d

0 0 1 0 0 1 1 0 0 0 0 0

5 9 6 8 7 7

                    

    .   

(3)

 a b c d e f g h i

               =                 

9 18 27 15 18 21 18 18 10 14 8 16

           .         

(4)

De projectiematrix A kan met de diagonale projecties nog steeds niet volgens de normale manier ge¨ınverteerd worden om een unieke oplossing te krijgen, want deze is 1. Overbepaald (beschikbare projecties > n2 ) Als we naar Tabel 2 kijken zien we dat voor de systemen (3) en (4) geldt: rang= n2 , dus deze hebben wel een unieke oplossing. Ook zien we dat bij een grotere fotomatrix F het aantal pixels sneller toeneemt dan de rang. De rang van de systemen n ≥ 4 zullen dus nooit in de buurt komen van de rang n2 om een unieke oplossing te krijgen voor F. Foto’s met n ≥ 4 pixels hebben geen unieke oplossing, maar oneindig veel oplossingen. 2. Niet vierkant. Projectiematrix A heeft rang n2 in systemen (3) en (4), dus deze twee systemen hebben een unieke oplossing.

19

beschikbare projecties: n rij-sommen, n kolom-sommen en 4n-6 diagonaal-sommen systeem Foto pixels projecties rang A (3) 2×2 4 6 4 overbepaald en niet vierkant (4) 3×3 9 12 9 overbepaald en niet vierkant 4×4 16 18 15 over- en onderbepaald en niet vierkant 5×5 25 24 20 onderbepaald en niet vierkant 6×6 36 30 25 onderbepaald en niet vierkant 7×7 49 36 30 onderbepaald en niet vierkant 8×8 64 42 35 onderbepaald en niet vierkant Tabel 2: Eigenschappen systeem Ax=b voor verschillende afmetingen foto’s.

20

5.2

Kleinste kwadraten methode en de pseudo-inverse

Om een oplossing van niet vierkante stelsels te bepalen gebruiken we de kleinste kwadraten methode. ¯ is degene met minimale norm van het De standaard oplossing onder alle oplossingen van A¯ x=b residu:

krkminimaal = ¯ b − b = kA¯ x − bk ¯ die r miniWe zijn op zoek naar een oplossing die voldoet aan de normaalvergelijking, dus voor x maliseert: AT A¯ x = AT b. ¯ is niet uniek als AT A niet inverteerbaar is. We zoeken dus de oplossing x ¯, De oplossing voor x waarbij bovendien geldt kxkminimaal . We moeten onthouden dat de rij-ruimte en null-ruimte van de 2 ¯ kan worden opgedeeld projectiematrix A loodrecht op elkaar staan in Rn . Dit betekent dat elke x ¯ k ∈ Kol(AT )en w ∈ Null(A). in twee loodrecht op elkaar staande stukken, x ¯ ¯ 0 een van de kleinste kwadraten oplossingen is voor de vergelijking A¯ We nemen aan dat x x = b. ¯ k in de rij-ruimte en w in de null-ruimte dan krijgen we Nemen we x ¯0 = x ¯ k + w. x Nu hebben we drie belangrijke punten: ¯ ¯ k op zichzelf een oplossing van A¯ 1. Omdat Aw = 0 is component x x = b. ¯ A¯ xk = A (¯ xk + w) = A¯ x0 = b ¯ delen dezelfde component x ¯ k in de rij-ruimte en verschillen 2. Alle oplossingen van A¯ x = b alleen in de null-ruimte component w. ¯ k + w gebruiken we de wet van Pythagoras, omdat de twee 3. Voor de lengte van een oplossing x componenten loodrecht op elkaar staan: k¯ xk + wk2 = k¯ xk k2 + kwk2 . ¯ gelijk is aan x ¯ k . Als we de null-ruimte component Hieruit concluderen we dat de kxkminimaal = x ¯ k wordt de kleinste gelijk stellen aan 0, blijft er een oplossing over die geheel in de rij-ruimte ligt. x kwadraten minimum norm oplossing genoemd. De matrix die de kleinste kwadraten minimum norm oplossing van het systeem Ax = b optimaal bepaalt is de pseudo-inverse 1 A+ [8]: ¯ k = A+ b. x

1

A+ wordt ook wel de Moore-Penrose inverse genoemd, naar zijn ontdekkers, maar staat bekender onder de naam gegeneraliseerde inverse. Er zijn matrices die een aantal van dezelfde eigenschappen delen als die voor A+ , die ook gegeneraliseerde inverse heten. Om verwarring te voorkomen gebruiken we de term pseudo-inverse.

21

Voor elke A is er een unieke pseudo-inverse A+ , die voldoet aan de volgende 4 voorwaarden: 1. AA+

T

= AA+

2. A+ A

T

= A+ A

3. A+ AA+ = A+ 4. AA+ A = A We kunnen de volgende gevallen voor A onderscheiden: • Als A inverteerbaar dan AA−1 = A−1 A = I (=eenheidsmatrix) waardoor A−1 = A+ [10]. • Als A een p × n2 −matrix met p > n2 en AT A is inverteerbaar dan A+ = (AT A)−1 AT [3]. ¯ = Pb met P de orthogonale • Als A en AT A geen van beide inverteerbaar dan AA+ b = A¯ x=b 2 T + n projectie van R op de Kol(A ), waaruit volgt P=AA = A(AT A)+ AT . Uit de 4 voorwaarden kunnen we o.a. afleiden dat: + • A+ =A [10]. • Kol(A+ )=Rij(A) en Rij(A+ )=Kol(A), waaruit volgt rang(A+ )=rang(A). Aldus AA+ en A+ A zijn symmetrisch (AA+ 6=A+ A) en rang(A) = rang(AA+ ) = rang(A+ A). • Als A een p × n2 matrix, waarvan p de dimensie van b en n de orde van de fotomatrix F, dan is A+ een n2 × p matrix. En dus is A+ b een n2 × 1 vector [10].

22

5.3

ART

Om medisch bruikbare informatie te krijgen moeten de gereconstrueerde foto’s gedetailleerd zijn. De precisie gebeurt ongeveer in 1 pixel per millimeter. Als we bijvoorbeeld menselijke hersens in een doos doen van 20 × 20 × 10 cm3 , dan moeten we de grijswaarden van 4 × 106 pixels berekenen. Om dit te doen hebben we dus 4 × 106 onafhankelijke projecties nodig en moeten we een aantal miljoen data punten berekenen uit een aantal miljoen projecties. Voor een stationair object als de hersens kunnen we het probleem versimpelen door data te verzamelen per plak waardoor de projectiematrix 20.000 data punten krijgt met minstens 20.000 projecties (A is 20.000 × 20.000) [6]. Vanwege de grote afmetingen van de projectiematrix A werken bijna alle algoritmen, voor een oplossing van een lineair systeem van vergelijkingen, met iteratieve technieken die geen pseudoinverse van matrix A nodig hebben. Dit zou anders een te grote rekenkracht van de computer eisen [3]. De Algebra¨ısche Reconstructie Techniek (ART) is zo’n iteratieve techniek. ART werkt door de informatie te gebruiken van 1 projectie-som bi (= 1 element uit de projectiesom-vector b). Dan worden de pixels gecorrigeerd die hebben bijgedragen aan de projectie-som bi . Hierna wordt dezelfde procedure toegepast op de volgenden projectie. ART wordt ook wel de ray-by-ray reconstructiemethode genoemd. ¯ de nulvector genomen. Voor de allereerste projectie-som wordt voor x ¯ oud worden gecorrigeerd per projectie: De oude, benaderde pixelwaarden x ¯ nieuw = x ¯ oud + x ¯ correctie . x ¯i en kolomvector x ¯; Voor het systeem A¯ x = b schrijven we per rijvector a ¯ = ai x ¯ oud + ai x ¯ correctie = bi . ¯ oud + x ¯ correctie = ai x ai x Dit levert, met behulp van de pseudo-inverse van één projectie ai , voor de correctie op de oude pixelwaarden: ¯ correctie x

¯ oud = a+ i bi − ai x + T ¯ oud = ai aTi ai bi − ai x −1 T ¯ oud . = ai aTi ai bi − ai x

=

aTi ¯ oud ) (bi − ai x ai aTi

¯ correctie invullen in x ¯ nieuw krijgen we het ART algoritme: Als we x ¯ nieuw = x ¯ oud + x

aTi ¯ oud ). (bi − ai x T ai ai

Het nemen van ART iteraties breken we af als de reconstructie op het oog acceptabel is.

23

6

Numerieke experimenten

Bij de reconstructie van een dwarsdoorsnede van een lichaam hebben we de pseudo-inverse en het ART algoritme nodig. Om een zo scherp mogelijke reconstructie te maken is het ook belangrijk dat we weten hoe de projecties genomen moeten worden zodat ze, samengevoegd tot een projectiematrix A, deze projectiematrix een rang van n2 krijgt. Als doel hierbij hebben we de volgende onderzoeksvraag. Onderzoeksvraag: Hoe moeten de projecties worden gekozen om met zo min mogelijk projecties een reconstructie van voldoende kwaliteit te maken?

6.1

Testproblemen

In de experimenten drukken we de kwaliteit van een reconstructie uit in de fout: Pn2 nieuw x ¯ − xi i=1

i

n2

.

¯ nieuw De fout gebruiken we om te meten hoe ver de benaderde pixelwaarden x zijn geconvergeerd i naar de echte pixelwaarden xi . (In de praktijk weten we de echte pixelwaarden niet, maar de projectiestrategie staat dan al vast. We willen nu juist onderzoeken welke projectiestrategie we moeten nemen.) Om onze projectiestrategie te testen definiëren we een aantal testproblemen (zie Figuur 13): • Kleine testproblemen: – Foto van random kleuren: 2 × 2, 3 × 3, 4 × 4 en 11 × 11 pixels: ∗ De kleuren zijn in RGB grijswaarden met een precisie van 10−4 ; ∗ Voor een foto van n × n zijn minstens n2 projecties nodig. • Grote testproblemen: – Foto van een gezicht: 100 × 100, Foto van een hersenscan: 390 × 390; ∗ De kleuren zijn in RGB grijswaarden met een precisie van 10−4 ; ∗ De berekening van de pseudo-inverses neemt te veel tijd bij het bepalen van de optimale hoek voor een minimum aan projecties; ∗ Voor een acceptabele reconstructie voor een foto van n × n hebben we minder dan n2 projecties nodig.

24

Figuur 13: Illustraties van de testproblemen. Hierbij onderzoeken we de fout bij het nemen van gerichte projecties en random projecties voor: • Horizontale projecties (voor verticale projecties geldt hetzelfde); • Horizontale en verticale projecties; • Waaierstraal projecties; – Pseudo-inverse van de projectiematrix A; – ART.

25

6.2

Horizontale projecties

We onderzoeken met ART wat er gebeurt met het residu bij het nemen van één projectie en of het mogelijk is om met alleen horizontale projecties een foto te reconstrueren.

Figuur 14: De originele foto van 2 × 2 Bij iedere ART iteratie wordt het residu van de projectie 0 gemaakt. Dit betekent dat het residu evenredig over de pixels wordt verdeeld. Als we twee keer achter elkaar dezelfde projectie nemen kan het residu van deze projectie niet minder worden en zal de fout dus ook niet veranderen (minder worden). De fout van de pixels kunnen zowel in de projectie als in de hele foto nog heel erg groot zijn. Bij het nemen van alleen horizontale projecties (maakt niet uit in welke volgorde) kan het residu per projectie, na alle n rijen te hebben gehad, ook niet minder worden, omdat het residu van elke rij op 0 is gesteld en verdeeld is over de pixels. (Hetzelfde geldt voor verticale projecties.) In Figuur 15 zien we een foto en grafieken van foto’s van 2×2 waarbij 50 keer 2 horizontale projecties (= 100 ART iteraties) zijn genomen. De grijswaarden worden evenredig over de horizontale pixels verdeeld waardoor er maar twee kleuren te zien zijn. Er zijn dus maar twee grijswaarden, ongeacht hoeveel projecties er worden genomen. De fout gaat niet naar 0, maar blijft na de tweede projectie constant.

Figuur 15: Bij het nemen van 50 keer 2 horizontale projecties (= 100 ART iteraties) zien we in de reconstructie van de foto 2 strepen. De pixelwaarden convergeren dus niet naar de originele pixelwaarden. Dit kunnen we ook zien in de grafiek van de fout doordat de fout niet naar 0 gaat. Conclusies: Door met ART alleen horizontale projecties te nemen voor een foto van n×n zullen we n verschillende strepen van grijswaarden zien in de foto en kunnen we dus geen foto reconstureren. In Hoofdstuk 5.1 hebben we gezien dat als de projectiematrix A onderbepaald is we geen foto kunnen reconstrueren met de pseudo-inverse van A. Hier hebben we twee projecties bij 4 pixels, dus de projectiematrix is 26

zeker onderbepaald. Het aantal projecties is n voor een foto van n × n en dus altijd onderbepaald.

27

6.3

Horizontale en verticale projecties

We onderzoeken met ART of we een foto kunnen reconstrueren met alleen horizontale en verticale projecties.

Figuur 16: Originele foto’s van 2 × 2 (links) en 3 × 3 (rechts). Als we nu alleen horizontale en verticale projecties nemen krijgen we de reconstructie van de foto, de grafieken van de grijswaarden per pixel en de grafiek van de fout in Figuur 17. De fout blijft constant na twee horizontale en twee verticale projecties. In Figuur 18 zien we dat voor een foto van 3 × 3 na 6 projecties (3 horizontale en 3 verticale projecties) de fout ook weer constant blijft.

Figuur 17: 25 keer 4 projecties levert 100 ART iteraties.

Figuur 18: 17 keer 6 projecties levert 102 ART iteraties. Conclusies: Bij een foto van n × n zullen we n2 verschillende lijnen zien in de grafieken voor de grijswaarde per pixel. 28

In de tabel in Hoofdstuk 5.1 hebben we al gezien dat de matrix A voor alleen horizontale en verticale projecties onderbepaald is en er geen unieke oplossing is. 6.3.1

Random of gericht horizontale en verticale projecties

We onderzoeken of het nemen van random projecties sneller gaat dan het nemen van een van te voren vastgestelde volgorde van de projecties. Als we naar Figuur 19 en 20 kijken zien we dat er meer ART iteraties nodig zijn als we ze random nemen, dan wanneer we ze gericht nemen om tot een minimale fout te komen en de pixelwaarden niet verder convergeren. Maar het aantal te nemen projecties blijft 4 (en 6 voor een foto van 3 × 3) voor gericht en random genomen projecties.

Figuur 19: 100 ART iteraties van de 4 projecties gericht (linker twee Figuren) en random (rechter twee Figuren).

Figuur 20: 100 ART iteraties van de 6 projecties gericht (linker twee Figuren) en random (rechter twee Figuren). Conclusies: Als we de horizontale en verticale projecties random nemen zijn er meer ART iteraties nodig voordat de fout minimaal is. We kunnen de projecties dus beter in een vastgestelde volgorde nemen en deze volgorde herhalen.

29

6.4

Waaierstralen

We hebben in Paragraaf 4.1 gezien dat de projectiestrategiën, voor één plak, de waaierstraal als basis hebben. De waaierstraal neemt vanuit een vast punt opeenvolgende projecties met een vaste hoek tussen elke projectie. We passen de waaierstraal toe op de pseudo-inverse van de projectiematrix A (zie Hoofdstuk 5.2) en op de methode van reeksonwikkeling met ART (zie Hoofdstuk 5.3). In Hoofdstuk 5.1 hebben we gezien dat een foto van n × n een rang van n2 moet hebben om de foto te reconstrueren. Om de rang te verhogen moet elke, aan de projectiematrix A toegevoegde projectie, 1 pixel verschil hebben met alle projecties die al zijn genomen. Dit is wat de waaierstraal doet bij een juiste hoek tussen de projecties. Voor een foto van n × n kunnen we een waaierstraal nemen vanuit de punten van de linker, boven, rechter en onder ribbe. Omdat de pixels van een foto discrete waarden hebben en de projecties continue functies zijn maken we deze functies discreet door alle reëele waarden tot en met het eerstvolgende gehele getal dit eerstvolgende gehele getal te maken. Als we een projectie uitgedrukt als functie y(x) hebben en een foto als y 0 (x0 ) dan krijgen we voor x ∈ (0, ..., 1] ⇒ x0 = 1, x ∈ (1, ..., 2] ⇒ x0 = 2, ... , x ∈ (n − 1, ..., n] ⇒ x0 = n. Hetzelfde geldt voor y. Zie Figuur 21. Voor de boven, rechter en onder ribbe geldt dezelfde discretisatie van functie y(x) naar foto y 0 (x0 ) bij een rotatie van 90 graden rond het middelpunt (0.5n, 0.5n). In Figuur 23, 24 en 25 zijn de twee verschillende schalen samengevoegd tot de schaal van de foto en zien we welke delen van de projecties in welke pixel terecht komen.

Figuur 21: Links: Waaierstralen kunnen we nemen vanuit de linker, boven, rechter en onder ribbe van een foto. rechts: Drie waaierstralen met 30 graden tussen de projecties. De schaal van de foto (links: y 0 (x0 )) is anders dan de schaal van de projecties (rechts: y(x)).

30

6.5

Waaierstraal met de pseudo-inverse van de projectiematrix A voor een foto van 3 × 3

Om te kijken hoe we het beste een waaierstraal met zo min mogelijk projecties kunnen nemen gebruiken we een foto van 3 × 3. Om de foto te reconstrueren nemen we de pseudo-inverse van de projectiematrix A en kijken we bij welk aantal waaierstraalprojecties de fout van de reconstructie voldoende klein wordt. Dit doen we met 3 verschillende varianten bij de experimenten voor het toepassen van een waaierstraal. We nemen een waaierstraal vanuit elke pixel y 0 ∈ [1, n], vanuit y = 1 of vanuit y = 0.5n. Ook kijken we of we het aantal projecties nog verder kunnen verkleinen door het beginpunt x van de waaierstraal te variëren van x = −3 tot en met x = 0. Figuur 6.5 is de originele foto waarmee we experimenteren.

Figuur 22: De originele foto.

31

In Figuur 23, 24 en 25 is als voorbeeld een foto van 3 × 3 genomen met 30 graden als hoek tussen de projecties. Voor y ∈ [1, 2, 3]

Figuur 23: Links: x = 0, midden: x = −1, rechts: x = −2. Voor y = 1

Figuur 24: Links: x = 0, midden: x = −1.5, rechts: x = −2.5. Voor y = 0.5n

Figuur 25: Links: x = 0, midden: x = −0.5, rechts: x = −2. We presenteren hier de resultaten voor y = 0.5n. In de Appendix A presenteren we de resultaten voor y 0 ∈ [1, n] en y = 1.

32

6.5.1

Waaierstraal vanuit beginpunt x = 0, y = 0.5n

We onderzoeken wat het minimum aantal projecties is bij het nemen van een waaierstraal vanuit punt (0, 0.5n). Deze waaierstraal nemen we vanuit één ribbe, vanuit twee opeenvolgende ribben, vanuit twee tegenover elkaar staande ribben en vanuit alle ribben. We zien per genomen ribbe in de middelste grafiek het aantal projecties per hoek en in de rechter grafiek de log van de fout per hoek. Hierdoor kunnen we per hoek aflezen bij welk aantal projecties er een minimum fout is. Er zijn in de grafiek van de log van de fout per hoek pieken te zien. Dit komt doordat de projecties bij verandering van de hoek verschuiven. Bij zo een verschuiving kan de manier waarop de pixels per projectie worden meegenomen veranderen en dit kan net het verschil uit maken tussen bijvoorbeeld een rang van 8 en een rang van 9 voor de projectiematrix. Bij een rang van 9 hebben we een minmale fout en dus een unieke oplossing. Bij een rang van 8 niet. Het minimum aantal projecties is 10 (zie Figuur 26, 27 en 28). Het maakt niet uit of we 1 of 2 ribben nemen. Bij 4 ribben hebben we een minimum aantal projecties van 12 (zie Figuur 29). E´ en ribbe We zien in de grafiek dat er bij een hoek van 16 graden en groter het aantal projecties kleiner wordt dan 10, maar de fout is hier niet minimaal. Voor een hoek groter dan 14 graden is er dus geen unieke oplossing.

Figuur 26: Minimum aantal projecties = 10, hoek = 14 graden.

33

Linker en boven ribbe We zien in de grafieken dat er voor 28 en 32 graden een minimum aantal projecties van 10 is. Tussen 28 en 32 graden is er ook een minimum aantal projecties van 10, maar hier is de fout te groot voor een unieke oplossing. We zien ook dat het minimum aantal projecties naar 8 gaat voor 33, 34 en 35 graden, maar hier is de fout ook te groot. De rang van de projectiematrix kan hier nooit 9 worden en er is dus geen unieke oplossing mogelijk.

Figuur 27: Minimum aantal projecties = 10, hoek = 32 graden. Linker en rechter ribbe Voor een hoek van 29 graden en groter zien we dat het minimum aantal projecties van 10 naar 8 gaat, maar de fout is voor een hoek groter dan 29 graden te groot om een unieke oplossing te krijgen.


34

Alle ribben In de grafiek van de fout zien we ook weer pieken met een hele grote fout. Ook hier komt dit weer door de verschuiving van de projecties waardoor de pixels per projectie zodanig verschillen dat de rang kleiner wordt dan 9 en er geen unieke oplossing mogelijk is. We zien een gat in de grafiek van het aantal projecties per hoek, omdat het aantal projecties hier kleiner is dan 10.

Figuur 29: Minimum aantal projecties = 12, hoek = 50 graden. In de grafieken met het aantal projecties hebben we de dubbele projecties niet meegeteld. Als er meerdere hoeken mogelijk zijn hebben we de grootste hoek genomen om zo min mogelijk dubbele projecties door een foto te krijgen. Conclusie ribbe één ribbe twee opeenvolgende ribben tegenover elkaar liggende ribben vier ribben

minimum projecties 10 10 10 12

hoek (graden) 14 32 29 50

Tabel 3: Conclusie voor waaierstralen vanuit beginpunt (0, 0.5n).

We kunnen het aantal projecties nog verkleinen door de x-coördinaat van de waaierstraal verder van de foto te plaatsen.

35

6.5.2

Waaierstraal vanuit beginpunten x ∈ [−3, ..., 0], y = 0.5n

In Paragraaf 6.5.1 hebben we gezien wat het minimum aantal projecties met de waaierstraal is (met bijbehorende hoek) vanuit één ribbe, vanuit twee opeenvolgende ribben, vanuit twee tegenover elkaar staande ribben en vanuit alle ribben. Het aantal projecties proberen we nu nog verder te verkleinen door de afstand van de waaierstraal tot de foto te vergroten. Dit doen we door de x-waarde van het beginpunt van de waaierstraal te variëren van x = −3 tot en met x = 0. E´ en ribbe met een hoek van 14 graden voor de waaierstraal. De hoek voor minimum aantal projecties per waaierstraal is 14 graden voor x = 0. Variëren we het beginpunt van de waaier van x = −3 tot en met x = 0, dan zien we in Figuur 30 dat het minimum aantal projecties 10 blijft. Als we de hoek 13 graden maken (zie Figuur 31) in plaats van 14 graden kunnen we het aantal projecties wel verminderen.

Figuur 30: Minimum aantal projecties = 10, hoek = 14 graden, x = 0.

Figuur 31: Minimum aantal projecties = 9, hoek = 13 graden, x = [−0.4, −0.3]. Linker en boven ribbe met een hoek van 32 graden voor de waaierstraal. De hoek voor minimum aantal projecties per waaierstraal is 32 graden voor x = 0. Variëren we het beginpunt van de waaier van x = −3 tot en met x = 0, dan zien we in Figuur 32 dat het minimum aantal projecties 10 blijft. Het minimum aantal projecties in Figuur 32 is 10 voor x ∈ [−0.1, 0].

36

Figuur 32: Minimum aantal projecties = 10, hoek = 32 graden, x ∈ [−0.1, 0]. Linker en rechter ribbe met een hoek van 29 graden voor de waaierstraal. De hoek voor minimum aantal projecties per waaierstraal is 29 graden voor x = 0. Variëren we het beginpunt van de waaier van x = −3 tot en met x = 0, dan zien we in Figuur 33 dat het minimum aantal projecties 10 blijft. Het minimum aantal projecties in Figuur 33 is 10 voor x ∈ [−0.1, 0].

Figuur 33: Minimum aantal projecties = 10, hoek = 29 graden, x ∈ [−0.1, 0]. Alle ribben met een hoek van 50 graden voor de waaierstraal. De hoek voor minimum aantal projecties per waaierstraal is 50 graden voor x = 0. Variëren we het beginpunt van de waaier van x = −3 tot en met x = 0, dan zien we in Figuur 34 dat het minimum aantal projecties 12 blijft. Het minimum aantal projecties in Figuur 34 is 12 voor x ∈ [−0.7, 0].

Figuur 34: minimum projecties = 12, hoek = 50 graden, x ∈ [−0.7, 0].

37

Conclusie ribbe één ribbe twee opeenvolgende ribben tegenover elkaar liggende ribben vier ribben



x∈ [−0.4, −0.3] [−0.1, 0] [−0.1, 0] [−0.7, 0]

Tabel 4: Conclusie voor waaierstralen vanuit beginpunten met y = 0.5n. Conclusies voor de waaierstraal vanuit beginpunten x ∈ [−3, ..., 0], y = 1 en x ∈ [−3, ..., 0],y ∈ [1, 2, 3] De conclusies voor y = 1, 2, 3 en y = 1 staan in onderstaande tabellen. De uitwerkingen gaan op dezelfde manier als voor y = 0.5n en deze zijn te vinden in de Appendix A. Conclusie ribbe één ribbe twee opeenvolgende ribben tegenover elkaar liggende ribben vier ribben



x∈ [−0.6, −0.5] [−0.6, 0] −0.7 ∪ [−0.6, 0] [−2.8, −1.9]

Tabel 5: Conclusie voor waaierstralen vanuit beginpunten met y = 1, 2, 3.




x∈ [−0.7, −0.6] ∪ 0 −0.4 [−0.4, 0] [−0.4, 0]

Tabel 6: Conclusie voor waaierstralen vanuit beginpunten met y = 1.

38

6.5.3

Conclusie voor de pseudo-inverse met de waaierstraal.

Om een minimum aantal projecties voor een foto van 3 × 3 te krijgen (zie Hoofdstuk 5.1) hebben we in Tabel 7 de projectiestrategiën met bijbehorende hoek en preciese afstand tot de foto voor het beginpunt van de waaierstralen gegeven. Restricties voor de waaierstraal bij een foto van 3 × 3 ribbe y∈ x∈ hoek (graden) één ribbe [1, 2, 3] [0.5, 0.6] 40 één ribbe 1 0 ∪ [0.6, 0.7] 10 twee opeenvolgende ribben 1 0.4 31 één ribbe 1.5 [0.3, 0.4] 13 Tabel 7: Minimum aantal projecties van 9 voor het gebruik van waaierstralen met de 3 verschillende projectiestrategiën. Voor de drie projectiestrategiën krijgen we als conclusie voor een foto van n × n Tabel 8. De coördinaten voor de beginpunten van de waaierstralen kunnen we omrekenen naar het algemene geval voor een foto van n × n. Maar als we dit toepassen op foto’s van n × n dienen deze waarden eerder als richtlijnen dan als kwantitatieve conclusies. Wat we weten is dat het aantal projecties voor een vaste hoek minder wordt als het beginpunt van de waaierstraal verder van de foto wordt gezet. Ook weten we dat de hoek voor de waaierstraal kleiner moet worden genomen als de afmeting van de foto groter wordt om per projectie een pixel toe te kunnen voegen zodat er een projectiematrix A van rang n2 gevormd kan worden. Hierdoor kunnen we wel gerichter zoeken naar een juiste hoek bij een van te voren bepaalde afstand tot de foto. Restricties voor de waaierstraal bij een foto van n × n ribbe y∈ x∈ 1 één ribbe y ∈ [1, 2, ..., n − 1, n] [− 5 n, − 16 n] 7 één ribbe 1 [− 30 n, − 51 n] ∪ 0 4 twee opeenvolgende ribben 1 − 30 n 2 1 één ribbe 0.5n [− 15 n, − 10 n] Tabel 8: Minimum aantal projecties van n × n voor het gebruik van waaierstralen met de 3 verschillende projectiestrategiën.

39

6.6 6.6.1

Waaierstraal met ART Volgorde van de projecties

We onderzoeken met ART of de volgorde van de projecties uit maakt bij het nemen van onafhankelijke projecties. Nemen we achtereenvolgens een projectie van de eerste pixel uit 1 rij en vullen de volgende projectie aan met één pixel en doen dit tot we de hele rij met pixels hebben gehad dan krijgen we voor een foto van 3 × 3 bijvoorbeeld de projectiematrix (zie bijbehorende Figuur 35) :  A=

 0 0 0 0 0 0 1 0 0  0 0 0 0 0 0 1 1 0 . 0 0 0 0 0 0 1 1 1

(5)

Van de projectiematrix A nemen we met ART eerst de eerste rij, dan de tweede rij en dan de derde rij. Hierna beginnen we weer met de eerste rij en gaan zo elke rij van de projectiematrix A af tot we 100 ART iteraties hebben gehad met de drie projecties. Dat de pixelwaarden, waar een projectie doorheen gaat, convergeren naar de originele pixelwaarden is te zien door de foto’s met elkaar te vergelijken. De pixels waar geen projectie doorheen gaat blijven zwart. Dit is te zien in Figuur 35.

Figuur 35: 3 projecties, achtereenvolgens: origineel, 3 ART iteraties, 15 ART iteraties, 99 ART iteraties. In Figuur 36 zien we dat bij genoeg herhaling van de drie projecties (100 ART iteraties delen door 3 projecties), de pixelwaarden waar de projecties doorheen gaan, convergeren naar de originele pixelwaarden. Als er geen projectie door een pixel gaat zal deze dus geen pixelwaarde krijgen en niet worden meegenomen in de foutcorrectie bij toepassing van ART. Dit is in Figuur 36 ook te zien doordat in de grafieken van de grijswaarde per pixel maar 3 grafieken te zien zijn. De andere zes blijven op waarde 0. De fout convergeert niet naar 0, omdat er nog een fout in de 6 pixels zit waar geen projectie doorheen gaat.

Figuur 36: Twaalf keer 3 onafhankelijke projecties is 99 ART iteraties. De projecties gaan door één rij pixels in de foto. De gebroken lijnen zijn de originele pixelwaarden van de foto.

40

Vergelijken we in Figuur 37 de fout van de projecties vanuit links en vanuit rechts, dan zien we dat voor een minimale fout het aantal ART iteraties hetzelfde blijft. Het maakt dus niet uit of we vanuit links of rechts beginnen met de projecties en binnen een set projecties hangt de convergentie dus niet af van de kleur.

Figuur 37: 1e en 2e grafiek: projecties van links naar rechts, 3e en 4e grafiek: projecties van rechts naar links. We vergelijken de minimale fout en de grijswaarden per pixel bij gericht en random genomen projecties. In Figuur 38 is te zien dat voor minimale fout meer ART iteraties nodig zijn bij het nemen van random projecties.

Figuur 38: 3 projecties van links naar rechts, 1e en 2e grafiek: gericht, 3e en 4e grafiek: random. Conclusie: Elke pixel moet in minstens één van de n2 onafhankelijke projecties zitten. Het maakt niet uit in welke volgorde deze groep projecties wordt genomen als het totaal aantal ART iteraties maar niet random wordt genomen. Binnen 1 projectie maakt het niet uit in welke volgorde de pixels worden genomen voor spreiding van het residu over de pixels.

41

6.6.2

Nauwkeurigheid met ART

We onderzoeken wat er gebeurt met het minimum benodigde aantal projecties ten opzichte van de fout bij het vermeerderen van het aantal ART iteraties om een zo nauwkeurig mogelijke reconstructie van de foto te krijgen. We nemen een waaierstraal vanuit vier verschillende x-waarden (0, −0.5, −1 en −2) en y = 0.5n = 1.5. Als we de hoek van de waaierstraal van 0 graden laten toenemen tot en met 90 graden, dan krijgen we de grafieken in onderstaande 3 Figuren (Figuur 39, 40 en 41). Hoe meer ART iteraties (30, 100, 300) er worden gebruikt hoe nauwkeuriger de grijswaarden van de pixels worden berekend. Bij 300 ART iteraties zien we dat het minimum aantal projecties 11 is. Voor x = −1 zien we dat bij 9 projecties de fout ook afneemt. Nemen we 600 ART iteraties voor x = −1, dan zien we in Figuur 42 dat bij 9 projecties de fout ook naar de maximale machineprecisie van 10−16 gaat.

Figuur 39: 30 keer aantal projecties, achtereenvolgens: x = 0, x = −0.5, x = −1, x = −2.



42

Figuur 42: 600 keer aantal projecties, x = −1. Conclusie: Bij minstens n2 projecties, die samen een projectiematrix van rang n2 vormen, convergeert ART naar de exacte oplossing.

43

6.6.3

Benodigde precisie van de fout voor een scherpe foto.

We onderzoeken wat de precisie van de fout moet zijn om een zo scherp mogelijke foto te krijgen. Nemen we 11 projecties met x = 0 en y = 1.5 voor een foto van 3 × 3 om een maximale convergentie van de grijswaarden per pixel te bepalen dan zien we in Figuur 43 dat hiervoor ongeveer 450 ART iteraties (41 keer 11 projecties) nodig zijn. Ook zien we dat er twee dubbele projecties worden genomen aan het begin van de waaier en aan het eind van de waaier. Dit maakt het minimum aantal projecties 11 − 2 = 9 projecties.

Figuur 43: 300 keer 11 projecties met een waaierstraal vanaf punt x = 0 en y = 1.5. Herhalen we de 11 projecties 41 keer (= 450 ART iteraties), dan krijgen we een precisie van de fout van bijna 10−4 (zie Figuur 44).

Figuur 44: Een precisie van de fout in de buurt van 10−4 is voldoende voor een scherpe foto van 3 × 3. In Figuur 45 zien we dat de pixelwaarden voor een foto van 4 × 4 maximale convergentie hebben voor ongeveer 2800 ART iteraties (97 keer 29 projecties). Dit levert een precisie van bijna 10−4 .

Figuur 45: 100 keer 29 projecties met een waaierstraal vanaf punt x = 0 en y = 1.5. Een precisie van de fout in de buurt van 10−4 is voldoende voor een scherpe foto van 4 × 4. Conclusie: Een precisie van de fout in de buurt van 10−4 is voldoende om een scherpe foto te krijgen. De scherpte van de foto hangt af van de precisie van de RGB-waarden van de kleuren die de pixels hebben die ook een precisie van 10−4 hebben. 44

6.6.4

Overeenkomsten tussen de pseudo-inverse en ART.

We onderzoeken de overeenkomsten tussen de pseudo-inverse en het gebruik van ART. 2 Voor een scherpe reconstructie van een foto van 3 × 3 nemen we (x, y) = (− 15 n, 0.5n) met een −4 foutprecisie in de buurt van 10 (zie Paragraaf 6.6.2). In Figuur 46 is te zien dat de fout per hoek met de pseudo-inverse en ART hetzelfde is als voor de waaierstraal hetzelfde (x, y) startpunt en dezelfde hoek is gebruikt.

Figuur 46: midden: Foutberekening per hoek met de pseudo-inverse, rechts: Foutberekening per hoek met ART (1000 keer per per hoek van de waaierstraal). In Figuur 47 en 48 zien we dat de reconstructie voor de pseudo-inverse en ART hetzelfde is als voor de waaierstraal hetzelfde (x, y) startpunt en dezelfde hoek is gebruikt. Bij een precisie van de fout, voor een hoek van 13 graden (9 projecties), in de buurt van 10−4 moeten we ongeveer 200 ART iteraties nemen (Zie Figuur 47). De reconstructie van de foto is scherp genoeg.

Figuur 47: Links: Reconstructie met de pseudo-inverse bij een hoek van 13 graden met waaierstraalprojecties, midden: Foutberekening per hoek met 200 ART iteraties, rechts: Reconstructie met ART bij een hoek van 13 graden met waaierstraalprojecties.

45

Figuur 48: Links: Reconstructie met de pseudo-inverse bij een hoek van 16 graden met waaierstraalprojecties, midden: Foutberekening per hoek met 200 ART iteraties, rechts: Reconstructie met ART bij een hoek van 16 graden met waaierstraalprojecties. Conclusie: Voor de pseudo-inverse en ART geldt dezelfde scherpte voor de reconstructie van een foto als hetzelfde startpunt (x, y) met dezelfde tussenliggende hoek is gebruikt voor de waaierstraal. Met andere woorden: De fout bij reeksontwikkeling met ART convergeert naar de oplossing met de pseudo-inverse van de projectiematrix A.

46

6.6.5

Conclusie voor waaierstraal met ART

De fout bij reeksontwikkeling met ART, convergeert naar de oplossing met de pseudo-inverse van de projectiematrix A, als hetzelfde startpunt (x, y) met dezelfde tussenliggende hoek is gebruikt voor de waaierstraal. Hiervoor moet, bij het gebruik van de waaierstraal om een foto van n × n te reconstrueren, elke pixel in minstens één van de n2 onafhankelijke projecties zitten. Anders zal de pixel geen correctiewaarde van ART krijgen en zal de ontbrekende pixel waarde 0 blijven (= zwart). Binnen één projectie maakt het niet uit in welke volgorde de pixels worden genomen voor spreiding van het residu over de pixels. Het maakt ook niet uit in welke volgorde de groep n2 onafhankelijke projecties, die samen dus een matrix van rang n2 vormen, wordt genomen. De nauwkeurigheid van de grijswaarden van de pixels hangt af van het aantal keren dat de groep onafhankelijke projecties wordt herhaald met ART. De precisie van de fout om een zo scherp mogelijke foto te krijgen hangt hierbij af van de precisie van de RGB-waarden van de grijswaarden. Of het totaal aantal ART iteraties wel of niet random genomen kunnen worden hangt af van de afmeting van de foto die gereconstrueerd wordt. Hoe groter de foto hoe kleiner de kans dat één projectie twee keer na elkaar genomen wordt. En ook wordt de kans groter dat alle projecties even vaak worden genomen. Bij een vierkante foto maakt het niet uit vanuit welke ribbe de projectiestrategie genomen wordt.

47

7

Validatie

7.1

Strategie om de parameters van te voren te kiezen.

We passen de projectiestrategie van de waaierstraal met y = 0.5n toe. Voor de afstand tot de foto 2 nemen we als richtlijn 15 n. De hoek voor een zo scherp mogelijke reconstructie kunnen we alleen voor een foto van 11 × 11, met behulp van de pseudo-inverse, optimaal berekenen. De foto van 100 × 100 en 390 × 390 nemen in matrixvorm te veel ruimte en rekentijd om een optimale hoek te berekenen. Hierdoor hebben we twee verschillende projectiestrategiën: • Een ideale projectiestrategie met de pseudo-inverse, waarbij de hoek exact wordt berekend. • Een practische projectiestrategie met ART, waarbij de hoek wordt geschat. Voor foto’s waarbij geen optimale hoek bepaald kan worden bepalen we eerst de kleinst mogelijke 2 hoek (zonder 0 graden) tussen twee opeenvolgende projecties, vanuit het punt (− 15 n, 0.5n), waarbij één pixel wordt toegevoegd. In Figuur 49 is dit hoek α. Voor een foto van 11 × 11 is deze hoek te klein om aan te geven. Met hoek α als referentie proberen we een steeds kleinere hoek uit om met ART te kijken welke een scherpere reconstructie geeft. Bij de ART iteraties zijn de dubbele projecties ook meegeteld. Het werkelijk aantal projecties kan e soms wel tot 31 deel van het getelde aantal projecties terug gebracht worden.

Figuur 49: Links: α is de kleinste hoek tussen twee opeenvolgende projecties waarbij één pixel wordt toegevoegd. rechts: De kleinste hoek tussen twee opeenvolgende projecties voor een foto van 11 × 11 is duidelijk kleiner dan die voor de foto van 3 × 3 (links).

48

7.2

Projectiestrategie voor een foto van 11x11

We vergelijken de ideale met de practische projectiestrategie voor een aanvaardbare reconstructie. Om een minimum aantal projecties te krijgen voor een foto van 11 × 11 nemen we een waaierstraal 2 vanaf punt (− 15 · 11, 0.5 · 11). In de grafiek van Figuur 50 is te zien dat er een optimale hoek is van 0.29 graden. Omdat de rang bij een hoek van 0.29 graden 121 is geeft de pseudo-inverse van de projectiematrix A een exacte reconstructie en dus ART ook. Bij de practische strategie komen we op een hoek van 0.34 graden. In Tabel 9 is te zien dat bij de hoek van 0.34 graden de rang 119 is en de fout met de pseudo-inverse een precisie heeft van 10−2 . Omdat, bij genoeg iteraties met ART, we dezelfde oplossing krijgen zullen we nooit in de buurt van de benodigde precisie van 10−4 komen. De optimale berekening met behulp van Figuur 50 doet 1 e het dus beter, omdat het resultaat makkelijk af te lezen is. Een verschil in hoek van 10 kan een verschil maken in de fout van 10−12 . De optimale hoek is dus bijna niet te schatten. Dit maakt van het bepalen van de hoek met de practische strategie erg veel werk.

Figuur 50: Er is een minimum aantal projecties bij een hoek van 0.29 graden (172 projecties met een rang van 121).

(x,y) 2 (− 15 · 11, 0.5 · 11)

11×11, pinv(A) hoek projecties 0.34 167 0.29 172

rang 119 121

fout 0.0280 5.93−15

Tabel 9: De pseudo-inverse voor de hoek van 0.34 graden met de practische en de hoek van 0.29 graden met de ideale projectiestrategie. In de tabel hieronder kunnen we zien dat voor een precisie in de buurt van 10−4 het aantal ART iteraties meer dan 500 · 172 moet zijn. Dan hebben we een exacte reconstructie zoals te zien in Figuur 52. Als we op het oog een acceptabele reconstructie kiezen, dan is er na 300 · 172 ART iteraties ook een goede reconstructie. De precisie van 10−4 is dus niet nodig.

2 (− 15

(x,y) · 11, 0.5 · 11)

hoek 0.29 0.29 0.29

11×11, ART projecties ART iteraties 172 100·172 172 300·172 172 500·172

fout 0.0214 0.0041 0.0011

tijd 5 sec 10 sec 11 sec

Tabel 10: 100, 300 en 500 ·172 ART iteraties bij een hoek van 0.29 graden.

49

Figuur 51: Pinv(A), links: origineel, midden: hoek = 0.34 graden, rechts: hoek = 0.29 graden.

Figuur 52: ART iteraties, 0.29 graden, links: 100 · 172, midden: 300 · 172, rechts: 500 · 172.

7.3

Projectiestrategie voor een foto van 100x100

Om een minimum aan projecties te krijgen maken we gebruik van de practische projectiestrategie. We komen op een hoek van 0.005 graden voor een waaierstraal vanuit punt (− 200 15 , 50). Bij een grotere hoek van 0.01 graden zien we in Figuur 53 dat er met 14.979 projecties ook al een duidelijke reconstructie is . Op het oog maakt het niet uit of we 1.000 · 14.979 of 10.000 · 14.979 ART iteraties nemen. Voor een kleinere hoek van 0.001 graden (zie Figuur 54) is de korreligheid tegenover het beginpunt van de waaierstraal bijna verdwenen. Maar het aantal projecties is 10 keer meer (149.781).

(x,y) 2 (− 15 · 100, 0.5 · 100)

hoek 0.01

0.001

100×100, ART projecties ART iteraties 14.979 1·14.979 1.000·14.979 10.000·14.979 149.781 1·149.781 100·149.781 200·149.781

fout 0.305070 0.064629 0.058657 0.306799 0.070845 0.034146

tijd 4 sec 44 min, 55 sec 6 uur, 48 min 41 min, 26 sec 1 uur, 19 min, 53 sec

Tabel 11: ART iteraties bij een hoek van 0.01 graden en een hoek van 0.001 graden.

50

Figuur 53: Hoek = 0.01, ART iteraties = links: 10.000·14.979.

1·14.979, midden:

1.000·14.979, rechts:

Figuur 54: Hoek = 0.001, ART iteraties = links: 1·149.781, midden: 100·149.781, rechts: 200·149.781.

7.4

Projectiestrategie voor een foto van 390 × 390

Om een minimum aan projecties te krijgen maken we gebruik van de practische projectiestrategie. We komen op een hoek van 0.0003 graden voor een waaierstraal vanuit punt (− 780 15 , 195). Deze hoek hebben we voor alle reconstructies toegepast.

2 (− 15

(x,y) · 390, 0.5 · 390)

hoek 0.0003

390×390, ART projecties ART iteraties 500.091 1·500.091 18·500.091

fout 0.2607 0.25520

tijd 20 min, 4 sec 6 uur, 1 minuut, 26 sec

Tabel 12: ART iteratgies bij een hoek van 0.0003 graden.

Figuur 55: Hoek = 0.0003, ART iteraties = links: origineel, midden: 1·500.091, rechts: 18·500.091. 51

Voor de reconstructie van Figuur 56 hebben we, vanwege het vele zwart in de hoeken van de foto, de totale hoek van alle projecties in de waaierstraal van 180 graden verkleind naar 80 graden. Afgezien van de korreligheid is er duidelijk te zien om welke reconstructie het gaat. Er zijn, inclusief dubbele projecties, 266.667 projecties gebruikt. Bij een rang van n2 = 152.100 is dit erg laag. Dit is een goeie reconstructie.

(x,y) 2 (− 15 · 390, 0.5 · 390)

hoek 0.0003

projecties 266.667

390×390, ART ART iteraties 1·266.667 190·266.667

fout 0.258712 0.169281

tijd 10 min, 22 sec 2 dagen, 1 uur, 34 min, 43 sec

Tabel 13: ART iteraties bij een hoek van 0.0003 graden met een totale projectiehoek voor de waaierstraal van 80 graden.

Figuur 56: Hoek = 0.0003, ART iteraties = links: 1·266.667, rechts: 190·266.667.

52

8

Conclusie

Het reconstrueren van een dwarsdoorsnede van een lichaam met een CT-scan wordt ook wel binnen de toegepaste wiskunde een inverse probleem genoemd. Als onderzoeksvraag hebben we: Hoe moeten de projecties worden gekozen om met zo min mogelijk projecties een reconstructie van voldoende kwaliteit te maken? Strategie Om tot een oplossing te komen hebben we gebruik gemaakt van de pseudo-inverse en het ART algoritme. We hebben geen rekening gehouden met ruis. Voor het stelsel Ax = b moet de rang van de projectiematrix A voor een foto van n × n n2 zijn om een unieke oplossing te krijgen. Omdat de projectiematrix A overbepaald en niet vierkant is, maken we gebruik van de pseudo-inverse van A. We hebben drie verschillende manieren gebruikt voor het nemen van projecties. Horizontale projecties, horizontale en verticale projecties en waaierstraalprojecties. We hebben voor de waaierstraal gekozen, omdat we hierbij ook voor grote reconstructies de exacte oplossing kunnen krijgen. De waaierstraal hebben we toegepast bij drie verschillende experimenten. Voor een beginpunt vanuit y = 0.5n, y = 1 en y ∈ [1, ..., n]. Het punt voor de waaierstraal vanuit y = 0.5n bij een foto van 3 × 3 waarbij de fout minimaal is bij een zo groot mogelijke hoek hebben we als referentie gebruikt voor een foto van n × n. Dit is 2 het punt (− 15 n, 0.5n) geworden voor een waaierstraal vanuit één ribbe. 2 Voor de projectiestrategie vanuit het punt (− 15 n, 0.5n) hebben we een ideale en een practische strategie. De ideale strategie maakt gebruik van de pseudo-inverse om de optimale hoek te bepalen en bij de practische strategie nemen we de kleinst mogelijke hoek tussen twee projecties, die één pixel toevoegt. Deze hoek kunnen we verkleinen tot er op het oog met ART een acceptabele reconstructie uit komt. De hoek van de practische strategie is een schatting en geen kwantitatieve conclusie, maar de reconstructies zijn op het oog wel aanvaardbaar. De practische strategie is nodig omdat de berekening van de pseudo-inverses, bij het bepalen van een optimale hoek, te veel tijd neemt (bij een reconstructie van n ≥ 100). Als dit gebeurt passen we in plaats van de ideale strategie, de practische strategie toe. De precisie van de fout van de reconstructie hangt samen met de nauwkeurigheid van de gebruikte RGB waarden van 10−4 . Een precisie van de fout van 10−4 is niet nodig. Als we op het oog een acceptabele reconstructie kiezen komen we op een precisie die zeker 102 kleiner is. Bij het reconstrueren van een foto gebruiken we het ART algoritme omdat deze geschikt is voor grote problemen. De manier waarop we met de pseudo-inverse de juiste positie voor de waaierstraal en de hoek bepalen geeft voor een foto van n × n een rang van n2 en de fout convergeert dus naar 0. Dit is dus een goeie manier voor het bepalen van het minst mogelijk aantal projecties. De afstand van 2 (− 15 n, 0.5n) voor foto’s n × n waarvan n groot (bijvoorbeeld n ≥ 100) is meer een richtlijn en de bijbehorende hoek wordt geschat. Over de coördinaten van de waaierstraal en de bijbehorende hoek kunnen we dus geen kwantitatieve conclusie geven voor een minimale fout.

53

Mogelijke verbeteringen De hoek voor een aanvaardbare reconstructie van grote afmeting wordt geschat en de aanvaardbaarheid bepalen we op het oog. Maar een aanvaardbare reconstructie kunnen we ook bepalen door naar de snelheid van de afname van de fout te kijken per projectie: q 2 n2 X ¯ nieuw ¯ oud x −x i i . n2 i=1

De (x, y) coördinaten van het beginpunt van de waaierstraal en de hoek, zouden we met de pseudoinverse kunnen berekenen voor grote reconstructies, door het matlab programma efficiënter te maken. Hierdoor kan ook de snelheid van het reconstrueren van een foto worden verhoogd. Het gebruik van een krachtiger computer helpt ook. We kunnen projecties maken die allemaal van pixel naar pixel gaan en samen precies n2 onafhankelijke projecties vormen (Zie Figuur 57). Hierbij hebben we geen pseudo-inverse meer nodig en dus ook geen krachtiger computer.

Figuur 57: Projecties van pixel naar pixel. Omdat we hierbij al een vast aantal projecties hebben die samen rang n2 vormen kunnen we de projecties bij een reconstructie van grote afmeting ook random nemen. Hoe groter de foto hoe minder kans dat twee dezelfde projecties na elkaar worden genomen en hoe groter de kans dat alle projecties even vaak voorkomen. Omdat grote reconstructies op het oog acceptabel zijn bij gebruik van minder dan deze n2 projecties pleit dit voor het gebruik van random projecties. Convergentie van de fout naar de originele pixelwaarden gebeurt, bij een van te voren bepaalde volgorde van de projecties, lineair. Dit is te zien in Hoofdstuk 6.6.3. Omdat de convergentie voorspelbaar is kan hier wellicht gebruik van worden gemaakt bijvoorbeeld met Richardson extrapolatie.

54

A A.1

Waaierstraal met de pseudo-inverse Waaierstraal vanuit beginpunten x = 0, y ∈ [1, 2, 3]

We onderzoeken wat het minimum aantal projecties is bij waaierstralen vanuit beginpunten (0, 1), (0, 2) en (0, 3). Deze waaierstralen nemen we vanuit één ribbe, vanuit twee opeenvolgende ribben, vanuit twee tegenover elkaar staande ribben en vanuit alle ribben. In Figuur 58, 59 en 60 zien we dat het minimum aantal projecties met een minimale fout 10 is. Het maakt niet uit of er één of twee ribben worden genomen. Als we alle vier de ribben nemen (zie Figuur 61) hebben we een minimum van 20 projecties. E´ en ribbe

Figuur 58: Minimum aantal projecties = 10, hoek = 40 graden. In Figuur 58 zien we dat er bij meerdere hoeken 10 projecties door de foto gaan. De 10 projecties kunnen per hoek van elkaar verschillen doordat de projecties niet altijd door dezelfde pixels gaan. Daarom zal niet bij elke waaierstraal (met bijbehorende hoek) met 10 projecties de fout naar de maximale machineprecisie 10−16 gaan. Linker en boven ribbe

Figuur 59: minimum projecties = 10, hoek = 72 graden.

55

Linker en rechter ribbe

Figuur 60: Minimum aantal projecties = 10, hoek = 70 graden. Alle ribben


In de grafieken met het aantal projecties hebben we de dubbele projecties niet meegeteld. Als er meerdere hoeken mogelijk zijn hebben we de grootste hoek genomen om zo min mogelijk dubbele projecties door een foto te krijgen. Conclusie ribbe één ribbe twee opeenvolgende ribben tegenover elkaar liggende ribben vier ribben



Tabel 14: Conclusie voor waaierstralen vanuit beginpunten (0, 1), (0, 2) en (0, 3).

We kunnen het aantal projecties verkleinen door de x-coördinaat van de waaierstraal verder van de foto te plaatsen. Dit doen we in Paragraaf A.2

56

A.2

Waaierstraal vanuit beginpunten x ∈ [−3, ..., 0], y ∈ [1, 2, 3]

In Paragraaf A hebben we gezien wat het minimum aantal projecties met de waaierstraal is (met bijbehorende hoek) vanuit één ribbe, vanuit twee opeenvolgende ribben, vanuit twee tegenover elkaar staande ribben en vanuit alle ribben. Het aantal projecties proberen we nu nog verder te verkleinen door de afstand van de waaierstraal tot de foto te vergroten. E´ en ribbe met een hoek van 40 graden voor de waaierstraal. De hoek voor minimum aantal projecties per waaierstraal is 40 graden voor x = 0. Variëren we het beginpunt van de waaier van x = −3 tot en met x = 0, dan zien we in Figuur 62 dat het minimum aantal projecties van 10 naar 9 gaat voor x ∈ [−0.6, −0.5].

Figuur 62: Minimum aantal projecties = 9, x ∈ [−0.6, −0.5]. Linker en boven ribbe met een hoek van 72 graden voor de waaierstraal. De hoek voor minimum aantal projecties per waaierstraal is 72 graden voor x = 0. Als we hiervoor de afstand tot de foto variëren krijgen we Figuur 63. Het aantal projecties kunnen we niet kleiner maken dan 10.

Figuur 63: Minimum aantal projecties = 10, x ∈ [−0.6, 0].

57

Linker en rechter ribbe met een hoek van 70 graden voor de waaierstraal. De hoek voor minimum aantal projecties per waaierstraal is 70 graden voor x = 0. Als we hiervoor de afstand tot de foto variëren krijgen we Figuur 64. Het aantal projecties wordt niet kleiner dan 10.

Figuur 64: Minimum aantal projecties = 10, x ∈ −0.7 ∪ [−0.6, 0]. Alle ribben met een hoek van 71 graden voor de waaierstraal. De hoek voor minimum aantal projecties per waaierstraal is 71 graden voor x = 0. Variëren we het beginpunt van de waaier van x = −3 tot en met x = 0, dan zien we in Figuur 65 dat het minimum aantal projecties van 20 naar 12 gaat voor x ∈ [−2.8, −1.9].

Figuur 65: Minimum aantal projecties = 12, x ∈ [−2.8, −1.9].




x∈ [−0.6, −0.5] [−0.6, 0] −0.7 ∪ [−0.6, 0] [−2.8, −1.9]

Tabel 15: Conclusie voor waaierstralen vanuit beginpunten met y = 1, 2, 3.

58

A.3

Waaierstraal vanuit beginpunt x = 0, y = 1

We onderzoeken wat het minimum aantal projecties is bij het nemen van een waaierstraal vanuit punt (0, 1). Deze waaierstraal nemen we vanuit één ribbe, vanuit twee opeenvolgende ribben, vanuit twee tegenover elkaar staande ribben en vanuit alle ribben. Bij één ribbe hebben we een minimum van 9 projecties. Bij twee ribben hebben we een minimum van 10 projecties (zie Figuur 67 en 68) en bij vier ribben is het minimum 12 projecties (zie Figuur 69). E´ en ribbe


Linker en boven ribbe


59

Linker en rechter ribbe


Alle ribben

Figuur 69: Minimum aantal projecties = 12, hoek = 50 graden. In de grafieken met het aantal projecties hebben we de dubbele projecties niet meegeteld. Als er meerdere hoeken mogelijk zijn hebben we de grootste hoek genomen om zo min mogelijk dubbele projecties door een foto te krijgen. Conclusie ribbe één ribbe twee opeenvolgende ribben tegenover elkaar liggende ribben vier ribben



Tabel 16: Conclusie voor waaierstralen vanuit beginpunt (0, 1).

We kunnen het aantal projecties verkleinen door de x-coördinaat van de waaierstraal verder van de foto te plaatsen.

60

A.4

Waaierstraal vanuit beginpunten x ∈ [−3, ..., 0], y = 1

In Paragraaf A.3 hebben we gezien wat het minimum aantal projecties met de waaierstraal is (met bijbehorende hoek) vanuit één ribbe, vanuit twee opeenvolgende ribben, vanuit twee tegenover elkaar staande ribben en vanuit alle ribben. Het aantal projecties proberen we nu nog verder te verkleinen door de afstand van de waaierstraal tot de foto te vergroten.

E´ en ribbe met een hoek van 10 graden voor de waaierstraal. De hoek voor minimum aantal projecties per waaierstraal is 10 graden voor x = 0. Variëren we het beginpunt van de waaier van x = −3 tot en met x = 0, dan zien we in Figuur ?? dat het minimum aantal projecties 9 blijft voor x ∈ [−0.7, −0.6] ∪ 0. In Hoofdstuk 5.1 hadden we al gezien dat het aantal projecties niet kleiner kan worden dan 9.

Figuur 70: Minimum aantal projecties = 9, hoek = 10 graden, x ∈ [−0.7, −0.6] ∪ 0. Linker en boven ribbe met een hoek van 31 graden voor de waaierstraal. De hoek voor minimum aantal projecties per waaierstraal is 31 graden voor x = 0. Variëren we het beginpunt van de waaier van x = −3 tot en met x = 0, dan zien we in Figuur 71 dat het minimum aantal projecties van 10 naar 9 gaat voor x = −0.4.

Figuur 71: Minimum aantal projecties = 9, hoek = 31 graden, x = −0.4.

61

Linker en rechter ribbe met een hoek van 31 graden voor de waaierstraal. De hoek voor minimum aantal projecties per waaierstraal is 31 graden voor x = 0. Variëren we het beginpunt van de waaier van x = −3 tot en met x = 0, dan zien we in Figuur 72 dat het minimum aantal projecties 10 blijft. De afstand tot de foto is 0 tot en met 0.4.

Figuur 72: Minimum aantal projecties = 10, hoek = 31 graden, x ∈ [−0.4, 0]. Alle ribben met een hoek van 50 graden voor de waaierstraal. De hoek voor minimum aantal projecties per waaierstraal is 50 graden voor x = 0. Variëren we het beginpunt van de waaier van x = −3 tot en met x = 0, dan zien we in Figuur 73 dat het minimum aantal projecties 12 blijft. De afstand tot de foto is 0 tot en met 0.4.

Figuur 73: Minimum aantal projecties = 12, hoek = 50 graden, x ∈ [−0.4, 0].




x∈ [−0.7, −0.6] ∪ 0 −0.4 [−0.4, 0] [−0.4, 0]

Tabel 17: Conclusie voor waaierstralen vanuit beginpunten met y = 1.

62

B B.1

Matlab code Grafieken met de pseudo-inverse: 1

%Pseudo-inverse: log van fout en aantal projecties t.o.v. de hoek. %Het beginpunt (x,y) voor de projecties, de ribbe en de afmeting van de %reconstructie kunnen we aangeven. clc; clear all; close all; %Afmeting van een reconstructie van mxm in pixels. m = input(’pixels foto= ’); %We gebruiken elke keer dezelfde vierkante foto. ff= rand(’state’); f= rand([m m]); rand(’state’,ff); p= reshape(f,m*m,1); %1= linker ribbe; 2= linker en boven ribbe; 3= linker, boven en rechter %ribbe; 4= linker, boven, rechter en onder ribbe. mmm= input(’ribbe= ’); %Beginpunt projecties x-coordinaat. xbegin= input(’beginpunt x-coordinaat= ’); %Grootte van 1 stap. h= 0.1; %Aantal stappen n= (m-xbegin)/ h; step= h*[0:1:n]; %X-coordinaat van een projectie. x= xbegin + step; xr= ceil(x); tic %Berekening van de fout per hoek for projhoek= 90 :-1 :1 mm= round(180/projhoek); xx= zeros(m*m,1); aantalprojecties= 0; %Kiezen van de ribben.

63

for K= 1: mmm %Beginpunt projecties y-coordinaat. Dit hebben we ook voor %ybegin= 1 en ybegin= 0.5*m gedaan. for ybegin= 1: 1: m %Maken van projectiematrix met de gegeven projecties for k= 1: mm %rc van een projectie in radialen gradenhoek= -90+(k-1)*projhoek; radhoek= 2*pi*gradenhoek/360; rc= tan(radhoek); %Y-coordinaat van een projectie. y= ybegin + rc * step; yr= ceil(y); %alle 0-en weglaten uit de vector a a= zeros(m); %De (x,y) coordinaten van de projectie omzetten naar %indexnummers voor de foto infoto= find ((xr<=m & xr>0) & (yr<=m & yr>0)); lengte= length(infoto); switch K case {1} %Linker ribbe for i= 2: lengte a(m+1-yr(infoto(i)),xr(infoto(i)))=1; end case{2} %Boven ribbe for i= 2: lengte a(xr(infoto(i)),yr(infoto(i)))=1; end case{3} %Rechter ribbe for i= 2: lengte a(yr(infoto(i)),m+1-xr(infoto(i)))=1; end otherwise %Onder ribbe for i= 1: lengte a(m+1-xr(infoto(i)),m+1-yr(infoto(i)))=1; end end

64

%ART algoritme aa= reshape(a,1,m*m); b= aa * p; a1=aa*aa’; if a1 > 0 aantalprojecties = aantalprojecties + 1; aa1(aantalprojecties,:) = aa; if (aantalprojecties > 1) & ((aa1(aantalprojecties,:) == aa1(aantalprojecties-1,:))) aantalprojecties = aantalprojecties - 1; aa1(aantalprojecties,:) = aa; end end end end end MM=aa1*p; M=pinv(aa1)*MM; hoek(projhoek)=projhoek; fout(projhoek)= ones(1,m*m) * abs(M-p) / (m*m); ap(projhoek)=aantalprojecties; end toc figure(6) imagesc(f); colormap(gray); caxis([0 1]) title(’Origineel’); figure(9) semilogy(hoek,fout); xlabel(’hoek’); ylabel(’log van fout’); figure(11) plot(hoek,ap,’o’); xlabel(’hoek’); ylabel(’aantal projecties’);

65

B.2

Grafieken met de pseudo-inverse: 2

%Pseudo-inverse: log van fout en aantal projecties t.o.v. beginpunt waaierstraal op x-as. %De hoek, de ribbe, de afmeting van de reconstructie en de y-coordinaat kunnen %we aangeven. clc; clear all; close all; %Afmeting van een reconstructie van mxm in pixels. m = input(’pixels foto= ’); %We gebruiken elke keer dezelfde vierkante foto. ff= rand(’state’); f= rand([m m]); rand(’state’,ff); p= reshape(f,m*m,1); %1= linker ribbe; 2= linker en boven ribbe; 3= linker, boven en rechter %ribbe; 4= linker, boven, rechter en onder ribbe. mmm= input(’ribbe= ’); %Grootte van 1 stap. h= 0.1; %rc van een projectie in graden projhoek= input(’hoek tussen projecties= ’); mm= round(180/projhoek); tic %Berekening van de fout per afstand (x-coordinaat) tot de foto. for xb = 10:1:40 xbegin= (xb-40)/10; %Aantal stappen n= (m-xbegin)/ h; step= h*[0:1:n]; %X-coordinaat van een projectie. x= xbegin + step; xr= ceil(x); xx= zeros(m*m,1); aantalprojecties= 0; %Kiezen van de ribben. for K= 1: mmm

66

%Beginpunt projecties y-coordinaat. Dit hebben we ook voor %ybegin= 1 en ybegin= 0.5*m gedaan. for ybegin= 1: 1: m %Maken van projectiematrix met de gegeven projecties for k= 1: mm %rc van een projectie in radialen gradenhoek= -90+(k-1)*projhoek; radhoek= 2*pi*gradenhoek/360; rc= tan(radhoek); %Y-coordinaat van een projectie. y= ybegin + rc * step; yr= ceil(y); %alle 0-en weglaten uit de vector a a= zeros(m); %De (x,y) coordinaten van de projectie omzetten naar %indexnummers voor de foto infoto= find ((xr<=m & xr>0) & (yr<=m & yr>0)); lengte= length(infoto); switch K case {1} %Linker ribbe for i= 2: lengte a(m+1-yr(infoto(i)),xr(infoto(i)))=1; end case{2} %Boven ribbe for i= 2: lengte a(xr(infoto(i)),yr(infoto(i)))=1; end case{3} %Rechter ribbe for i= 2: lengte a(yr(infoto(i)),m+1-xr(infoto(i)))=1; end otherwise %Onder ribbe for i= 1: lengte a(m+1-xr(infoto(i)),m+1-yr(infoto(i)))=1; end end %ART algoritme 67

aa= reshape(a,1,m*m); b= aa * p; a1=aa*aa’; if a1 > 0 aantalprojecties = aantalprojecties + 1; aa1(aantalprojecties,:) = aa; if (aantalprojecties > 1) & ((aa1(aantalprojecties,:) == aa1(aantalprojecties-1,:))) aantalprojecties = aantalprojecties - 1; aa1(aantalprojecties,:) = aa; end end end end end MM=aa1*p; M=pinv(aa1)*MM; xbb(xb)=xbegin; fout(xb)= ones(1,m*m) * abs(M-p) / (m*m); ap(xb)=aantalprojecties; end toc figure(6) imagesc(f); colormap(gray); caxis([0 1]) title(’Origineel’); figure(10) semilogy(xbb,fout1); xlabel(’beginpunt waaierstraal op x-as’); ylabel(’log van fout’); figure(12) plot(xbb,ap,’o’); xlabel(’beginpunt waaierstraal op x-as’); ylabel(’aantal projecties’);

68

B.3

Grafieken met ART: 1

%ART: log van fout t.o.v. het aantal projecties (en t.o.v de hoek). %Het beginpunt (x,y), de ribbe, de afmeting van de reconstructie en het aantal %herhalingen van de projecties kunnen we aangeven. clc; clear all; close all; %Afmeting van een reconstructie van mxm in pixels. m = input(’pixels foto= ’); %We gebruiken elke keer dezelfde vierkante foto. ff= rand(’state’); f= rand([m m]); rand(’state’,ff); p= reshape(f,1,m*m); mmmm= input(’aantal herhalingen van de projecties= ’); %1= linker ribbe; 2= linker en boven ribbe; 3= linker, boven en rechter %ribbe; 4= linker, boven, rechter en onder ribbe. mmm= input(’ribbe= ’); %Beginpunt projecties x-coordinaat. xbegin= input(’beginpunt x-coordinaat= ’); %Grootte van 1 stap. h= 0.1; %Aantal stappen n= (m-xbegin)/ h; step= h*[0:1:n]; %X-coordinaat van een projectie. x= xbegin + step; xr= ceil(x); tic %Berekening van de fout per hoek for projhoek= 90 :-1 :1 mm= round(180/projhoek); xx= zeros(m*m,1); ARTiteraties= 0; %Aantal herhalingen van de projecties vanaf de gekozen ribben.

69

for T= 1: mmmm %Kiezen van de ribben. for K= 1: mmm %Beginpunt projecties y-coordinaat. Dit hebben we ook voor %ybegin= 1 en ybegin= 0.5*m gedaan. for ybegin= 1: 1: m for k= 1: mm %rc van een projectie in radialen gradenhoek= -90+(k-1)*projhoek; radhoek= 2*pi*gradenhoek/360; rc= tan(radhoek); %Y-coordinaat van een projectie. y= ybegin + rc * step; yr= ceil(y); %alle 0-en weglaten uit de vector a a= sparse(m,m); %De (x,y) coordinaten van de projectie omzetten naar %indexnummers voor de foto infoto= find ((xr<=m & xr>0) & (yr<=m & yr>0)); lengte= length(infoto); switch K case {1} %Linker ribbe for i= 2: lengte a(m+1-yr(infoto(i)),xr(infoto(i)))=1; end case{2} %Boven ribbe for i= 2: lengte a(xr(infoto(i)),yr(infoto(i)))=1; end case{3} %Rechter ribbe for i= 2: lengte a(yr(infoto(i)),m+1-xr(infoto(i)))=1; end otherwise %Onder ribbe for i= 1: lengte a(m+1-xr(infoto(i)),m+1-yr(infoto(i)))=1; end 70

end %ART algoritme aa= reshape(a,m*m,1); b= p * aa; a1=aa’*aa; if a1 > 0 ARTiteraties = ARTiteraties + 1; r = b-aa’*xx; xx = xx + (r/a1) * aa; end end end end ap(projhoek)= ARTiteraties/mmmm; hoek(projhoek)= projhoek; fout(projhoek)= ones(1,m*m) * abs(p’-xx)/ (m*m); end end toc figure(6) imagesc(f); colormap(gray); caxis([0 1]) title(’Origineel’); figure(8) semilogy(ap,fout); xlabel(’aantal projecties’); ylabel(’log van fout’); figure(9) semilogy(hoek,fout); xlabel(’hoek’); ylabel(’log van fout’);

71

B.4

Grafieken met ART: 2

%ART: log van fout en grijswaarde per pixel t.o.v het aantal ART iteraties. %Het beginpunt (x,y), de ribbe, de afmeting van de reconstructie, de hoek %en het aantal herhalingen van de projecties kunnen we aangeven. clc; clear all; close all; %Afmeting van een reconstructie van mxm in pixels. m = input(’pixels foto= ’); %We gebruiken elke keer dezelfde vierkante foto. ff= rand(’state’); f= rand([m m]); rand(’state’,ff); p= reshape(f,1,m*m); mmmm= input(’aantal herhalingen van de projecties= ’); %1= linker ribbe; 2= linker en boven ribbe; 3= linker, boven en rechter %ribbe; 4= linker, boven, rechter en onder ribbe. mmm= input(’ribbe= ’); %Beginpunt projecties x-coordinaat. xbegin= input(’beginpunt x-coordinaat= ’); %Grootte van 1 stap. h= 0.1; %Aantal stappen n= (m-xbegin)/ h; step= h*[0:1:n]; %X-coordinaat van een projectie. x= xbegin + step; xr= ceil(x); %rc van een projectie in graden projhoek= input(’hoek tussen projecties= ’); mm= round(180/projhoek); xx= zeros(m*m,1); ARTiteraties= 0; tic %Aantal herhalingen van de projecties vanaf de gekozen ribben.

72

for T= 1: mmmm %Kiezen van de ribben. for K= 1: mmm %Beginpunt projecties y-coordinaat. Dit hebben we ook voor %ybegin= 1 en ybegin= 0.5*m gedaan. for ybegin= 1: 1: m %Maken van een projectievector for k= 1: mm %rc van een projectie in radialen gradenhoek= -90+(k-1)*projhoek; radhoek= 2*pi*gradenhoek/360; rc= tan(radhoek); %Y-coordinaat van een projectie. y= ybegin + rc * step; yr= ceil(y); %alle 0-en weglaten uit de projectie a= sparse(m,m); %De (x,y) coordinaten van de projectie omzetten naar %indexnummers voor de foto infoto= find ((xr<=m & xr>0) & (yr<=m & yr>0)); lengte= length(infoto); switch K case {1} %Linker ribbe for i= 2: lengte a(m+1-yr(infoto(i)),xr(infoto(i)))=1; end case{2} %Boven ribbe for i= 2: lengte a(xr(infoto(i)),yr(infoto(i)))=1; end case{3} %Rechter ribbe for i= 2: lengte a(yr(infoto(i)),m+1-xr(infoto(i)))=1; end otherwise %Onder ribbe for i= 1: lengte a(m+1-xr(infoto(i)),m+1-yr(infoto(i)))=1; 73

end end %ART algoritme aa= reshape(a,m*m,1); b= p * aa; a1=aa’*aa; if a1 > 0 ARTiteraties = ARTiteraties + 1; r = b-aa’*xx; xx = xx + (r/a1) * aa; end %De fout per pixel voor een reconstructie van 3x3. %Als ybegin = 0.5*m of ybegin = 1 dan moeten we van de index %hieronder m = 1 en ybegin = 1 maken! nrproj((T-1)*mm*m*mmm+(K-1)*m*mm+(ybegin-1)*mm+k)=ARTiteraties; fout((T-1)*mm*m*mmm+(K-1)*m*mm+(ybegin-1)*mm+k) = ones(1,m*m) * abs(p’-xx) / (m*m); fout1((T-1)*mm*m*mmm+(K-1)*m*mm+(ybegin-1)*mm+k)=xx(1); fout2((T-1)*mm*m*mmm+(K-1)*m*mm+(ybegin-1)*mm+k)=xx(2); fout3((T-1)*mm*m*mmm+(K-1)*m*mm+(ybegin-1)*mm+k)=xx(3); fout4((T-1)*mm*m*mmm+(K-1)*m*mm+(ybegin-1)*mm+k)=xx(4); fout5((T-1)*mm*m*mmm+(K-1)*m*mm+(ybegin-1)*mm+k)=xx(5); fout6((T-1)*mm*m*mmm+(K-1)*m*mm+(ybegin-1)*mm+k)=xx(6); fout7((T-1)*mm*m*mmm+(K-1)*m*mm+(ybegin-1)*mm+k)=xx(7); fout8((T-1)*mm*m*mmm+(K-1)*m*mm+(ybegin-1)*mm+k)=xx(8); fout9((T-1)*mm*m*mmm+(K-1)*m*mm+(ybegin-1)*mm+k)=xx(9); orig1((T-1)*mm*m*mmm+(K-1)*m*mm+(ybegin-1)*mm+k)=f(1); orig2((T-1)*mm*m*mmm+(K-1)*m*mm+(ybegin-1)*mm+k)=f(2); orig3((T-1)*mm*m*mmm+(K-1)*m*mm+(ybegin-1)*mm+k)=f(3); orig4((T-1)*mm*m*mmm+(K-1)*m*mm+(ybegin-1)*mm+k)=f(4); orig5((T-1)*mm*m*mmm+(K-1)*m*mm+(ybegin-1)*mm+k)=f(5); orig6((T-1)*mm*m*mmm+(K-1)*m*mm+(ybegin-1)*mm+k)=f(6); orig7((T-1)*mm*m*mmm+(K-1)*m*mm+(ybegin-1)*mm+k)=f(7); orig8((T-1)*mm*m*mmm+(K-1)*m*mm+(ybegin-1)*mm+k)=f(8); orig9((T-1)*mm*m*mmm+(K-1)*m*mm+(ybegin-1)*mm+k)=f(9); end end end end toc figure(6) imagesc(f); 74

colormap(gray); caxis([0 1]) title(’Origineel’); X = reshape(xx,m,m); figure(7) imagesc(X); colormap(gray); caxis([0 1]) title(’Reconstructie’); figure(3) semilogy(nrproj,fout); xlabel(’ART iteraties’); ylabel(’log van fout’); figure(2) plot(nrproj,fout1,nrproj,fout2,nrproj,fout3,nrproj,fout4,nrproj,fout5,nrproj, fout6,nrproj,fout7,nrproj,fout8,nrproj,fout9,nrproj,orig1,’-.’, nrproj,orig2,’-.’,nrproj,orig3,’-.’,nrproj,orig4,’-.’,nrproj, orig5,’-.’,nrproj,orig6,’-.’,nrproj,orig7,’-.’,nrproj,orig8,’-.’, nrproj,orig9,’-.’); xlabel(’ART iteraties’); ylabel(’grijswaarde per pixel’); fout= ones(1,m*m) * abs(p’-xx) / (m*m) ARTiteraties= ARTiteraties aantalprojecties= ARTiteraties/mmmm

75

B.5

Reconstructie met de pseudo-inverse

%Pseudo-inverse: reconstructie van hersenscan met beginpunt (-2*m/15,0.5*m) %vanuit 1 ribbe. De hoek kunnen we aangeven. clc; clear all; close all; %We maken een reconstructie van hersenscan390x390.png. We kunnen dit ook %doen voor gezicht.png. IU=imread(’hersenscan390x390.png’); f=im2double(IU); S=size(f); m=S(1); p= reshape(f,m*m,1); %Beginpunt projecties x-coordinaat. xbegin= -2*m/15; %Beginpunt projecties y-coordinaat. ybegin= 0.5*m; %Grootte van 1 stap. h= 0.1; %Aantal stappen n= (m-xbegin)/ h; step= h*[0:1:n]; %X-coordinaat van een projectie. x= xbegin + step; xr= ceil(x); %rc van een projectie in graden projhoek= input(’hoek tussen projecties= ’); mm= round(180/projhoek); xx= zeros(m*m,1); aantalprojecties= 0; tic %Maken van projectiematrix met de gegeven projecties for k= 1: mm %rc van een projectie in radialen gradenhoek= -90+(k-1)*projhoek; radhoek= 2*pi*gradenhoek/360;

76

rc= tan(radhoek); %Y-coordinaat van een projectie. y= ybegin + rc * step; yr= ceil(y); %alle 0-en weglaten uit de vector a a= zeros(m); %De (x,y) coordinaten van de projectie omzetten naar %indexnummers voor de foto infoto= find ((xr<=m & xr>0) & (yr<=m & yr>0)); lengte= length(infoto); %Linker ribbe for i= 2: lengte a(m+1-yr(infoto(i)),xr(infoto(i)))=1; end %ART algoritme aa= reshape(a,1,m*m); b= aa * p; a1=aa*aa’; if a1 > 0 aantalprojecties = aantalprojecties + 1; aa1(aantalprojecties,:) = aa; if (aantalprojecties > 1) & ((aa1(aantalprojecties,:) == aa1(aantalprojecties-1,:))) aantalprojecties = aantalprojecties - 1; aa1(aantalprojecties,:) = aa; end end figure(1) plot (x,y); hold on; axis ([0 m 0 m]) title(’Projecties’); end toc figure(6) imagesc(f); colormap(gray); caxis([0 1]) title(’Origineel’);

77

MM=aa1*p; M=pinv(aa1)*MM; X=reshape(M,m,m); figure(7) imagesc(X); colormap(gray); caxis([0 1]) title(’Reconstructie’); fout= ones(1,m*m) * abs(M-p) / (m*m) ap=aantalprojecties

78

B.6

Reconstructie met ART

%ART: reconstructie van hersenscan met beginpunt (-2*m/15,0.5*m) %vanuit 1 ribbe. De hoek en het aantal herhalingen van de projecties %kunnen we aangeven. clc; clear all; close all; %We maken een reconstructie van hersenscan390x390.png. We kunnen dit ook %doen voor gezicht.png. IU=imread(’hersenscan390x390.png’); f=im2double(IU); S=size(f); m=S(1); p= reshape(f,1,m*m); mmmm= input(’aantal herhalingen van de projecties= ’); %Beginpunt projecties x-coordinaat. xbegin= -2*m/15; %Beginpunt projecties y-coordinaat. ybegin= 0.5*m; %Grootte van 1 stap. h= 0.1; %Aantal stappen n= (m-xbegin)/ h; step= h*[0:1:n]; %X-coordinaat van een projectie. x= xbegin + step; xr= ceil(x); %rc van een projectie in graden projhoek= input(’hoek tussen projecties= ’); mm= round(180/projhoek); xx= zeros(m*m,1); ARTiteraties= 0; tic %Aantal herhalingen van de projecties vanaf de gekozen ribben. for T= 1: mmmm

79

%Maken van projectiematrix met de gegeven projecties for k= 1: mm %rc van een projectie in radialen gradenhoek= -90+(k-1)*projhoek; radhoek= 2*pi*gradenhoek/360; rc= tan(radhoek); %Y-coordinaat van een projectie. y= ybegin + rc * step; yr= ceil(y); %alle 0-en weglaten uit de vector a a= zeros(m); %De (x,y) coordinaten van de projectie omzetten naar %indexnummers voor de foto infoto= find ((xr<=m & xr>0) & (yr<=m & yr>0)); lengte= length(infoto); %Linker ribbe for i= 2: lengte a(m+1-yr(infoto(i)),xr(infoto(i)))=1; end %ART algoritme aa= reshape(a,m*m,1); b= p * aa; a1=aa’*aa; if a1 > 0 ARTiteraties = ARTiteraties + 1; r = b-aa’*xx; xx = xx + (r/a1) * aa; end end end toc figure(6) imagesc(f); colormap(gray); caxis([0 1]) title(’Origineel’); X = reshape(xx,m,m); figure(7) 80

imagesc(X); colormap(gray); caxis([0 1]) title(’Reconstructie’); fout= ones(1,m*m) * abs(p’-xx) / (m*m) ARTiteraties= ARTiteraties aantalprojecties= ARTiteraties/mmmm

81

Referenties [1] Herman, Gabor T., Fundamentals of Computerized Tomography: Image Reconstruction from Projections, Springer-Verlag London, second edition, 2009. [2] Katenburg, K.J., Sals, S., Sijbers, J., K¨ ubel, C., Midgley, P.A., Hernandez, J.C., Kaiser, U., Encina, E.R., Coronado, E.A., Tendeloo, G. van, 3D imaging of nanomaterials by discrete tomography, Ultramicroscopy, Article accepted 20 January 2009. [3] Buzug, Thorsten, Computed tomography: from photon statistics to modern cone-beam CT (Hst 6 algebraic and statistical reconstruction methods), Springer Berlin, 2008. [4] Chan, R., Thiagalingam, A., dAvila, A., Ho, I., Reddy, V., Manzke, R., Intraprocedural Imaging of Left Atrial and Pulmonary Vein Anatomy for Atrial Fibrillation Ablation, Philips Research North America, Briarcliff Manor, NY, USA, 2007. [5] Beckmann, E.C., President’s conference paper: CT scanning the early days, Lanmark, Beaconsfield, Bucks, UK, The British journal of radiology, 79(2006). [6] Herman, G.T., Image Reconstruction From Projections, Medical Image Processing Group, Department of Radiology, University of Pennsylvania, Philadelphia, PA 19104, USA, 1995. [7] Laurent, C., Efficient Implementation of Parallel Image Reconstruction Algorithms for 3D X-ray Tomography, 1992. [8] Strang, Gilbert, Linear Algebra and its Applications, Academic Press, Inc. New York, second edition, 1980. [9] Hounsfield, Godfrey N., Computed medical imaging, the medical system departmetnt of central research laboratories EMI, London, England, 8 december, 1979. [10] Krishnamurthy, E.V., Mahadeva Rao, T., Subramanian, K., Prabhu, S.S., short note: Reconstruction of Objects from Their Projections Using Generalized Inverses, Department of Applied Mathematics, Indian Institute of Science, Bangalore 560012, India, Communicated by A. Rosenfeld, Received 9 April 1974. [11] Gordon, Richard, Bender, Robert, Herman, Gabor T., Algebra¨ıc Reconstruction Techniques (ART) for Three-dimensional Electron Microscopy and X-ray Photography, Center for Theoretical Biology and Department of Computer Science, State University of New York at Buffalo, Amherst, N.Y. 14266, U.S.A., Received 12 August 1970.

82

CT-scan met zo min mogelijk projecties

Recommend Documents