Tomografi

Röntgen-skannings princippet


I 1970-erne var tiden moden til at skabe et apparatur som man længe havde været sikker på var en teknisk mulighed: et röntgenapparat som gør det muligt at se rigtigt ind i legemer som er uigennemtrængelige for lys - i modsætning til det gamle röntgenbillede hvor det kun er en projektion af det indre man kan se. Problemet er af to-dimensional natur, siden man kun behøver at danne billedet af et plant snit igennem legemet. Og dette må nødvendigvis ske ved at sende röntgen- eller lignende stråler igennem legemet fra et stort antal retninger beliggende i den pågældende plan, og for hver retning måle hvormeget strålens intensitet aftager, og så udfra alle disse informationer konstruere et billede. Men hvordan? Ved det sædvanlige röntgenapparat udsendes stråler fra röntgenrøret i alle retninger indenfor en kegle, og disse kan opfanges på en fotografisk plade eller på en fluorescerende glasplade. Men nu skal der tilsyneladende også matematik til for at danne billedet: informationerne skal omdannes til tal som skal bearbejdes i formler, og udfra de bearbejdede tal skal man få et billede frem på en skærm. Dette sidste var ikke noget problem, det havde man tv-skærme til, men det første problem kunne muligvis være formidabelt: der skal ideelt set benyttes informationer fra "alle tænkelige" linier indenfor det plane snit, og da en linie er karakteriseret ved en vinkel og en forskydning, er der måske i praksis brug for 300x300 = 90.000 linier, og for enhver af disse fås et tal (intensitetstabet) som skal opbevares da det skal bruges for alle punkterne på skærmen. Og denne skærm er måske på 400x400 = 160.000 punkter. Altså: 90.000 tal skal oplagres i en hukommelse og bearbejdes for hvert af de 160.000 punkter.

Og dog var det ikke dette oplagrings- og regneproblem man skulle slås mest med i 70-erne. Man rådede allerede dengang over computerkraft som var stor nok, for man kan nemlig lade nogle af operationerne foregå parallelt. Nej, det var matematikken: man kendte ikke de formler der skulle bruges - har jeg hørt. To fysikere fik i 1979 Nobelprisen for deres arbejde med denne nye teknik, og det var ikke mindst for at have opstillet formlerne. Men kan dette virkelig være sandt? Matematikerne har alle dage hævdet, at de fra deres elfenbenstårn har et glimrende udsyn over verden og at de er langt forud for denne i udvikling. Og når nogen har sat spørgsmål ved nytten af deres arbejde, har de altid forsvaret dette med at det engang vil vise sig nyttigt. Og her var endda tale om et problem hvis løsning måtte ligge langt under deres niveau.

Og selvfølgelig var et så elementært problem som der her er tale om, løst forlængst. Det er et centralt problem indenfor det fagområde som hedder integralgeometri, og denne disciplin havde eksisteret i over hundrede år. Og i 1917 offentliggjorde den østriske matematiker J.K.A. Radon (1887-1956) en afhandling, hvori han fremlagde og beviste de formler der skal til for at rekonstruere en funktion defineret på et område i planen udfra dens integral langs alle tænkelige linier igennem området. Og hvis denne funktion er tætheden (i en eller anden forstand) af et materiale fordelt over det plane område, så må integralet langs en ret linie være et mål for den modstand som en stråle der er følsom overfor materialet møder, når den trænger igennem materialet langs linien, det vil sige den aftagen i intensitet der sker i en strålen.

Radon har dog næppe tænkt den tanke at hans teori i praksis ville kunne bruges til at se ind i legemer. Röntgenbilleder havde man ganske vist haft i tyve år, men fra denne matematikløse teknik til hans regnekrævende formler er der lang vej. Radon havde udtænkt sin formel fordi den på hans tid var en oplagt udfordring og et naturligt trin i integralgeometriens udvikling - og måske viste formlen og dens bevis sig ligefrem at besidde en stor skønhed. Om formlen må man sige, at den ser ud som enhver erfaren matematiker vil forvente: den er enkel og ganske køn. Det er mere beviset som den matematikinteresserede hungrer efter at se, og her vil han blive glædelig overrasket. Men inden vi skitserer teorien skal vi formulere Radons formel og efterprøve den på computeren.

Hvis vi lader f(p) være materialets tæthed i punktet p og D være områdets diameter, lyder Radons formel således:

      For ethvert p defineres en funktion M(p, r) af de positive reelle tal r ved:

M(p, r) = middelværdien af modstanden i alle de rette linier som har afstand r fra p.

      Da gælder:

f(p) = -1/π ∫ (dM(p, r)/dr) dr/r (r fra 0 til D).

Der skal laves to programmer: ét der skanner og ét der gendanner. Vi vil forestille os at vi indenfor et rektangulært område har en fordeling af masse, men dette svarer til et billede i gråtoner, idet enhver tæthed i massen tildeles en gråtone. Og skanningen foregår ved at man for "enhver" tænkelig ret linie igennem billedet, udregner modstanden (intensitetstabet) langs linien, hvilket er det samme som at udregne integralet af gråtonetallene langs linien. Da det er utænkeligt at tage enhver tænkelig ret linie i betragtning, må vi udvælge et antal linier som i praksis repræsenterer alle tænkelige linier. Og dette kan foregå ved, at man ud fra centrum i billedet lader udgå rette linier ligeligt i et "stort" antal retninger, og for enhver sådan retning med "små" ligeligt beliggende mellemrum, tager normalen til retningen - altså forskyder normalen langs retningen. Resultatet er en matrix af reelle tal, hvis dimension er antallet af vinkler gange antallet af forskydninger. Gendannelsen af gråtonebilledet foregår ved Radons formel for f(p).

Men hvorfor ikke, i stedet for et trist gråtonebillede, benytte et farvebillede? En farve på en computerskærm er sammensat af tre farver, nemlig rød, grøn og blå, og hver af disse har en styrke som er et helt tal gående fra 0 til 255 - altså et tal af typen byte. Dette betyder at skærmen har 2563 = næsten 17 millioner forskellige farver. Farvebilledet er altså sammensat af tre billeder, som er "gråtonebilleder" i henh. rød, grøn og blå. For hver af disse tre "gråtonebilleder" kan vi danne en matrix af relle tal, og ved gendannelsen får vi for et punkt p i billedet tre tal, som vi nu lader være en farves bestandel af rød, grøn og blå, således at det vi gendanner er at farvebillede - som vi skal sammenligne med det originale farvebillede.

Der bliver nok nogle fejl i det gendannede billede, da vi jo har med integrationer og differentationer at gøre, og disse kun kan udføres med tilnærmelse. Og skal vi have en stor præcision vil programmet nok køre ret langsomt. På den anden side, da Radons formel jo er helt eksakt og da en moderne computer er hurtig, kan vi godt tillade os at stille store krav til billedet. At det måske vil tage en god tid at genskabe billedet, er irriterende, med det har som sagt ingen praktisk betydning, da flere af operationerne kan foregå parallelt, således at man kan lade flere operationsenheder køre samtidig.

Vi vil først udtænke et program som er simplest muligt, idet vi udfører integrationerne og differentationerne lige efter lærebogen, det vil sige tager sådanne tilnærmelser som optræder i deres definition. Så vil vi undervejs opdage, at vi med enkle midler kan forbedre programmet.

Resultatet af en skanning (for en farve) er en matrix af tal bestemt ved en retning (vinkel) og en afstand fra billedets centrum, og disse tal er hver især et integral langs normalen til retningen i afstanden fra centrum:



Denne to-dimensionale funktion vil vi kalde modstand(vinkel, afstand) - vinkel opfattes altså som x-koordinaten og afstand opfattes som y-koordinaten. For normalen bestemt ved (vinkel, afstand) skal vi udføre en integration af talværdier langs denne rette linie, derfor må vi kende koordinatsættene (i, j) for punkterne langs normalen. For et givet i svarer et j - og omvendt. Om det er i eller j der skal vælges som uafhængig variabel, kommer an på hvordan normalen ligger i forhold til koordinatakserne: Vi vælger i (x-koordinaten) for de normaler hvis vinkel er nærmest y-aksen (|cos(vinkel)| < |sin(vinkel)|), og j (y-koordinaten) for de normaler hvis vinkel er nærmest x-aksen (|cos(vinkel)| > |sin(vinkel)|). Men afstanden fra et koordinatsæt til det næste langs normalen er ikke éns for alle vinkler: jo mere skrå normalen ligger i forhold til henh. x-aksen og y-aksen, jo længere afstand er der imellem de udvalgte punkter på den. Derfor skal vi i det første tilfælde dividere med |sin(vinkel)|, og i det andet tilfælde dividere med |cos(vinkel)|.

Hvor mange vinkler og afstande vi skal medtage når vi skanner, kommer an på billedets størrelse og den kvalitet af det gendannede billedet vi ønsker. Men vi kan roligt vælge mange vinkler og afstande, da skanningen er ret hurtig - mindst 400 vinkler og afstande.

Ved gendannelsen skal vi først finde funktionen M(p, r) for ethvert punkt p (pixel) i billedet og "enhver" afstand r fra p - en sådan afstand kan højest være rektanglets diameter D = 2R. Og udfra M(p, r) skal vi udregne integralet

∫ (dM(p, r)/dr) dr/r (r fra 0 til D).

For et givet r ≥ 0 er M(p, r) summen af modstanden i alle de normaler til de udvalgte retninger som har afstanden r fra punktet p divideret med antallet af disse linier, altså antallet af vinkler. For enhver udvalgt vinkel v skal vi integrere tætheden langs den normal til retningen v, hvis afstand fra p er r. Hvis punktet p har koordinatsættet (a, b) er denne normal den samme som den normal (bestemt ved vinklen v) som har afstanden r1 = r + a·cos(v) + b·sin(v) fra origo (centrum), idet a·cos(v) + b·sin(v) er projektionen af linien fra origo til p ind på retningen v:



- med mindre at r1 bliver negativ, for så skal vi lægge π til v, da vi har valgt at operere med vinkler hele vejen rundt, men kun med positive afstande.

Herefter skal vi udregne Radons integral. I integranden ophæver de to dr-er hinanden, så vi behøver blot dele intervallet fra 0 til D op i et stort antal lige store stykker, og hvis vi kalder længden af hvert stykke d, skal vi for stykket der går fra r til r + d udregne tallet

(M(p, r + d) - M(p, r))/(r + d/2),

og disse tal skal adderes. Hvis summen kaldes s, skal vi i programmet først sætte s = 0 og r = 0, og succesivt gøre følgende indtil r > D: udregne (M(p, r + d) - M(p, r))/(r + d/2), addere dette tal til s, og derefter erstatte r med r + d. Til sidst skal vi gange s med -1/π og afrunde det tal vi får til et helt tal: dette tal er gråtonen i pixlen p - eller en af styrkegraderne for rød, grøn og blå hvis billedet er i farver.

Da dette primitive program kan forbedres betydeligt ved beskedne ændringer, vil vi straks gå igang hermed.

Den inddeling vi i udregningen af integralet i Radons formel foretager af intervallet [0, D], bør være finere jo nærmere vi kommer til 0, da integranden må være størst her, da der er en division med r. Vi vil begynde med et meget lille d, og succesivt gøre d større ved at erstatte d med kd, hvor k er et fast tal lidt større end 1 (f.eks. k = 1.1). For at danne M(p, r) skal vi for enhver udvalgt vinkel v bruge modstanden i den udvalgte afstand fra origo som er nærmest r1. Og her bør vi indføre en interpolation: hvis r1 og r2 er de udvalgte afstande på hver side af r1, bør vi lade modstand(v, r1) være

modstand(v, r1) + (r1 - r1)(modstand(v, r2) - modstand(v, r1))/(r2 - r1).

Dannelsen af en differenskvotient der er en tilnærmelse til dM(p, r)/dr i punktet r, foretages udfra de successivt forøgede d-værdier, som er dr. Men istedet for at lade bidraget til summationen fra r til r + d være d gange (dM (p, r')/dr')/r' taget i et punkt r' imellem r og r + d, vil vi benytte tilnærmelsesformlen:

∫ f(s) ds/s (s fra r til r + d) ≈ (f(r)/(r + d/4) + f(r + d)/(r + 3d/4)) d/2,

som bygger på værdierne af f(s) i de to endepunkter af det lille interval [r, r + d] (den er ikke svær at eftervise).

Desuden behøver vi i integrationen ikke at integrere helt til D, som er rektanglets diameter, vi kan nøjes med at integrere til R + afstanden fra origo til p.



Programmernes virkemåde


Download disse to programmer:     Scan     Recreate


   


Vi forudsætter at billedet er i BMP-format og ikke for stort (under 800 pixels i bredde og højde), og det gives navnet "pict". I programmet "Scan" skal indtastes antallet af vinkler og forskydninger (f.eks. 400-800 af hver). Resultatet af skanningen er en matrix af reelle tal, hvis første variabel går igennem antallet af vinkler og hvis anden variabel går igennem antallet af forskydninger - og tre sådanne matricer hvis billedet er i farver. Men et billede jo er en bredde gange højde-matrix af bytes, hvis billedet er i gråtoner, og tre sådanne, hvis billedet er i farver. Så derfor kan vi omdanne resultatet af skanningen til et billede, ved justering og afrunding af de reelle tal således at tallene bliver bytes. Dette billede, hvis bredde er antallet af vinkler og hvis højde er antallet af forskydninger, er givet navnet "trans_pict". Programmet "Recreate" tegner nu, udfra informationerne i dette billede, et billede som gerne skulle være meget tæt på "pict", og som er givet navnet "pict0". Dette billede behøver ikke nødvendigvis at have samme størrelse som det oprindelige billede ("pict"), så derfor indtastes bredde og højde af det ønskede billede ("pict0"). Desuden må indtastes en "lysfaktor", da resultatet af skanningen er relative tal. Den er et reelt tal omkring 1 - lav først et lille billede med lysfaktor 1, og juster derefter tallet til for eksempel 0.7. Et vindue viser gendannelsen som den skrider frem, trykkes på en tast lukker programmet og efterlader et billede ("pict0") med det hidtil tegnede.


   


Billedet til venstre er det "trans_pict" (vist i halv størrelse, da det er på 800x800 pixels) som er brugt i gendannelsen ovenfor (det højre billede). Billedet til højre er fået fra det venstre ved at forskyde farveværdierne ret kraftigt bort fra den middelgrå farve - der er kommet mere lyst og mere mørkt. Vi vil nu bruge dette i gendannelsen, men da der ved en sådan "ødelæggelse" af det skannede billede er en tendens til at det gendannede billede bliver mørkere i hjørnerne, har vi ladet billedet blive halvanden gange så stort og beskåret:



Vi ser at farveforskydningen i det skannede billede ikke har stor indflydelse på det gendannede billede - blot lidt øget kontrast.

Ved at tage udgangspunkt i et nogenlunde enkelt og interessant computerfremstillet billede, kan man ved skanning lave et "abstrakt kunstværk":


   


Et billede med bratte overgange og ensfarvede områder kan bruges til at teste programmet:


   


Hvis vi skanner et billede igen og igen, må vi tilsidst komme frem til et billede hvis skanning er identisk med billedet selv. For et helt hvidt billede ser skanningen ud som billedet til venstre (hvor kommer de fire kegler fra?), og det billede vi når frem ved at skanne igen og igen, ses til højre:


   



Bevis for Radons transformationsformel


Her er en skitse af Radons bevis for sin formel. Denne formel kaldes også for Radon-transformationen, og den hører blandt en type af matematiske formler som er særligt fascinerende, nemlig inversionsformlerne (som for eksempel Möbius- og Fourier-inversionen).

Inversionsprincippet kommer ind i billedet i forbindelse med en vis integralligning som vi får brug for at løse, det vil sige en ligning hvor den ubekendte er en funktion som indgår i en integraldannelse. Og den integralligning som skal løses lyder: givet en funktion g(y), find en funktion f(x) således at

g(y) = ∫ f(x)/√(x - y) dx (x fra y til ∞).

Men denne integralligning blev løst af Niels Henrik Abel (1802-29), og løsningsformlen er ganske tankevækkende:

f(x) = -1/π ∫ g'(y)/√(y - x) dy (y fra x til ∞)

- vi skal altså blot differentiere g(y) og bruge den samme formel (og dividere med -π). En så smuk inversionsformel må kunne ses i et højere perspektiv. Hvis faktoren 1/√(x - y) ikke havde været der, havde formlen været

g(y) = ∫ f(x) dx (x fra y til ∞).

Og denne integralligning kan løses ved blot at differentiere begge sider: f(y) = -g'(y). Disse to integraldannelser er specialtilfælde af følgende integraldannelse: Lad α være et positivt reelt tal, da kan vi for en funktion f(x) definere en ny funktion Iαf(y) ved integralet

Iαf(y) = 1/Γ(α) ∫ f(x) (x - y)α-1 dx (x fra y til ∞),

hvor Γ(x) er den såkaldte gammafunktion - her skal blot siges at for et naturligt tal n er Γ(n) = (n-1)! og Γ(½) = √π. Denne integraltransformation har følgende tre elegante egenskaber:

Iα(Iβf) = Iα+βf      (Iαf)' = Iα(f')      (I1f)' = -f

Den midterste siger at operationerne Iα og differentation kan ombyttes, og at den sidste (hvor α = 1) er den ovennævnte formel g'(y) = -f(y) (idet Γ(1) = 1).

Vi får Abels inversionsformel ved at sætte α = ½, så har vi nemlig at Abels integralligning kan skrives g = √π I1/2f, og så følger af de tre formler at I1/2(g') = (I1/2g)' = √π (I1/2(I1/2f))' = √π (I1f)' = -√π f, og dette er Abels løsningsformel.

Og nu til hvordan Radon omformer Abels inversionsformel til sin egen inversionsformel. Vi lader f(p) være materialets tæthedsfunktionen. Lad L(θ, ρ) være den rette linie som er normal til vinklen θ (fra x-aksen) og som har afstanden ρ fra origo, og lad f^(θ, ρ) være integralet af f(p) langs L(θ, ρ) (det tal vi har kaldt "modstanden" i linien):

f^(θ, ρ) = ∫ f(q) ds (s langs L(θ, ρ)),

idet q er det til parameteren s svarende punkt på L(θ, ρ). Og lad for et punkt p og et tal r ≥ 0, Mpf^(r) være middelværdien af f^(θ, ρ) for alle tangenterne til cirklen om p med radius r, det vil sige:

Mpf^(r) = 1/(2π) ∫ f^(θ, ρ) dθ (θ fra 0 til 2π),

hvor ρ er afstanden fra origo af den tangent hvis normal har vinkel θ, det vil sige: ρ = r + projektionen af vektoren fra origo til p på retningen θ. For ethvert θ er der to tangenter til cirklen som er normal til θ: den ene knyttet til θ, den anden knyttet til θ + π, da ρ skal være positiv.

Radons inversionsformel siger at vi fra denne funktion Mpf^(r) kommer tilbage til funktionen f(p) ved:

f(p) = -1/π ∫ (dMpf^(r)/dr) dr/r (r fra 0 til D).

I definitionen af Mpf^(r) indgår to integrationer - én langs en ret linie og én langs en cirkel - og Radon viser at disse i en vis forstand kan ombyttes, idet

Mpf^(r) = 2 ∫ Mpf(r*)/√(1 - (r/r*)2) dr* (r* fra r til ∞),

hvor Mpf(r*) er middelværdien af f(q) langs cirklen C(p, r*) med centrum i p med radius r*, det vil sige

Mpf(r*) = 1/(2πr*) ∫ f(q) ds (s langs C(p, r*)).

Går vi stykket s ud af en tangent til cirklen med centrum i p og radius r - fra skæringspunktet - kommer vi til et punkt med afstanden r* fra p - altså beliggende på cirklen C(p, r*). Da r* er hypotenuse i den retvinklede trekant med kateterne r og s, er s = √(r*2 - r2), og heraf følger at ds/dr* = r*/√(r*2 - r2) = 1/√(1 - (r/r*)2), og dermed at det langs linien L(θ, ρ) gælder at

∫ f(q) ds (s fra 0 til ∞) = ∫ f(q)/√(1 - (r/r*)2) dr* (r* fra r til ∞)

(q er punktet svarende til s). Men det er ∫ f(q) ds (s fra -∞ til ∞) vi skal bruge, og da der for et givet r* er to q - q+ og q- - som har samme afstand s fra skæringspunktet, skal vi i det andet integral erstatte f(q) med f(q+) + f(q-), så når vi ombytter integrationerne har vi gennemløbet C(p, r*) to gange, derfor 2-tallet i udtrykket for Mpf^(r).



Men denne sammenhæng mellem Mpf^(r) og Mpf(r*) er Abels integralligning: for hvis vi erstatter r med √y og r* med √x har vi at

Mpf^(√y) = ∫ Mpf(√x)/√(x - y) dx (x fra y til ∞)

(idet vi har benyttet at d(√x)/dx = 1/(2√x)). Så derfor har vi ifølge Abel at

Mpf(√x) = -1/π ∫ (d(Mpf^)(√y)/dy)/√(y - x) dy (y fra x til ∞).

Og benytter vi at d(Mpf^)(√y)/dy = (Mpf^)'(√y) /(2√y), erstatter y med r2 og sætter x = 0, får vi Radons formel (da Mpf(0) = f(p)):

f(p) = -1/π ∫ (dMpf^(r)/dr) dr/r (r fra 0 til ∞).



Programmerne

(skrevet i Pascal)


Skannings-programmet


ww er den halve bredde og wh er den halve højde af "pict"-billedet - og ra er dets radius. For et punkt på en normal i dette billede er (i, j) dets koordinatsæt set fra origo (centrum af billedet). na er det indtastede antal af vinkler og nd er det indtastede antal af forskydninger - altså bredde og højde af "trans_pict"-billedet. For p (vinkler) gående fra 0 til na - 1, og q (forskydninger) gående fra 0 til nd - 1 gøres det nedenstående. s-værdierne (som til at begynde med er sat til 0) er de søgte farveværdier. For et p er vv = (p + 0.5) * 2 * pi / na den tilknyttede vinkel i radianer, og h = ra * q / nd er forskydningen i forhold til radius.


      vv := (p + 0.5) * 2 * pi / na;

      u := cos(vv);

      v := sin(vv);

      h := ra * q / nd;

      s1 := 0;

      s2 := 0;

      s3 := 0;

      if abs(u) > abs(v) then

            begin

                  for j := -wh to wh - 1 do

                        begin

                              i := trunc((h - j * v) / u) + 1;

                              if (i >= -ww) and (i < ww) then

                                   begin

                                         s1 := s1 + mat1[ww + i, wh + j] / abs(u);

                                         s2 := s2 + mat2[ww + i, wh + j] / abs(u);

                                         s3 := s3 + mat3[ww + i, wh + j] / abs(u);

                                    end;

                        end;

            end

      else

            begin

                  for i := -ww to ww - 1 do

                        begin

                              j := trunc((h - i * u) / v) + 1;

                              if (j >= -wh) and (j < wh) then

                                    begin

                                          s1 := s1 + mat1[ww + i, wh + j] / abs(v);

                                          s2 := s2 + mat2[ww + i, wh + j] / abs(v);

                                          s3 := s3 + mat3[ww + i, wh + j] / abs(v);

                                    end;

                        end;

            end;

      mem1[p, q] := s1;

      mem2[p, q] := s2;

      mem3[p, q] := s3;


mat[ww + i, wh + j] er "pict"-billedets farveværdi i punktet (ww + i, wh + j) og mem[p, q] er "trans_pict"-billedets farveværdi i punktet (p, q). Men mem[p, q]-tallene er umiddelbart reelle tal, og de skal omdannes til bytes. Vi gør det ved at finde den størst forekommende s-værdi, smax, og lade byte-værdierne være round(255 * mem[p,q] / smax).



Gendannelses-programmet


For i gående fra 0 til bredde - 1 og for j gående fra 0 til højde - 1 gøres det nedenstående. ww er den halve bredde og wh er den halve højde og ra er radius af billedet. (a, b) er koordinatsættet til punktet (hvis farveværdier skal findes) set fra origo. na er bredden (antallet af vinkler) og nd er højden (antallet af forskydninger) af "trans_pict"-billedet. p går fra 0 til na - 1 og bestemmer vinklen v = (p + 0.5) * 2 * pi / na i radianer. r er det r som optræder i Radons integral. For en vinkel v er r afstanden fra (a, b) til normalen, og r1 (= r + a * cos(v) + b * sin(v)) er denne normals afstand fra origo. d-værdierne er de løbende dr-værdier i Radons integral. k (som er sat til 1.06) er det tal som de følgende d-værdier skal ganges med. mem[p, q] er farveværdien i punktet (p, q) i "trans_pict"-billedet. m-værdierne (som til at begynde med er sat til 0) finder M((a, b), r). m0-værdierne er de foregående m-værdier. De to første operationer må behandles for sig, ved z henh. 0 og 1. s-værdierne (som til at begynde med er sat til 0) finder Radons integral. g-værdierne (som er differenskvotienter) bruges i denne udregning, h-værdierne er de foregående g-værdier. n-værdierne - bestemt ved s-værdierne - er de søgte farveværdier i punktet (i, j).


      a := i - ww;

      b := j - wh;

      r0 := ra + sqrt(sqr(a) + sqr(b));

      d := 1 / 10;

      r := d / 2;

      s1 := 0;

      s2 := 0;

      s3 := 0;

      z := 0;

0:

      p := 0;

      m1 := 0;

      m2 := 0;

      m3 := 0;

1:

      v := (p + 0.5) * 2 * pi / na;

      r1 := r + a * cos(v) + b * sin(v);

      p1 := p;

      if r1 < 0 then

            begin

                  r1 := -r1;

                  if p1 < na div 2 then

                        p1 := p1 + na div 2

                  else

                        p1 := p1 - na div 2;

            end;

      q := trunc(nd * r1 / ra);

      ip := nd * r1 / ra - q;

      if r1 < ra then

            begin

                  m1 := m1 + mem1[p1, q] + ip * (mem1[p1, q + 1] - mem1[p1, q]);

                  m2 := m2 + mem2[p1, q] + ip * (mem2[p1, q + 1] - mem2[p1, q]);

                  m3 := m3 + mem3[p1, q] + ip * (mem3[p1, q + 1] - mem3[p1, q]);

            end;

      p := p + 1;

      if p < na then

            goto 1;

      if z = 0 then

            begin

                  m10 := m1;

                  m20 := m2;

                  m30 := m3;

                  r := r + d;

                  z := 1;

                  goto 0;

            end;

      g1 := (m10 - m1) / d;

      g2 := (m20 - m2) / d;

      g3 := (m30 - m3) / d;

      m10 := m1;

      m20 := m2;

      m30 := m3;

      if z = 1 then

            begin

                  h1 := g1;

                  h2 := g2;

                  h3 := g3;

                  r := r + d;

                  d := k * d;

                  z := 2;

                  goto 0;

            end;

      s1 := s1 + (h1 / (r + d / 4) + g1 / (r + 3 * d / 4)) * d / 2;

      s2 := s2 + (h2 / (r + d / 4) + g2 / (r + 3 * d / 4)) * d / 2;

      s3 := s3 + (h3 / (r + d / 4) + g3 / (r + 3 * d / 4)) * d / 2;

      h1 := g1;

      h2 := g2;

      h3 := g3;

      r := r + d;

      d := k * d;

      if r < r0 then

            goto 0;

      n1 := round(u * s1);

      n2 := round(u * s2);

      n3 := round(u * s3);


s-værdierne må korrigeres for at få n-værdierne. For det første skal der divideres med antallet af vinkler na og med π-et fra Radons formel. Og da vi i Scan-programmet har ganget tallene med 255/smax, skal vi nu gange med smax/255. Da smax højest kan være 255 * "diameter af "pict"", skal det tal u vi korrigerer med, være mindre end "diameter af "pict""/(na * π), men her er "pict" det oprindelige billede, og vi kender ikke smax, så derfor må vi gange med et tal li, som vi har kaldt "lysstyrke". Altså u = li * 2 * ra / (na * π).