Denne side har til formål at give et indblik i den matematiske teori bag teknikken brugt i CT-skannere, der i dag bruges på hospitaler rundt omkring i hele verden. Den er tiltænkt som en del af undervisningmaterialet til SRP projektet Hvordan ser vi det skjulte, der udbydes fra DTU Compute. Materialet her på siden er et sammenkog hentet fra forskellige resources og forsøgt genformidlet på et niveau forståeligt for en gymnasieelev på tredje år. Nederst på siden vil være links til yderligere resourcer, men vær obmærksom på, at de kan være på engelsk og højt niveau.
Lidt Historie
Johan Radon
Det er nok de færreste, der går og tænker over at for bare 50 år siden, eksisterede CT-skannere slet ikke, og mulighederne for at se ind i kroppen uden at åbne op, var ganske begrænsede. Man kunne tage røntgenbilleder, men den slags billeder giver slet ikke den samme detaljegrad som CT. I dag bruges røntgen primært til at observere knoglestrukturer og sygdom i blødt væv i kroppen.
Røntgenstråler har været kendt siden slutningen af 1800-tallet, hvor de blev opdaget og navngivet efter den tyske fysiker Willhelm Röntgen, der var den første til at systematisk studere strålerne.
Allerede for omkring 100 år siden i 1917 lagde den østrigske matematiker Johann Radon det matematiske grundlag for CT-skanning, som vi kender det i dag. Han studerede projektioner, det at tage et objekt og afbillede i en lavere dimension, og mulighederne for at rekonstruere objektet, hvis man kender nok afbilledeninger af det. For et simpelt eksempel forestil dig en figur i LEGO, og hvor mange billeder fra forskellige vinkler du ville have brug for, hvis du skulle bygge den igen kun ud fra billederne.
Der skulle dog gå lang tid før CT-skannerne kom. Først over 50 år senere i 1971, blev den første praktiske brugbare CT-skanner markedsført. Det var Allan McLeod Cormack og Godfrey Hounsfields fælles indsats op gennem 50'erne og 60'erne, der gjorde maskinen til virkelighed og de modtog for deres arbejde Nobelprisen i Medicin i 1979.
Sidenhen er CT-skannere blevet langt mere sofistikerede, som du også kan se herunder. Hounsfields prototype og den første serie af CT-skannere skannede langs parallelle linier, mens nyere modeller skanner langs en vifte af linier udspringende fra det samme punkt. Det har en række tekniske fordele, som vi ikke vil gå i detaljer med her. Den teori vi senere vil gennemgå beskæftiger sig med metoden med parallelle linier, men kan generalisers til vifte situationen. Det kunne måske være et spændende projekt i sig selv.
Nu vil vi begynde at se nærmere på den matematik, der ligger bag CT-skanning, som Johann Radon udledte tilbage i 1917. For at arbejde med dette får vi brug for en række værktøjer. Noget vil måske virke lidt svært, men tager vi det et skridt af gangen, går det forhåbentlig alligevel.
Før vi går i dybden med matematikken, så lad os lige snakke lidt om, hvad der fysisk foregår når vi skanner et object. For det første har vi et afgrænset objekt, som vi ønsker at skanne. Typisk vil vi placere objektet i et koordinatsystem og beskrive det med en funktion $ f(x,y) $ (Se Slide 1 herunder). Løst formulere beskriver funktionen $ f $, hvor tæt vores objekt er forskellige steder. Stråling passerer f.eks. ikke lige nemt gennem knogle- og vævmateriale.
Dernæst placerer vi en strålingskilde på en ene side af objektet og en detektor på den anden side. Strålingens bane gennem objektet danner en lige linie fra kilden til detektoren (Slide 2). Ved at bevæge detektoren og kilden på tværs af liniens retning opnår vi en hel skare af parallelle linier gennem objektet. Målingen langs hver af disse linier noteres og danner en profil (Slide 3 og 4).
Ved at måle fra en masse forskellige vinkler opnår vi forskellige profiler, som vi kan prøve at konstruere billedet ud fra.
Det her vil formentlig være relativt simpelt for de flestes, men vi vil alligevel gå lidt systematisk til værks for at komme rundt om, hvad vi får brug for. Lad os først tale kort om forskellen på funktioner, ligninger og identiteter. Rent teknisk vil vi nok sige at funktioner er identiteter, så lad os vende tilbage til dem om lidt.
Så godt som alle udtryk der involverer et lighedstegn ($ = $) er enten en ligning eller en identitet. Den basale forskel er, at ligniger løser vi (typisk ved at finde en ubekendt $ x $), mens identiter er de værktøjer, vi bruger til at løse ligninger. Identiteter er udsagn der altid er sande, mens ligninger kun er opfyldt for bestemte indsatte værdier.
Eksempel
Er følgende en ligning eller en identitet? $$[ \cos^2 x + \sin^2 x = 1 ]$$ Man kunne måske forledes til at tro, at vi kan løse udtrykket og finde $ x $, men det viser sig, at udtrykket er sandt for alle reelle tal $ x $. Udtrykket er altså en identitet.
Identiteter kommer ofte med en række præmisser som skal være opfyldt for, at den gælder. Pythagoras sætning $$[a^2 + b^2 = c^2 ]$$ er fx en identitet, men holder kun for retvinklede trekanter hvor kateterne har længderne $ a $ og $ b $ og hypotenusen længden $ c $.
Et problem kan ofte bruge en identitet til at opstille ligningen, der ved løsning giver løsningen på problemet. Hvis du er givet en retvinklet trekant med hypotenuse $ 5 $ og katete $ 4 $, kan du opstille ligningen ved brug af Pythagoras sætning $$[ a^2 + 4^2 = 5^2, ]$$ hvor $ a $ er vores ubekendte katete længde. Svaret her er $ a = 3 $, men hvis man ikke kender identiteten i Pythagoras sætning er det ikke så simpelt et problem at løse.
Så hvordan passer identiteter og ligninger ind med funktioner? En funktion er grundlæggende en identitet vi selv definerer. Et simpelt eksempel kunne være $ f(x) = ax + b $. Vi vælger at $ f(x) $ skal være lig $ ax + b $ for alle (eller måske bare nogle, i så fald specificerer vi hvilke,) værdier af $ x $. Det er en identitet, men der er ikke langt til en ligning. Vi kan spørge om for hvilke værdier af $ x $, det gælder, at $ f(x) = c $? Så skal vi bruge identiteten vi valgte for $ f(x) $. $$[ f(x) = c ]$$ $$[ ax + b = c ]$$ $$[ x = \frac{c-b}{a} ]$$ så længde $ a \neq 0 $.
En vigtig detalje ved funktioner er at en værdi $ x $ aldrig kan give anledning til mere end en værdi $ f(x) $. Du tænker måske, "Hør, hov?! Hvad med funktionen $ f(x) = \sqrt{x} $? Vi ved da, at kvadratroden har både en postiv og en negativ løsning!" Det er sandt, men set på netop den måde er kvadratroden ikke en funktion. Når vi taler om kvadratroden som en funktion begrænser vi os altid til kun den ene af løsningerne (typisk den positive).
Men nok om $ f(x) $, denne er kun en funktion i een variabel. Vi vil have flere, vi vil have to variable: $ g(x,y) $. De fleste principper for funktioner af en variabel falder også naturligt for funktioner af to eller flere variable. Hvis du kender til vektorer, kan du faktisk tænke på $ g(x,y) $ som en funktion af kun en variabel, nemlig en vektor $ \mathbf{w} = (x,y) $: $$[ g(\mathbf{w}) = g(x,y). ]$$ Meget af det man ser for funktioner af en variabel som fx differentiation og kædereglen herfra kan se overvældende ud når man ser der for funktioner med flere variable, men i vektornotationen ligner de udtrykkene vi kender. Men nok om vektorerne.
Et værktøj vi ofte bruger til bedre at forstå funktioner som fx $ f(x) $ er visualisering. Vi tegner grafen for vores funktion i et $ xy $-koordinatsystem, hvor $ f(x) $ afsættes op af $ y $-aksen og $ x $ naturligt ud af $ x $-aksen. For funktioner som $ g(x,y) $ der er i to variable, kan vi gøre noget lignende, men vi har brug for en akse mere, så vi har et $ xyz $-koordinatsystem. Vi laver den naturlige kobling og afsætte $ x $ og $ y $ ud af deres respektive akser, $ x $-aksen og $ y $-aksen, og afsætter således værdien $ g(x,y) $ op af $ z $-aksen. I stedet for en linie eller kurve er grafen for $ g(x,y) $ en flade.
Funktionen i vores skanningsproblem vil beskrive absorbansen i objektet vi skanner. Vi indlægger simpelthen et koordinatsystem over vores objekt, og så beskriver værdien $ g(x,y) $ hvor meget materialet absorberer/dæmper en stråle i netop det punkt. På grund af de meget bratte grænser for et fysisk objekt vil vores funktion $ g(x,y) $ være ikke kontinuert/sammenhængende, se fx billedet af Shepp-Logan fantomet.
Hvis du gerne vil lære endnu mere om funktioner med flere variable kan du læse i eNoterne for det introducerende matematik kursus her på DTU. De kan findes frit tilgængelige her: eNote 15: Funktioner af to variable.
Polære Koordinater
Polære koordinater er et rimelig simpelt koncept. Vi har allerede i folkeskolen stiftet bekendtskab med de cartesiske koordinater. De cartesiske koordinater er simpelthen bare $ xy $-koordinaterne som vi kender dem. Det cartesiske koordinatsæt $ (3,4) $ er 3 ud af $ x $-aksen og 4 op af $ y $-aksen. Men hvad er så polære koordinater? Ideen er ganske simpel. Hvis du står et sted $ A $ og noget spørger dig, hvor stedet $ B $ er, vil du formentlig oftest svare ved at angive i hvilken retning $ B $ er, og hvor langt væk det er. Polære koordinater er det samme. Her står vi i origo ($(0,0)$), og et punkt angives med en retning (vinkel med $ x $-aksen) og afstand.
I formler får vi følgende relation: Hvis $ (x,y) $ er et punkt angivet ved cartesiske koordinater, da er punktet polære koordinater $ (r, \theta) $ givet ved $$[ r = \sqrt{ x^2 + y^2 } ]$$ $$[ \theta = \tan^{-1}\left(\frac{y}{x}\right) + \chi_{x< 0}\pi ]$$ hvor $ \chi_{x<0} $ er 0, hvis $ x $ er positiv og 1 hvis $ x $ er negativ. Omvendt kan vi også gå den anden vej, hvilket er en lille smule nemmere: $$[ x = r\cos\theta ]$$ $$[ y = r\sin\theta. ]$$
Så hvornår bruger vi polære koordinater? Ofte når vi skal arbejde med cirkulære områder er de smarte at bruge. Men som vi skal se lige om lidt, kan de også være yderst brugbare når vi skal beskrive linier i vores koordinatsystem.
Vi lærer om planens ligning i gymnasiet:
$$[ a(x-x_0) + b(y-y_0) + c(z-z_0) = 0, ]$$
hvor $ \mathbf{n} = (a,b,c) $ er normalvektor til planen. På lignende vis har du nok stiftet bekendtskab med den simplere liniens ligning:
$$[ a(x-x_0) + b(y-y_0) = 0, ]$$
hvor $ \mathbf{n} = (a,b) $ er en normalvektor til linien. Vi kan således beskrive en linie, hvis vi kender en normalvektor $ \mathbf{n} =(a,b) $ og et punkt $ P = (x_0,y_0) $. Dette er en måde at beskrive linien på, men man kan gøre det på andre måder også.
Vi kan istedet for at kigge på punkter overveje 1) hvad er afstanden $ \rho $ fra linien til origo? og 2) i hvilken vinkel $ \theta $ ligger den korteste afstand til linien? Faktisk de to størrelse $ (\rho, \theta) $ nok til at udfylde liniens ligning. Vi sætter simpelthen $ P = (\rho\cos\theta,\rho\sin\theta) $ og $ \mathbf{n} = (\cos\theta,\sin\theta) $.
Noget særlig smart ved at beskrive vores linier på denne måde er, at det nemt beskriver alle de linier strålerne i vores skanner følger. For en bestemt vinkel $ \theta $ kan vi få alle parallelle linier bare ved at variere på $ \rho $.
Siden vi gerne vil kunne regne på hvor meget dæmpning en stråle har langs en linie, vil det være praktisk at kunne beskrive punkterne på linien. Til dette benytter vi en parameterfremstilling af for linien. En parameterfremstilling er i princippet bare to funktioner $ x(t) $ og $ y(t) $, der beskriver en bevægelse langs henholdsvis $ x $-aksen og $ y $-aksen. Parameterfremstilling at en ret linie er simpelt, vi beskriver først et punkt på linien, og derpå liniens retning.
$$[ \mathbf{r}(t) = \begin{bmatrix} x(t) \\ y(t) \end{bmatrix} = \begin{bmatrix} \rho\cos\theta \\ \rho\sin\theta \end{bmatrix} + \begin{bmatrix} -t\sin\theta \\ t\cos\theta \end{bmatrix}. ]$$
Læg mærke til hvordan første led bare er stedvektoren til vores punkt $ P $, og andet led er stedvektoren (normaliseret til længden 1) hattet og så ganget med $ t $. Ethvert punkt $ P_0 $ på vores linie kan nu beskrives som $ \mathbf{r}(t_0) = (x(t_0),y(t_0)) $ for et passende valg af $ t_0 $.
Fra vores parameterfremstilling kan vi også nemt opstille liniens ligning for vores linie. Der gælder om prikproduktet at $ a \bullet b $ præcis er længden af vektor $ a $ projekteret vinkelret ned på vektor $ b $ skaleret med længden af $ b $.
For en stedvektor til et punkt $ (x,y) $ på vores linie ser vi ret tydeligt at den vinkelrette projektering af $ (x,y) $ ned på vores stedvektor $ P $ har længden $ \rho $. Ergo må enhver stedvektor til et punkt $ (x,y) $ på vores linie prikket med vores stedvektor $ P $ være $ \rho \cdot \|P\| = \rho^2 $. Vi har altså
$$[ \begin{bmatrix} x \\ y \end{bmatrix} \bullet\begin{bmatrix} \rho\cos\theta \\ \rho\sin\theta \end{bmatrix} = x\rho\cos\theta + y\rho\sin\theta = \rho^2 ]$$
og dermed liniens ligning
$$[ x\cos\theta + y\sin\theta = \rho. ]$$
På gymnasiet lærer vi integralregning for funktioner af en variabel. For en funktion $ f(x) $ er integralet af $ f(x) $ stamfunktionen $ F(x) $, der opfylder at $ F'(x) = f(x) $. For bestemte integraler over et interval $ I = [a,b] $ har integralet en anden fortolkning som arealet mellem $ f(x) $ og $ x $-aksen, der ligger inden for $ I $. Vi vil i det følgende afsnit se hvordan vi kan integrere funktioner af to variable langs linier, på næsten præcis samme måde som vi ellers kender det.
Så lad os se på en funktion af to variable $ g(x,y) $ og en linie, $ \ell $, i $ xy $-planen givet ved $ x\cos\theta + y\sin\theta = \rho $. Vi skriver nu integralet af $ g(x,y) $ langs $ \ell $ som
$$[ \int_\ell g(x,y)\,ds. ]$$
Det ser måske lidt vildt ud, men lad os tage det et skridt af gangen. Vi tager udgangspunkt i parameterfremstillingen $ \mathbf{r}(t) $ i sidste afsnit. Fra den har vi følgende to ligninger: $ x = \rho\cos\theta - t\sin\theta $ og $ y = \rho\sin\theta + t\cos\theta $. Vi kan indsætte disse ind i vores integrale herover:
$$[ \int_\ell g(\rho\cos\theta - t\sin\theta, \rho\sin\theta + t\cos\theta)\,ds. ]$$
Det måske stadig lidt vildt ud, men husk at $ \rho $ og $ \theta $ begge bare er nogle tal der beskriver vores linie, så når vi har valgt linien er de konstante. Vi har nu sikret os at alle punkterne $(x,y) $ vi sætter ind i vores funktion $ g $ ligger på linien, vi skal bare variere på $ t $. Hvis vi skriver
$$[ h(t) = g(\rho\cos\theta - t\sin\theta, \rho\sin\theta + t\cos\theta) ]$$
er $ h(t) $ jo bare en helt normal funktion af een variabel.
I det "normale" integrale er $ dx $ på en måde længden det infinitisimale liniestykke vi bevæger os langs under integrationen, $ ds $ har samme rolle her. Længden $ s $ af et stykke af vores linie kan beskrives som $ s = t\|(-\sin\theta,\cos\theta)\| $. Vi har således
$$[ \frac{ds}{dt} = \left\|\begin{bmatrix} -\sin\theta \\ \cos\theta\end{bmatrix}\right\| ]$$
og dermed $ ds = \|(-\sin\theta,\cos\theta)\|\cdot dt = dt $. Så vores integrale bliver altså
$$[ \int_{t_0}^{t_1} g(\rho\cos\theta - t\sin\theta, \rho\sin\theta + t\cos\theta)\,dt = \int_{t_0}^{t_1} h(t)\,dt, ]$$
hvor $ t_0 $ og $ t_1 $ bestemmer punkterne på linien vi integrerer imellem. Når vi har bestemt $ h(t) $, kan vi altså bestemme vores integrale på præcis samme manér, som vi plejer.
Læg mærke til hvordan længden af vores retningsvektor kom ind i bestemmelsen af $ ds $. Havde vores retningsvækter haft længden 2 havde vi haft $ ds = 2dt $. Mere generelt har vi at $ ds = \|\mathbf{r}'(t)\|dt $, hvilket vi har brug for, hvis vi vil bruge andre kurver end rette linier. Helt generelt kan vi for en vilkårlig kurve $ k $ parametriseret med $ \mathbf{r}_k(t) = (x(t),y(t)) $ bestemme
$$[ \int_k g(x,y) \,ds = \int_{t_0}^{t_1} g(\mathbf{r}_k(t))\,\|\mathbf{r}_k'(t)\|dt = \int_{t_0}^{t_1} g(x(t),y(t))\,\|(x'(t),y'(t))\|dt. ]$$
Eksempel
Lad os prøve at regne et eksempel. Givet er funktionen $ g(x,y) = (1-x^2)(1-y^2) $ i området $ x,y \in [0,1] $ og linien $ \ell $ være givet ved $ y\sin(\pi/3) + x\cos(\pi/3) = -1/3 $, dvs. $ \rho = -1/3 $ og $ \theta = \pi/3 $ radianer. Fladen $ g $ kan ses på første slide i figuren nederst i eksemplet og linien $ \ell $ er angivet med blå. $ h(t) $ kan ses som den røde kurve. På slide to er arealet vi finder fremhævet.
Vi kan parameterisere linien som
$$[ \mathbf{r}_\ell(t) = \begin{bmatrix}x(t) \\ y(t)\end{bmatrix} = -\frac{1}{3}\begin{bmatrix} \cos(\pi/3) \\ \sin(\pi/3) \end{bmatrix} + t\begin{bmatrix} -\sin(\pi/3) \\ \cos(\pi/3) \end{bmatrix} = -\frac{1}{3}\begin{bmatrix} \frac{1}{2} \\ \frac{\sqrt{3}}{2} \end{bmatrix} + t\begin{bmatrix} -\frac{\sqrt{3}}{2} \\ \frac{1}{2} \end{bmatrix}. ]$$
Så er $ h(t) $ givet ved
$$[ h(t) = g\left(-\frac{1}{6} - t\frac{\sqrt{3}}{2}, -\frac{\sqrt{3}}{6} + t\frac{1}{2}\right) = \frac{385}{432}+\frac{\sqrt{3}}{108}t-\frac{73}{72}t^2-
\frac{\sqrt{3}}{12}t^3+\frac{3}{16}t^4.
]$$
Ved at løse $ x(t) = \pm 1 $ og $ y(t) = \pm 1 $ kan vi se at $ t_0 = -7/(3\sqrt{3}) $ og $ t_1 = 5/(3\sqrt{3}) $. Vi kan nu regne integralet
$$[ \int_\ell g(x,y)\,ds = \int_{-\frac{7}{3\sqrt{3}}}^{\frac{5}{3\sqrt{3}}} \frac{385}{432}+\frac{\sqrt{3}}{108}t-\frac{73}{72}t^2-
\frac{\sqrt{3}}{12}t^3+\frac{3}{16}t^4\,dt ]$$
$$[ = \left[\frac{385}{432}t+\frac{\sqrt{3}}{216}t^2-\frac{73}{216}t^3-
\frac{\sqrt{3}}{48}t^4+\frac{3}{80}t^5\right]_{-\frac{7}{3\sqrt{3}}}^{\frac{5}{3\sqrt{3}}} = \frac{848\sqrt{3}}{1215}. ]$$
På tredje slide (er ikke med på print) kan figuren ses med en række markerede slices, samt profilkurven for den valgte retning $ \theta $.
Første slide: Funktionen $ g $, linien $ \ell $ og kurven $ h(t) $.
Så hvad bruger vi kurveintegralerne til? Hvis objektet vi gerne vil skanne er placeret i et $ xy $-koordinatsystem og absorbtionsevnen ("tætheden") af objektet i punktet $ (x,y) $ er beskrevet ved $ g(x,y) $, kan vi nu bestemme dæmpningen af en stråle der passerer igennem vores objekt som integralet langs linien strålen følger.
Radontransformationen
Radontransformationen, $ \mathcal{R} $, er en smule langhåret. Den er en funktion, der tager en funktion af 2 variable som input og giver en ny funktion af 2 variable som output. Så hvis $ g(x,y) $ er en funktion af 2 variable er radontransformationen af $ g $, $ \mathcal{R}g(\rho,\theta) $, også en funktion af 2 variable (her $ \rho $ og $ \theta $). Det er ikke uden grund, at vi genbruger $ \rho $ og $ \theta $ her, hvilket vi vil se.
Lad $ \ell_{\rho,\theta} $ være linien givet ved $ x\cos\theta + y\sin\theta = \rho $. Hvis vi regner kurveintegralet af en funktion $ g(x,y) $ langs denne linie giver det os et tal, $ I $.
$$[ I = \int_{\ell_{\rho,\theta}} g(x,y)\,ds ]$$
Men $ I $ varierer når vi ændrer på $ \rho $ og $ \theta $. Hvis vi fastholder $ \theta $, får vi en funktion $ p_\theta(\rho) $, som vi kalder profilkurven for $ \theta $.
$$[ p_\theta(\rho) = \int_{\ell_{\rho,\theta}} g(x,y)\,ds ]$$
Den siger noget om vores funktion set præcis fra vinklen $ \theta $. Radontransformationen er samlingen af alle disse profilkurver.
Radontransformationen af $ g(x,y) $ er givet ved
$$[ r(\rho,\theta) = \mathcal{R}g(\rho,\theta) = \int_{\ell_{\rho,\theta}} g(x,y)\,ds = \int_{-\infty}^{\infty} g(\rho\cos\theta - t\sin\theta, \rho\sin\theta + t\cos\theta) \,dt, ]$$
Det bør nævnes at for at Radontransformationen eksisterer har vi brug for at $ g $ opfylder visse egenskaber. Hvis vi fx kræver at det $ g(x,y) = 0 $ når $ \|(x,y)\| > K $ (dvs. uden for en cirkel med radius $ K $) for et eller anden tal $ K > 0 $ er det klart at vi kan begrænse integralet fra $ -\infty \rightarrow \infty $ til $ -K \rightarrow K $ (formentlig kan vi begrænse mere, men det er ikke til at sige uden at kende $ g $).
Så hvad er Radontransformationen af $ g $ så? Hvis vi forstår $ g $ som absorbansen af et objekt og kigger på formlen, kan vi se at $ \mathcal{R}g $ i et punkt $ (\rho,\theta) $ er præcis den samlede dæmpning observeret langs linien beskrevet ved $ \rho $ og $ \theta $. Herunder kan ses et par animerede billeder og deres korresponderende Radontransformation. Læg mærke til at på Radontransformationen (højre billede) er billedet på højre og venstre side af den røde linie faktisk det samme blot spejlet. Det er fordi tætheden af billedet er ens om du kommer fra den ene eller den anden side.
Med Radontransformationen er vi nu i stand til at projektere vores objekt fra en masse vinkler. Men hvad med den anden vej? Hvordan kan vi få processen til at gå baglens så vi kan finde vores oprindelige objekt igen? En naiv men simpel måde, der faktisk godt kan virke lidt er en metode der kaldes back projection. Essentielt set går back projection ud på at for hver vinkel $ \theta $ tage profilkurven $ p_\theta(\rho) $ og smøre den ud over vores område. Formuleret som et integrale, ser det sådan her ud
$$[ \widehat{g}(x,y) = \mathcal{B}r = \int_0^\pi r(x\cos\theta + y\sin\theta, \theta) d\theta. ]$$
Her er $ \mathcal{B} $ back projection-funktionen, der tager Radontransformationen $ r(\rho,\theta) $ og giver funktionen $ \widehat{g}(x,y) $. Vi bruger her $ \widehat{g} $ i stedet for $ g $ for at gøre opmærksom på, at det altså ikke er helt den rigtige løsning vi finder her. Faktisk kan det vises at
$$[ \widehat{g}(x,y) = h\ast g(x,y) = \int h(s,t)g(x-s,y-t)\,ds\,dt ,\ \text{hvor}\ h(x,y) = \frac{1}{\sqrt{x^2 + y^2}}. ]$$
$ h \ast g $ kaldes foldningen af $ h $ og $ g $ (convolution på engelsk). Det betyder at $ \widehat{g} $ er en mere udglattet version af $ g $. $ \widehat{g}(x,y) $ inderholder værdien af $ g $ i $ (x,y) $, men også nogle bidrag fra andre punkter tæt på, $ (x \pm \delta x, y \pm \delta y) $.
En mere avanceret metode, der giver meget skarpere billeder, er filtered back projection. Filtered back projection gør brug af et teknisk redskab kaldet Fouriertransformationen, der viser sig at være nært beslægtet med Radontransformationen. Formlen ser således ud
$$[ \widehat{g}(x,y) = \mathcal{FB}r = \frac{1}{4\pi^2}\iiint r(\rho,\theta) e^{i\omega(x\cos\theta + y\sin\theta - \rho)}\omega\,d\rho\,d\omega\,d\theta. ]$$
hvor $ \mathcal{FB} $ er vores filtered back projection-funktion og $ i $ er det komplekse tal med egenskaben $ i^2 = -1 $. Groft sagt er det, der foregår i filtered back projection, at vi i stedet for at kigge på vores radon transformation, kigger på, hvilke bøljer den består af. Hvor hurtige de er (frekvensen) og hvor store de er (amplituden). Fouriertransformationen af en funktion fortæller os noget om frekvenserne og amplituden af bøljerne i funktionen. For at lave spring (svarende til kanter i et billede) i en funktion kræves høje frekvenser, så det vi gør for at få dem til at stå tydeligere frem, er simpelthen at vægte de høje frekvenser mere før vi rekonstruerer vores billede. Det er det $ \omega $'et i formlen laver.
For mere information om Fourier transformationen og en udledning af filtered back projection formlen herover, se Advanced sektionen på WhyDoMath.org: Tomography.
På billedet her ser du Shepp-Logan fantomet (øverst til venstre) og radontransformationen/sinogrammet (nederst til venstre), samt rekonstruktion med back projection (øverst til højre) og med filtered back projection (nederst til højre).
Octave
Octave er et gratis program til numeriske beregninger. Det ligner i høj grad det kommericelle produkt MATLAB, der ofte benyttes akademisk til computerdrevne beregninger. Vi vil i denne sektion gennemgå nogle eksempler på brug af Octave.
% =======================================
% PROFILE CURVE EDITOR
% =======================================
%
clear all % Delete All Variables
close all % Close all figure windows
clc % Clear the command window
% Check if image package is loaded
try %
radon(0,0); % (A): If this fails it runs (B)
catch, %
pkg load image; % (B): Loads image package
end %
% Set global variables
global ax_im ax_rad dx xn yn;
global selected patch_handles;
% Setup grid
nGrid = 3;
dx = 1/nGrid;
selected = false(nGrid);
patch_handles = zeros(nGrid);
% Draw figure and axes
fig = figure(1, 'Position',[100,100,1000,500]);
ax_im = axes('Position',[1/50,1/25,23/50,23/25]);
for i = 1:4
ax_rad(i) = axes('Position',[55/100,(i-1)*23/100+1/25,43/100,22/100]);
end
hold(ax_im, 'on')
% Help lines
for i = 1:(nGrid-1)
plot(ax_im, [0,1], [i*dx,i*dx], '--', 'Color',[.5,.5,.5])
plot(ax_im, [i*dx,i*dx], [0,1], '--', 'Color',[.5,.5,.5])
end
% Axis limits
xlim(ax_im, [0,1])
ylim(ax_im, [0,1])
set(ax_im, 'xtick',[], 'ytick',[])
set(ax_rad, 'xtick',[], 'ytick',[])
% Labels
labels = 'ABCD';
for i = 1:4
uicontrol('Style','Text', 'String',labels(i), 'Units','normalized',...
'Position',[51/100,(i-1)*23/100+1/25,2/100,22/100]);
end
% Grid Centers
xn = dx/2 + dx*[0:nGrid-1];
yn = dx/2 + dx*[0:nGrid-1];
% Functions
function p = drawPatch(xn, yn)
global ax_im dx;
p = patch(ax_im, xn + dx/2*[-1,1,1,-1], yn + dx/2*[-1,-1,1,1],[.5,.5,.5]);
end
% Callback
function select(s, e)
global ax_im xn yn;
global selected patch_handles;
p = get(ax_im,'currentpoint');
p = p(1,1:2);
if p(1) < 0 || p(1) > 1 || p(2) < 0 || p(2) > 1
return
end
[val,Xidx] = min(abs(xn-p(1)));
[val,Yidx] = min(abs(yn-p(2)));
if selected(Xidx,Yidx)
selected(Xidx,Yidx) = false;
delete(patch_handles(Xidx,Yidx));
else
selected(Xidx,Yidx) = true;
patch_handles(Xidx,Yidx) = drawPatch(xn(Xidx),yn(Yidx));
end
drawRadon();
end
function drawRadon()
global selected ax_rad;
E = ones(25);
im = kron(selected,E);
[RT,XT] = radon(im,[0:45:135]);
c = [.6,0,0;
0,.6,0;
0,0,.6;
.6,.6,0];
for i = 1:4
plot(ax_rad(5-i),XT,RT(:,i),'-','Color',c(i,:),'LineWidth',2);
set(ax_rad,'xtick',[],'ytick',[])
end
end
% Set callbacks
set(fig,'windowbuttondownfcn',@select);
Om mig
Mit navn er Bjørn Jensen, og jeg er uddannet Civilingeniør inden for matematik fra DTU Compute med specialisering i funktionalanalyse og optimering.
Johan Radon
Slide 3
Shepp-Logan fantomet
Første slide: Funktionen $ g $, linien $ \ell $ og kurven $ h(t) $.
Andet slide: Det integrerede område.
