Lineaire Algebra: Speciaal Onderwerp 1 Tomografie in Flatland Jan Brandts 31 augustus 2009
1
Tomografie in Flatland
Flatland is een tweedimensionale platte wereld waarin bewoners leven die driehoekig of vierkant of zelfs meerhoekig zijn. Als ultieme veelhoek schuifelen er zelfs cirkels in rond, zoals je hebt kunnen zien in Flatland The Movie. Deze tweedimensionale wereld leent zich prima voor het bestuderen van een vereenvoudigd model voor tomografie, al was het alleen maar omdat we van boven kunnen meekijken naar het inwendige van de bewoners. Dit is eigenlijk een beetje een merkwaardig idee, vooral omdat de Flatlanders zelf alleen maar elkaars buitenkant zien, oftewel, elkaars omtrek. Ze gruwelen bij het idee dat wij drie-dimensionale wezens zomaar hun inwendige zou kunnen bekijken. Je wilt waarschijnlijk zelf ook niet weten wat een vierdimensionaal wezen zou zien als hij op je “neer”kijkt vanuit de vierde dimensie! Toch zou dit voor medische toepassingen natuurlijk fantastisch zijn. Immers, er zou nooit meer iemand opengesneden hoeven te worden, noch om te kijken wat er aan de hand is, noch om bepaalde reparaties uit te voeren. Ook buiten de medische hoek zijn natuurlijk toepassingen te vinden, zoals bagage-controles op lucht- en andere havens. De perfecte bankoverval wordt ineens een re¨ele optie. Niets, maar dan ook helemaal niets, kan nog verborgen blijven.
Afbeelding 1. Doorsnedes van het inwendige verkregen middels medische tomografie. Helaas kunnen we geen gebruik maken van de vierde dimensie, en misschien is dat maar goed ook. Toch kunnen we middels wat lineaire algebra een heel eind komen om ons doel om een object van binnen te bekijken zonder het open te snijden, waar te maken. Achterliggend idee is dat je misschien wel heel vakkundig je handen zodanig kunt vouwen dat de schaduw ervan op de muur net een konijn is, maar als men met meerdere lampen vanuit verschillende richtingen 1
op die handen gaat schijnen en alle schaduwen worden bekeken, zal men toch concluderen dat het je handen zijn die worden beschenen, en zeker geen konijn.
1.1
Modelaannamen over anatomie en rookgedrag van Flatlanders
We gaan ervan uit dat ook Flatlanders elkaar in hun donkere en barbaarse verleden al veelvuldig hebben opengesneden, zodat de anatomie van de diverse soorten Flatlanders goed bekend is. Daarnaast nemen we aan dat behoorlijk wat Flatlanders de dubieuze gewoonte hebben om (platte) sigaretten te roken. Zoals ook in andere dimensies, resulteert dit in verhoogde concentraties van roet in de diverse organen, en in de longen natuurlijk het meest. Tot slot nemen we aan dat de roetverzadiging binnen ´e´en en hetzelfde orgaan van een Flatlander gelijk is. Dus, ieder orgaan heeft zijn eigen, homogene roetverzadiging. Om deze concentraties te meten bestaat er in Flatland een soort straling met de volgende eigenschappen. Als je deze straling door een stukje weefsel van een Flatlander stuurt waar geen roet in zit, dan gaat het er ongehinderd doorheen. Is het stukje weefsel echter volledig verzadigd van roet, dan gaat er helemaal niets doorheen. Ligt de roetverzadiging echter tussen de nul en volledig in, dan hangt de hoeveelheid straling die door het weefsel heen gaat af van de lengte van de weg die de straling hierdoorheen moet afleggen.
1.2
Een wiskundig model
We gieten dit als volgt in een wiskundig model. Eerste stap is dat we de roetverzadiging aan zullen duiden met een getal, bijvoorbeeld C, dat tussen de nul en de ´e´en ligt. Hierbij spreken we het volgende af. • C = 0 is volledige verzadiging: zwart van het roet, er gaat niets doorheen; • C = 1 staat voor helemaal geen roet, alles gaat erdoorheen; • 0 < C < 1 staat voor tussenliggende waarden van roetverzadiging. Voor het gemak zullen we het hebben over zwart als C = 0, over wit als C = 1, en grijs als 0 < C < 1. Laat nu ` de afstand zijn waarover de straling door een stukje weefsel gaat. Je kan je voorstellen dat hoe langer deze afstand, hoe minder straling er overblijft, tenzij C = 0. We beschrijven dit met de volgende formule, waarbij A de hoeveelheid straling is bij aanvang, en M de hoeveelheid die overblijft na over een afstand ` door weefsel met roetverzadiging C is gegaan. We zullen ervan uit gaan dat dan M = C ` A.
(1)
Deze formule heeft een aantal realistische eigenschappen. Stel bijvoorbeeld dat je een hoeveelheid A aan straling over een afstand m door een stukje tissue met grijstint C stuurt, dan blijft er C m A. Stuur je deze overgebleven hoeveelheid vervolgens over een afstand k door nog een stukje weefsel met dezelfde waarde van C, dan blijft er C k (C m A) over. Omdat geldt dat C k (C m A) = C k+m A
(2)
is hetzelfde resultaat als wanneer je formule (1) gebruikt om te zien wat je overhoudt als je diezelfde straling in ´e´en keer over een afstand k + m door dit weefsel stuurt. E´en en ander is ge¨ıllustreerd in Afbeelding 2. 2
Afbeelding 2. Schematische voorstelling van sommige aspecten van het model Andere realistische eigenschappen van formule (1) zijn • Als C = 0 dan C ` = 0 voor alle positieve `. Dun of dik zwart maakt niet uit. • Als C = 0 en ` = 0, dan is C ` = 00 = 1. Zwart met dikte nul is wit! • Als C = 1 dan C ` = 1 voor alle positieve `, ´en voor ` = 0. Wit is wit, voor alle diktes.
1.3
Een uitgewerkt voorbeeld
Zoals gezegd is de anatomie van Flatlanders bekend. Hiermee bedoelen we dat we precies weten waar de diverse organen zich in het lichaam bevinden, en welke ruimte ze innemen. Hieronder staan van een driehoekige en een vierkante Flatlander de indeling in drie en vier organen getekend.
Afbeelding 3. Anatomie van een driehoek en een vierkant. We nemen de driehoek als voorbeeld. Deze is onderverdeeld in drie driehoeken, de drie organen, die alledrie een bepaalde grijstint hebben. Het is mogelijk om deze grijstinten te berekenen door op drie manieren een straal met intensiteit A = 1 door deze Flatlander te sturen, en te meten wat ervan overblijft. Dit is weergegeven in Afbeelding 4.
3
Afbeelding 4. Het sturen van drie stralen door een driehoekige Flatlander. De informatie die deze procedure oplevert is de volgende. Er gaat een straal met intensiteit A = 1 over een afstand van k door het orgaan met grijstint C die vervolgens over een afstand k door het orgaan met grijstint D gaat. Stel dat je met de daarvoor benodigde apparatuur meet dat van die A = 1 er M1 aan intensiteit overblijft. Je weet dan op grond van het model dat blijkbaar C k D k = M1 . (3) Ondanks dat k bekend is, is dit niet voldoende om C en D afzonderlijk uit te rekenen. Echter, de andere twee stralen geven eveneens aanleiding tot meetgegevens M2 en M3 en we weten dat Dk E k = M2 en C k E k = M3 . (4) De laatste drie vergelijkingen bij elkaar zijn wel voldoende om te bepalen wat C, D en E zijn. Om dat in te zien is het nodig om je de volgende rekenregels voor logaritmes te herinneren, log(ab) = log(a) + log(b) en
log(ap ) = p log(a), voor alle a, b > 0 en p ∈ R.
(5)
Door nu van de bovenstaande drie vergelijkingen links en rechts de logaritme te nemen, en de rekenregels correct toe te passen, volgt dat = log(M1 ) k log(C) + k log(D) k log(D) + k log(E) = log(M2 ) (6) k log(C) + k log(E) = log(M3 ) De onbekenden in deze vergelijkingen zijn log(C), log(D), en log(E). De meetgegevens en dus ook hun logaritmes, en de afstanden k zijn bekend. Als we voor het gemak x1 = log(C), x2 = log(D), x3 = log(E) schrijven, staat in (6) niets anders dan = log(M1 ) kx1 + kx2 kx2 + kx3 = log(M2 ) (7) kx1 + kx3 = log(M3 ) Dit is een stelsel van drie lineaire vergelijkingen voor drie onbekenden x1 , x2 en x3 . We kunnen x1 , x2 en x3 hieruit oplossen en daarna C, D en E bepalen, in termen van M1 , M2 , M3 en k.
4
1.4
Discussie en interpretatie van de lineaire algebra
Dit uitgewerkte voorbeeld laat zien hoe een concreet probleem resulteert in lineaire algebra. In dit probleem staan de drie onbekende roetverzadigingen C, D, E centraal. Een manier om deze drie getallen te weten te komen, is ervoor te zorgen dat je er lineaire vergelijkingen voor afleidt. Iedere straal die je door de Flatlander stuurt zorgt voor zo’n vergelijking. Opmerking 1: De logaritme van een meetwaarde is het inproduct tussen de vector van onbekende logaritmes van roetconcentraties in, en de vector van afgelegde afstanden door de drie organen. Bijvoorbeeld, de eerste vergelijking in (6) laat zich herschrijven als * k log(C) + log(M1 ) = k log(C) + k log(D) + 0 log(E) = k , log(D) . (8) 0 log(E) Ook bij andere modellen komt het voor dat je inproducten kunt uitrekenen met een onbekende vector, om op die manier de vergelijkingen op te stellen die de onbekende ontmaskeren. Opmerking 2: Het liefst zou je natuurlijk het inproduct nemen van bijvoorbeeld de vector met getallen 1, 0, 0 en de vector met onbekenden, omdat dat de direct oplosbare vergelijking C = M1 zou opleveren. Helaas staat de geometrie (vorm) van de driehoekige Flatlander dat niet toe: het is onmogelijk om een rechte lijn door hem heen te tekenen die door precies ´e´en van zijn organen gaat. Bij de vierkante Flatlander uit Afbeelding 3 gaat dit overigens wel! Ieder orgaan moet minstens ´e´en straal door zich heen krijgen, omdat anders de corresponderende onbekende roetverzadiging in geen enkele vergelijking voorkomt, en je hem dus ook niet kunt uitrekenen. Echter, niet ieder drietal stralen is geschikt: Opmerking 3: Als je in Afbeelding 4 de straal door C en D een halve centimeter naar links of rechts zou verplaatsen evenwijdig aan de getekende straal, resulteert dit in vergelijking (??) vermenigvuldigd met een scalair, respectievelijk kleiner en groter dan ´e´en. Deze drie vergelijkingen corresponderen met ´e´en en hetzelfde vlak in R3 en zijn dus niet voldoende om C, D en E te bepalen. Immers, de doorsnede van die drie vlakken is het vlak zelf. Gerelateerd hieraan zien we nu ook een alternatief voor de procedure gevolgd in Sectie ??. Opmerking 4: Stuur behalve de eerste straal door C en D er een tweede doorheen die niet evenwijdig is aan de eerste, en die door C een lengte p aflegt ongelijk aan de lengte q door D. Als M de meetwaarde is bij deze straal, vinden we dat C p Dq = M en dus p log(C) + q log(D) = log(M ).
(9)
Door de bovenste vergelijking uit (6), behorend bij de eerste straal, hier p/k maal vanaf te trekken kunnen we log(D) en dus ook log(C) berekenen, al voordat we een straal door E hebben gestuurd! Zulke strategie¨en zijn in de tomografische praktijk helaas vrijwel onmogelijk te bepalen omdat ze voor ieder object weer anders zijn, zoals blijkt uit de verschillen tussen de twee Flatlanders uit Afbeelding 3, en al helemaal niet door het apparaat dat de straling uitzendt zelf.
1.5
Tot slot: de praktijk in medische tomografie
In de praktijk is het meestal niet zo dat er van tevoren al bekend is in welke gebieden de roetverzadiging hetzelfde is, zoals we hier aannamen voor de organen van de Flatlanders. Om 5
dit probleem te omzeilen, wordt het object in denkbeeldige kleine kubusjes verdeeld en nemen we aan dat binnen ieder van die kubusjes het concentratie homogeen is. Dit is in de realiteit natuurlijk niet helemaal waar, dus we introduceren hier een modelleerfout. Door de kubusjes maar klein genoeg te kiezen zal deze fout kleiner en kleiner worden. Als je een menselijk hoofd van grofweg 30 bij 30 bij 30 centimeter echter in kubusjes van 1 bij 1 bij 1 millimeter verdeelt, heb je wel te maken met zo’n 27 miljoen uit te rekenen onbekenden, waarvoor je dus ook minstens 27 miljoen vergelijkingen nodig hebt.
Afbeelding 5. Opdelen in vakjes, een straal door ieder vakje. Dit verklaart waarom er wereldwijd nog steeds veel onderzoek plaatsvindt in het snel oplossen van grote stelsels vergelijkingen. Daarnaast zijn er methoden die genoegen nemen met een goede benadering van de oplossing, en hierdoor veel rekentijd uitsparen. Ook in Nederland wordt hier veel aan gewerkt. De bijdrage H.A. van der Vorst (1992). Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems. SIAM J. Sci. Statist. Comput. 13(2):631–644. van Prof. dr. Henk van der Vorst van de Universiteit Utrecht (inmiddels met emeritaat) over dit onderwerp is bekroond als het wereldwijd vaakst geciteerde wiskundige onderzoeksartikel gemeten over het decennium 1990-2000! Merk op dat de snelheid van het oplossen van stelsels vergelijkingen vooral heel belangrijk wordt als je real-time wilt meekijken naar bepaalde processen binnen iemands schedel.
Afbeelding 6. Real-time studie van de effecten van mobiele telefoons op de hersenen.
6
Je zou dan de oplossing van zo’n stelsel met 27 miljoen onbekenden in een fractie van een seconde uit willen rekenen. Ondanks dat computers steeds sterker worden en de oplosmethoden steeds ingenieuzer, zijn dit soort uitdagingen nog steeds niet doeltreffend afgehandeld. Opmerking: Als je een computer leert om Gauss-eliminatie toe te passen op stelsels vergelijkingen, neem de benodigde rekentijd toe met de derde macht van het aantal onbekenden. Dus duurt het acht keer zolang om een stelsel met twee keer zoveel onbekenden op te lossen. Ga na: om een stelsel van 27 miljoen onbekenden in ´e´en seconde op te kunnen lossen zou je een stelsel van twee onbekenden in een onmogelijk korte tijd op moeten kunnen lossen. Natuurlijk is er ook aandacht voor het reduceren van het aantal benodigde stralen door met betere modellen te werken. Immers, niet iedere patient verdraagt tientallen miljoenen verschillend gerichte doses straling. Duidelijk is dat tomografie vele interessante multi-disciplinaire vragen opwerpt, waarvan de wiskundige vragen zeker niet het onbelangrijkst zijn!
2
Opgaven en vooruitblik
Met onderstaande opgaven oefen je je vaardigheden in het opstellen van vergelijkingen voor een gegeven aantal onbekenden. Ze laten zien dat er verschillende vergelijkingen mogelijk zijn die allemaal dezelfde oplossing hebben. Sommige vergelijkingen laten zich echter veel gemakkelijker oplossen dan andere, en het kan enorm lonen hier goed over na te denken. Opgave 1: Een vierkante Flatlander met afwijkende anatomie In Afbeelding 3 zagen we de anatomie van een vierkante Flatlander. Omdat het mogelijk is om door ieder orgaan precies ´e´en straal te sturen, zijn de vergelijkingen voor de vier onbekenden wel heel gemakkelijk op te lossen. Ze zijn namelijk ontkoppeld. Echter, er zijn ook vierkante Flatlanders met de volgende anatomie:
Afbeelding 7. Vier metingen aan een vierkante Flatlander met afwijkende anatomie. Neem voor het gemak aan dat de zijden van het vierkant lengte twee hebben. (a) Stel de vergelijkingen op die resulteren uit de stralen zoals aangegeven in Afbeelding 7. (b) Los deze vergelijkingen op: druk iedere Cj uit in termen van M1 t/m M4 . (c) Geef een alternatief voor de vier getekende stralen gebaseerd op Sectie 1.4, Opmerking 4.
7
Opgave 2: De n-hoekige Flatlander met n driehoekige organen Bekijk nu een n-hoekige Flatlander met n gelijke hoeken, en met n driehoekige organen. De drie hoekpunten van zo’n orgaan zijn gelijk aan het middelpunt van de Flatlander, en twee op´e´envolgende hoekpunten van de Flatlander zelf. De driehoek uit Afbeelding 3 en het vierkant uit Afbeelding 7 zijn hier voorbeelden van.
Afbeelding 8. n-hoekige Flatlanders met een eenvoudige, regelmatige anatomie. Ontwikkel voor algemene waarden van n een strategie om de onbekende roetgehaltes in deze organen zo gemakkelijk mogelijk op te kunnen lossen. Opgave 3: De 5-hoekige Flatlander met 11 organen Onderstaande vijfhoekige Flatlander heeft een wat ingewikkeldere anatomie, hij bestaat uit maar liefst elf organen.
Afbeelding 9. Vijfhoekige Flatlander met een ingewikkeldere anatomie Stel, je bent aleen ge¨ınteresseerd in het roetgehalte van de inwendige vijfhoek. Bedenk een geschikte strategie om deze zo eenvoudig mogelijk te bepalen.
2.1
Aansluitend projectwerk bij Simuleren en Programmeren
Bij het vak Simuleren en Programmeren in het voorjaarssemester is het mogelijk om dit onderwerp te kiezen voor projectwerk. Hierbij zal je dan in een groep van drie of vier personen strategie¨en ontwikkelen voor het algemene geval van een in n × n hokjes opgedeeld vierkant, waarbij n erg groot kan worden. De berekening van bijvoorbeeld de afstand die een straal door een hokje aflegt (en dus van de co¨effici¨enten van de lineaire vergelijkingen ) en de berekening van de oplossing, alsmede de visualisering ervan in zwart-wit of zelfs in kleur (rood-groenblauw tinten) laat je dan zoveel mogelijk over aan de computer middels een zelfgeschreven computerprogramma. Hierbij leer je onder andere over de zgn. kleinste kwadraten methode. 8
