Flerdimensionell analys i bildbehandling

Flerdimensionell analys i bildbehandling

Erik Melin 27 november 2006

1. Förord

Målet med den här lilla uppsatsen är att ge några exempel på hur idéer från kursen flerdimensionell analys kan användas i bildbehandling på datorer, ett område som i och med digitalkamerans utbredning fått stor betydelse för många människor.

Uppsatsen är inte menad som en kurs i hur man konkret programmerar en dator till att behandla digitala bilder, med de numeriska approximationer och den digitala geometri man då måste ha förståelse för. Istället gör vi en matematisk modell för bilder som direkt kan behandlas med analysens verktyg.

2. En matematisk modell för datorbilder

En bild på datorskärmen är uppbyggd av bildelement som kan anta olika färger och som kallas pixlar. Det vanliga (troligtvis av tekniska orsaker) är att pixlarna är kvadratiska eller rektangulära. Om skärmen har upplösningen 1024×768 punkter betyder det att en pixel kan adresseras av två koordinater (x, y), där x och y är heltal och 0 6 x 6 1023 och 0 6 y 6 767.

Vi behöver en matematisk modell för datorskärmen. Låt som vanligt Z beteckna heltalen och R beteckna de reella talen. Ett slutet intervall av reella tal betecknas [a, b] då a, b ∈ R. Om a och b är heltal låter vi [a, b] Z = [a, b]∩Z, dvs mängden {a, a + 1, . . . , b − 1, b} förutsatt att a < b.

Med dessa beteckningar kan datorskärmen, D, i första stycket skrivas som D = [0, 1023] _Z × [0, 767] _Z , alltså mängden av ordnade par (x, y) där x och y hämtas ur respektive intervall.

Vad är då en bild på D? Jo, en bild får vi om vi tillordnar varje pixel en färg. Om C betecknar mängden av möjliga färger, är en bild en funktion f : D → C.

För att komma vidare behöver vi titta närmare på C, mängden av färger.

För en svartvit bild vi C = {svart, vitt}, som i allmänhet kodas om till

1

2 2 EN MATEMATISK MODELL FÖR DATORBILDER

mängden {0, 1} där 0 motsvarar svart och 1 vitt. Jag vet inte om det idag finns någon som minns grafikstandardena CGA, EGA och VGA, med 4, 16, respektive hela 256 färger, men färgmodellen där är densamma som för svartvita bilder – fast med något större palett.

Mer intressant blir det för oss med gråskalebilder. Det är ganska vanligt att tillåta 256 ljusintensiteter, där 0 motsvarar svart och 255 är vitt. Alltså kan vi låta C = [0, 255] _Z . Nu har vi möjlighet att matematiskt definiera rimliga operationer med bilder. Om f, g : D → C är två gråskalebilder, är funktionen (x, y) 7→ f (x, y)/2 en mörkare version av f (alla ljusintensiteter har delats med två), medan negativet till f ges av funktionen

(x, y) 7→ 255 − f (x, y).

Om vi låter 0 6 t 6 1 representera tiden så är avbildningen (x, y, t) 7→ (1 − t)f (x, y) + tg(x, y)

en animation; en mjuk övergång från bilden f vid tiden t = 0 till bilden g då t = 1.

Men det finns ett problem med konstruktionerna ovan. Det är ju inte säkert att resultatet av operationen är ett heltal mellan 0 och 255, alltså att resultatet ligger i C. Om det resultatet inte är ett heltal kan man naturligtvis avrunda. Om resultatet är för stort eller för litet är det kanske rimligt att begränsa det genom en min/max operation, t.ex.

max(0, min(f (x, y), 255))).

För att slippa den här typen av komplikationer, ska vi låta C = R. Det visar sig också praktiskt att skala om värdemängden så att 0 motsvarar svart och 1 motsvarar vitt. Värden över 1 kan då tolkas som punkter som är ljusare än vad som kan återges, medan negativa värden är fysikaliskt något svårtolkade.

För att vi ska kunna behandla bilderna med verktyg från analysen, måste vi också göra något med definitionsområdet. Istället för att betrakta heltals- området D = [0, 1023] _Z × [0, 767] _Z ska vi utöka definitionsområdet till det reella området R = [0, 1023] × [0, 767] ⊆ R ² . Notera att D ⊆ R.

Nu uppstår följande problem: Om vi har en funktion definierad på D, hur ska vi utöka den till det större området R? Naturligtvis ska den utökade funktionen överensstämma med orginalfunktionen i heltalspunkterna. Men mellan heltalspunkterna finns uppenbarligen väldigt många val. Kanske finns det någon funktion med särskilt bra egenskaper?

Vi ska inte närmare gå in på vilka egenskaper som är önskvärda, men vi formulerar en sats nedan, som visar att vi i alla fall kan anta att den utökade funktionen är kontinuerligt deriverbar.

Sats 2.1. Givet en funktion q : Z ² → R finns en kontinuerligt deriverbar

funktion Q : R ² → R, så att Q(x, y) = q(x, y) för alla heltal x och y.

2.1 Färgbilder 3

Vi ger inget av bevis av satsen, men den som genomför övning 2.1 nedan får ett bevis genom att explicit konstruera en sådan funktion Q.

För att göra situationen ännu enklare för oss, ska vi anta att funktionen f är definierad på hela R ² och inte bara på rektangeln (skärmen) R. Men vi kräver att f (x, y) = 0 om vi kommer tillräckligt långt bort från skärmen.

Med andra ord, att f (x, y) = 0 utanför någon sluten begränsad mängd. Man säger då att f har kompakt stöd.

Definition 2.2. En gråskalebild är en kontinuerligt deriverbar funktion f : R ² → R med kompakt stöd.

Övning 2.1. a) Låt a och b vara två reella tal. Konstruera en deriverbar funktion f _ab : R → R, så att f ab (0) = a, f ab (1) = b, och f _ab ⁰ (0) = f _ab ⁰ (1) = 0.

Ledning: Det finns ett polynom i grad 3 med de angivna egenskaperna.

b) Låt F : R ³ → R definieras genom F (a, b, x) = f ab (x). Om f _ab är polynomet från ledningen till a), är F ett polynom i variablerna a, b, och x, och därför en deriverbar funktion. Sätt

g(x, y) = g _abcd (x, y) = F (F (a, b, x), F (c, d, x), y)

Visa att g(0, 0) = a, g(1, 0) = b, g(0, 1) = c och g(1, 1) = d. Visa vi- dare, genom att använda kedjeregeln, att ^∂g _∂x (0, y) = ^∂g _∂x (1, y) = 0 och att

∂g

∂y (x, 0) = ^∂g _∂y (x, 1) = 0.

c) Antag att q : Z ² → R är en godtycklig funktion. Definiera en funktion Q : R ² → R på följande sätt. Låt

Q(x, y) = g abcd (x − bxc , y − byc),

där funktionen g _abcd beror på x och y genom sambanden a = q(bxc , byc), b = q(bxc + 1, byc), c = q(bxc , byc + 1) och d = q(bxc + 1, byc + 1). (Beteckningen bxc betyder att x avrundas till närmaste mindre heltal, sålunda är b2.5c = 2, b3c = 3 och b−2.5c = −3.)

Visa att Q(x, y) är kontinuerligt deriverbar och att för alla heltal x och y gäller Q(x, y) = q(x, y). Funktionen Q är alltså ett exempel på en funktion av den typ som omtalas i Sats 2.1. Med andra ord är satsen nu bevisad.

2.1. Färgbilder

En vanlig modell för färgbilder är den så kallade RGB-modellen, där man anger en färg genom intensiteten för rött, grönt och blått. Alltså kan man betrakta en färgbild (i RGB-modellen) som en funktion f : R ² → R ³ . I vissa sammanhang kan man behandla en färgbild som tre separata gråskalebilder.

I andra fall är det nödvändigt att ta hänsyn till den sammanlagda färginfor- mationen.

Det finns också många andra modeller för färg. Ett viktigt exempel är

HSL, där bokstäverna står för hue, saturation och luminance, alltså nyans,

4 3 DERIVATOR OCH KANTDETEKTERING

mättnad och ljusintensitet. Här blir funktionsavbildningarna lite mer kom- plicerade, eftersom nyanskomponenten avbildas till en cirkel (färgcirkeln) och inte till en tallinje.

Färgmodeller är ett intressant område i sig med en hel del matematik involverad, men vi ska inte fördjupa oss mer i ämnet.

2.2. Något om digital geometri

Digitala bilder är ändliga, diskreta, objekt, inte oändliga och kontinuerliga som i vår modell. En dator kan dessutom bara hantera ändliga mängder, så för att få effektiva algoritmer för bilder arbetar datorn direkt med den ursprungliga digitala bilden.

Istället för att beräkna derivata genom en gränsvärdesprocess där, låt oss säga, h → 0, beräknar man en differens då exempelvis h = 1; differential- kalkylen ersätts av en differenskalkyl. Att man slipper gränsvärdesprocesser kan ju tyckas vara en fördel, eftersom gränsvärden kan vara besvärliga att beräkna. Men det finns nackdelar också. Undersöker man bevisen för t.ex.

kedjeregeln eller produktregeln ser man att det finns resttermer som försvin- ner när h → 0. Om vi tar h = 1 försvinner inte resttermerna och kalkylen med differenser blir därför mer komplicerad än kalkylen med differentialer.

I en digital bild är det inte heller klart vad man ska mena med en kurva eller en rät linje. Hur kan man se på en mängd pixlar att de utgör en rät linje? Finns begreppet kontinuitet för funktioner definierade på en digital bild? Frågar som dessa försöker man besvara inom digital geometri, som är ett område av matematiken.

3. Derivator och kantdetektering

Ett standardproblem inom bildanalys är att hitta kanter i en bild för att kunna dela upp bilden i olika komponenter. Vad utmärker då en kant i en gråskalebild? Jo, att ljusstyrkan ändras snabbt i någon riktning. Vi kan an- vända gradienten för att undersöka kanter.

Vi börjar med att undersöka problemet i en dimension,. Låt f : R → R vara en deriverbar funktion i en variabel, där värdemängden tolkas som lju- sintensitet. När ljusintensiteten ändras snabbt är derivatan stor (till belop- pet), och vi säger att bilden har en kant där |f ⁰ (x)| > M för någon positiv konstant M (se Figur 1). En svårighet är att välja en bra konstant. Ju mind- re M är, desto fler (men svagare) kanter hittas. I praktiken måste M , som ofta kallas tröskelvärde, anpassas för varje tillämpning.

Låt nu f : R ² → R vara en gråskalebild. Gradienten, ∇f = ( ^∂f _∂x , ^∂f _∂y ), i

punkten (x ₀ , y 0 ), är en vektor vars längd ger den maximala förändringen av

funktionen f och som pekar i den riktning som den maximala förändringen

sker.

5

Figur 1: En endimensionell funktion och dess derivata. Kant har vi då funk- tionen ändrar sig snabbt, det vill säga då derivatan är stor till beloppet.

På samma sätt som i det endimensionella fallet låter vi en punkt, (x ₀ , y 0 ) i bilden räknas som en kant om |∇f (x ₀ , y ₀ )| > M för någon konstant M . Med andra ord, för en gråskalebild f , kan kantmängden med tröskelvärde M definieras genom

Kant _M = {(x, y) ∈ R ² ; |∇f (x, y)| > M }.

Man kan visa att Kant _M är en sluten och begränsad mängd. Däremot finns det inga garantier att Kant _M är en kurva. Tvärtom består Kant _M i allmänhet av områden med positiv area (jämför Figur 1 där kantområdena är intervall), och kanske enstaka punkter. Att tolka den information som kommer från gradienten och sätta ihop kanten till sammanhängande, tunna, kurvor är ett problem som man försöker lösa inom bildanalysen. Problemet är svårt.

Vi nöjer oss med att ge en illustration. Den vänstra bilden i Figur 2 föreställer Orienteringsklubben Linnés klubbgård i västra Uppsala. Den hög- ra bilden är kantmängden (det svarta), för något visst tröskelvärde. Vissa kanter syns tydligt, medan träden i bakgrunden bara blir meningslöst brus.

4. Integraler och oskärpa

Låt f vara en gråskalebild och antag att bildområdet är rektangeln R som ges av 0 6 x 6 a och 0 6 y 6 b för några positiva tal a och b. Speciellt kräver vi att f = 0 utanför R. Det medför att RR

R f (x, y)dA = RR

R

f (x, y)dA.

Ibland har man behov av att jämna ut en bild, att göra den suddigare.

Vi har redan sett exempel på hur brus ställer till problem när man ska hitta kanter i bilden.

Den totala ljusintensiteten som sänds ut av bilden, kan man få genom att

addera ljusintensiteterna för alla pixlar. I vår kontinuerliga modell motsvarar

6 4 INTEGRALER OCH OSKÄRPA

Figur 2: En gråskalebild och dess kantmängd. Bilden ska ses principiellt; de (numeriska) metoder som används för att faktiskt beräkna den tas inte upp i den här artikeln.

det integration. Den totala ljusintensiteten, I, är alltså I =

Z Z

R

f (x, y)dA.

En mycket utsuddad version av bilden är den man får om man tilldelar alla bildpunkter den genomsnittliga ljusintensiteten, alltså den konstanta bilden g(x, y) = c, för (x, y) ∈ R, där c ges av

c = 1

Area(R) Z Z

R

f (x, y)dA = 1 ab

Z Z

R

f (x, y)dA.

(Ett litet problem med den här definitionen är att g gör ett språng ner till 0 utanför R. Det går naturligtvis att lösa genom att runda till funktionen på kanten.)

Den konstanta bilden är naturligtvis alldeles för utsuddad för att vara av någon större praktisk betydelse. En kompromiss är att istället beräkna den genomsnittliga ljusintensiteten i en omgivning av varje punkt, t.ex. i en skiva med radie r.

Låt r > 0 och låt D _r = {(x, y); x ² +y ² 6 r ² } vara den slutna cirkelskivan med radie r. Om f är en kontinuerlig funktion i planet, så är

A r = 1

Area(D _r ) Z Z

D

f (x, y)dA = 1 πr ²

Z Z

D

f (x, y)dA (1) det genomsnittliga funktionsvärdet i D _r . Notera att

r→0 lim A r = f (0, 0),

vilket kan bevisas t.ex. med hjälp av medelvärdessatsen för integraler.

4.1 Mer avancerad oskärpa 7

Vi ska skriva om integralen (1) med hjälp av funktionen H _r (x, y) som definieras genom

H r (x, y) =

 1/πr ² om x ² + y ² 6 r ² , d.v.s. om (x, y) ∈ D _r

0 annars.

Det gäller nu att A _r = 1

πr ² Z Z

D

f (x, y)dA = Z Z

R

H _r (x, y)f (x, y)dA.

Med hjälp av integralen ovan kan vi alltså beräkna medelvärdet av f i en skiva med radie r kring punkten (0, 0) . För att beräkna medelvärdet kring en annan godtycklig punkt (x ₀ , y 0 ) behöver vi bara förflytta skivan på vanligt sätt. Definiera

g(x ₀ , y ₀ ) = Z Z

R

H _r (x − x ₀ , y − y ₀ )f (x, y)dA. (2) Funktionen g i en punkt (x ₀ , y ₀ ) är alltså medelvärdet av f i en skiva med radie r kring denna punkt. Vidare gäller att

r→0 lim g(x 0 , y 0 ) = f (x 0 , y 0 ).

Ju större r man väljer, desto suddigare blir bilden, och låter man r → 0 får man tillbaka orginalbilden. Ett problem som vi har sopat under mattan är att bilden är en rektangel. Går skivan, D _r , utanför bilden (där f = 0) blir resultatet för mörkt. Detta måste i en verklig implementation hanteras på något sätt (Övning: Hur kan man lösa det?), men vi ignorerar det problemet.

Som förberedelse inför nästa avsnitt noterar vi att RR

R

H _r (x, y) = 1.

Om integralen antog något annat värde, skulle vi ändra den genomsnittliga ljusintensiteten i bilden. Det är inte önskvärt.

4.1. Mer avancerad oskärpa

Istället för att bara använda medelvärdet, kan det vara en bra idé att använda ett vägt medelvärde, där närliggande punkter får större betydelse än punkter långt borta. Det finns flera möjligheter, men ett alternativ som ofta används är ett så kallat Gaussiskt filter.

I en variabel har vi funktionen G _σ (x) = e ⁻

för en parameter σ > 0.

Figur 3 föreställer funktionskurvan för σ = 1. Ju större σ desto långsammare går funktionen mot 0 när x → ±∞.

I två dimensioner låter vi G _σ (x, y) = _2πσ ¹

e ⁻

; funktionsytan är allt-

så i princip den endimensionella grafen roterad ett varv kring den beroende

8 4 INTEGRALER OCH OSKÄRPA

x

3 2 1 0

−1

−2

−3

1,0

0,75

0,5

0,25

0,0

Figur 3: Kurvan y = e ^−x

^/2 . Om kurvan roteras runt y-axeln får man en klockliknande yta.

axeln. Av den anledningen kallas grafen ibland för en Gaussisk klocka. Skal- faktorn framför motsvarar faktorn _πr ¹

i definitionen av H _r (x, y) i föregående avsnitt. Vi har nämligen, genom att byta till polära koordinater,

Z Z

R

G σ (x, y)dA = 1 2πσ ²

Z Z

R

e ⁻

dA

= 1

2πσ ² Z _2π

0

dθ Z ∞

0

e ^−r

^/2σ

r dr = 1.

Vi applicerar den Gaussiska funktionen på en bild på samma sätt som vi i föregående avsnitt applicerade H _r (x, y). Om f är en gråskalebild ges en bild med Gaussisk oskärpa (Gaussian blur) av funktionen

g(x 0 , y 0 ) = Z Z

R

G σ (x − x 0 , y − y 0 )f (x, y)dA.

Parametern σ brukar kallas för standardavvikelse, av orsaker som kommer att få sin förklaring i någon kurs om sannolikhetsteori. I Figur 4 illustreras Gaussisk oskärpa och kantmängden för den oskarpa bilden.

4.2. Allmänna filter

Operationen vi har utfört i det här och föregående avsnitt kan generaliseras.

Om H : R ² → R är en integrerbar funktion och f är en gråskalebild, definierar H ett filter på f genom operationen

(H ∗ f )(x ₀ , y ₀ ) = Z Z

R

H(x ₀ − x, y ₀ − y)f (x, y) dA (3)

Uttrycket H ∗ f är alltså en ny funktion i två variabler. Av orsaker vi inte

går in på brukar man, som ovan, integrera med H(x ₀ − x, y ₀ − y) istället för

4.2 Allmänna filter 9

Figur 4: Gaussisk oskärpa (σ = 3 pixlar), och kantmängden för den oskarpa bilden. Kantmängden för den oskarpa bilden innehåller betydligt mindre brus än kantmängden för orginalbilden.

H(x − x ₀ , y − y ₀ ) som vi gjort tidigare. Är H jämn i x och y (vilket både H _r och G _σ är) gör det ingen skillnad.

Den filtrerade bilden (f ∗ H) brukar kallas en faltning av funktionerna f

och H (På engelska heter faltning convolution). Faltning är en mycket bety-

delsefull matematisk operation, och förekommer ofta inom sannolikhetsteori

och transformteori.