Einführung Bildverarbeitung

(1)

Einführung Bildverarbeitung

D. Schlesinger () Einführung Bildverarbeitung 1 / 8

(2)

Was sind Bilder?

(3)

Was sind Bilder?

Bilder sindVektoren.

„Ähnliche“ Vektoren entsprechen nicht unbedingt „ähnlichen“ Bildern.

„Ähnliche“ Bilder entsprechen nicht unbedingt „ähnlichen“ Vektoren.

Beispiele: Verschiebung, andere geometrische Transformationen, Farbtransformationen ...

⇒Bei Bildern ist die Räumliche Anordnung der Elemente extrem wichtig – geht aber meist verloren.

(4)

Was sind Bilder?

Bilder sindMatrizen.

Was bedeuten z.B. die Eigenwerte, die Determinante?

Kann man Bilder wie Matrizen multiplizieren?

⇒typische Operationen und Eigenschaften sind nutzlos.

Bilder sindGraphen.

Pixel sind Knoten, 4- (oder 8- oder komplizierter) Nachbarschaftstruktur.

Ohne Weiteres zu eingeschränkt (nur Gitter)→attributierte, gelabelte Graphen etc.

Wird eher für höhere Stufen der Verarbeitung (bei Computer Vision) verwendet Bilder sindFunktionen(Abbildungen).

Sind sie z.B. konvex, stetig, differenzierbar? Wie sind sie zu spezifizieren?

Bilder sindBilder

Bildverarbeitung – wie geht man damit um.

(5)

Themen

1. Menschliches Sehen (0.5 DS – eher zum Spaß) 2. Lineare Filterung, Morphologische Operationen (2 DS)

Effiziente Algorithmen (Übungsaufgaben) 3. Diffusion Filters: „Bilder sind Funktionen“

4. Fourier-Transformation (Klassik) 5. Besondere Punkte, Bildmerkmale

mit Hinblick auf Erkennung (Computer Vision) 6. Momente, Hauptkomponentenanalyse, Clustering 7. Diskrete und Kontinuierliche Energieminimierung (2 DS)

Am Beispiel der Entrauschung 8. Kameramodelle, 3D-Geometrie 9. Anwendungsbeispiel: Stereo

2, 5, 7 und 8 zusammen in einer Anwendung

(6)

Voraussetzungen

„Solides mathematisches Grundlagenwissen“. Was heißt das?

Beispiele – man muss zumindest verstehen, worum es geht:

ln

Y

i

f(xi) =

X

i

lnf(xi)

min

x f(x) =−max

x

−f(x)

arg min

x

f(x) = arg min

x

lnf(x) min

x

X

y

f(x,y)≥

X

y

min

x f(x,y)

n

X

i=1

ailnxi→max

x

s.t.xi≥0,

X

i

xi= 1

xi∼ai

Konkret: Geometrie (sin, cos), Lineare Algebra (Vektoren, Matrizen, SVD), Funktionen (Ableitungen, Gradienten, Integrale, Reihen) ...

(7)

Übungen

Freier Besuch – Konsultationen

Übungsaufgaben: alles, was mit BV/ME zu tun hat, Aufgabenkatalog (später, im Laufe), eigene Aufgaben sind willkommen

Abschluss: Punktesystem – 1-3 Punkte/Aufgabe, 4 Punkte insgesamt, Gruppen bis max. 3 Leute möglich

Umgebungen: Linux(Windows), C/C++, Qt, OpenCV (siehe Beispiele aufwww. ...) Abgabe: Kommentierte Quelltexte, ggf. Dokumentation und/oder Ausarbeitung, Einleitung zum Kompilieren/Aufrufen/Ein- und Ausgabe

Per E-Mail anDmytro.Shlezinger@tu-dresden.de

(8)

Organisatorisches

Vorlesungsscripte, Folien, Aufgaben, Info:

http://www1.inf.tu-dresden.de/∼ds24/lehre/bvme_ss_2012/bvme_ss_2012.html Abschluss:

mündliche Prüfung (Fachgebiete 8 SWS, Vertiefungsgebiete 12 SWS) Mit Übungen – 4 SWS, ohne – 2 SWS

Literatur:

– David A. Forsyth, Jean Ponce: ”Computer vision – A modern approach”

– Klaus D. Tönnies: ”Grundlagen der Bildverarbeitung”

– Internet, Google, Wikipedia ...

– Papers ...

Kommentare, Wünsche, Kritik sind willkommen (auch anonym via Mail-Formular).

(9)

Bildverarbeitung: Das menschliche Sehen

D. Schlesinger () BV: Das menschliche Sehen 1 / 8

(10)

Hierarchie der Verarbeitung

Die Verarbeitung ist hierarchisch aufgebaut Die Sehnerven tauschen in der Sehkreuzung (chiasma opticum) die Hälfte der Nervenfasern, so dass die linken Gesichtshälften zur rechten Hemisphäre projizieren und umgekehrt Erhaltung der Topologie in früheren Schichten Rückkopplung!!!

(11)

Das Auge

Die Netzhaut enthält unterschiedliche Sensortypen:

Zäpfchen (Farben, 6 Mio.) und Stäbchen (Grauwerte, 120 Mio.)

Licht muss zuerst die Neuronenschicht durchqueren bevor es auf die Fotosensoren trifft (Glättung), nur in der fovea centralis werden diese unmittelbar von den Lichtstrahlen getroffen

Retina→Ganglienzellen (1 Mio.)→Sehnerv→...

Wahrnehmung der Farben Räumliche Auflösung Sakkaden

(12)

Wahrnehmung der Farben

Was ist Licht? – Spektrum, d.h. eine Funktion der Wellenlänge Menge aller Spektren ist unendlichdimensional!!!

Spektrale Auflösung des Auges ist relativ schlecht, weil∞ →3 (Projektion) – nureineFarbe an einem bestimmten Ort.

(13)

Räumliche Auflösung

.

(14)

Räumliche Auflösung

In Fovea ist die Dichte viel höher

Die Information ist im Auge durch Ganglienzellen bereits vorverarbeitet (Vergleiche: 3072×2304=7MPixel→2.4 MB RGB JPEG-Bild verlustfrei) Kein Standbild, sondern „Video“

Abtasten der Szene mittels Sakkaden

(15)

Sakkaden

Augenbewegungen wechseln sich mit Fixationsperioden ab. Bei konjugierten Augenbewe- gungen bewegen sich die Augen nie gleichmäßig, sondern springen in Sakkaden (von ca.

15-100 ms Dauer) von einem Fixationspunkt zum nächsten, große Sakkaden werden oft von Kopfbewegungen begleitet. Sakkaden werden durch Konturen oder Stellen geführt, die für die Bedeutung des Gesehenen besonders wichtig sind (z.B. Augen, Mund).

(16)

Optische Täuschungen

Gibt es optische Täuschungen in Computer Vision?

Ja – schlecht gestellte Probleme.

(17)

Bildverarbeitung: Filterung

D. Schlesinger () Bildverarbeitung: Filterung 1 / 17

(18)

Allgemeines

Klassische Anwendung: Entrauschung

(Fast) jeder Filter basiert auf einem Modell (Annahme): Signal + Rauschen

Pipeline: Modell→Aufgabe→Lösung→Algorithmus (Programm) – Mittelwertfilter: Modell→...→Lösung

– Medianfilter: Aufgabe→Lösung – Paar andere Filter, Beispiele – Lineare Filterung: Algorithmen

(19)

Mittelwertfilter

D⊂Z²– Definitionsbereich (Gitter),r – Pixelr= (i,j),r∈D Bild ist eine Funktionx:D→C (Farbraum),xr – Farbe des Pixelsr Seiydas „Ideale“ Signal (nicht verrauscht) und

xdas verrauschte Signal (Beobachtung) Aufgabe: Man siehtx, man berechney.

Rauschmodell – Gaussche Wahrscheinlichkeitsverteilung für Abweichungen der Farben p(xr|y_r) =N(xr;yr, σ)∼exp

h

−kxr−yrk² 2σ²

i

Weitere Annahme – Pixel sind von einander unabhängig verrauscht:

p(x|y)∼

Y

r

exp

h

−kxr−yrk² 2σ²

i

Aufgabe nach dem Maximum Likelihood Prinzip:p(y|x) =p(y)p(x|y)/p(x)→maxy

Ohne zusätzlicher Annahmen überp(y) ist die Lösung trivial:yr=xr für aller Eine Annahme über das Signalyist notwendig!!!

(20)

Mittelwertfilter

Annahme:

In einer kleinen UmgebungW(r)⊂Deines Pixelsr isty(fast) konstant –yr

ML:

Y

r⁰∈W(r)

exp

h

−kx_r0−yrk² 2σ²

i

→max

yr

ln ... : F(yr) =

X

r⁰

kx_r0−yrk²→min

yr

Ableiten: ∂F

∂yr

=

X

r⁰

(x_r0−yr) =

X

r⁰

x_r0− |W| ·yr= 0

⇒ yr= 1

|W|

X

r⁰

xr⁰ (Mittelwert)

(21)

Median Filter

Andere Zielfunktion – entspricht einem anderen Rauschmodell (Einfachheit halberC=R– Grauwertbild):

F(yr) =

X

r⁰

|x_r0−yr| →min

yr

Problem: nicht differenzierbar. Gute Neuigkeit – konvex.

(Ableitung an der Tafel).

Lösung – Medianfilter.

(22)

Vergleich: Gaussches Rauschen

(a) Original (b) Verrauscht

(c) Gaussche Glättung (d) Medianfilter

(23)

Vergleich: „Salz und Pfeffer“ Rauschen

In manchen Pixeln ist der Farbwert zufällig gesetzt (gleichwahrscheinlich)

(e) Original (f) Verrauscht

(g) Gaussche Glättung (h) Medianfilter

(24)

Eine „andere“ Aufgabe

Glättung des Hintergrunds verbessert die räumliche Wahrnehmung der Szene

(Xue Bai, Guillermo Sapiro)

Umgekehrte Aufgabe: Rekonstruiere 3D aus Unschärfe.

Konsequenz: Die Filter liefern keine “Antworten“,

sie geben ein lokales Maß für die höheren Stufen der Verarbeitung.

(25)

Noch paar andere Aufgaben

Faltung:

yr =

X

r⁰

x_r0·g_r0−r

mit der Maskeg:Z²→R Beispiel – Mittelwert:

gr =

n

_1/|W_| _wenn _r_∈_W(0)

0 sonst.

Welche weitere Masken sind denkbar (nützlich, wofür etc.)?

– Etwas „zartere“ Glättungen, z.B. mit dem Gausschen Kern:gr ∼exp

−krk²/2σ²

– Kontrasterhöhung:gr∼α·1I(0)−β·exp

−krk²/2σ²

(Unsharp Mask)

– Kantendetektoren, etc.

(26)

Eine ganz andere Aufgabe

„Salz und Pfeffer“ Rauschen mit 90%

Lokale Filterung ist offensichtlich kaum hilfreich

Explizite Modellierung des Signals ist notwendig (in diesem Fall – MRF)

Das Modell entspricht dem Bild „exakt“. Dies macht genaue Rekonstruktion möglich.

(27)

Fazit

Explizite Annahmen über das Signal (lokale Eigenschaften) – lokale Filterung Signal wird durch seine Autokorrelationen beschrieben – Wiener Filter

Signal ist die Überlagerung der Signale bestimmter Frequenzen – Fourier Analyse Signal wird durch eine Differentialgleichung beschrieben – Variationelle Ansätze Explizite Modellierung statistischer Abhängigkeiten zwischen den Pixelwerten – MRF u.s.w.

Die Eigenschaften des Signals

(das Modell, die Annahmen, die Art der zu lösenden Aufgabe etc.) sind ausschlaggebend für die Wahl des Filters.

(28)

Lineare Filterung – Algorithmen

Eindimensional:r∈N

Beispiel: Mittelwertfilteryr=_|W|¹

r+w

P

r⁰=r−w

x_r0 (|W|= 2w+ 1)

Naiver Algorithmus (direkt nach der Formel):

for r= 0 bis n sum= 0

for r⁰=r−w bis r+w sum=sum+x_r0

yr =sum/|W| Zeitkomplexität:n· |W|

(29)

Lineare Filterung – Algorithmen

Idee: (Bildchen)

r+w

X

r⁰=r−w

xr⁰=

r+w

X

r⁰=0

xr⁰−

r−w−1

X

r⁰=0

xr⁰ = ˜xr+w−x˜r−w−1

Besserer Algorithmus:

Berechne ˜xr für alle r:

for r= 0 bis n

˜

xr = ˜xr−1+xr

Berechne yr: for r= 0 bis n

yr = (˜xr+w−x˜r−w−1)/|W| Zeitkomplexität:n

(30)

Lineare Filterung – Algorithmen

Verallgemeinerung auf 2D:

Verwendung des „Integralbildes“ ˜x Berechne ˜xr:

for (i,j) = (0,0) bis (m,n) (zeilenweise)

˜

xr = ˜xi,j−1+ ˜xi−1,j−x˜i−1,j−1+xr

Berechne yr:

for (i,j) = (0,0) bis (m,n) yr = (˜xi+w,j+w−˜xi−w,j+w−

˜

xi+w,j−w+ ˜xi−w,j−w)/|W|

+ -

- +

(31)

Lineare Filterung – Algorithmen

Ein etwas komplizierteres Beispiel (wieder 1D, Bildchen):

yr =

r

X

r⁰=r−w

(w−r+r⁰)·x_r0

Man überlege, ob sichyr+1aus demyr effizient berechnen lässt:

yr+1 =

r+1

X

r⁰=r+1−w

(w−r−1 +r⁰)·x_r0

=

r+1

X

r⁰=r+1−w

(w−r+r⁰)·x_r0−

r+1

X

r⁰=r+1−w

x_r0

=

r

X

r⁰=r−w

(w−r+r⁰)·x_r0+w·xr+1−

r

X

r⁰=r−w

x_r0

= yr+w·xr+1−¯xr

¯

xr kann mit linearer Zeitkomplexität berechnet werden (Mittelwertfilter)

⇒die gesamte Zeitkomplexität ist linear.

(32)

Lineare Filterung – Algorithmen

Faltungen allgemein:

y=x∗g yi=

∞

X

j=−∞

xi−j·gj

Eigenschaften:

– sind kommutativ, d.h.x∗g=g∗x;

– sind assoziativ, d.h. (x∗g¹)∗g²=x∗(g¹∗g²);

– sind distributiv mit „+“, d.h.x∗(g¹+g²) =x∗g¹+x∗g². Identische Faltung (ändert das Signal nicht):g_j^I= 1I(j= 0) Inverse Faltungen:g∗g⁻¹=g^I

Beispiel (j= 0 ist fett gekennzeichnet):

g^diff = [. . . ,0,0,0,1,−1,0, . . .] Differential Operator g^int = [. . . ,0,0,0,1, 1,1, . . .] Integral Operator

(33)

Lineare Filterung – Algorithmen

Der Trick zur effizienten Berechnung basiert auf der folgenden Umwandlung:

x∗g=x∗gÎ∗g=x∗gînt∗g^diff∗g= (x∗gînt)∗(g∗g^diff) oder sogar mehr

x∗g= (x∗gînt∗. . .∗gînt)∗(g∗g^diff∗. . .∗g^diff) Die Faltungx∗gîntbraucht lineare Zeitkomplexität,

mit ˜g=g∗g^diff ∗. . .∗g^diff (Vorberechnung) wird erreicht, dass ˜gschwach besetzt ist.

(34)

Bildverarbeitung: Morphologische Filterung

D. Schlesinger () BV: Morphologische Filterung 1 / 8

(35)

Dilation und Erosion

Zunächst für binäre Bilderx:D→ {0,1}

UND (Erosion): yr =

V

r⁰∈W(r)

x_r0 y=x W

ODER (Dilation): yr =

W

r⁰∈W(r)

x_r0 y=x⊕W Wwird „Strukturelement“ genannt (Kreis, Rechteck usw.)

(a) Original (b) Verrauscht

(c) Erosion (d) Dilation

(36)

Opening und Closing

Opening:x◦W= (x W)⊕W

Closing:x•W= (x⊕W) W

(nicht kommutativ)

(37)

Erweiterung auf Grau(Farb)wertbilder

Erosion:yr= min

r⁰∈W(r)x_r0

Dilation:yr= max

r⁰∈W(r)

x_r0

RGB – zum Beispiel Getrennt für R, G und B (andere Varianten auch möglich)

(m) Original (n) Verrauscht

(o) min (p) max (q) max (r) min

(38)

1D-min-Algorithmus

Aufgabe in 1D:yr= ^r+Wmin

r⁰=r−Wx_r0

Naiver Algorithmus: laut der Formel – probiere alleWElemente und wähle das minimale.

Zeitkomplexität –O(nW)

Die Idee zur Beschleunigung:

– Aktualisiere die „geordnete“ Menge der Werte [xr⁰...].

– In jedem Schritt soll ein Element aus der Menge entfernt werden und ein Element hinzugefügt werden.

– Nutze die Datenstrukturen, die es erlauben diese zwei Operationen inO(lnW) durchzuführen.

⇒Zeitkomplexität –O(nlnW)

(Erklärung an der Tafel)

(39)

Verallgemeinerung

StrukturelementW⊂R² wird durch Strukturfunktionw:R²→Rersetzt:

Erosion:yr= min

r⁰∈Ω x_r0+w(r−r⁰)

Das Vorige ist ein Spezialfall:

w(x) =

n

₀ _wenn _x_∈_W

∞ sonst.

Vergleiche mit der linearen Filterungyr =

P

r⁰

x_r0·g(r−r⁰) Ist dasselbe bis auf die verwendeten Operationen.

⇒Morphologische Filterung ist „lineare Filterung im (min,+) Semiring“.

(z.B. die Erosion entspricht dem Mittelwertfilter)

(40)

Separierbare (lineare) Filter

(Achtung!!! – etwas andere Bezeichnungen)

Sei eine Maskeg(i,j) als Produktg1(i)·g2(j) darstellbar.

Beispiel: Gaussche Glättung

g(i,j)∼exp[−(i²+j²)/2σ²] = exp[−i²/2σ²]·exp[−j²/2σ²] =g1(i)·g2(j) Die Faltungx∗g:

X

ij

x(i,j)·g1(i)·g2(j) =

X

i

g1(i)·

hX

j

x(i,j)·g2(j)

i

=

X

i

g1(i)·˜x(i)

x∗g=x∗(g1∗g2) = (x∗g1)∗g2

SeiW ein Fenster der Größem×m. Zeitkomplexität:O(m²)→O(m).

Mittelwertfilter ist separierbar, viele andere auch, z.B. Gaussche Glättung.

min-Filter ist auch „separierbar“, d.h. die Zeitkomplexität in 2D istO(nlnW), Median – leider nicht.

(41)

Box-Filter

Die sequenzielle Anwendung der Mittelfilter approximiert die Gaussche Glättung.

(basiert auf Central Limit Theorem)

x∗gGauss≈(((x∗gM)∗gM). . .)∗gM

Zeitkomplexität:

O(nW)→O(nm) mitm die Anzahl der Anwendung von Mittelwertfilter (5-6 Mal).

(42)

Bildverarbeitung: Fourier-Transformation

D. Schlesinger () BV: Fourier-Transformation 1 / 16

(43)

Allgemeines

Bilder sind keine Vektoren.

Bilder sind Funktionenx:D→C(Menge der Pixel in die Menge der Farbwerte).

Allerdings kann man eine Funktion auch als ein Vektor verstehen (darstellen).

Zum Beispiel:

eine Funktionf(x),x∈R,f(x)∈R, d.h.f :R→Rist ein Vektor imR^∞

Somit kann man mit den Funktionen alles machen, was man mit Vektoren machen kann:

addieren, multiplizieren, Skalarprodukt etc.

Beispiel:

Skalarprodukt zweier Funktionen ist (bis auf paar Details) der Korellationskoeffizient.

Man kann über Funktionsräume sprechen, Basisfunktionen, zu einander orthogonale Funktionen, linear unabhängige Funktionen etc.

Konsequenz: Bilder sind mehr als Vektoren, Vektoren sind sie aber auch.

(44)

Allgemeines

Die Aufgabe – zerlege einen Vektorx∈Rⁿ auf seine „Komponenten“ in einem Basis x=

X

i

vi·λi,

mit den Basisvektorenvi∈Rⁿ und den Koeffizientenλi∈R. Äquivalent – löse ein lineares Gleichungssystemx=V·λmitλ∈Rⁿ Wichtig (Eigenschaften des Basis):

– Die Basisvektoren sollen den Raum aufspannen, d.h. eine solche Zerlegung existiert für allex.

– Die Vektoren sind linear unabhängig, d.h. ein Vektorvi lässt sich nicht als eine lineare Kombination anderer Vektorenvj darstellen – die Zerlegung einesxist dann eindeutig.

Spezialfall – orthonormierter Basis:

– allevi sind zu einander orthogonal, d.h.hvi,vji= 0 füri6=j.

– allevi haben dieselbe Länge (= 1), d.h.hv_i,vii= 1.

Dann gilt:λi=hx,vii.

(45)

Allgemeines

Übergang zu Funktionen:

Der Funktionsraum ist unendlichdimensional→

→unendlich viele Basisfunktionenvi(x), d.h.v(x,y) (yersetzti) sowie

→eine kontinuierliche Funktionλ(y).

Die Aufgabe ist, eine gegeben Funktionf(x) auf Basisfunktionen zu zerlegen:

f(x) =

Z

y

v(x,y)λ(y)dy

Orthonormierter Basis heißt:

– Orthogonal:

Z

x

v(x,y⁰)v(x,y⁰⁰)dx= 0 für alley⁰6=y⁰⁰. – Normiert:

Z

x

v(x,y)v(x,y)dx=const für alley.

Dann gilt:

λ(y) = „hi“ =

Z

x

f(x)v(x,y)dx

(46)

Fourier-Transformation

Funktionsraum:

alle periodische Funktionen mit der Periode 2π, d.hf(x) =f(x+k·2π),k∈Z.

Basisfunktionen: sin(kx) und cos(kx),k= 0, . . . ,∞ Eigenschaften:

– orthonormiert (trivial),

– spannen den Funktionsraum auf (Jean Baptiste Joseph Fourier, 1822) Zerlegung:

f(x) = a0

2 +

∞

X

k=1

akcos(kx) +bksin(kx)

mit

a0= 1/π

Z

2π

0

f(x)dx

ak= 1/π

Z

2π

0

f(x) cos(kx)dx

bk= 1/π

Z

2π

0

f(x) sin(kx)dx

(47)

Fourier-Transformation

Ersetzt man

∞

P

1

durch

k_max

P

1

, so wird die Ausgangsfunktionf(x) approximiert.

(48)

Fourier-Transformation

Komplexe Schreibweise:

Grundlage – Euler Identität:

e^ikx = cos(kx) +i·sin(kx) e^−ikx = cos(kx)−i·sin(kx) Zerlegung:

f(x) =

∞

X

k=−∞

cke^ikx

Koeffizienten:

ck= 1 2π

Z

2π

0

f(x)e^−ikxdx=

(

a0/2 k= 0 1/2·(ak−ibk) k>0 1/2·(a−k+ib−k) k<0

(49)

Fourier-Transformation

Erweiterung auf beliebige periodische Signale: cos(kx)→cos(^2πkx_T ) Erweiterung auf nichtperiodische Signale (GrenzübergangT→ ∞)

→die Koeffizienten werden kontinuierlich

→die Reihec0,c1, . . .wird zur komplexen Funktion reellwertiges Argumentes F(u) =R(u) +I(u)

Amplitudenspektrum:

|F(u)|=

p

R²(u) +I²(u) Phasenspektrum:

φ(u) = tan⁻¹ I(u) R(u)

(50)

Fourier-Transformation

2D Diskrete Fourier-Transformation:

F(u,v) = 1 MN ·

M−1

X

x=0 N−1

X

y=0

f(x,y)e−i2π(xu/M+yv/N)

mitM undN – Breite und Höhe,xundy– Bildkoordinaten,uundv– Frequenzen.

Inverse dazu:

f(x,y) =

M−1

X

u=0 N−1

X

v=0

F(u,v)ei2π(xu/M+yv/N)

(51)

Fourier-Transformation

Beispiele – charakteristische Amplitudenspektren:

(52)

Fourier-Transformation

Beispiele – Amplitude vs. Phase:

(53)

Fourier-Transformation

Beispiele – Richtungen:

(54)

Fourier-Transformation

Beispiele – Erkennung der Richtung:

(55)

Faltungstheorem

F[f∗g] =F[f]· F[g]

F– Operator (Fourier-Transformation),

F(v) =F[f] – das Abbild der Funktionf im Fourier-Raum.

Beweis:

f(x) =F⁻¹[F(v)] =

Z

^∞

−∞

F(v)e^2πivxdv, g(x) = ... analog

f∗g=

Z

^∞

−∞

g(x⁰)f(x−x⁰)dx⁰=

=

Z

∞

−∞

g(x⁰)·

h Z

∞

−∞

F(v)e^2πiv(x−x⁰⁾dv

i

dx⁰=

=

Z

∞

−∞

F(v)·

h Z

∞

−∞

g(x⁰)e^−2πivx⁰dx⁰

i

e^2πivxdv=

=

Z

^∞

−∞

F(v)·G(v)e^2πivxdv=F⁻¹[F(v)·G(v)]

⇒ F[f∗g] =F[f]· F[g]

(56)

Faltungstheorem

Konsequenz 1:

Eine Faltungf∗gkann durchf∗g=F⁻¹[F(v)·G(v)] implementiert werden.

Zeitkomplexität –O(nlogn).

Konsequenz 2:

Jeder Filter hat seine Spektralcharakteristika im Fourier-Raum

→Spektralanalyse,

→Entwicklung der Filter mit bestimmten Spektralcharakteristika.

(57)

Andere Transformationen

Die Bilder sagen„Wo“, aber nicht„Was“.

Die Spektren sagen„Was“, aber nicht„Wo“.

Windowed Fourier-Transformation, Schnelle Fourier-Transformation (FFT – nur schnell) Cosine Transformation (1D, Diskret, DCT-II):

F(u) =

N−1

X

x=0

f(x)·cos

h

_π

N x+1 2

u

i

Wavelet Transformation (1D, Kontinuierlich):

F(a,b)∼

Z

^∞

−∞

f(x)·ψ x−b a

dx

mit der „Mutter“-Funktionψ(·)

(z.B. „Complex mexican hat wavelet“ etc.).

(58)

Bildverarbeitung: Diffusion Filters

D. Schlesinger () Bildverarbeitung: Diffusion Filters 1 / 10

(59)

Diffusion – Idee

Motiviert durch physikalische Prozesse – Ausgleich der Konzentration eines Stoffes.

Konzentration ist eine Funktion im Raum d.h.

u:Rⁿ→R, oft zum Beispielu:R³→Rin der Physik.

Räumlicher Gradientder Konzentration∇u= (_∂x^∂u

1,_∂x^∂u

2, . . .) verursacht

„flux“j:Rⁿ→Rⁿ (Massbewegung, Vektorfeld) – Ficksches Gesetz:

j=−D· ∇u,

Dist eine positiv definite symmetrische Matrix –Diffusion Tensor.

Aus der Erhaltung der Masse folgt (tist die Zeit)

∂u

∂t =−divj= div(D· ∇u) mitDivergenz: divj(x) =^∂j_∂x¹^(x)

1 +^∂j_∂x²^(x)

2 +. . .

(60)

Diffusion für Bilder

Das Bild wird als initiale Verteilung der Konzentration interpretiert:

u(x,y,t= 0) =I(x,y)

Das „Bild“ wird entsprechend ∂u/∂t= div(D∇u) mit der Zeit geändert.

Diffusion TensorDsteuert die Entwicklung der Verteilung der Konzentration in Zeit.

Fälle nachD:

skalar → isotropisch

allgemein → anisotropisch

unabhängig vonu → linear abhängig vonu → nichtlinear

⇒Alle vier Varianten möglich.

(61)

Lineare Isotropische Diffusion

Diffusion Tensor ist eine Konstante, d.h.D=c·1I (1I ist die Einheitsmatrix):

u(x,0) =I(x), ∂u

∂t = div(c· ∇u) =c· 4u mit demLaplace Operator4u= div(∇u) = ^∂u²

∂x∂x+_∂y∂y^∂u². Fragen:u(x,∞) = ? Stationär?

4u= 0 beit→ ∞ ⇒u(x,∞) = ? (Linear, Bilinear ... )

Das hängt voncund von Regularisierung am Rand ab. In meist der Fälle – Glättung.

Der einfachste Fall (Homogene Diffusion):Dist die Einheitsmatrix (oderc= 1).

Es existiert die analytische Lösung:

u(x,t) = (G^√_2t∗I)(x)

d.h. Faltung des AusgangsbildesI mit dem Gausschen Glättungskern mitσ=√ 2t

(62)

Finite Differenzen

Diskretisiere (Homogene Diffusion)

∂u(x,u,t)

∂t = ∂u(x,y,t)²

∂x∂x +∂u(x,y,t)²

∂x∂x Ableitungen (kontinuierlich)→Differenzen (diskret):

∂u(x,u,t)

∂t =u(x,y,t+τ)−u(x,y,t)

τ +O(τ)

∂u(x,y,t)²

∂x∂x =u(x+h,y,t)−2u(x,y,t) +u(x−h,y,t)

h² +O(h²)

∂u(x,y,t)²

∂y∂y =u(x,y+h,t)−2u(x,y,t) +u(x,y−h,t)

h² +O(h²)

τist der Zeitschritt,h– räumliche Auflösung.

O(·) werden vernachlässigt, alles wird zusammengesetzt und nachu(x,y,t+τ) umgestellt:

u(x,y,t+τ) =

1−4τ h²

u(x,y,t) + τ

h² u(x+1,y,t) +u(x−1,y,t) +u(x,y+1,t) +u(x,y−1,t)

Explizites Schema: die neuen Werte werden direkt aus den alten berechnet.

Stabil, wenn alle „Gewichte“ nichtnegativ, d.h. _h^τ₂ ≤¹

4.

(63)

Finite Differenzen

Implizites Schema: Divergenzen für dennächstenZeitpunkt werden verwendet u(x,y,t+τ) =

1−4τ h²

u(x,y,t) + τ

h² u(x+1,y,t) +u(x−1,y,t) +u(x,y+1,t) +u(x,y−1,t)

wird zu

u(x,y,t+τ) =u(x,y,t)+

τ

h² u(x+1,y,t+τ) +u(x−1,y,t+τ) +u(x,y+1,t) +u(x,y−1,t+τ)−4u(x,y,t+τ)

Die neuen Werte können nicht direkt berechnet werden, sondern es entsteht ein System linearer Gleichungen.

Sehr groß – so viel Gleichungen wie viel Pixel, dafür aber dünn besetzt

⇒spezielle iterative Methoden (Jakobi ...).

Explizit: instabil, schnell

Implizit: stabil, langsam (in jedem Zeitpunkt ein System lösen)

(64)

Lineare Isotropische Diffusion

Die Idee – bei Anwesenheit der Kanten weniger glätten c· 4u≡c(x,y,I)· 4u mitc(x,y,I) vorberechnet aus dem Bild.

Sehr oftc(x,y,I) = g(|∇I(x,y)|²) – eine positive fallende Funktion (Diffusivität) der quadratischen Länge des Bildgradienten.

(65)

Nichtlineare Isotropische Diffusion

Die Idee – Kanten sind besser im entrauschten Bild (unbekannt)

∂u

∂t = div(g(|∇I|²)∇u) wird zu ∂u

∂t = div(g(|∇u|²)∇u) Spezialfall – TV-flow:∂u/∂t= div(_|∇u|^∇u)

– Keine weitere Kontrastparameter

– Stückweise konstante Grauwertverläufe – Segmentierung ähnlich Problem:∞bei|∇u|= 0→Regularisierungg(s²) = √¹

s²+ε

(a) Verrauschtes Ausgangsbild (b) Gaussche Glättung (c) Nichtlineare Diffusion

Implizites Schema führt zum SystemnichtlinearerGleichungen in jedem Zeitpunkt.

(66)

Shock Filter

Die Idee: „dilation“ in der Nähe des Maximums und „erosion“ in der Nähe des Minimums

∂u

∂t =−sign(4u)· |∇u|

(67)

Shock Filter

(68)

Bildverarbeitung: Kontinuierliche Energieminimierung

D. Schlesinger BV: Kontinuierliche Energieminimierung() 1 / 9

(69)

Idee

Statt zu sagen, wie die Lösung geändert werden muss (explizite Algorithmus, Diffusion), werden die erwünschten Eigenschaften der Lösung explizit formuliert.

Die Ausprägungen eines Objektes werden durch Abbildungen repräsentiert.

Beispiele:

„Menge der Pixel“→„Menge der Farben“(alles diskret).

„Menge der Pixel“→„Kontinuierlicher Grauwertbereich“.

Eigenschaften des Modells werden mittels „Energie“ dargestellt – Funktion, die „ungünstige“ Abbildungen bestraft.

Die Aufgebe wird zu einem Optimierungsproblem – suche nach der günstigsten Abbildung.

Fälle:

Definitionsbereich: kontinuierlich –R² diskret – Menge der PixelR Wertebereich: kontinuierlich –R diskret – z.B. [0. . .255]

Heute: Wertebereich kontinuierlich, Beispiel – Entrauschen (denoising).

(70)

Diskreter Definitionsbereich

R∈Z²– die Pixelmenge,E⊂R²– die Nachbarschaftstruktur (z.B. 4-Nachbarschaft) x:R→Z– das Ausgangsbild,y:R→R– die gesuchte Abbildung (das restaurierte Bild).

EnergieE:R^|R|→Rbesteht (normalerweise) aus zwei Teilen:

1) Der Daten-Term:

Ed(y) =

X

r∈R

(xr−yr)² (entspricht Gausschem Rauschen).

2) Der Modell-Term:

Em(y) =

X

rr⁰∈E

(yr−y_r0)²

Annahme: der rekonstruierte Grauwertverlauf (die Abbildungy) soll glatt sein Die Optimierungsaufgabe:

y^∗= arg min

y

Ed(y) +αEm(y)

(71)

Diskreter Definitionsbereich

Lösungsweg – Ableitungen Nullsetzen:

∂

∂yr^∗

hX

r∈R

(xr−yr)²+α

X

rr⁰∈E

(yr−y_r0)²

i

=

yr^∗−xr^∗+α

X

r⁰:r^∗r⁰∈E

(yr^∗−y⁰_r) = 0

⇓

(1 + 4α)yij−αyij−1−αyij+1−αyi−1j−αyi+1j=xij ∀i,j System linearer Gleichungen mirn=|R|Variablen undnGleichungen:

A·y=x mit

y= (y1,y2, . . . ,yn)∈Rⁿ– die Lösung, x= (x1,x2, . . . ,xn)∈Zⁿ– das Ausgangsbild,

ak= 1 + 4α,akl=−αwenn die entsprechenden Pixel benachbart sind, sonst 0.

Das System kann bezüglichymithilfe Standardmethoden

(Gaussche Eliminierung,LU-Dekomposition,Ainvertieren usw.) gelöst werden – das ist aber leider sehr Zeitaufwendig (nur im 1D-Fall effizient).

(Merke: online/offline Varianten).

(72)

Diskreter Definitionsbereich

Die MatrixAist schwach besetzt→iterative Methoden.

Jacobi Methode:

Man zerlegeA=D+M mit einer diagonalen MatrixM:

Ay=x ⇔ (D+M)y=x ⇔ Dy=x−My ⇔ y=D⁻¹(x−My) y^(k+1)=D⁻¹(x−My^(k))

Vorteile: extrem einfach, parallelisierbar

Nachteile: immer noch zu langsam, konvergiert nur beik→ ∞, konvergiert nur wenn die Matrix streng diagonal dominant ist, d.h.|a_ii|>

P

j6=i|a_ij|, was glücklicherweise für das Beispiel gerade der Fall ist.

Andere Algorithmen:

Gauss-Seidel, Successive Over-relaxation (schneller), Konjugierte Gradienten (bessere Konvergenz), Multigrid Methoden (viel schneller aber komplizierter) etc.

Das obige Beispiel ist sehr einfach – quadratische Energie⇒lineare Gleichungssystem.

Probleme wenn nicht differenzierbar, nicht konvex etc.

(73)

Kontinuierlicher Definitionsbereich

Definitionsbereich wird zuR⊂R²,

die Abbildungy:R²→Rist somit eine Funktion, Die Energie wird zum EnergiefunktionalE:R^∞→R.

„Calculus of Variations“, Variationelle Ansätze.

Beispiel – das Entrauschen:

E(y) =

Z

R

h

y(r)−x(r)

2

+α|∇y(r)|²

i

dr→min

y

Gâteaux Ableitungen entlang „Richtungen“h:R²→R (Richtungsableitungen im Funktionsraum):

∂Eh(y)

∂y = lim

ε→0

E(y+εh)−E(y)

ε = dE(y+εh) dε

ε=0

Euler-Lagrange Gleichungen:

im Optimum sind alle Gâteaux Ableitungen (d.h. für alleh) Null.

(74)

Kontinuierlicher Definitionsbereich

d dε

Z

R

h

(y+εh−x)²+α|∇(y+εh)|²

i

dr

ε=0

= // koordinatenweise inR²

d dε

Z

R

h

(y+εh−x)²+α ∂

∂r1

(y+εh)

²

+α ∂

∂r2

(y+εh)

²

i

dr

ε=0

=

2

Z

R

h

(y+εh−x)h+α ∂

∂r1

(y+εh)∂h

∂r1

+α ∂

∂r2

(y+εh)∂h

∂r2

i

dr

ε=0

=

2

Z

R

h

(y−x)h+α ∂y

∂r1

∂h

∂r1

+α ∂y

∂r2

∂h

∂r2

i

dr= // partielle Integration

2

Z

R

h

(y−x)h−α ∂²y

∂r₁²h

−α ∂²y

∂r₂²h

i

dr+. . .(Grenzeffekte) =

2

Z

R

(y−x−α4y)h dr+. . .(Grenzeffekte) = 0 ∀h

⇓

(75)

Kontinuierlicher Definitionsbereich

⇒ y−x−α4y= 0 ∀r∈R , und für die Grenzen ∂hn,∇yi

∂R = 0

Relation zum Fall diskretes Definitionsbereiches:

Diskretisiert man die Bedingungen und schreibt sie für alle Pixel (i,j) auf, d.h.

yi,j−xi,j−α (yi−1,j−2yi,j+yi+1,j) + (yi,j−1−2yi,j+yi,j+1)

= 0 so entsteht dasselbe lineare Gleichungssystem wie beim diskreten Definitionsbereich:

yi,j(1 + 4α)−αyi−1,j−αyi+1,j−αyi,j−1−αyi,j+1=xi,j ∀(i,j).

Relation zur Diffusion:

(Anti)Gradient Verfahren zur Minimierung der EnergieE(y):

y^(t+1)=y^(t)−∂E(y)

∂y =y^(t)+α4y+ (x−y) Vergleiche mit der linearen isotropischen Diffusion:

u^(t+1)=u^(t)+∂u

∂t =u^(t)+c4u

(76)

Kontinuierlicher Definitionsbereich

Erweiterungen (kompakte Schreibweise):

E(y) =

Z

R

h

(y−x)²+αΨ(|∇y|²)

i

dr→min

y

mit einem Regularisator Ψ:

Ψ(s²) =s² – Tikhonov

Ψ(s²) =

√

s² – Total Variation

Ψ(s²) = 1−λ²exp(−^s²

λ²) – Perona-Malik Ψ(s²) =

n

₀ _wenn _s2= 0

1 else – Potts-Modell

Euler-Lagrange Gleichungen:

div Ψ⁰(|∇y|²)∇y

−y−x α = 0

(77)

Bildverarbeitung: Diskrete Energieminimierung

D. Schlesinger () BV: Diskrete Energieminimierung 1 / 11

(78)

Entrauschung → Segmentierung

Beide Definitionsbereich und Wertebereich sind diskret.

R∈Z²– die Pixelmenge,E⊂R²– die Nachbarschaftstruktur (z.B. 4-Nachbarschaft) x :R →Z– das Ausgangsbild,y :R→K – die gesuchte Abbildung (das restaurierte Bild).k∈Krepräsentiert den „wahren“ Grauwert (Label).

Die Energieminimierung:

y^∗= arg min

y

Ed(y) +αEm(y)

z.B. arg min

y

hX

r∈R

(xr−yr)²+α

X

rr⁰∈E

(yr−yr⁰)²

i

Farben→semantische Bedeutungen (Werte eines Merkmals)...

Die Menge der Pixel ist auf „sinnvolle“ Teilmengen zu partitionieren.

(79)

Segmentierung

Original

A possible segmentation

r r r r

r r

Data terms Compactness terms

Penalty Zero

k= 3: Shadow k= 2: Forest k= 1: Field

Dissimilarity measure Observed features

y^∗= arg min

y

hX

r

qr(yr) +

X

rr⁰

g_rr0(yr,y_r0)

i

(80)

Iterated Conditional Modes

y^∗= arg min

y

hX

r

qr(yr) +

X

rr⁰

g_rr0(yr,y_r0)

i

Die Idee: wähle (lokal) immer wieder das energetisch günstigste Label bei fixiertem Rest [Besag, 1986].

Wiederhole oft für aller: yr= arg min

k

h

qr(k) +

X

r⁰:rr⁰∈E

g_rr0(k,y_r0)

i

(ME: synchrone Dynamik in Hopfield-Netzen)

+ Extrem einfach, parallelisierbar.

− „Koordinatenweise“ Optimierung

→konvergiert nicht zum globalen Optimum selbst bei einfachen Modellen.

(81)

Iterated Conditional Modes

Erweiterung: fixiere nicht alle Variablen bis auf eine, sondern nur eine Teilmenge so, das der Rest einfach optimierbar ist (zum Beispiel eine Kette oder ein Baum).

Für Bilder – Zeilenweise/Spaltenweise Optimierung.

→durch Dynamische Programmierung exakt und effizient lösbar.

(82)

Äquivalente Transformationen (Reparametrisierung)

Zwei AufgabenA= (q,g) undA⁰= (q⁰,g⁰) sind zu einanderäquivalent, wenn

hX

r

qr(yr) +

X

rr⁰

g_rr0(yr,y_r0)

i

=

hX

r

q_r⁰(yr) +

X

rr⁰

g⁰_rr0(yr,y_r0)

i

für alle Labellingsygilt.

A(A) – Äquivalenzklasse (alle zuAäquivalenten Aufgaben).

ÄquivalenteTransformationen:

Φ = ϕr(k)∀r,k, ϕrr⁰(k), ∀rr⁰,k

ϕr(k) +

X

r⁰:rr⁰∈E

ϕ_rr0(k) = 0 ∀r,k