Variationell Regularisering av Inversproblem med
Tillämpning inom Skiktröntgen
Variational Regularization of Inverse Problems with Appli-
cation in Computed Tomography
Examensarbete för kandidatexamen i matematik vid Göteborgs universitet
Kandidatarbete inom civilingenjörsutbildningen vid Chalmers
Louise Christensen
Sara Nilsson
Oskar Pauli
Gustav Svensson
Institutionen för Matematiska vetenskaper
CHALMERS TEKNISKA HÖGSKOLA
GÖTEBORGS UNIVERSITET
Göteborg, Sverige 2022

Variationell Regularisering av Inversproblem med Tillämpning
inom Skiktröntgen
Examensarbete för kandidatexamen i matematik vid Göteborgs universitet
Louise Christensen
Examensarbete för kandidatexamen i matematik inom Matematikprogrammet, in-
riktning Tillämpad matematik, vid Göteborgs universitet
Gustav Svensson Sara Nilsson
Kandidatarbete i matematik inom civilingenjörsprogrammet Teknisk matematik vid
Chalmers
Oskar Pauli
Handledare: Axel Ringh
Institutionen för Matematiska vetenskaper
CHALMERS TEKNISKA HÖGSKOLA
GÖTEBORGS UNIVERSITET
Göteborg, Sverige 2022

Förord
Detta kandidatarbete är skrivet vid Göteborgs Universitet och Chalmers Tekniska högskola v̊aren
2022. Arbetet är författat av Louise Christensen, Sara Nilsson, Oskar Pauli och Gustav Svensson.
De individuella bidragen är dokumenterade i gruppens dagbok, samt i en individuell tidslogg och
sammanfattas i bidragsrapporten nedan.
Louise och Sara har lagt störst fokus p̊a teoridelen och har samarbetat genom hela processen.
Oskar och Gustav har haft störst fokus p̊a programmeringen och att ta fram resultat.
Slutligen vill vi tacka v̊ar handledare Axel Ringh för all hjälp och stöd.
Bidragsrapport
I bidragsrapporten är huvudförfattare angivet med fet stil, medan resterande namn är medförfat-
tare. Bidragsrapporten ses p̊a nästa sida.
Bidragsrapport
Rapportskrivning
Avsnitt Titel Namn
Populärvetenskaplig text Gustav
Sammanfattning Sara
Abstract Oskar
1 Inledning Louise, Sara
2.1 Bakgrund för skiktröntgen Louise, Sara
2.1.1 Modellering av absorption Louise
2.2.1 Fram̊atproblem Louise, Sara
2.2.2 Inversproblem Louise, Sara
3.1 Optimeringsteori Louise, Sara
3.1.1 Steepest descent-metoden Louise
3.1.2 Gradientprojektionsmetoden Louise
3.2 Variationell regularisering Louise, Sara, Gus-
tav
3.2.1 Tikhonovregularisering Louise, Gustav
3.2.2 Total variationsregularisering Louise, Sara, Gus-
tav
3.3 Optimering av Tikhonovregularisering Louise
3.4 Optimering av total variationsregularisering Louise, Sara
4.1 Implementering av optimeringslösare Louise, Sara
4.2 Analys av regulariseringsparametrar med Sara, Gustav
hjälp av L-kurvmetoden
4.2.1 L-kurvillustrationer Oskar
5.1 Resultat Louise, Oskar
5.2 Rekonstruktion genom bak̊atprojektion Sara, Oskar
5.2.1 Jämförelse av resultat med bak̊atprojektion Oskar
6.1 Utveckling Louise, Sara, Oskar
6.2 Samhälleliga och etiska aspekter Louise, Sara, Oskar
6.3 Slutsats Louise, Gustav
MATLAB-kod och övrigt
Tikhonovregularisering Louise, Oskar, Gus-
tav
Total variationsregularisering Louise, Sara,Oskar,
Gustav
L-kurvmetoden Louise, Oskar, Gus-
tav
Sinogram Sara
Visning av weakly coercive Louise, Gustav
Populärvetenskaplig presentation
Vad är det som egentligen händer när vi tar röntgenbilder? Hur kan vi f̊a fram bilder av skelett,
muskler, fett med mera? Bakom detta finns det en rad olika matematiska verktyg och resultat som
tillsammans gör det möjligt att f̊a fram bilder p̊a exempelvis en hjärna eller ett ben. I denna uppsats
kommer vi beskriva matematikens roll inom röntgenbildtagning - mer specifikt inom skiktröntgen.
Röntgenstr̊alning upptäckes redan 1895 av Wilhelm Röntgen. Han var ute efter att studera
n̊agot som kallas elektriska urladdningsfenomen i lufttomt rum. Hans idé var att skjuta elektrisk
ström genom ett rör som var täckt med en svart kartong s̊a att inget ljus kunde ta sig igenom. I
Röntgens labb fanns det en pappskiva med fosfor p̊a. Fosfor har egenskapen att sända ut synligt
ljus vid kontakt med vissa typer av str̊alning [ne.se1][humanprogress.org]. Efter att ha skjutit
den elektriska strömmen genom röret märkte han att pappskivan glimrade till och började genast
studera detta fenomen. Han l̊aste till och med in sig själv i sitt laboratorium för att studera
detta i sju veckor utan att prata med n̊agon. Det var efter dessa veckor han tog världens första
röntgenbild av handen p̊a sin fru. P̊a bilden kunde hennes vigselring ses hänga p̊a skelettet [Kje17].
Denna nya upptäckt blev fort populär runt om i världen, och självklart även i Sverige. Här togs
första röntgenbilden tre veckor efter att W. Röntgen hade föreläst om denna nya metod för första
g̊angen [ne.se3]. Han tog inget patent p̊a tekniken d̊a han ins̊ag potentialen inom sjukv̊ard och ville
f̊a ut den s̊a fort som möjligt till allmänheten [Kje17].
Sedan 70-talet har det även tillkommit skiktröntgen, där str̊alningen skjuts in fr̊an flera olika
vinklar och mäter absorptionen i dessa. Syftet med skiktröntgen är att kunna urskilja olika vävnader
i kroppen utan att behöva göra medicinska ingrepp. Detta är en väldigt effektiv metod för att
upptäcka exempelvis tumörer eller annan sjuk vävnad. [Bom08, s. 177-178]
D̊a återst̊ar fr̊agan, vad har matematik med detta att göra? Är det inte bara att skjuta str̊al-
ningen genom kroppen och se vad som kommer ut p̊a andra sidan? Jo, p̊a sätt och vis. Men faktum
är att matematiken spelar en stor roll i processen av att f̊a fram en bild. Röntgenbilder kan ses
som skuggbilder där skuggorna representerar mängden absorption av röntgenstr̊alningen när den-
na rör sig genom vävnaden. Men att endast veta absorptionen i en vävnad säger inte s̊a mycket
om vad som kan gömma sig, exempelvis tumörer. Därför tas skuggbilder fr̊an flera olika vinklar
inom skiktröntgen för att sedan kunna kombinera dessa för att se hur vävnaden ser ut. Det in-
nebär att vi känner till utfallet, allts̊a absorptionen, men inte orsaken. Detta är vad vi kallar för
ett inversproblem inom matematiken. Det kan uppst̊a inversproblem som är illaställda, vilket kan
beskrivas som ’väldigt sv̊ara att lösa’, dock inte omöjliga. Med hjälp av matematiska verktyg som
optimering och variationell regularisering kan vi hitta en lösning till v̊art illaställda inversproblem,
allts̊a återskapa röntgenbilden.
Optimering innebär att försöka hitta den bästa lösningen givet vissa förutsättningar. Det finns
olika sätt att g̊a tillväga gällande optimering. I v̊art arbete använder vi oss av en startgissning av
röntgenbilden för att sedan förbättra denna i flera omg̊angar. Denna process upprepas till dess att
vi inte längre ser n̊agon förbättring.
Variationell regularisering använder sig av tv̊a termer - ett huvuduttryck och ett feluttryck. Vi
vill ligga n̊agorlunda vid samma ställe, allts̊a h̊alla oss s̊a snarlik originalbilden som möjligt. Därför
best̊ar feluttrycket oftast av n̊agon typ av bestraffning, som i v̊art fall beror p̊a storleken av v̊ar
beräknade lösning. Den straffar oss om vi hittar lösningar som rör p̊a sig för mycket. Lösningen blir
därmed unik eftersom vissa lösningar till problemet förkastas när de rör sig för l̊angt fr̊an tidigare
observationer. [ne.se2]
Vi kan allts̊a med hjälp av dessa matematiska verktyg lösa det illaställda inversproblemet. Mer
specifikt; vi kan med hjälp av matematiken ta fram exempelvis röntgenbilder av ett barns brutna
arm eller skanningsbilder av ett kranium.
Denna uppsats använder en förbestämd bild som vi vill återskapa med matematiska verktyg. För
att återskapa denna s̊a bra s̊a möjligt kommer vi att undersöka olika typer av regulariseringsmetoder
och ta fram fyra olika röntgenbilder. Vi kommer att använda dessa för att kunna jämföra metoderna
och göra en visuell jämförelse. Uppsatsens utvecklingspotential kommer även att undersökas.
Sammanfattning
Detta kandidatarbete har fokuserat p̊a lösning av illaställda inversproblem vid rekonstruk-
tion av data med hjälp av variationell regularisering. All programmering har skett i MATLAB.
Data har simulerats fr̊an MATLAB:s Shepp-Loganfantom med upplösning 256×256 pixlar.
Fantomen har rekonstruerats med hjälp av optimering av olika regulariseringsmetoder. Opti-
meringsalgoritmer har implementerats för varje regularisering. Vi har använt Tikhonov- och
total variationsregularisering, b̊ada med och utan ett bivillkor för icke-negativitet, och jäm-
fört dessa med varandra. Regulariseringarna har berott p̊a olika parametrar, där en av dessa
är regulariseringsparametern λ. Regulariseringsparametern har utvärderats genom testning
och L-kurvmetoden. Resultaten har presenterats i form av fyra rekonstruerade bilder med re-
spektive värden p̊a regulariseringsparametern och jämförts med varandra. Det fanns tydliga
skillnader mellan Tikhonov- och total variationsregularisering där vi särskilt ser en skillnad i
hur de hanterar brus. Möjlig vidareutveckling av arbetet diskuterades.
Abstract
The aim of this study was to solve ill-posed inverse problems when reconstructing data by
using variational regularization theory. This problem was solved using MATLAB. Data was
simulated from MATLAB’s 256×256 Shepp-Logan phantom, which also acted as our reference
image. This phantom was reconstructed using Tikhonov regularization and total variation
regularization, both with and without a non-negativity constraint. The regularization methods
rely on a regularization parameter λ for which we used L-curves and testing to find appropriate
values. Results was presented in the form of the reconstructions together with their respective
values of the regularization parameter and compared against each other. There were clear
differences between Tikhonov regularization and total variation regularization, mainly in how
they handled noise. Further development was discussed.
Inneh̊all
1 Inledning 1
2 Skiktröntgen och inversproblem 2
2.1 Bakgrund kring skiktröntgen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
2.1.1 Modellering av absorption . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
2.2 Introduktion till inversproblem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.2.1 Fram̊atproblem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.2.2 Inversproblem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
3 Optimerings- och variationell regulariseringsteori 5
3.1 Optimeringsteori . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
3.1.1 Steepest descent-metoden . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
3.1.2 Gradientprojektionsmetoden . . . . . . . . . . . . . . . . . . . . . . . . . . 6
3.2 Variationell regularisering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
3.2.1 Tikhonovregularisering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
3.2.2 Total variationsregularisering . . . . . . . . . . . . . . . . . . . . . . . . . . 8
3.3 Optimering av Tikhonovregularisering . . . . . . . . . . . . . . . . . . . . . . . . . 10
3.4 Optimering av total variationsregularisering . . . . . . . . . . . . . . . . . . . . . . 10
4 Implementering 11
4.1 Implementering av optimeringslösare . . . . . . . . . . . . . . . . . . . . . . . . . . 12
4.2 Analys av regulariseringsparametrar med hjälp av L-kurvmetoden . . . . . . . . . 13
4.2.1 L-kurvillustrationer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
5 Resultat och jämförelse 14
5.1 Resultat . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
5.2 Rekonstruktion genom bak̊atprojektion . . . . . . . . . . . . . . . . . . . . . . . . . 16
5.2.1 Jämförelse av resultat med bak̊atprojektion . . . . . . . . . . . . . . . . . . 16
6 Utveckling och slutsats 17
6.1 Utveckling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
6.2 Samhälleliga och etiska aspekter . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
6.3 Slutsats . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
Referenser 19
A Appendix 20
A.1 Bevis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
A.1.1 Sats 3.4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
A.1.2 Sats 3.5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
A.1.3 Normens konvexitet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
A.2 MATLAB-kod . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
A.2.1 Sinogram . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
A.2.2 Tikhonovregularisering med steepest descent, L-kurva . . . . . . . . . . . . 22
A.2.3 Tikhonovregularisering med gradientprojektion, L-kurva . . . . . . . . . . . 24
A.2.4 Total variationsregularisering med steepest descent, L-kurva . . . . . . . . . 26
A.2.5 Total variationsregularisering med gradientprojektion, L-kurva . . . . . . . 28
1 Inledning
Inom m̊anga olika omr̊aden finns det ett behov av att ta fram bilder av otillgänglig information.
Detta behov kan matematiken hjälpa till med. Baserat p̊a mätdata är det önskvärt att ta fram det
okända, det vill säga, baserat p̊a en känd effekt ta fram en okänd orsak. Detta kan beskrivas som ett
inversproblem. Det är viktigt för samhället att kunna lösa s̊adana problem d̊a inversproblem dyker
upp p̊a många ställen, som exempelvis röntgen, ultraljud, ekolokalisering etc. Att lösa inversproblem
är därför väldigt relevant.
Skiktröntgen är ett s̊adant problem, där vi vill ta fram en bild av insidan av kroppen. Det är
ett viktigt redskap för att diagnostisera olika sjukdomar, skador och kan minska andelen invasiva
ingrepp. Inom skiktröntgen är det allts̊a relevant att kunna rekonstruera s̊a bra och precisa bilder
som möjligt. Vid röntgenprocessen skjuts röntgenstr̊alar genom en kropp med en viss ing̊angsinten-
sitet, där vävnaden absorberar delar av röntgenstr̊alarna. Efter att röntgenstr̊alarna har passerat
genom kroppen kan utg̊angsintensiteten mätas. Detta genererar en mängd data. Modelleringen av
denna absorptionen ger ett förh̊allande mellan mätdata och absorptionen av data. Detta förh̊allan-
de kan beskrivas som Radontransformen. Denna kommer vi att använda som v̊ar fram̊atoperator
under rekonstruktionen.
Målet är att rekonstruera en bild p̊a insidan av kroppen som är s̊a bra som möjlig. Dock är
inversproblemet i detta fall illaställt, det vill säga att vi inte kan ta fram inversen med direkta be-
räkningar och därmed rekonstruera bilden. Istället måste vi använda oss av s̊a kallad regularisering.
Mer specifikt variationell regularisering som möjliggör lösning av det illaställda inversproblemet.
Genom att använda olika regulariseringsmetoder kan vi göra det illaställda problemet välställt.
Detta arbete kommer att fokusera p̊a fyra olika typer av regularisering - Tikhonovregularisering
och total variationsregularisering, b̊ada med och utan ett bivillkor för icke-negativitet. Vi kan med
hjälp av optimering av kostnadsfunktionen för respektive regulariseringsmetod rekonstruera bilden.
Regulariseringarna kommer undersökas och deras egenskaper kommer jämföras, dock vill vi enbart
fokusera p̊a att återskapa Shepp-Loganfantomen som finns inbyggd i MATLAB. Vi kommer därmed
inte kunna generalisera v̊ar optimeringslösare till andra typer av inversproblem. Arbetet fokuserar
p̊a en studie av dessa specifika metoder och inte p̊a att göra en fullskalig utvärdering av alla
typer av regulariseringar som lösning till inversproblem. Vidare innebär detta att utvärderingen
av lösningarna blir begränsad till de framtagna regulariseringsmetoderna. Arbetet kommer inte att
bidra med ny forskning, utan är en fokusstudie av dessa specifika metoder.
Under optimering av rekonstruktionsmetoder kan det uppst̊a olika samhälleliga och etiska aspek-
ter som behöver tas i beaktande. Exempel p̊a detta är integritetsaspekter, balansen mellan kvalitet
p̊a rekonstruktion och tids̊atg̊ang och kostnadsaspekter. Detta diskuteras i avsnitt 6.2.
Detta kandidatarbete är skrivet till matematik- och ingenjörstuderande som har matematiska
förkunskaper inom linjär algebra och matematisk analys. Arbetets teoridel redogör för de teoretiska
kunskaper som är nödvändiga för att först̊a projektet. Arbetet är uppbyggt av en introduktion till
inversproblem, grundläggande optimeringsteori och variationell regularisering som är nödvändig
för först̊aelsen av implementering av optimeringslösaren. Avslutningsvis presenteras resultatet, en
diskussion kring detta och utvecklingspotentialen för arbetet.
1
2 Skiktröntgen och inversproblem
Vi inleder med en genomg̊ang av hur skiktröntgen fungerar och hur absorption av röntgenstr̊alning i
ett material kan modelleras. Vidare beskrivs generell teori bakom inversproblem samt optimerings-
och regulariseringsteori. Detta introduceras för att sedan tillämpa denna teori vid rekonstruktion
av MATLAB:s Shepp-Loganfantom, se Figur 4 i kapitel 4.
2.1 Bakgrund kring skiktröntgen
Skiktröntgen används bland annat inom sjukv̊arden för att ta detaljerade bilder av kroppens organ.
Denna typ av röntgen utförs genom att röntgenstr̊alar med en viss intensitet skjuts genom en
kropp fr̊an flera olika vinklar och sedan mäts av en detektor p̊a motsatt sida kroppen[1177.se].
När röntgenstr̊alar rör sig genom kroppen absorberas delar av str̊alningen vilket resulterar i en
reducerad intensitet när str̊alarna har passerat genom kroppen. Detta generar mätdata vilken kan
användas till att återskapa bilden.
I följande sektion introducerar vi modelleringen av absorption, där vi f̊ar fram relationen mellan
absorptionsintensiteten och mätdata. Denna relation visar sig vara Radontransformen.
2.1.1 Modellering av absorption
För att modellera absorptionen vid skiktröngten behöver vi ta fram relationen mellan intensiteten
av str̊alningen I och absorptionsfunktionen f . Vid skiktröntgen skjuts det, som tidigare nämnts, ett
antal röntgenstr̊alar linjärt genom objektet fr̊an olika vinklar. Vi kan parametrisera varje enskild
str̊ale med enhetsvektorn, u ∈ R2, dess ortogonal, û ∈ R2, och avst̊andet, s, fr̊an origo till en punkt
x p̊a röntgenstr̊alen [Ker16, s. 21]. Detta illustreras i Figur 1.
Figur 1: Illustration av parametriseringen av en str̊ale genom ett objekt. u är enhetsvektorn, û
dess ortogonal, s avst̊and och t längd.
Värdet av f(x) i en given punkt, x, beror p̊a vävnaden och dess absorption. Olika vävnader
absorberar olika mycket str̊alning. Betrakta exempelvis en kropp - denna best̊ar av olika vävna-
der med olika absorptionskoefficienter. Skelett har en högre absorption än exempelvis organ och
framst̊ar därmed tydligare p̊a en röntgenbild [ne.se3]. Vi kan allts̊a återskapa en bild av objek-
tet genom att bestämma absorptionsfunktionen f . Fr̊an [Ker16, s. 22] f̊ar vi följande ekvation för
intensitetsdifferensen vid avst̊andet s fr̊an origo och vid längden t längs linjen i punkten su+ tû
∆I(su+ tû) = −I(su+ tû)f(su+ tû)∆t.
L̊at ∆t → 0, vilket ger differentialekvationen
d
I(su+ tû) = −I(su+ tû)f(su+ tû).
dt
2
Vi delar med intensiteten p̊a b̊ada sidor och löser differentialekvationen. Ing̊angs- och utg̊angsin-
tensiteten betecknas med I0 respek(tive Il, vilket ger oss relationen mellan intensiteten och absorp-tionsfunktionen ) ∫
I0(s,u)
ln = f(su+ tû)dt, (1)
Il(s,u) R
som är givet av integration längs str̊alningslinjerna. Allts̊a är logaritmen av ing̊angsintensiteten
dividerat med utg̊angsintensiteten lika med linjeintegralen av absorptionsfunktionen längs str̊a-
lens linje. Ing̊angsintensiteten, I0, är känd och utg̊ansintensiteten, Il, mäts av en detektor och
är därmed ocks̊a känd, se Figur 1. Därför har vi de facto mätt linjeintegralen av den kvantitet
vi är intresserade av. Vi har nu ett uttryck som förbinder mätdata med en absorptionsfunktion
vilket ger oss möjligheten att ta fram en bild av insidan av objektet. Högerledet i ekvation (1)
är precis Radontransformen [Fee15, s. 14]. All data vi f̊ar är genererad fr̊an integration längs al-
la individuella röntgenstr̊alar vi skjuter genom objektet. Vi kan allts̊a mäta värdet av integralen
av absorptionsfunktionen, men inte av själva absorptionsfunktionen f . Detta är ett inversproblem
vilket g̊as igenom i nästa avsnitt.
Figur 2: Illustration av en geometri för skiktröntgen där parallella str̊alar används.
Som indikerat p̊a Figur 2 ovan kommer vi att använda oss av parallella str̊alar när vi simu-
lerar datan. I kapitel 6 kommer vi även att nämna möjligheten att använda fan-beam str̊alar till
simuleringen av mätdata.
2.2 Introduktion till inversproblem
För att kunna rekonstruera en bild av ett givet objekt måste vi först̊a den bakomliggande teorin
kring inversproblem och den möjliga problematiken vid lösning av dessa. Som namnet antyder är
inversproblemet motsatsen till ett annat problem. Det senare kallas ofta för fram̊atproblem, och
innan vi kan ge oss p̊a inversproblemet m̊aste vi först först̊a detta.
2.2.1 Fram̊atproblem
Informellt handlar fram̊atproblemet om att beräkna effekten baserat p̊a en given och känd orsak.
Mer formellt: fram̊atproblemet är att givet fram̊atoperatorn A : X → Y och vektorn f ∈ X,
beräkna b ∈ Y , [LLM16, s. 52], det vill säga
Af = b.
Fram̊atproblemet inom skiktröntgen är att givet fram̊atoperatorn, A, och absorptionsfunktio-
nen, f , beräkna data b. I v̊art arbete är A och b kända, men f saknas. Allts̊a har vi ett s̊a kallat
inversproblem.
3
2.2.2 Inversproblem
Vid inversproblem, motsatt fram̊atproblem, är effekten given, men orsaken är okänd. Det vill säga,
givet operatorn A : X → Y och vektorn b ∈ Y vill vi bestämma f ∈ X s̊adan att Af = b [LLM16,
s. 122].
För att bestämma en exakt lösning av inversproblemet behöver problemet vara välställt. Enligt
[Had1923] är ett problem välställt om följande krav är uppfyllda;
a) en lösning existerar för alla b ∈ Y ;
b) lösningen är unik;
c) lösningen är kontinuerligt beroende av data.
Den första punkten, att en lösning existerar, innebär att operatorn A är surjektiv. Andra punkten
innebär att A är injektiv, allts̊a att det enbart finns en lösning till problemet. Om b̊ade krav a)
och b) är uppfyllda s̊a är A bijektiv, och inversoperatorn, A−1, kan definieras p̊a ett entydigt sätt
[LLM16, s. 130]. Det tredje kravet, att lösningen är kontinuerlig, innebär att A−1 är kontinuerlig.
Om problemet är välställt kan vi bestämma f genom f = A−1b. Ett problem som inte uppfyller
dessa krav sägs vara illaställt.
Exempel 1
Betrakta följande ekvationssystem, Af = b, med A : R3 → R3;    1 2 3 2 4 6
x1
  
1
x2 = 2 .
0 0 0 x3 0
Här ser vi tydligt att kolumnvektorerna i A är linjärt beroende d̊a exempelvis tredje kolumnen kan
skrivas som en linjär kombination av de tv̊a första och andra kolumnen som en multipel av första,
         1 2 3 1      ·   2 + 4 = 6 och 2 2 = 
2
 4 .
0 0 0 0 0
Detta innebär att rang(A) = 1, och att lösningen ej är unik. Allts̊a är problemet illaställt.
Inom skiktröntgen är Radontransformen fram̊atoperatorn A, medan matrisen b är den mätdata
som skanningen producerar. I v̊art arbete har vi ett inversproblem som är illaställt d̊a Radon-
transformens invers inte är kontinuerlig, vilket inte uppfyller det tredje kravet för ett välställt
problem enligt [Had1923]. Att denna inte är kontinuerlig innebär att det inte är möjligt att hitta
en tillfredsställande lösning p̊a inversproblemet, oavsett numerisk metod [Ker16, s. 5].
Datan b beror p̊a absorptionen genom objektet, vilket ger oss ett antal mätvärden för varje
inskjutningsvinkel. När vi skjuter röntgenstr̊alar genom ett objekt är mätdatan p̊averkad av en
viss mängd brus, ϵ. Vi beskriver därför problemet som
b = Af + ϵ,
där ϵ ∼ N(0, δ). D̊a röntgenstr̊alar best̊ar av fotoner, [arpansa.gov.au], kan en diskret brusfördel-
ning som Poissonfördelningen vara ett annat alternativ för att modellera bruset. Vi väljer dock att
använda oss av normalfördelningen d̊a v̊art fokus i tillämpningen ligger p̊a jämförelse av regularise-
ringsmetoder. Bruset används ocks̊a för att säkerställa att vi vid lösningen av inversproblemet inte
beg̊ar ett s̊a kallat inversbrott. Ett inversbrott innebär att den teoretiska och syntetiska datan är
för lika varandra, vilket kan medföra att rekonstruktionerna har för hög artificiell kvalitet [MS12,
s. 7]. Denna kvalitet kan vara bättre än vad som f̊as vid användning av riktig mätdata.
4
För att lösa v̊art illaställda inversproblem behöver vi göra det välställt. En metod för detta är
variationell regularisering som med hjälp av optimering säkerställer att en lösning uppfyller de tre
tidigare nämnda kraven. I följande kapitel redogör vi först för optimeringsteori, för att sedan g̊a
vidare med teorin av olika former av variationell regularisering som vi kommer att använda oss av
för rekonstruktion av fantombilden.
3 Optimerings- och variationell regulariseringsteori
För att kunna göra lösa det illaställda problemet kommer vi att använda oss av optimerings- och
variationell regulariseringsteori. I de följande avsnitten g̊ar vi igenom teori kring optimering och
variationell regularisering, för att sedan redogöra för hur vi implementerar och tillämpar denna
teori.
3.1 Optimeringsteori
Optimering är ett omr̊ade inom matematiken som handlar om att hitta en funktions maximala
eller minimala värde, allts̊a att hitta ett optimum till en målfunktion. Optimeringen av en s̊adan
m̊alfunktion kan vara med eller utan bivillkor. Ett standardproblem inom optimering är p̊a formen
min F (x)
(2)
s.a. x ∈ S,
där S ⊆ Rn är en icke-tom mängd och F : Rn → R. [AEP17, s. 3-18]
Inom optimeringsteori existerar ett viktigt begrepp som vi kommer att använda oss av, nämligen
konvexitet. En konvex funktion kan definieras enligt följande.
Definition 3.1. [AEP17, Definition 3.39, s. 65] (Konvex funktion) Antag att S ⊆ Rn. En funktion
F : Rn → R ∪ {+∞} är konvex i y∈ S om
x ∈ S 
λ ∈ (0, 1)  =⇒ F (λy + (1− λ)x) ≤ λF (y) + (1− λ)F (x).
λy + (1− λ)x ∈ S
Funktionen F är konvex p̊a S om den är konvex för varje y ∈ S.
En funktion är allts̊a konvex om en linjär interpolation av funktionsvärdena aldrig är mindre
än själva funktionen mellan punkterna. Vi kan d̊a, för tv̊a godtyckliga punkter i den konvexa
mängden, alltid rita en linjär interpolation mellan tv̊a punkter i mängden S som är större eller lika
med funktionen F (x).
För att g̊a vidare behöver vi fler verktyg och definierar begreppen globalt minimum och lokalt
minimum. För den senare definitionen använder vi den öppna Euklidiska bollen Bϵ(x
∗) := {y ∈
Rn; ‖y − x∗‖2 < ϵ}, med radie ϵ och med centrum i x∗ [AEP17, s. 82].
Definition 3.2. [AEP17, Definition 4.1, s. 82] (Globalt minimum) Betrakta standardproblemet (2).
L̊at x∗ ∈ S. Vi säger att x∗ är ett globalt minimum av F p̊a S om F uppn̊ar sitt lägsta värde p̊a
S i x∗. Med andra ord, x∗ ∈ S är ett globalt minimum av F p̊a S om
F (x∗) ≤ F (x), x ∈ S
h̊aller.
Definition 3.3. [AEP17, Definition 4.2, s. 82] (Lokalt minimum) Betrakta standardproblemet (2).
L̊at x∗ ∈ S.
a) Vi säger att x∗ är ett lokalt minimum av F p̊a S om det existerar en tillräckligt liten boll
som skär S runt x∗ s̊a att x∗ är en globalt optimal lösning i den mindre mängden. Med andra ord,
x∗ ∈ S är ett lokalt minimum av F p̊a S om
∃ϵ > 0 s.a. F (x∗) ≤ F (x), x ∈ S ∩Bϵ(x∗). (3)
5
b) Vi säger att x∗ ∈ S är ett strikt lokalt minimum av F p̊a S om, i (3), det r̊ader strikt olikhet
d̊a x =6 x∗.
Dessa definitioner använder vi i följande satser.
Sats 3.4. [AEP17, Theorem 4.3, s. 82] (Fundamentalsatsen för Global Optimalitet) Betrakta stan-
dardproblemet (2) och anta att S är en konvex mängd och att F är konvex p̊a S. D̊a är varje lokalt
minimum av F p̊a S även ett globalt minimum.
Sats 3.4 innebär att om vi hittar ett lokalt minimum av F p̊a den konvexa mängden S har vi
även hittat funktionens globala minimum. För att hitta minimum kan vi använda oss av följande
sats om nödvändiga och tillräckliga krav för global optimalitet.
Sats 3.5. [AEP17, Theorem 4.18, s. 94] (Nödvändiga och Tillräckliga Krav för Global Optimalitet)
Antag att F : Rn → R är konvex och är C1 p̊a Rn. D̊a gäller
x∗är ett globalt minimum av F p̊a Rn ⇔ ∇F (x∗) = 0n.
Denna sats säger att om gradienten∇F (x∗) = 0 i en viss punkt x∗ s̊a är x∗ ett globalt minimum
av F , och vice versa. Vi kan med hjälp av Sats 3.4 och Sats 3.5 hitta den optimala lösningen till
standardproblemet (2). Bevis för Sats 3.4 och Sats 3.5 kan ses i appendix A.1.1 respektive A.1.2.
Det finns många olika typer av optimeringsalgoritmer som kan användas till att hitta en opti-
mal lösning. Vi kommer att använda oss av tv̊a iterativa optimeringsalgoritmer, steepest descent-
metoden1 och gradientprojektionsmetoden, för att ta fram minimum. Steepest descent-metoden
används för optimering utan bivillkor, medan gradientprojektionsmetoden hanterar bivillkor p̊a
lösningen.
3.1.1 Steepest descent-metoden
För att minimera standardproblemet (2) utan bivillkor använder vi oss av den iterativa metoden
steepest descent. Denna metod bygger p̊a en algoritm som söker efter en nedstigningsriktning2 pk
fr̊an en given punkt xk. Fr̊an denna punkt tar vi ett steg, αk > 0, i nedstigningsriktningen och
uppdaterar positionen,
xk+1 = xk + αpk.
För att säkerställa att en riktning är nedstigande m̊aste F (xk+αkpk) < F (xk) gälla. Det finns
flera olika metoder för att bestämma nedstigningsriktningen. Steepest descent-metoden definierar
pk = −∇F (xk). Efter uppdateringen av positionen undersöker vi om ett termineringskrav är
uppfyllt. Termineringskravet är oftast definierat som en minsta toleransgräns mellan position xk
och xk+1 [AEP17, s. 291]. Om termineringskravet inte är uppfyllt börjar algoritmen om fr̊an
position xk+1 med att söka efter en ny nedstigningsriktning.
Med steepest descent-metoden är v̊ar till̊atna mängd Rn×n. Vi visar sedan i avsnitt 3.2.2 att
det finns ett optimum i denna mängd och vi vet därmed att algoritmen konvergerar.
3.1.2 Gradientprojektionsmetoden
Vid optimering har vi behov av en algoritm som kan hantera bivillkor, mer specifikt bivillkoret
x ≥ 0. Vi använder oss därför av optimeringsmetoden gradientprojektion som kan hantera detta
villkor. Gradientprojektionsmetoden liknar steepest descent-metoden och använder motsvarande
den negativa gradienten, −∇F (xk), som nedstigningsriktning. Vid uppdatering av lösningen an-
vänder gradientprojektion sig av projektionen av xk − αk∇F (xk), enligt [AEP17, s. 332]
xk+1 = ProjX [xk − αk∇F (xk)], k = 1, ...,
där α > 0. ProjX(y) = argmin‖y − x‖, det vill säga punkten x ∈ X som ligger närmast punkten
x∈X
y. Vid en position utanför den till̊atna mängden projiceras punkten tillbaka in p̊a mängden. För-
tydligat innebär detta att om xk −∇F (xk) ovan är en till̊aten lösning uppdateras positionen p̊a
1Direkt översättning är ’brantaste nedstigningen’. Vi kommer att använda oss av den engelska termen steepest
descent.
2Direkt översättning av den engelska termen, descent direction. Vi kommer att använda oss av nedstigningsrikt-
ning.
6
samma sätt som steepest descent-metoden. Om xk − α∇F (xk) ligger utanför mängden projiceras
punkten till närmaste punkt i den till̊atna mängden.
Betrakta det 1-dimensionella fallet med R+ som den till̊atna mängden. Om x ligger utanför R+
projiceras denna till närmaste punkt i den till̊atna mängden. Av Figur 3 framg̊ar det tydligt att
närmaste punkten är exakt 0.
I v̊art fall är den till̊atna mängden Rn×n+ . Detta innebär att inga element i v̊ar lösning f̊ar vara
mindre än 0. Ett element i x som är mindre än 0 projiceras till närmaste punkten i den till̊atna
mängden R+. Som med steepest descent-metoden visar vi i avsnitt 3.2.2 att ett optimum i Rn×n+
existerar och att gradientprojektionsmetoden därmed konvergerar.
Figur 3: Illustration av gradientprojektion i 1-dimension med R+ som till̊aten mängd.
3.2 Variationell regularisering
Variationell regularisering är en metod för att göra ett illaställt problem välställt och bygger p̊a
minsta kvadratmetoden. Kort beskrivet formuleras rekonstruktionsproblemet som ett minimerings-
problem som sedan optimeras. Mer specifikt, målet är att hitta matrisen g ∈ Rn×n som ger den
optimala lösningen f i minimeringsproblemet
f = argmin‖Ag − b‖22 + λG(g).
g
Första termen är en datamatchningskomponent som matchar given data b mot data genererad av
fram̊atoperatorn, A, och ett givet g. G(g) är en regulariseringsfunktional och 0 < λ < ∞ är en
regulariseringsparameter. Regulariseringsfunktionalen och regulariseringsparametern säkerställer
att problemet är bijektivt och att inversen är kontinuerlig. Vi har därmed ett välställt problem
enligt [Had1923]. Vi definierar kostnadsfunktionen som
H(g) := ‖Ag − b‖22 + λG(g). (4)
Vi kommer att använda oss av tv̊a typer av variationell regularisering - Tikhonovregularisering
och total variationsregularisering. Dessa best̊ar av likadana datamatchningskomponenter, men med
olika regulariseringsfunktionaler G(g). Fram̊atoperatorn A är Radontransformen, medan AT är
bak̊atprojektion utan filter. Det senare finns inbyggt i MATLAB som iradon.m [Niev1986, s. 80]
och används när vi g̊ar igenom optimering av de olika regulariseringsmetoderna i avsnitt 3.3,
respektive avsnitt 3.4.
För att säkerställa att den till̊atna mängden för matrisen g i regulariseringsmetoderna har ett
optimum måste vi visa att kostnadsfunktionen (4) är weakly coercive3 enligt definition 3.6. Om
kostnadsfunktionen är weakly coercive vet vi att det finns en kompakt mängd i Rn av globala
optimum fr̊an Weierstrass Sats 3.7.
Definition 3.6. [AEP17, Definition 4.4, s. 84] (Weakly coercive) Antag att S ⊆ Rn är en icke-tom
och sluten mängd, och F : S → R en given funktion.
Vi säger att F är weakly coercive om
lim F (x) = ∞
∥x∥→∞
x∈S
h̊aller.
Sats 3.7. [AEP17, Theorem 4.6, s. 86] (Weierstrass Sats)
Antag att S ⊆ Rn är en icke-tom och sluten mängd, och F : S → R en nedre semi-kontinuerlig
funktion p̊a S. Om F är weakly coercive p̊a S existerar det en icke-tom, sluten och begränsad
(kompakt) mängd av globala optimum till standardproblemet (2).
3Vi har ej hittat n̊agot etablerat svenskt uttryck och därför kommer vi att använda det engelska.
7
För att kunna tillämpa optimeringsteorin och säkerställa att det finns en lösning till v̊ara regula-
riseringsmetoder måste vi visa att det finns optimum i de till̊atna mängderna. Mer specifikt m̊aste
vi visa att b̊ade Tikhonov- och total variationsregulariseringen är weakly coercive och konvexa.
Först redogör vi för Tikhonovregularisering för att sedan g̊a igenom total variationsregularisering.
3.2.1 Tikhonovregularisering
Vanligtvis är Tikhonovregularisering första valet för illaställda linjära problem. Denna regularise-
ring ger en viss utjämning [MS12, s. 63]. Tikhonovregulariseringen ges av
fλ,ϵ = argmin‖Ag − b‖2 22 + λ‖g‖2. (5)
Först observerar vi att uttrycket (5) best̊ar av Euklidiska normer vilket innebär att funktio-
nen är konvex d̊a normen alltid är konvex. Det senare kan visas med hjälp av triangelolikheten
‖x + y‖ ≤ ‖x‖ + ‖y‖, se appendix A.1.3. Vidare är Tikhonovregulariseringen weakly coercive d̊a
datamatchningstermen är större eller lika med noll för alla g ∈ Rn×n och regulariseringsfunktio-
nalen g̊ar mot oändligheten d̊a ‖g‖2 → ∞. Det finns allts̊a en kompakt delmängd i den till̊atna
mängden som inneh̊aller åtminstone ett globalt optimum enligt Sats 3.7.
Vi kan nu applicera optimeringsteorin och hitta en optimal lösning f . Lösningen av denna kan
dock inneh̊alla negativa element, vilket fysikaliskt skulle innebära att vävnaden i dessa positio-
ner avger str̊alning och inte absorberar. För att skapa en fysikaliskt korrekt lösning lägger vi till
bivillkoret g ≥ 0
fλ,ϵ = argmin‖Ag − b‖22 + λ‖g‖22. (6)
g≥0
P̊a samma sätt som för uttrycket i (5) kan vi analogt visa med hjälp av Sats 3.7 att det existerar
ett optimum för (6).
Trots att regularisering utan bivillkor inte är fysikaliskt möjligt löser vi problemet med och utan
bivillkoret g ≥ 0 för att kunna jämföra bivillkorets p̊averkan och jämföra kvaliteten av resultaten.
3.2.2 Total variationsregularisering
En annan form av regularisering är total variationsregularisering (TV-regularisering), där regula-
riseringsfunktionalen inneh̊aller den diskretiserade differentieringsoperatorn L och definieras med
2,1-normen. 2,1-normen är en 1-norm av en 2-dimensionell vektor och beräknas genom att ta 1-
normen av 2-normen av vektorn. Till skillnad fr̊an Tikhonovregularisering är denna regularisering
kantbevarande, vilket innebär att den filtrerar bort brus samtidigt som den beh̊aller skarpa kanter
[MS12, s. 83]. Följande uttryck för TV-regulariseringen ges av [MS12, s. 90]
fλ,ϵ = argmin‖Ag − b‖22 + λ‖Lg‖2,1. (7)
Motsvarande Tikhonovregularisering kan vi även här lägga till bivillkoret g ≥ 0
fλ,ϵ = argmin‖Ag − b‖22 + λ‖Lg‖2,1. (8)
g≥0
Regulariseringsfunktionalen ‖Lg‖2,1 är definierad genom avst̊andet mellan pixlarna. [BLS16, s. 9]
ger oss att ‖Lg‖2,1 kan uttryckas
∑n ∑n
‖Lg‖2,1 = | ‖∇sgi,j‖2 |,
i=1 j=1
där ∇sgi,j = (gi,j+1 − gi,j , gi+1,j − gi,j) är en 2-dimensionell vektor och 1-normen är summan av
absolutbeloppen. Detta är den spatiella gradienten som uttrycks av spatiella derivator vilket ger
uttryck för de fysikaliska ändringarna i rummet. | ‖∇sgi,j‖2 | är 2,1-normen som uttrycks
| ‖∇sgi,j‖2 | = |√‖(gi,j+1 − gi,j , gi+1,j − g 2i,j)‖2 |
= | (gi,j+1 − gi,j)2 + (gi+1,j − gi,j)2|.
8
Detta ger oss att ∑n ∑n √
‖Lg‖2,1 = | (g 2i,j+1 − gi,j) + (gi+1,j − g )2i,j |. (9)
i=1 j=1
Denna funktional är dock inte differentierbar d̊a derivatan ej är kontinuerlig kring origo. För att
hantera detta approximerar vi absolutbeloppet i√ekvation (9) med
|t| 2β := t + β, (10)
där β > 0 är ett godtyckligt litet tal som säkerställer kontinuitet och differentierbarhet. När vi
använder approximation (10)√i ekvation (9) f̊ar vi∑n ∑n (√ )2
(g − g )2 2i,j+1 i,j + (gi+1,j − gi,j) + β =
∑i=1 ∑j=1n n √
(gi,j+1 − g 2 2i,j) + (gi+1,j − gi,j) + β.
i=1 j=1
TV-regulariseringen kan d̊a skrivas som
∑n ∑n √
fλ,ϵ = argmin‖Ag − b‖22 + λ (g 2i,j+1 − gi,j) + (g 2i+1,j − gi,j) + β. (11)
i=1 j=1
För att hantera randen kommer vi att använda oss av ett periodiskt randvillkor, det vill sä-
ga att vi p̊a randen definierar gi,n+1 = gi,1 respektive gn+1,j = g1,j . TV-regulariseringen med
approximationen av 1-normen kommer att implementeras med och utan bivillkoret g ≥ 0.
Sommed Tikhonovregularisering är kostnadsfunktionen för TV-regularisering konvex och weakly
coercive. Datamatchningstermen är en Euklidisk norm och regulariseringstermen, ‖Lg‖2,1, är en
2, 1-norm vilken är konvex. För att visa att kostnadsfunktionen till TV-regulariseringen (7) är
weakly coercive använder vi oss av definitionen av nollrummet Ker(C) = {x ∈ Rn|Cx = 0}
[LLM16, s. 217 ]. Vi vill visa att Ker(A) ∩Ker(L) = {0}. Om detta gäller s̊a g̊ar åtminstone en
av termerna i kostnadsfunktionen mot oändligheten d̊a ‖g‖2 → ∞ och är därmed weakly coercive.
Vi antar att det existerar ett g ∈ Ker(A) ∩Ker(L), g 6= 0. Vi vet att
Lg = [ (gi+1,j − gi,j , gi,j+1 − gi,j) ],
och eftersom L är en linjär operator
L : Rn×n → Rn×n×2
g −7 → [ (gi+1,j − gi,j , gi,j+1 − gi,j) ]ni,j=1,
vet vi att nollrummet till L existerar och är definierat som Ker(L) = {g ∈ Rn×n|Lg = 0}.
Nollrummet Ker(L) är när gi+1,j = gi,j och gi,j+1 = gi,j . Det vill säga att varje element i g har
samma värde för alla index (i, j) och eftersom att vi har antagit att g 6= 0, är gi,j 6= 0 för alla (i, j).
Nollrummet till L är allts̊a det linjära delrummet av alla bilder med konstant pixelvärde i hela
bilden. Detta vill säga att vi f̊ar en konstant n× n matris 
1 1 · · · · · · 1  . .1 1
. . ..
  . g = α ..
. . .. .
.
. .
.
. ..
 
 , (12)
 .. . . . . 1 1
1 · · · · · · 1 1
9
där α ∈ R\{0}. Vi har sedan tidigare antagit att g ∈ Ker(A), g 6= 0. D̊a fram̊atoperatorn A är
Radontransformen (1) och g definierat som (12) är Ag linjeintegralen över den konstanta matrisen
g. Diskretiseringen av Radontransformen ger en summation över alla str̊alar som skickas genom
pixlarna i bilden g. En summation av ett ändligt antal positiva tal är strikt positiva, analogt för
ett ändligt antal negativa tal. Allts̊a kan det g ∈ Ker(L) inte ligga i Ker(A), med andra ord
g ∈/ Ker(A). Därmed har vi en motsägelse och det gäller att Ker(A) ∩ Ker(L) = {0}, vilket
skulle visas. Eftersom snittet av nollrummen är {0}, innebär detta att b̊ada, eller åtminstone
en av termerna i TV-regulariseringen g̊ar mot oändligheten när ‖g‖2 → ∞. Detta betyder att
TV-regulariseringen (7) är weakly coercive och utifr̊an Sats 3.7 vet vi att det finns en kompakt
delmängd som best̊ar av globala optimum till funktionen.
Analogt med uttrycket i (7) kan vi med hjälp av Sats 3.7 visa att det existerar ett optimum för
uttryck (8).
I följande sektioner kommer vi att redogöra för optimeringen av de tv̊a regulariseringsmetoderna
och hur vi tar fram gradienterna av kostnadsfunktionerna.
3.3 Optimering av Tikhonovregularisering
Baserat p̊a teori kring de tv̊a optimeringsmetoderna steepest descent och gradientprojektion kan
vi använda Tikhonovregularisering för att rekonstruera fantombilden. För att kunna implementera
optimeringsmetoderna måste vi bestämma gradienten av kostnadsfunktionen för Tikhonovregula-
riseringen, HT (g). Genom att differentiera ekvation (5) med avseende p̊a g tar vi fram gradienten
till HT (g).
Vi gör en omskrivning av regulariseringsfunktionalen
∑n
‖g‖2 = g2i = gT g.
i=1
P̊a motsvarande sätt kan datamatchningstermen skrivas som
‖Ag − b‖2 = (Ag − b)T (Ag − b) = gTATAg − 2bTAg + bT b.
Kostnadsfunktionerna i (5) och (6) definieras, i enlighet med (4), som
H (g) = gTATT Ag − 2bTAg + bT b+ λ‖g‖2. (13)
För att f̊a fram gradienten differentierar vi ekvation (13) och f̊ar
∇HT (g) = 2ATAg − 2AT b+ 2λg
⇔
∇HT (g) = 2(ATA+ λI)g − 2AT b.
Med hjälp av denna gradient som nedstigningsriktning har vi nu möjligheten att hitta den optimala
lösningen och därmed absorptionsfunktionen f .
3.4 Optimering av total variationsregularisering
Analogt med Tikhonovregularisering definierar vi kostnadsfunktion av TV-regularisering som HTV (g),
vilken vi behöver bestämma gradienten av
∑n ∑n √
HTV (g) = ‖Ag − b‖22 + λ (g 2 2i,j+1 − gi,j) + (gi+1,j − gi,j) + β,
i=1 j=1
allts̊a ∑n ∑√ n
∇HTV (g) = ∇‖Ag − b‖22 + λ∇ s (gi,j+1 − gi,j)2 + (gi+1,j − g 2 i,j) + β . (14)
i=1 j=1
10
Gradienten av datamatchningstermen är likadan som tidigare, nämligen 2ATAg − 2AT g. Vi vill
beräkna gradienten i index(k, l) vilket motsvarar pixel (k, l)∑n ∑n √ ∂ (g 2 2 i,j+1 − gi,j) + (gi+1,j − gi,j) + β .
∂gk,l i=1 j=1
Vi f̊ar ett uttryck som best̊ar av tre termer när vi deriverar över en dubbelsumma med avseende
(k, l)
(√
∂
(g 2 2√ k,l+1 − gk,l) + (gk+1,l − gk,l) + β∂gk,l
+√(g 2 2k,l − gk,l−1) + (gk+1,l−1 − gk,l−1) + β)
+ (g 2k,l − gk−1,l) + (gk−1,l+1 − gk− )21,l + β .
Gradienten av regulariseringsfunktionalen beräknas genom att använda kedjeregeln. Derivatan
av första termen blir d̊a
√ 2gk,l − gk,l+1 − gk+1,l 2gk,l − gk,l+1 − gk+1,l= .
(g − g )2k,l+1 k,l + (g 2k+1,l − gk,l) + β | ‖∇sgk,l‖2 |β
Andra termen blir
√ gk,l − gk,l−1 gk,l − gk,l−1= ,
(g 2k,l − gk,l−1) + (gk+1,l−1 − g 2k,l−1) + β | ‖∇sgk,l−1‖2 |β
och tredje termen
√ gk,l − gk−1,l gk,l − gk−1,l= .
(g 2 2k,l − gk−1,l) + (gk−1,l+1 − gk−1,l) + β | ‖∇sgk−1,l‖2 |β
Vi f̊ar en gradient som inneh̊aller fyra förskjutningar där första termen representerar en förskjut-
ning ned̊at och till höger, andra termen representerar en förskjutning åt vänster och den tredje
termen representerar en förskjutning upp̊at. Som tidigare nämnts kommer vi att använda oss av
ett periodiskt randvillkor där sista indexet i en rad jämförs mot första indexet i samma rad. Vi
gör motsvarande för kolumnerna. TV-regulariseringens gradient, (14), med approximationen av
1-normen blir
[ ]n
∇ T − T 2gk,l − gk,l+1 − gk+1,l gk,l − gk,l−1 gk,l − gk−1,lHTV (g) = 2A Ag 2A b+λ + + .| ‖∇sgk,l‖2 |β | ‖∇sgk,l−1‖2 |β | ‖∇sgk−1,l‖2 |β k,l=1
Baserat p̊a gradienterna för Tikhonov- respektive TV-regularisering kan vi nu tillämpa opti-
meringsteorin för att hitta minimum för de fyra olika regulariseringsmetoderna. Vi har därmed
möjligheten att lösa v̊art inversproblem och ta fram en rekonstruktion av fantombilden. I nästa
kapitel g̊ar vi igenom hur vi har implementerat metoderna i MATLAB, och vilka antaganden om
parametrarna vi har gjort för att f̊a fram lösningar. Efter att vi har redogjort för implementeringen
kommer vi att presentera v̊ara resultat i kapitel 5.
4 Implementering
I detta kapitel g̊ar vi igenom implementeringen av optimeringslösarna, val av parametervärde och
presenterar en analys av regulariseringsparametrarna med hjälp av L-kurvmetoden. Alla parametrar
som har introducerats i teoridelen har implementerats under programmeringen.
11
4.1 Implementering av optimeringslösare
Implementering av optimeringslösarna har gjorts med hjälp av programmeringsspr̊aket MATLAB.
All kod kan ses i appendix. Vi har använt oss av Shepp-Loganfantomen som ursprungsfigur, se Figur
4. Denna är inbyggd i MATLAB och ger bättre förutsättningar för att kunna testa och jämföra
v̊ara resultat med den. Utifr̊an denna fantombild med upplösningen 256×256 pixlar simulerar vi
mätdata. Bilden har representerats av en matris med samma dimensioner som bildens upplösning,
nämligen en 256×256 matris.
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
Figur 4: MATLAB:s Shepp-Loganfantom. Fantombilden används för simulering av data och som
referensfigur.
Vi har programmerat optimeringslösare som använder sig av MATLAB:s befintliga funktio-
ner radon.m och iradon.m för rekonstruktion. Radontransformen har använts som fram̊atoperator,
medan AT är ofiltrerad bak̊atprojektion som tas fram med hjälp av iradon.m [mathworks.com].
Mätdatan har genererats utifr̊an fantombilden med funktionen radon.m som skickar 367 pa-
rallella str̊alar genom bilden. Detta görs för varje vinkel över 180◦. Startvinkeln är 1◦. Vidare är
vinklarna likformigt fördelade med 1◦ mellan varje par av vinklar. För att undvika inversbrott i
regulariseringen har vi förskjutit startpunkten med 0, 5◦, vilket innebär att 0, 5◦ är startpunkten
och vinklarna fördelar sig med 1◦ mellan varje par upp till 179, 5◦.
60
60
50
50
40 40
30 30
20 20
10
10
0
0
(a) Sinogram utan brus. (b) Sinogram med 5% brus.
Figur 5: En representation av simulerad data fr̊an Shepp-Loganfantomen, utan och med 5% brus.
Vi har lagt till normalfördelat brus, ϵ ∼ N(0, δ), där δ motsvarar 5% av standardavvikelsen fr̊an
den simulerade mätdatan. Mängden brus är valt baserat p̊a kvaliteten av rekonstruktionsbilderna.
12
Bruset läggs till för att imitera det faktiska brus som finns p̊a riktig mätdata. Mätdatan har
sammanställts i ett s̊a kallat sinogram, se Figur 5.
Fyra regulariseringar har tillämpats - Tikhonovregularisering och TV-regularisering, b̊ada med
och utan bivillkoret g ≥ 0. En approximation av ‖Lg‖2,1 och ett periodiskt randvillkor har
tillämpats för TV-regularisering. Optimeringsalgoritmerna är beroende av steglängden α. Vi har
valt steglängden α = 11000 . Denna är vald tillräckligt liten för att säkerställa konvergens. TV-
regulariseringens approximationsparameter β säkerställer differentierbarhet kring origo. Den är
vald s̊a att den inte p̊averkar negativt och är implementerad med värdet β = 12000 . Vi har imple-
menterat ett termineringsvillkor för optimeringslösarna i form av en toleransgräns p̊a 10−3. Detta
innebär att om ‖xk+1 − xk‖2 är mindre än 10−3 termineras programmet.
Alla regulariseringar är beroende av en regulariseringsparameter λ. Regulariseringsparametrar-
na har tagits fram med hjälp av en kombination av testning och L-kurvmetoden, som g̊as igenom
i nästa avsnitt.
4.2 Analys av regulariseringsparametrar med hjälp av L-kurvmetoden
Det finns flera olika metoder för att f̊a fram ett bra värde p̊a regulariseringsparameteren λ, men
ingen av dessa garanterar en optimal lösning [MS12, s. 72-76]. Vi har valt att utvärdera regu-
lariseringsparametrarna med hjälp av L-kurvmetoden som baseras p̊a att hitta en balans mellan
logaritmen av datamatchningstermen och logaritmen av regulariseringsfunktionalen. Denna me-
tod är ej konvergent utan ger endast en indikation för vad det optimala värdet av λ borde vara.
Med lämpliga parametervärden ger metoden en L-formad kurva med ett mjukt hörn där det op-
timala värdet bör ligga nära hörnet. Valet av detta parametervärde kan beskrivas som följer: om
λ > 0 är litet ges en lösning som mestadels är p̊averkad av logaritmen av datamatchningstermen,
log‖Ag − b‖22, eftersom lösningen blir straffad om den är för l̊angt ifr̊an den exakta lösningen. Vid
stora λ blir minimeringsproblemet straffat om logG(g) ökar. Vi kan därmed g̊a miste om en bra
lösning till optimeringsproblemet, d̊a regulariseringsparametern inte till̊ater en liten ökning av g.
[CMRS00, s. 425]
Figur 6: Illustration av en idealiserad L-kurva. Läge 1 representerar underregularisering, läge 3
representerar överregularisering och läge 2 indikerar en bra balans.
En idealiserad L-kurva är illusterad i Figur 6. Läge 1 i Figur 6 representerar lösningen när
λ är litet, allts̊a är lösningen underregulariserad vilket gör att regulariseringsfunktionalen blir
obetydlig, och vice versa för läge 3. Där är λ stort och lösningen blir mestadels p̊averkad av
regulariseringsfunktionalen medan datamatchningstermen istället blir obetydlig. Därav kan vi dra
slutsatsen att en bra balans av värdet p̊a regulariseringsparametern bör hittas vid läge 2 i Figur 6.
13
4.2.1 L-kurvillustrationer
Teorin om L-kurvor tillämpas här p̊a v̊ara regulariseringsmetoder. I Figur 7 sammanställs L-kurvor
för de olika regulariseringsmetoderna. Vid observation av de fyra L-kurvorna ser vi att Figur 7d,
som representerar TV-regularisering med gradientprojektionsmetoden, liknar mest en typisk L-
kurva. Därför tros en bra regulariseringsparameter för denna regularisering ligga nära hörnet p̊a
kurvan. Vi valde den fjärde punkten fr̊an höger p̊a Figur 7d, vilket motsvarar λ = 0.6. Även i Figur
7c, TV-regularisering med steepest descent-metoden, ser vi en L-kurvstruktur där ett värde runt
den sjunde punkten fr̊an höger, λ = 0.3, indikerar en bra balans mellan logaritmen av datamatch-
ningstermen och logaritmen av regulariseringsfunktionalen. L-kurvorna för Tikhonovregularisering
är presenterade i Figur 7a och 7b, där syns inte lika tydliga L-kurvor. Om vi istället observerar
värdena p̊a axlarna ser vi att datamatchningstermen förändras mer i proportion mot regularise-
ringsfunktionalen p̊a lösningen i Figur 7b. Av denna anledning väljer vi ett parametervärde som
ligger till vänster om mitten, vilket i detta fall motsvarar λ = 0.1. I Figur 7a resonerar vi p̊a
liknande sätt oss fram till punkt fyra fr̊an höger, λ = 0.15.
7.432 7.429
7.431 7.428
7.43 7.427
7.429 7.426
7.428 7.425
7.427 7.424
7.426 7.423
7.425 7.422
7.424 7.421
7.423 7.42
7.422 7.419
6.4 6.45 6.5 6.55 6.6 6.65 6.7 6.75 6.95 7 7.05 7.1 7.15 7.2 7.25 7.3 7.35 7.4 7.45
2 2
log||Ag-b|| log||Ag-b||
2 2
(a) Tikhonovregularisering utan bivillkor. (b) Tikhonovregularisering med bivillkor.
9.35 8.9952
8.995
9.3
8.9948
9.25
8.9946
9.2
8.9944
9.15
8.9942
9.1
8.994
9.05
8.9938
9 8.9936
6 6.5 7 7.5 8 8.5 9 9.5 6.5 7 7.5 8 8.5 9 9.5
2
log||Ag-b|| 2log||Ag-b||
2 2
(c) TV-regularisering utan bivillkor. (d) TV-regularisering med bivillkor.
Figur 7: L-kurvor för respektive regulariseringsmetod. Ett bra parametervärde bör ligga i hörnet
p̊a kurvan.
5 Resultat och jämförelse
I detta kapitel presenterar vi resultaten av rekonstruktionerna som tagits fram med de olika regula-
riseringsmetoderna. Alla bilder presenteras med varje metods individuella regulariseringsparameter
λ.
14
log||Lg|| 2
2,1 log||g||2
log||Lg|| 2
2,1 log||g||2
5.1 Resultat
Varje rekonstruerad bild presenteras med en färgskala fr̊an 0 till 1, vilket indikerar värdet p̊a
pixlarna, och därmed absorptionen. Fantombilden i Figur 4 representerar en bild av ett huvud.
Den vita delen representerar skelettet vilken har en hög absorptionskoefficient. De gr̊a och svarta
delarna indikerar mjukare vävnad med en lägre absorptionskoefficient. Färgskalan indikerar allt-
s̊a hur mycket en given vävnad absorberar. Det framg̊ar inte av bilderna, men b̊ade Tikhonov-
och TV-regularisering utan bivillkor inneh̊aller små negativa pixelvärden. De negativa värdena
framg̊ar därför inte av färgskalan som har ett minsta värde 0. I Figur 8 nedan ses de fyra olika
rekonstruktionerna av Shepp-Loganfantomen.
I översta raden i Figur 8 nedan presenteras rekonstruktionerna med Tikhonovregularisering.
Bilden till vänster, 8a, är rekonstruerad med hjälp av steepest descent-metoden och bilden till
höger, 8b, är rekonstruerad med hjälp av gradientprojektionsmetoden. Motsvarande uppdelning
gäller för rekonstruktionerna med TV-regularisering p̊a nedre raden.
1 1
0.8 0.8
0.6 0.6
0.4 0.4
0.2 0.2
0 0
(a) Tikhonovregularisering utan bivillkor, λ = (b) Tikhonovregularisering med bivillkor, λ =
0.15. 0.1.
1 1
0.8 0.8
0.6 0.6
0.4 0.4
0.2 0.2
0 0
(c) TV-regularisering utan bivillkor, λ = 0.3. (d) TV-regularisering med bivillkor, λ = 0.6.
Figur 8: Presentation av rekonstruktioner av Shepp-Loganfantomen framtaget med Tikhonov- och
TV-regularisering, b̊ada med och utan bivillkor, samt respektive metods regulariseringsparameter
λ.
I Figur 8a ses rekonstruktionen av bilden med Tikhonovregularisering utan bivillkor. Rekon-
struktionen har tydliga konturer och olika typer av vävnad syns. Bilden är kornig vilket gör det
sv̊art att tyda konturerna av de tre sm̊a ljusa ovalerna i nederdelen av fantomen. Bilden är tydligt
p̊averkad av brus, vilket vid närmare inspektion kan ses även i de svarta omr̊adena. Vidare s̊a
har Tikhonovregulariseringen med bivillkor i Figur 8b ett liknande beteende som Figur 8a. Tyd-
liga konturer ses, men bruset i bilden finns fortfarande kvar, även om det är n̊agot mildare än
15
rekonstruktionen utan bivillkor. I det svarta omr̊adet verkar bivillkoret ha sorterat bort en större
mängd brus. Denna förändring gör det lättare att avgöra form och storlek p̊a mindre omr̊aden. Om
vi kollar p̊a TV-regulariseringen utan bivillkor i Figur 8c skiljer sig denna jämfört med tidigare
rekonstruktioner. Konturerna framst̊ar tydligare och kanterna är skarpare och vitare. Brusigheten
i bilden har även reducerats s̊a att bilden upplevs slätare. Slutligen kan vi konstatera att TV-
regulariseringen med bivillkor i Figur 8d är mycket lik Figur 8c, framförallt om vi ser till skärpan
i kanterna. Brusniv̊an inuti fantomen är n̊agot lägre och upplevs som den bästa rekonstruktionen.
Sammanfattningsvis lyckas alla fyra rekonstruktioner återskapa fantombilden med olika kva-
litet av konturer och hantering av brus. De stora skillnaderna är skärpan i kanter och brusighet
inuti fantomen. Detta beteende kan förklaras med hur regulariseringarna är konstruerade. I TV-
regularisering används en spatiell gradient som jämför intilliggande pixlar med varandra. Denna
metod minimerar allts̊a med avseende p̊a avst̊anden till andra pixlar vilket p̊averkar den slutgiltiga
bilden där rekonstruktionen har slätat ut överg̊angen mellan olika vävnader.
5.2 Rekonstruktion genom bak̊atprojektion
En standardmetod för rekonstruktion av data fr̊an parallella str̊alar är bak̊atprojektion, [Her09,
s. 157], som antingen kan vara ofiltrerad eller filtrerad. Bak̊atprojektionens fördel är att den kan
beräknas fort, men nackdelen är att den inte kan hantera fall när datan inneh̊aller en stor mängd
brus. De tv̊a metoderna skiljer sig åt d̊a filtrerad bak̊atprojektion använder ett spatialfilter som
läggs p̊a sinogrammet innan bak̊atprojektionen genomförs [Sch 20, s. 222]. Detta filter filtrerar
bort höga frekvenser och beh̊aller l̊aga frekvenser, vilket gör rekonstruktionen mindre brusig. Vi
kommer inte att g̊a in djupare p̊a teorin bakom denna metod, utan använder den för jämförelse av
v̊ara rekonstruktionsresultat. I Figur 9 illustreras rekonstruktionen av Shepp-Loganfantomen med
5% brus genom ofiltrerad respektive filtrerad bak̊atprojektion med iradon.m.
1
70
0.9
65
0.8
60
0.7
55
0.6
50
0.5
45
0.4
40
0.3
35
0.2
30
0.1
25
0
(a) Ofiltrerad bak̊atprojektion. (b) Filtrerad bak̊atprojektion.
Figur 9: Ofiltrerad och filtrerad bak̊atprojektion vid rekonstruktion av Shepp-Loganfantom med
5% brus.
5.2.1 Jämförelse av resultat med bak̊atprojektion
Skalorna p̊a de tv̊a olika typerna av bak̊atprojektion är olika. Den ofiltrerade bak̊atprojektionen har
värden mellan 25 och 75 medan den filtrerade bak̊atprojektionens rekonstruktion har pixelvärden
mellan 0 och 1. Vit färg i fantomen och de fyra bilderna i Figur 8 motsvarar värde 1, vilket innebär
att Figur 9a hade varit helt vit i denna skalan.
Rekonstruktionen med ofiltrerad bak̊atprojektion i Figur 9a skiljer sig tydligt fr̊an rekonstruk-
tionerna i Figur 8 d̊a det inte g̊ar att urskilja konturer vilket gör den oanvändbar.
I Figur 9b, filtrerad bak̊atprojektion, ser vi ett resultat som omedelbart liknar rekonstruktio-
nen av Tikhonovregularisering utan bivillkor. Bilden är tydligt brusig med tydliga konturer kring
16
absorptionen av skelettet. Filtrerad bak̊atprojektion hanterar inte bruset lika bra som Tikhonovre-
gularisering med bivillkor, samt TV-regularisering med och utan bivillkor, d̊a det förekommer brus
i den svarta delen av rekonstruktionen.
6 Utveckling och slutsats
I detta kapitel diskuterar vi möjliga vidareutvecklingar av v̊art arbete och avslutar med en slutsats.
6.1 Utveckling
V̊art arbete har varit begränsat till fyra olika typer av variationell regularisering som vi har op-
timerat med tv̊a olika metoder. För vidare undersökning hade det varit intressant att ta in fler
typer av regularisering och jämföra dessa för att f̊a en mer övergripande först̊aelse för kvaliteten
och skillnaden hos de olika typerna av regularisering. Detta skulle medföra en bredare först̊aelse
kring regularisering vid rekonstruktion av bilder. En möjlighet är ocks̊a att lägga ett större fokus
p̊a valet av parametrar och känslighetsanalys av dessa.
En annan fr̊aga som kan vara intressant att undersöka är hur snabbt de olika lösningarna
konvergerar. Om vi skulle fokusera p̊a att optimera lösningarna med avseende p̊a tid kan kvaliteten
bli lidande. Det hade d̊a varit viktigt att hitta en bra balans mellan kvalitet och snabbhet.
Vi har ocks̊a i v̊art projekt valt att enbart fokusera p̊a data fr̊an en fantombild. Det kan därför
vara relevant att utöka arbetet till att inkludera flera och olika fantombilder för större komplexitet.
Dessutom undersöka hur v̊ar optimeringslösare skulle kunna hantera detta och undersöka vilken
p̊averkan detta skulle ha p̊a resultatet.
All v̊ar data har, som tidigare nämnts, simulerats utifr̊an parallella str̊alar. För framtida un-
dersökningar kan det vara relevant att skapa rekonstruktioner av bilder med hjälp av exempelvis
fan-beam str̊alar. Användning av fan-beam str̊alar kräver en annan uppställning av fram̊atopera-
torn d̊a str̊alningsgeometrin är annorlunda. I praktiken kräver skiktröntgen med fan-beam str̊alar
färre rotationer när röntgenstr̊alarna skjuts in. Detta innebär att processen g̊ar fortare, men är mer
kostsam än vid användning av parallella str̊alar. [Her09, s. 47-49]
Brus har lagts till p̊a all data, men vi kan inte garantera att det är representativt för riktig
data. Det kan därför vara intressant att undersöka andra former av brus för modelleringen av data,
exempelvis Poissonfördelat brus.
6.2 Samhälleliga och etiska aspekter
Som tidigare nämnts i inledningen kan det uppst̊a olika samhälleliga och etiska aspekter när det
kommer till optimering av rekonstruktionsmetoder, exempelvis kostnad för tagning av röntgenbil-
der, b̊ade för den enskilda personen men även för samhället. Vidare kan kontroll av sprickor i broar
och pelare ocks̊a ses som ett inversproblem och vi kan även där prata om för- och nackdelar med
detta. Ytterligare situationer där inversproblem kan uppst̊a är vid dekryptering där en fr̊aga om
integritet kan förekomma.
Metoderna behöver ha en bra balans mellan kvaliteten p̊a rekonstruktionen och tids̊atg̊angen.
En högre kvalitet kan medföra längre skannings- och rekonstruktionstid vilket kan leda till att färre
människor f̊ar tillg̊ang till denna teknik p̊a grund av kö eller kostnadsaspekter. Om skannings-
och rekonstruktionskvalitet istället väljs med fokus p̊a att minimera tids̊atg̊ang och kostnad kan
det innebära att fler människor har tillg̊ang, men att kvaliteten istället blir lidande. Arbetet har
använt simulerad data fr̊an MATLAB och har p̊a s̊a sätt ej använt n̊agon data som inneh̊aller
personuppgifter.
Optimeringsmr̊adet är i ständig utveckling, där metoderna för rekonstruktion förbättras och
skapar bilder av högre kvalitet. Exempelvis är ultraljud ocks̊a ett inversproblem där det kan uppst̊a
en relevant diskussion kring utvecklingspotentialen och konsekvenserna av denna utveckling. I takt
med att det finns möjlighet att generera ultraljudsbilder av bättre kvalitet kan större information
ges om exempelvis foster.
17
6.3 Slutsats
Vi har studerat omr̊aden inom inversproblem, optimering och regularisering. Utifr̊an de resultat
vi har f̊att fram kan vi dra slutsatsen att v̊ara optimeringslösningar konvergerar. Efter att ha un-
dersökt, med hjälp av L-kurvor, vilket λ som ger en bra balans mellan datamatchningstermen och
regulariseringsfunktionalen har vi kunnat ta fram rekonstruktioner av fantombilden. Rekonstruk-
tionerna lyckas till viss grad efterlikna fantombilden.
Vi har lyckats att ta fram fyra rekonstruktioner av fantombilden, vilket var syftet med detta
arbete. Det är därmed möjligt, baserat p̊a de uppställda regulariseringsmetoderna, att lösa invers-
problemet och ta fram absorptionsfunktionen f med hjälp av optimering.
18
Referenser
[1177.se] https://www.1177.se/Halland/behandling–hjalpmedel/undersokningar-och-
provtagning/bildundersokningar-och-rontgen/datortomografi/ , 2022-03-03, 14:50.
[AEP17] N. Andréasson, A. Evgrafov, M. Patriksson An Introduction to Continuous Optimization,
Studentlitteratur AB, 3rd edition, (2017)
[arpansa.gov.au] https://www.arpansa.gov.au/understanding-radiation/what-is-
radiation/ionising-radiation/x-ray, 2022-05-02, 10:24
[BLS16] B.Bajic, J. Lindblad, N. Sladoje, Restoration of images degraded by signal-dependent noise
based on energy minimization: an empirical study, Journal of Electronic Imaging, Vol. 25(4),
(Jul/Aug 2016)
[Bom08] J. Boman, Datortomografins matematik Om en matematisk teori med m̊anga nya tillämp-
ningar, Matematiska Institutionen Stockholms Universitet, (2008), s. 177-186
[CMRS00] D. Calvetti, S. Morigi, L. Reichel, F. Sgallari, Tikhonov regularization and the L-curve
for large discrete ill-posed problems, Journal of Computational and Applied Mathematics
123, (2000), s. 423446
[Fee15] T. G. Feeman, The Mathematics of Medical Imaging, Springer International Publishing,
(2015), s. 1-35
[Her09] G. T. Herman, Fundamentals of Computerized Tomography, Springer International
Publishing, (2009)
[Had1923] J. Hadamard, Lectures on Cauchy’s Problem in Linear Partial Differential Equations,
Dover Publications, New York, edition anno 1952, (1923)
[humanprogress.org] https://www.humanprogress.org/heroes-of-progress-pt-26-wilhelm-rontgen/
, 2022-05-09, 14:15
[Ker16] M. Kern, Numerical Methods for Inverse Problems, John Wiley & Sons, NY, 1st edition
(2016)
[Kje17] T. Kjellström, Medicinhistoria Röntgen genom tiderna,
https://www.doktorn.com/artikel/medicinhistoria-r%C3%B6ntgen-genom-tiderna/, 2022-
05-06, 12:30
[LLM16] D. C. Lay, S. R. Lay, J.J. McDonald, Linear Algebra and its Applications, Pearson Educa-
tion Limited, 5th edition, (2016)
[mathworks.com] https://se.mathworks.com/help/images/ref/radon.html , 2022-03-22, 13:07
[mathworks.com] https://se.mathworks.com/help/images/ref/iradon.html, 2022-05-02, 11:34
[MS12] J. L. Mueller, S. Siltanen, Linear an Nonlinear Inverse Problem with Practical Applications,
SIAM, (2012)
[ne.se1] Nationalencyklopedin, lysämne. https://www.ne.se/uppslagsverk/ordbok/svensk/lys%C3
%A4mne, 2022-05-09, 13:50
[ne.se2] Nationalencyklopedin, regularisering. http://www.ne.se/uppslagsverk/encyklopedi/l̊ang/reg
ularisering, 2022-05-03, 10:10
[ne.se3] Nationalencyklopedin, röntgenstr̊alning. http://www.ne.se/uppslagsverk/encyklopedi/l̊ang
/röntgenstr̊alning, 2022-05-02, 15:45
[Niev1986] Y. Nievergelt, Elementary inversion of Radon’s Transform, Society for Industrial and
Applied Mathematics, SIAM review vol. 28 No. 1, (1986), s. 79-84
[Sch 20] R. Schofield et. al, Image reconstruction: Part 1 understanding filtered back projection,
noise and image acquisition, Journal of Cardiovascular Computed Tomography 14, (2020),
s. 219-225
19
A Appendix
A.1 Bevis
A.1.1 Sats 3.4
Bevis av sats 3.4. [AEP17, s. 82]
Antag att x∗ är ett lokalt minimum, men inte ett globalt. Betrakta en vektor x ∈ S med egenskapen
att f(x) < f(x∗). L̊at λ ∈ (0, 1). D̊a S och f är konvexa gäller det att λx + (1 − λ)x∗ ∈ S, och
f(λx + (1 − λ)x∗) ≤ λf(x) + (1 − λ)f(x∗). Genom att välja λ > 0 tillräckligt liten f̊ar vi en
motsägelse om lokal optimalitet av x∗. □
A.1.2 Sats 3.5
Bevis av sats 3.5. , [AEP17, s.91-94]
⇒)
Antag att x∗ är ett lokalt minimum, men att ∇f(x∗) 6= 0n. L̊at p := −∇f(x∗), och undersök
Taylorexpansionen runt x = x∗ i p′s riktning:
f(x∗ + αp) = f(x∗) + α∇f(x∗)Tp+ o(α),
där 0 : R → R är s̊adan att o(s)/s → 0as när s → 0. Vi f̊ar
f(x∗ + αp) = f(x∗)− α||∇f(x∗)||2 + o(α) < f(x∗)
för alla α > 0 små nog, d̊a ||∇f(x∗)|| 6= 0. Detta avslutar första delen av beviset.
⇐)
D̊a f är konvex gäller det för varje y ∈ Rn att,
f(y) ≥ f(x∗) +∇f(x∗)T (y − x∗) = f(x∗),
där likheten följer fr̊an att ∇f(x∗) = 0n. Beviset klart. □
A.1.3 Normens konvexitet
Bevis normens konvexitet. Vi har fr̊an triangelolikheten att
||x+ y|| ≤ ||x||+ ||y||.
D̊a kan vi skriva
||λx+ (1− λ)y|| ≤ ||λx||+ ||(1− λ)y|| = λ||x||+ (1− λ)||y||.
□
20
A.2 MATLAB-kod
A.2.1 Sinogram
%% Sinogram utan och med 5% brus
clf ; clear ; clc ;
xd = 256; % Bildens upplosning.
obj = phantom(xd); % Fantombilden med dimension
256 x256
angle = 1:1:180; % Inskjutningsvinklar.
data = radon(obj , angle); % Radontransformen av obj.
v = 0.05* std(data (:)); % 5% brus.
datarad = data + v*randn (367, 180); % Tillagt 5% brus.
% Plot sinogram utan brus
subplot(1, 2, 1)
imshow(data ,[]), axis( ' square ' )
colorbar
% Plot sinogram med 5% brus
subplot(1, 2, 2)
imshow(datarad ,[]), axis( ' square ' )
colorbar
21
A.2.2 Tikhonovregularisering med steepest descent, L-kurva
%% Steepest Descent
clf;clear;clc;
xd = 256; % Dimensioner bild
alpha = 1/1000; % Steglangden
object = phantom(xd); % Fantombilden 256 x256
angle = 1:1:180; % Inskjutningsvinklar simulerad data
rec_angle = 0.5:1:179.5; % Inskjutningsvinklar rekonstruktion
data = radon(object ,angle); % Simulerad data fra fantom
noise = 0.05* std(data (:)); % Genererar brus
data = data + noise*randn (367 ,180); % Simulerad data med brus
tol = 1e-3; % Tolerans terminering
% Varden pa regulariseringsparametern som testas
lambda = [10^-3 3*10^ -3 10^-2 3*10^ -2 10^-1 1.5*10^ -1 2*10^ -1
2.5*10^ -1 3*10^ -1];
% Gradient som bestammer stegriktning
grad = @(x, angle , data , lambda) 2*( iradon(radon(x,angle)-data ,
angle , ' none ' ,xd) + lambda*x);
% Initierar tomma vektorer och celler som sedan lagrar resultat
Lg_norm = zeros(1,length(lambda));
Ag_b_norm = zeros(1,length(lambda));
saved_pics = cell(1,length(lambda));
for i = 1: length(lambda)
pic = zeros(xd); % Skapa bild med startpunkt 0
while true
% Stegriktning
df = grad(pic ,rec_angle ,data ,lambda(i));
tmp_pic = pic - alpha*df; % Ny losning
% Kollar om normen av skillnaden mellan x_k+1 - x_k
% ar mindre an toleransnivan
pic_diff = sqrt(sum(( tmp_pic (:)-pic(:)).^2));
if pic_diff < tol
pic = tmp_pic;
break
end
% Uppdaterar och sparar var nya losning
pic = tmp_pic;
saved_pics{i} = pic;
end
%Log av reg.funktional
Lg_norm(i) = log(norm(pic)^2);
%Log av datamatchning
22
Ag_b_norm(i) = log(norm(radon(pic , rec_angle)-data)^2);
end
plot(Ag_b_norm ' , Lg_norm ' , ' .- ' ) %Plot L-kurva
xlabel ("log||Ag-b||_2^2") %x-label
ylabel ("log||g||_2^2") %y-label
23
A.2.3 Tikhonovregularisering med gradientprojektion, L-kurva
%% Gradientprojektion
clf;clear;clc;
xd = 256; % Dimensioner bild
alpha = 1/1000; % Steglangden
object = phantom(xd); % Fantombilden 256 x256
angle = 1:1:180; % Inskjutningsvinklar simulerad data
rec_angle = 0.5:1:179.5; % Inskjutningsvinklar rekonstruktion
data = radon(object ,angle); % Simulerad data fra fantom
noise = 0.05* std(data (:)); % Genererar brus
data = data + noise*randn (367 ,180); % Simulerad data med brus
tol = 1e-3; % Tolerans terminering
% Varden pa regulariseringsparametern som testas
lambda = [10^-3 3*10^ -3 10^-2 3*10^ -2 10^-1 0.15 0.2 0.25 3*10^ -1];
% Gradient som bestammer stegriktning
grad = @(x, angle , data , lambda) 2*( iradon(radon(x,angle)-data ,
angle , ' none ' ,xd) + lambda*x);
% Initierar tomma vektorer och celler som sedan lagrar resultat
Lg_norm = zeros(1,length(lambda));
Ag_b_norm = zeros(1,length(lambda));
saved_pics = cell(1,length(lambda));
for i = 1: length(lambda)
pic = zeros(xd); % Skapa bild med startpunkt 0
while true
% Stegriktning
df = grad(pic ,rec_angle ,data ,lambda_test);
% Uppdaterar positionen ( bilden )
tmp_pic = pic - alpha*df;
% Projicerar alla bilder med negativt varde till 0
tmp_pic(tmp_pic < 0) = 0;
% Kollar om normen av skillnaden mellan x_k+1 - x_k
% ar mindre an toleransnivan
pic_diff = sqrt(sum(( tmp_pic (:)-pic(:)).^2));
if pic_diff < tol
pic = tmp_pic;
break
end
% Uppdaterar och sparar var nya losning
pic = tmp_pic;
saved_pics{i} = pic;
end
% Log av reg.funktional
24
Lg_norm(i) = log(norm(pic)^2);
% Log av datamatchning
Ag_b_norm(i) = log(norm(radon(pic , rec_angle)-data)^2);
end
plot(Ag_b_norm ,Lg_norm , ' .- ' ) % Plot L-kurva
xlabel ("log||Ag-b||_2^2") % x-label
ylabel ("log||g||_2^2") % y-label
25
A.2.4 Total variationsregularisering med steepest descent, L-kurva
%% TV Gradient Projection , L-curve
clf;clear;clc;
xd = 256; % Dimensioner bild
beta = 0.0005; % Approximeringsparameter
alpha = 1/1000; % Steglangden
object = phantom(xd); % Fantombilden 256 x256
angle = 1:1:180; % Inskjutningsvinklar simulerad data
rec_angle = 0.5:1:179.5; % Inskjutningsvinklar rekonstruktion
data = radon(object ,angle); % Simulerad data fra fantom
noise = 0.05* std(data (:)); % Genererar brus
data = data + noise*randn (367 ,180); % Simulerad data med brus
tol = 1e-3; % Tolerans terminering
% TV -regularisering gradient
grad_data = @(x, angle , data) 2* iradon(radon(x,angle)-data ,angle , '
none ' ,xd);
grad_spec_norm = @(x, y, beta) ((x.^2 + y.^2 + beta).^(1/2));
% Olika lambda till undersokning
lambda = [3*10^ -3 10^-2 3*10^ -2 10^-1 3*10^ -1 4*10^ -1 5*10^ -1
6*10^ -1 7*10^ -1 1 3];
% Tomma vektorer och cella for att spara undan information
Lg_norm = zeros(1,length(lambda));
Ag_b_norm = zeros(1,length(lambda));
saved_pics = cell(1,length(lambda));
tmp_left_shift = zeros(xd);
tmp_right_shift = zeros(xd);
tmp_up_shift = zeros(xd);
tmp_down_shift = zeros(xd);
tmp_up_right_shift = zeros(xd);
tmp_left_down_shift = zeros(xd);
for i = 1: length(lambda)
pic = zeros(xd); % Skapa bild med startpunkt 0
while true
% Forskjutning , left
tmp_left_shift (:,1:end -1) = pic(:,2:end);
tmp_left_shift (:,end) = pic(:,1);
% Forskjutning , right
tmp_right_shift (:,2:end) = pic(:,1:end -1);
tmp_right_shift (:,1) = pic(:,end);
% Forskjutning , up
tmp_up_shift (1:end -1,:) = pic(2:end ,:);
tmp_up_shift(end ,:) = pic(1,:);
% Forskjutning , down
26
tmp_down_shift (2:end ,:) = pic(1:end -1,:);
tmp_down_shift (1,:) = pic(end ,:);
% Forskjutning , left down
tmp_left_down_shift (2:end ,:) = tmp_left_shift (1:end -1,:);
tmp_left_down_shift (1,:) = tmp_left_shift(end ,:);
% Forskjutning , up right
tmp_up_right_shift (1:end -1,:) = tmp_right_shift (2:end ,:);
tmp_up_right_shift(end ,:)= tmp_right_shift (1,:);
% Implementerar forskjutningarna
left_shift = pic - tmp_left_shift;
right_shift = pic - tmp_right_shift;
up_shift = pic - tmp_up_shift;
down_shift = pic - tmp_down_shift;
left_down_shift= tmp_left_down_shift - tmp_left_shift;
up_right_shift = tmp_up_right_shift - tmp_up_shift;
% Berakna gradienten
grad_tot = grad_data(pic , rec_angle , data) + lambda(i)*((
right_shift + down_shift)./( grad_spec_norm(down_shift ,
right_shift , beta)) + up_shift ./( grad_spec_norm(up_shift
, up_right_shift , beta)) + left_shift ./( grad_spec_norm(
left_shift , left_down_shift , beta)));
% Uppdaterar positionen (bilden)
tmp_pic = pic - alpha*grad_tot;
% Beraknar normen av den nya och gamla position
pic_diff = sqrt(sum(( tmp_pic (:)-pic(:)).^2));
% Kollar toleransen , om mindre an toleransgrans terminer
if pic_diff < tol
pic = tmp_pic;
iterationer = iterationer + 1;
break
end
% Uppdaterar och sparar var nya losning
pic = tmp_pic;
saved_pics{i} = pic;
end
Lg_norm(i) = log(sum(sum(abs(pic)))); % Log av reg.funktional
Ag_b_norm(i) = log(norm(radon(pic , rec_angle)-data)^2); % Log
datamatchning
end
plot(Ag_b_norm ,Lg_norm , ' .- ' ) % Plot L-kurva
xlabel ("log||Ag-b||_2^2") % x-label
ylabel ("log||Lg||_{2 ,1}") % y-label
27
A.2.5 Total variationsregularisering med gradientprojektion, L-kurva
%% TV Gradient Projection , L-curve
clf;clear;clc;
xd = 256; % Dimensioner bild
beta = 0.0005; % Approximeringsparameter
alpha = 1/1000; % Steglangden
object = phantom(xd); % Fantombilden 256 x256
angle = 1:1:180; % Inskjutningsvinklar simulerad data
rec_angle = 0.5:1:179.5; % Inskjutningsvinklar rekonstruktion
data = radon(object ,angle); % Simulerad data fra fantom
noise = 0.05* std(data (:)); % Genererar brus
data = data + noise*randn (367 ,180); % Simulerad data med brus
tol = 1e-3; % Tolerans terminering
% TV -regularisering gradient
grad_data = @(x, angle , data) 2* iradon(radon(x,angle)-data ,angle , '
none ' ,xd);
grad_spec_norm = @(x, y, beta) ((x.^2 + y.^2 + beta).^(1/2));
% Olika lambda till undersokning
lambda = [3*10^ -3 10^-2 3*10^ -2 10^-1 3*10^ -1 4*10^ -1 5*10^ -1
6*10^ -1 7*10^ -1 1 3];
% Tomma vektorer och celler for att spara undan information
Lg_norm = zeros(1,length(lambda));
Ag_b_norm = zeros(1,length(lambda));
saved_pics = cell(1,length(lambda));
tmp_left_shift = zeros(xd);
tmp_right_shift = zeros(xd);
tmp_up_shift = zeros(xd);
tmp_down_shift = zeros(xd);
tmp_up_right_shift = zeros(xd);
tmp_left_down_shift = zeros(xd);
for i = 1: length(lambda)
iterationer = 0; % Antal iterationer = 0
pic = zeros(xd); % Skapa bild med startpunkt 0
while true
% Forskjutning , left
tmp_left_shift (:,1:end -1) = pic(:,2:end);
tmp_left_shift (:,end) = pic(:,1);
% Forskjutning , right
tmp_right_shift (:,2:end) = pic(:,1:end -1);
tmp_right_shift (:,1) = pic(:,end);
% Forskjutning , up
tmp_up_shift (1:end -1,:) = pic(2:end ,:);
tmp_up_shift(end ,:) = pic(1,:);
28
% Forskjutning , down
tmp_down_shift (2:end ,:) = pic(1:end -1,:);
tmp_down_shift (1,:) = pic(end ,:);
% Forskjutning , left down
tmp_left_down_shift (2:end ,:) = tmp_left_shift (1:end -1,:);
tmp_left_down_shift (1,:) = tmp_left_shift(end ,:);
% Forskjutning , up right
tmp_up_right_shift (1:end -1,:) = tmp_right_shift (2:end ,:);
tmp_up_right_shift(end ,:)= tmp_right_shift (1,:);
% Implementerar forskjutningarna
left_shift = pic - tmp_left_shift;
right_shift = pic - tmp_right_shift;
up_shift = pic - tmp_up_shift;
down_shift = pic - tmp_down_shift;
left_down_shift= tmp_left_down_shift - tmp_left_shift;
up_right_shift = tmp_up_right_shift - tmp_up_shift;
% Berakna gradienten
grad_tot = grad_data(pic , rec_angle , data) + lambda(i)*((
right_shift + down_shift)./( grad_spec_norm(down_shift ,
right_shift , beta)) + up_shift ./( grad_spec_norm(up_shift
, up_right_shift , beta)) + left_shift ./( grad_spec_norm(
left_shift , left_down_shift , beta)));
% Uppdaterar positionen (bilden)
tmp_pic = pic - alpha*grad_tot;
% Projicerar alla bilder med negativt varde till 0
tmp_pic(tmp_pic < 0) = 0;
% Beraknar normen av den nya och gamla position
pic_diff = sqrt(sum(( tmp_pic (:)-pic(:)).^2));
% Kollar toleransen , om mindre an toleransgrans terminer
if pic_diff < tol
pic = tmp_pic;
break
end
pic = tmp_pic;
iterationer = iterationer + 1; % Iterationer
saved_pics{i} = pic; % Sparar undan bild for
lambda(i)
end
Lg_norm(i) = log(sum(sum(abs(pic)))); % Log av reg.funktional
Ag_b_norm(i) = log(norm(radon(pic , rec_angle)-data)^2); % Log
datamatchning
end
plot(Ag_b_norm ,Lg_norm , ' .- ' ) % Plot L-kurva
xlabel ("log||Ag-b||_2^2") % x-label
ylabel ("log||Lg||_{2 ,1}") % y-label
29