Normat 56:4, 177–186 (2008) 177

Datortomograﬁns matematik –

Om en matematisk teori med många nya tillämpningar

Jan Boman

Matematiska Institutionen

Stockholms Universitet

106 91 Stockholm

jabo@math.su.se

Datortomograﬁ utvecklades under 1970-talet och kom i mycket allmänt bruk under

1980-talet, trots att en datortomograf redan då kostade någon miljon kronor. Att

datortomograﬁ är en soﬁstikerad och mycket eﬀektiv teknik för röntgenundersök-

ning för medicinsk diagnostik är i dag välkänt. Mindre välkänt är förmodligen

det faktum att avancerad matematik spelar en väsentlig roll för datortomografens

funktion.

Vid datortomograﬁ avbildar man ett plant snitt genom det område som man

vill undersöka. Metoden kallas därför ibland för skiktröntgen, och ordet tomograﬁ

är bildat av det grekiska ordet tomos, som betyder sektion eller snitt. Anledningen

till att matematiken behövs kan sägas vara att det är omöjligt att direkt avbilda

eller mäta de storheter som man är intresserad av, och att apparaten därför i stället

mäter helt andra storheter. Matematiken utgör förbindelselänken mellan det man

mäter och det man vill veta. Det som man vill veta är värdet på absorptionskoef-

ﬁcienten för röntgenstrålning i varje punkt i ett plant snitt genom det undersökta

området. Olika slags vävnad har olika täthet och därmed något olika absorptions-

koeﬃcient. Således kan sjuk vävnad, exempelvis tumörvävnad, ofta kännas igen

genom något avvikande absorptionskoeﬃcient. Om vi känner ett gott närmevärde

på absorptionskoeﬃcienten i ett tillräckligt tätt gitter av punkter i det plana snit-

tet, så är det mycket lätt att med datorns hjälp rita upp en bild av snittet genom

att tätare medium ges högre svärtning på bilden eller tvärtom. Låt oss tänka oss

att vi söker värdet på absorptionskoeﬃcienten i varje ruta i ett rutnät bestående av

kvadrater med 0, 5 millimeters sida. Om det undersökta området är en cirkelskiva

med 25 centimeters diameter, så har vi alltså bortåt 200 000 kvadrater och lika

många tal som ska bestämmas.

Problemet är att vi inte kan mäta dessa värden direkt, utom möjligen genom

att såga ut den tunna skiva som vi vill studera och genomlysa den vinkelrätt mot

sitt eget plan. Men detta vill vi förstås inte göra med en levande människa och inte

heller med en dyrbar maskin, exempelvis en rymdraket. Vad man kan göra är att

genomlysa skivan i dess eget plan, i många olika riktningar.

Låt oss försöka uttrycka matematiskt den information som en sådan undersök-

ning kan ge. Antag att vi skickar en röntgenstråle med känd intensitet I

genom

området så som i ﬁguren nedan. Med en detektor mäter vi den utgående strålens

178 Jan Boman Normat 4/2008

intensitet I

. Kvoten I

är ett mått på den totala absorptionen längs strålens

hela väg. Om strålen passerar genom n stycken rutor med absorptionskoeﬃcient

, . . . , a

, där talet n är exempelvis 500, så kan man säga att vi härigenom har

mätt värdet av

(1) a

+ a

+ . . . + a

i=1

Närmare bestämt gäller att log(I

) är

proportionell mot kvantiteten (1). Vi har

därmed fått ett känt samband mellan ett

antal av alla de obekanta talen. Genom att

genomlysa i samma riktning men med pa-

rallellförﬂyttat läge (se ﬁguren) kan vi få

ytterligare 500 samband. Lika många till

kan vi få genom att genomlysa i den vin-

kelräta riktningen.

Poängen är nu att man kan mäta upp ytterligare samband mellan de obekanta

storheterna genom att genomlysa området i ﬂer riktningar, exempelvis parallellt

med rutornas diagonaler. En liten komplikation blir då att det inte längre är så

lätt att säga hur lång väg vår röntgenstråle passerar genom var och en av de små

rutorna. Men detta går att ta hänsyn till; jag går inte in på dessa detaljer här.

Det ﬁnns inte heller någon anledning att inskränka sig till rutnätets diagonalrikt-

ningar. För varje ny riktning som vi genomlyser området får vi ett stort antal nya

linjära samband mellan de sökta talen. Sammanfattningsvis, om man genomlyser

den tunna plana skivan i låt oss säga 500 olika riktningar och registrerar cirka 500

mätvärden för varje riktning, så får man totalt 250 000 mätvärden som vart och

ett måste vara en vägd summa av ett stort antal av de okända talen; vikterna i

summan kommer sig som sagt av att strålen i fråga har passerat olika långt genom

de olika rutorna längs dess väg. Därmed kan problemet beskrivas som ett gigan-

tiskt linjärt ekvationssystem med (minst) cirka 200 000 obekanta. Sådana har vi ju

goda och välkända metoder för att lösa, och vi har utmärkta snabba datorer. Så

problemet förefaller att vara löst, åtminstone i princip.

Emellertid vore ett sådant ekvationssystem på tok för stort även för dagens da-

torer, åtminstone för att behandlas med standardmetoder. Ty antalet operationer

som krävs för att lösa ett linjärt ekvationssystem med vanlig Gausselimination är

proportionellt mot kuben på antalet obekanta, och 200 000

= 8 · 10

. Om da-

torn gör en miljard operationer i sekunden så skulle kalkylen alltså ta nästan 10

sekunder, det vill säga några månader!

Vi behöver därför en bättre matematisk analys av problemet. Det visar sig att

det är bättre att formulera och analysera problemet i kontinuerliga termer. Vi frågar

därför åter efter värdet av absorptionskoeﬃcienten i varje punkt i den (tunna) plana

skivan. Absorptionskoeﬃcienten kan då betraktas som en funktion på skivan, det

vill säga som en funktion på ett plant område. Om vi inför ett koordinatsystem i

skivans plan, kan vi alltså betrakta den sökta funktionen som en funktion f(x, y)

av två variabler. Nästa steg är att uttrycka våra mätvärden i termer av denna

funktion.

Normat 4/2008 Jan Boman 179

Antag att vi gör en mätning med hjälp av en röntgenstråle parallellt med x-axeln

längs linjen y = b. Denna stråles dämpning längs hela sin väg utgörs av summan av

all dämpning längs dess väg genom området, och det är lätt att se att i detta fall

(logaritmen för) den totala dämpningen måste vara proportionell mot integralen

(2)

f(x, b)dx.

Obervera att den ändliga summan (1) skulle kunna skrivas som

i=1

f(x

, b), där

= f(x

, b) är värdet av absorptionskoeﬃcienten i punkten (x

, b), en punkt i

den i:te rutan längs strålens väg. Om vi betecknar den i:te kvadratiska rutans sida

med ∆x

, så är denna summa lika med en konstant gånger

i=1

f(x

, b)∆x

, ett

uttryck som vi ju känner igen som ett närmevärde till integralen (2). Jag har inte

satt ut några gränser, eftersom det kan underförstås att vi ska integrera över den

mängd av x för vilka (x, y) ligger inom området, vilket vi kan anta är ett intervall.

Alternativt kan vi deﬁniera funktionen f (x, y) som noll överallt utanför området

och integrera från −∞ til ∞.

Uttrycket (2) brukar kallas för linjeintegralen eller kurvintegralen av funktionen

f över den räta linjen L, och betecknas

(3)

fds,

där ds betecknar längdelementet längs linjen L och anger att integrationen sker

med avseende på båglängden. Med hjälp av en parameterframställning av den räta

linjen kan man lätt deﬁniera integralen (3) för vilken rät linje L som helst. Därmed

har vi gjort oss oberoende av koordinatsystemet och kan nu beskriva de mätvärden

som erhålls med vilken strålriktning som helst. Det matematiska problem som

datortomografen måste lösa kan därmed, i sin kontinuerliga form, uttryckas på

följande sätt:

(P)

En obekant funktion f är deﬁnierad på ett område i planet, exempelvis

en cirkelskiva. Antag att man känner värdet av linjeintegralen

f ds

för varje rät linje L genom området. Beräkna funktionen f .

I en berömd artikel från 1917 löstes detta problem av den österrikiske matemati-

kern Johann Radon. I artikeln gavs en explicit formel för värdet av funktionen f i

punkten (x, y) uttryckt i värdena av alla linjeintegralerna

f ds. När Radon löste

detta problem hade han givetvis ingen aning om att lösningen skulle få viktiga tek-

niska tilllämpningar, utan han studerade problemet uteslutande därför att frågan

var matematiskt intressant och visade sig ha en vacker lösning. Efter utvecklingen

av datortomograﬁ har Radons problem och många av dess varianter och generali-

seringar studerats intensivt av matematiker, och avbildningen L 7→

f ds kallas

numera för Radontransformen av funktionen f .

I praktiken undersöker man vanligen ett tredimensionellt område genom att

avbilda en hel serie av parallella plana snitt på så sätt som ovan beskrivits. Praktiskt

sker detta oftast genom att patienten successivt förﬂyttas i riktning vinkelrätt

mot de plana snitten. Detta innebär förstås att man för undersökning av ett tre-

dimensionellt område använder sig endast av röntgenstrålar som är parallella med

ett ﬁxt plan. För att eﬀektivisera datainsamlingen och minimera stråldosen har

180 Jan Boman Normat 4/2008

man i vissa nyare datortomografer övergivit principen att analysera en tunn skiva

i taget och går i stället direkt på problemet att rekonstruera funktionen f i ett tre-

dimensionellt område. Mätdata baseras i dessa fall vanligen på en större mängd av

strålriktningar, inte alla parallella med ett och samma plan (multislice tomography,

helical tomography etc.). Att på ett eﬀektivt sätt utnyttja denna större mängd av

mätdata leder till mycket svårare, genuint tredimensionella, matematiska problem,

som ännu är långt ifrån lösta.

Den stora framgången med datortomograﬁ för medicinsk diagnostik har lett

till att samma eller liknande mättekniker har utvecklats för många andra ända-

mål. Förutom röntgenstrålning används bland annat ultraljud, gammastrålning

och elektronstrålar. Mätmetoder baserade på kärnspinnresonans (MRI, Magnetic

Resonance Imaging) leder till liknande matematiska problem som datortomograﬁ.

En av de tidigaste tillämpningarna av tomograﬁska metoder uppkom inom geo-

den seismiska vågens ankomst till en annan punkt B för ett stort antal par av

punkter A och B beräknar man utbredningshastigheten för seismiska vågor i olika

punkter i jordens inre och drar därav slutsatser om den fysiska beskaﬀenheten hos

jordklotets inre delar.

En snarlik teknik används för malm- och oljeprospektering. En explosion arran-

geras vid jordytan, och tidpunkt och styrka hos reﬂexerna mäts med detektorer vid

ett stort antal punkter på jordytan. Experimentet upprepas med ett antal andra

lägen för explosionen. Med utgångspunkt från dessa mätdata försöker man beräkna

beskaﬀenheten hos jordlagren under mätområdet.

Inom industrin används tomograﬁska metoder för produktkontroll och process-

kontroll.

Elektronmikroskopi kan framgångsrikt kombineras med tomograﬁ genom att

provet genomlyses med elektronstrålar i ett stort antal riktningar.

Många ﬂer exempel skulle kunna ges. Området utvecklas mycket snabbt. Sökning

med Google på ”tomography” ger 10 miljoner träﬀar!

De olika tillämpningarna ger upphov till liknande, men ändå olika, matematiska

problem. Analys av problem med anknytning till Radontransform och tomograﬁ är

därför ett aktivt forskningsområde inom matematik.

För utvecklingen av datortomograﬁn belönades fysikern Allan Cormack, upp-

vuxen i Sydafrika och senare verksam vid Tufts University nära Boston, och den

brittiske ingenjören Godfrey Hounsﬁeld med nobelpriset i fysiologi eller medicin

1979. Under en kort tid som sjukhusfysiker på 1950-talet sysslade Cormack med

strålterapi mot cancer. Vid strålterapi önskar man planera bestrålningen så att

stråldosen mot tumören maximeras medan övrig vävnad skadas så litet som möj-

ligt. Cormack insåg att man skulle kunna förbättra strålningsplaneringen om man

kände den variabla absorptionskoeﬃcienten för strålningen i det bestrålade områ-

det. Med syfte att bestämma denna formulerade Cormack, ovetande om Radons

arbete 1917, datortomograﬁns fundamentala problem (P) och löste det. Resultaten

publicerades i två artiklar i början av 1960-talet. Dessa artiklar rönte till en bör-

jan föga uppmärksamhet. Oberoende av Cormack byggde Hounsﬁeld i början av

1970-talet den första datortomografen. Därefter gick utvecklingen mycket snabbt.

På återstoden av dessa sidor ska jag säga några ord om hur man analyserar det

matematiska problemet och hur beräkningen görs. Jag måste börja med att införa

några beteckningar.

Normat 4/2008 Jan Boman 181

Om α är en vinkel och p ett reellt tal är

x = −t sin α + p cos α, y = t cos α + p sin α, t ∈ R,

en parameterframställning med båglängden som parameter av en linje L(α, p) med

vinkeln α mot y-axeln och avståndet |p| från origo. Om L = L(α, p) så kan inte-

gralen (3) därmed skrivas

(4) g

(p) =

L(α,p)

f ds =

f(−t sin α + p cos α, t cos α + p sin α)dt.

Index R vid integraltecknet betyder att integrationen sker över hela reella axeln R.

Våra mätdata kan alltså beskrivas av ett ändligt antal funktioner g

(p) av en reell

variabel p för olika värden på α, α = α

, . . . , α

, där 0 ≤ α

< π. I datortomografer

kan man välja exempelvis m ≈ 500. Naturligtvis är i praktiken även p en diskret

variabel, men det bortser vi från här. Problemet är som sagt att beräkna funktionen

f(x, y) utgående från funktionerna g

(p).

Vi kommer att behöva använda några enkla egenskaper hos Fouriertransforma-

tionen. För en funktion u(t) av en variabel deﬁnieras Fouriertransformen bu genom

bu(τ) =

u(t)e

−itτ

dt, τ ∈ R.

Funktionen u(t) kan (under lämpliga förutsättningar på funktionen u) återvinnas

från Fouriertransformen genom den så kallade inversionsformeln

(5) u(t) =

2π

bu(τ)e

itτ

dτ, t ∈ R.

Formeln (5) är den matematiska formuleringen av det välkända faktum att ett god-

tyckligt tidsförlopp kan beskrivas som en summa av sinusformade svängningar med

olika frekvenser; det komplexa talet bu(τ) anger amplitud och fas för svängningen

med frekvensen τ. För en funktion f(x, y) av två variabler deﬁnieras Fouriertrans-

formen genom en motsvarande dubbelintegral

(6)

f(ξ, η) =

f(x, y)e

−i(xξ+yη)

dxdy,

och inversionformeln lyder i detta fall

(7) f(x, y) =

(2π)

f(x, y)e

i(xξ+yη)

dξdη.

Nu ﬁnns ett enkelt samband mellan den 1-dimensionella Fouriertransformen av

funktionen g

(p) i (4) och den 2-dimensionella Fouriertransformen av f. Om vi

nämligen sätter in η = 0 i deﬁnitionen (6) av den 2-dimensionella Fouriertransfor-

men

f och skriver dubbelintegralen som en upprepad integral, så får vi

f(ξ, 0) =

(

f(x, y)dy)e

−ixξ

(x)e

−ixξ

dx = bg

(ξ),

(8)

182 Jan Boman Normat 4/2008

där g

(x) är uttrycket (4) svarande mot α = 0, det vill säga

(p) =

f(p, y)dy.

(Det är förstås likgiltigt om vi kallar variabeln för x eller p och om integrations-

variabeln heter y eller t.) Formeln (8) innebär därmed att den 1-dimensionella

Fouriertansformen av mätdata svarande mot strålriktning parallellt med y-axeln

är lika med restriktionen av den 2-dimensionella Fouriertransformen av den sökta

funktionen f till en linje genom origo i ξη-planet, nämligen linjen η = 0. Genom en

lätt räkning inses att samma utsaga gäller för godtycklig strålriktning, nämligen i

formelspråk

(9) cg

(τ) =

f(τ cos α, τ sin α)

för alla τ ∈ R och alla vinklar α. Det betyder att om vi känner g

(p) för alla p och

ett visst α, så kan vi räkna ut värdet av den 2-dimensionella Fouriertransformen

f(ξ, η) för alla (ξ, η) av formen (ξ, η) = (τ cos α, τ sin α) för τ ∈ R. Då α genomlöper

ett stort antal värden α = α

så kommer dessa linjer att ganska väl täcka en

stor del av ξη-planet, varmed den 2-dimensionella Fouriertransformen

f(ξ, η) är

approximativt känd. Det återstår därmed endast att invertera den 2-dimensionella

Fouriertransformen för att bestämma f (x, y).

Detta resonemang visar att den sökta funktionen måste vara entydigt bestämd

av alla funktionerna g

(p) — i varje fall om g

(p) vore känd för alla α. Resonemang-

et ger också följande metod för att beräkna f(x, y) utifrån mätdata g

(p). Beräkna

Fouriertransformen av var och en av de uppmätta funktionerna g

(p) = g

(p) för

k = 1, 2, . . . , m. Använd formeln (9) för att beskriva den 2-dimensionella Fourier-

transformen av f på ett gitter bestående av ett stort antal räta linjer genom origo.

Beräkna till sist f genom att invertera den 2-dimensionella Fouriertransformen.

För numerisk beräkning av Fouriertransformer är mycket eﬀektiva algoritmer

kända sedan 1960-talet (Fast Fourier Transform, FFT). För det sista steget i den

nyssnämnda proceduren, nämligen invertering av den 2-dimensionella Fouriertrans-

formen

f(ξ, η), behöver man emellertid känna

f(ξ, η) på ett rektangulärt gitter sna-

rare än ett radiellt. Problemet att ﬁnna tillräckligt goda närmevärden för

f(ξ, η)

på ett rektangulärt gitter utgående från kända värden på ett radiellt gitter visade

sig dock vara svårare än väntat. Man utför därför vanligen i praktiken beräkningen

på ett annat sätt, som nu ska beskrivas.

För var och en av de m funktionerna g

(p) bildar vi en ny funktion h

(p) genom

formeln

(10) h

(p) =

(p − q)K(q)dq,

där K(p) är en viss ﬁx funktion, som strax ska beskrivas. Uttrycket i högra ledet

ovan kallas för faltningen av funktionerna g

och K. För var och en av funktionerna

(p) bildar vi sedan en funktion av två variabler (x, y) som är konstant på linjer

med normalriktningen (cos α

, sin α

) genom

(11) h

(x cos α

+ y sin α

Normat 4/2008 Jan Boman 183

Till sist beräknar man halva medelvärdet av alla funktionerna (11); resultatet är

(12) f(x, y) ≈

(x cos α

+ y sin α

Det kan tyckas anmärkningsvärt att en så enkel formel ger svaret på ett så kom-

plicerat problem. Noteras bör också att denna beräkningsprocedur lät sig beskrivas

i helt elementära termer; Fouriertransformation och annan avancerad matematik

behövdes inte.

Och hur stort är räknearbetet? Om vi antar att var och en av funktionerna

(p) representeras genom n stycken tal g

), i = 1, . . . , n, så kräver evalueringen

av integralen (10) ungefär n operationer för varje p, alltså ≈ n

operationer, och

för samtliga k får vi därmed totalt ≈ mn

operationer. Evaluering av formeln

(12) kräver ≈ m operationer för varje punkt (x, y), och om vi har ≈ n

punkter i

området, så kräver detta ytterligare ≈ mn

operationer. Om m ≈ n ger detta totalt

2mn

≈ n

operationer. Med exemplvis n ≈ 500 krävs således ≈ 10

operationer,

vilket datorn utför på bråkdelen av en sekund.

Jag ska skissera ett bevis för formeln (12) för x = y = 0, det vill säga

(13) f(0, 0) ≈

(0).

Härav följer (12) ganska lätt för godtyckligt (x

, y

) genom att man tillämpar (13)

på funktionen f(x

+ x, y

+ y). Vi börjar med att använda inversionsformeln (7)

i två variabler för (x, y) = (0, 0), det vill säga

f(0, 0) =

(2π)

f(ξ, η)dξdη.

Polära koordinater ξ = r cos α, η = r sin α, ger

f(0, 0) =

(2π)

2π

∞

f(r cos α, r sin α)r drdα.

Om vi gör variabeltransformationen α 7→ α + π i den yttre integralen och trans-

formationen r 7→ −r i den inre integralen, får vi samma uttryck med den inre

integralen ersatt av

−∞

f(r cos α, r sin α)|r| dr. Genom att addera den så erhållna

formeln till den föregående får vi

(14) 2f(0, 0) =

(2π)

2π

∞

−∞

f(r cos α, r sin α)|r| drdα.

Den 1-dimensionella inversionsformeln (5) använd för t = 0 innebär att

(15) g

(0) =

2π

(τ)dτ.

184 Jan Boman Normat 4/2008

Om vi kombinerar detta med formeln (9), ﬁnner vi att

(16) g

(0) =

2π

∞

−∞

f(τ cos α, τ sin α)dτ.

Om vi inte hade haft faktorn |r| i integralen (14), så hade vi nu kunnat dra slutsat-

sen av (14) att f(0, 0) är approximativt lika med en konstant gånger medelvärdet

av alla g

(0).

För att ta hänsyn till faktorn |r| väljer vi nu en funktion K(p) med egenskapen

att dess Fouriertransform uppfyller villkoret

K(τ ) ≈ |τ |

för alla τ i ett ganska stort intervall |τ | < A och

K(τ ) = 0 utanför ett något större

intervall. Deﬁnitionen av h

genom (10) tillsammans med formeln

h(τ) = bu(τ)bv(τ )

för Fouriertransformen av en faltning h(t) =

u(t − s)v(s)ds ger nu att

(τ) = bg

(τ)

K(τ ) ≈ bg

(τ)|τ| =

f(τ cos α

, τ sin α

)|τ| för |τ | < A.

Analogt med (16) ﬁnner vi därmed att

(17) h

(0) =

2π

∞

−∞

(τ)dτ ≈

2π

∞

−∞

f(τ cos α

, τ sin α

)|τ|dτ.

Om vi nu ersätter den yttre integralen i (14) med en Riemannsumma enligt for-

meln

2π

ϕ(α)dα ≈ (2π/m)

k=1

ϕ(α

) och använder uttrycket (17) för den inre

integralen, så får vi formeln (13).

Kalkylen illustreras i bilderna 1 - 5 på sidan 186. Bild 1 är en välkänd testbild,

den så kallade Shepp-Logan-fantomen, som har varit i bruk sedan 1970-talet och

som påminner om ett plant snitt genom en skalle. Bild 2 illustrerar samtliga mätta

data, och kallas i tomograﬁkretsar för sinogrammet. Om L(α, p) är linjen med

ekvation

(18) x cos α + y sin α = p,

så är den vertikala axeln α-axeln (med ett halvt varv graderat från 0 till 500) och

den horisontella axeln är p-axeln. Gråskalevärdet hos punkten (p, α) i sinogrammet

anger alltså det uppmätta värdet av integralen av f över linjen L(α, p). Namnet

sinogram kommer sig av att ett föremål som beﬁnner sig i punkten (x, y) ger upphov

till sinuskurvan (18) i p α-planet i sinogrammet.

För att underlätta tolkningen av Bild 2 har jag gjort Bild 3, som är ett horisonellt

utsnitt svarande mot ett ﬁxt α i sinogrammet, nämligen α = 400 i den skala som

används. Observera att grafen i Bild 3 säger precis samma sak som gråskaleremsan.

Normat 4/2008 Jan Boman 185

Indata till beräkningen, det vill säga mätdata, består alltså precis av 500 sådana

funktioner som anges av grafen i Bild 3.

Bild 4 illustrerar rekonstruktionen, det vill säga funktionen (12). Funktionen

K(p) har deﬁnierats genom att

K(τ ) har satts lika med |τ| cos(τ cπ/2) för |τ| < 1/c

och lika med noll för övriga τ, och konstanten c har valts lika med 0, 003. I praktiken

använder man sig av ett antal ﬁnesser för att få bästa möjliga rekonstruktion, men

här har beräkningen gjorts exakt enligt beskrivningen i texten ovan. Vid noggrann

jämförelse av Bild 4 och Bild 1 ser man att rekonstruktionen ser helt korrekt ut,

förutom att den är något oskarpare än testbilden. Oskärpan beror huvudsakligen på

att vi genom att ersätta funktionen τ 7→ |τ | i Fourierrummet med den nyssnämnda

funktionen

K(τ ) har försummat höga frekvenser.

Vid vissa tillämpningar av den tomograﬁska metoden (exempelvis inom elektron-

mikroskopi) ﬁnns det ett intervall av riktningar α inom vilket man inte kan mäta

integralerna

L(α,p)

f ds. I sådana situationer är det mycket svårt eller omöjligt

att rekonstruera den sökta funktionen noggrant. Bild 5 illustrerar detta. Bilden

visar en delsumma av summan (12) där talet m = 500 har ersatts med 350. Detta

innebär att mätdata har utnyttjats endast från strålriktningar svarande mot sju

tiondelar av ett halvt varv i stället för ett halvt varv. Följden har blivit att alla

konturer som är parallella med de saknade strålriktningarna är både suddiga och

felaktigt återgivna.

Sammanfattningsvis, det som man har kunnat se eller mäta direkt med röntgen

i detta fall är 500 gråskaleremsor, alternativt grafer, liknande dem i Bild 3. Det är

ju uppenbart att ﬁnare strukturer hos fantombilden, exempelvis de små cirkelski-

vorna i mitten av bilden, är omöjliga att se på någon av de 500 röntgenbilderna,

och det samma gäller i själva verket även de grövre strukturerna, möjligen med

undantag av fantombildens yttre skal, ”skallbenet". Men de små cirkelskivorna har

tydligen ändå lämnat — för våra ögon osynliga — spår i de 500 graferna, för med

matematikens hjälp kan vi se dem och exakt ange deras läge.

Ett stort tack riktas till Hans Rullgård, som har skrivit Matlabprogrammet och givit

värdefulla kommentarer till tidigare versioner av manuskriptet.

Litteratur

Frank Natterer: The mathematics of computerized tomography, Wiley & Sons, 1986.

Nytryck utgivet av Society for Industrial and Applied Mathematics 2001.

Frank Natterer och Frank Wübbeling: Mathematical methods in image reconstruc-

tion, Society for Industrial and Applied Mathematics, 2001.

Charles Epstein: Introduction to the mathematics of medical imaging, Pearson

Education/Prentice Hall, 2003.

Wikipedia: ”Tomography” (http://en.wikipedia.org/wiki/Tomography), eller

”Computed tomography” (http://en.wikipedia.org/wiki/Computed_tomography).

186 Jan Boman Normat 4/2008

50 100 150 200 250 300 350 400 450 500

100

150

200

250

300

350

400

450

500

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Bild 1

50 100 150 200 250 300 350 400 450 500

100

150

200

250

300

350

400

450

500

0.1

0.2

0.3

0.4

0.5

0 50 100 150 200 250 300 350 400 450 500

0.2

0.4

0.6

0.8

50 100 150 200 250 300 350 400 450 500

Bild 2 Bild 3

50 100 150 200 250 300 350 400 450 500

100

150

200

250

300

350

400

450

500

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

50 100 150 200 250 300 350 400 450 500

100

150

200

250

300

350

400

450

500

−0.2

0.2

0.4

0.6

0.8

Bild 4 Bild 5