Skip to content

Zpracování obrazu pomocí minimalizace energie

octicon-noteNote

Variační filtrování. Aktivní křivky a plochy (geodetický model, Chan-Vese model). Minimalizace pomocí grafových řezů. Variační optický tok.
PA166

Variační filtrování neřeší lokálně jednotlivé pixely, ale snaží se minimalizovat enetgii celého obrázku.

Typickým předpisem takovéhle funkce je

Ef(u)=Ω((uf)2similarity+αΨ(u2)smoothness)dxE_f(u) = \int_{\Omega} ( \underbrace{(u - f)^2}_\text{similarity} + \alpha \underbrace{\Psi (| \nabla u |^2)}_\text{smoothness} ) dx

Chceme, aby byl obrázek co nejpodobnější původnímu (data/similarity term), ale zároveň co nejvíc vyhlazený (smoothness term/regularizer/penaliser). Poměr těchto dvou členů určuje parametr α\alpha. Regularizer Ψ\Psi určuje “metodu” vyhlazení:

Tikhonov (linear) — identitaΨ(u2)=u2\Psi(\lvert \nabla u \rvert^2) = \lvert \nabla u \rvert^2
CharbonnierΨ(u2)=2λ21+u2/λ22λ2\Psi(\lvert \nabla u \rvert^2) = 2 \lambda^2 \sqrt{1 + \lvert \nabla u \rvert^2 / \lambda^2} - 2 \lambda^2
Total VariationΨ(u2)=2u2=2u\Psi(\lvert \nabla u \rvert^2) = 2\sqrt{\lvert \nabla u \rvert^2} = 2\lvert\nabla u\rvert

Tuto funkci můžeme minimalizovat dvěma způsoby:

Provedeme diskretizaci a poté parciálně zderivujeme pro jednotlivé pixely. Tím dostáváme sadu lineárních funkcí položených nule. Tyto poté položíme do matice, kterou řešíme. Matice má hodnoty pouze kolem diagonály, což znamená, že je velmi rychle řešitelná např. pomocí Thomas method.

Pokud je obraz 2D, bereme všechny pixely jako dlouhou 1D posloupnost a řešíme obdobně.

Matice je ale obrovská, proto můžeme využít např. Gauss-Seidel metodu, s pomocí které matici řešíme iterativně po řádcích.

Použijeme Variační kalkulus, který dokáže minimalizovat Funkcionály

  • Funkcionál
    Bere na vstupu funkci a přiřazuje jí nějakou hodnotu u(x)Ru(x) \Rightarrow \mathbb{R}.

Pro řešení můžeme nyní využít Euler-Lagrangeovy rovnice. Pro funkcionál Ef(u)=Ω((uf)2+αΨ(u2))dxdyE_f(u) = \int_{\Omega} ( (u - f)^2 + \alpha \Psi (| \nabla u |^2) ) dxdy získáme předpis pro minimalizaci:

0=ufαdiv(Ψ(u2)u)0 = u - f - \alpha \text{div}(\Psi'(\lvert \nabla u \rvert^2) \nabla u)

Nyní můžeme hezky vidět vztah k difuznímu filtrování (otázka PGV10):

Variační metodaufα=div(Ψ(u2)u)\frac {u - f}{\alpha} = \text{div}(\Psi'(\lvert \nabla u \rvert^2) \nabla u)
Difuzní filtrovánítu=div(Ψ(u2)u)\partial_t u = \text{div}(\Psi'(\lvert \nabla u \rvert^2) \nabla u)

Přičemž vidíme, že

  • Tikhonov regulariser ~ Linear diffusivity
    • Ψ(u2)=u2\Psi(\lvert \nabla u \rvert^2) = \lvert \nabla u \rvert^2
    • Ψ(u2)=1\Psi'(\lvert \nabla u \rvert^2) = 1
  • Charbonnier regulariser ~ Charbonnier diffusivity
    • Ψ(u2)=2λ21+u2/λ22λ2\Psi(\lvert \nabla u \rvert^2) = 2 \lambda^2 \sqrt{1 + \lvert \nabla u \rvert^2 / \lambda^2} - 2 \lambda^2
    • Ψ(u2)=11+u2/λ2\Psi'(\lvert \nabla u \rvert^2) = \frac{1}{\sqrt{1 + \lvert \nabla u \rvert^2 / \lambda^2}}
  • Total Variation regulariser ~ Total Variation diffusivity s difusivitou 1u\frac{1}{\lvert \nabla u \rvert}
    • Ψ(u2)=2u\Psi(\lvert \nabla u \rvert^2) = 2\lvert\nabla u\rvert
    • Ψ(u2)=1u\Psi'(\lvert \nabla u \rvert^2) = \frac{1}{\lvert \nabla u \rvert}
  • Segmentace
    Chceme rozdělit obrázek podle něčeho (region based, edge based, texture, …)
  • Segmentace založená na energii
    Hledáme konturu s minimální energií. Závisí na volbě iniciální křivky - hledá pouze lokální minimum.

Působí dvě energie:

  • Vnitřní - tvar a plynulost křivky
  • Vnější - omezení z obrázku

Kontura zadaná, jako parametrická křivka p(s):[0,1]R    p(s)=(x(s),y(s))Tp(s): [0, 1] \rightarrow \mathbb{R}\ \ \ \ p(s) = (x(s), y(s))^T. Podobně, jako u variačního filtrování, definujeme funkcionál energie, který se snažíme minimalizovat:

E(p)=01Eint(p(s))+Eext(p(s))dsE(p) = \int_0^1 E_{int}(p(s)) + E_{ext}(p(s)) ds

Kde EintE_{int} je vnitřní energie, kde α\alpha a β\beta jsou váhy pro délku (α\alpha) a křivost (β\beta) křivky:

Eint(p(s))=α(s)p(s)2+β(s)p(s)2E_{int}(p(s)) = \alpha(s) \lvert p'(s) \rvert^2 + \beta(s) \lvert p''(s) \rvert^2

A EextE_{ext} je vnější energie, kde PP je potenciální funkce:

Eext(p(s))=P(p(s))E_{ext}(p(s)) = P(p(s))
  • Potenciál
    Potenciál určuje vztah křivky k obrázku. Může být definován různými způsoby:

    • Hrany - P(x,y)=f(x,y)2P(x, y) = -\lvert \nabla f(x, y) \rvert^2
    • Hrany po rozmazání - P(x,y)=(Gσ(x,y)f(x,y))2P(x, y) = -\lvert \nabla (G_{\sigma}(x, y) * f(x, y) )\rvert^2
    • Linie - P(x,y)=f(x,y)P(x, y) = -f(x, y)
    • Linie po rozmazání - P(x,y)=Gσ(x,y)f(x,y)P(x, y) = -G_{\sigma}(x, y) * f(x, y)

Minimalizaci funkcionálu chceme opět udělat Euler-Lagraneovou rovnicí, tedy:

αp(s)+βp(s)P=0, kde p(0),p(0),p(1),p(1) jsou zadaneˊ.\alpha p''(s) + \beta p'''(s) - \nabla P = 0 \text{, kde } p(0), p'(0), p(1), p'(1) \text{ jsou zadané}.

Na derivace můžeme použít finite differences a přepíšeme výsledek do matice AX=FextAX = F_{ext}. Blbý ale je, že F*{ext} závisí na p a p závisí na F*{ext}. Proto se to řeší iterativně pt=Fint(p)+Fext(p)p(s,0)=p0(s)\frac{\partial p}{\partial t} = F_{int}(p) + F_{ext}(p)\\p(s, 0) = p_0(s)

  • Problémy

  • Reparametrizace

    • Body mají tendenci se slévat do míst s nízkou potenciální energií
    • Musíme jednou za čas rozmístit body pravidelně podél vzniklé křivky
  • Normalizace (Diskretizace času)

    • Velikost gradientu není ničím omezena a může nám “vystřelit” body úplně do kša
    • Místo klasického gradientu vezmeme gradient normalizovaný, tedy jen směr
  • Diskretizace prostoru

    • Potenciální funkce je spočítaná pouze v diskrétních bodech, ale křivka může procházet i mimo tyto body
    • Na dopočítání použijeme například bilineární interpolaci
  • Problémy předchozí detekce hran

    • Výsledkem detekce hran je často binární obrázek, kde by ale křivka nevěděla, kam se pohybovat - kuličky jsou na rovině

    • Můžeme obrázek Gaussovsky rozmazat, nebo dopočítat vzdálenost všech bodů od bodů hran (distance function from the edge pixels)

      width=400

  • Balonování

    • Kontury mají tendenci se slévat do středu do jednoho bodu, kde zaniknou

    • Často chceme detekovat něco, kde začneme uvnitř

    • K externím silám přidáme tzv. balonovou sílu ve směru normály křivky

      Fext=k1n(s)kPPF_{ext} = k_1 n(s) - k \frac{\nabla P}{|\nabla P|}

    • Parametr k1k_1 ovlivňuje, jestli se křivka “nafukuje”, nebo “vyfukuje”, můžeme měnit během iterování

  • GVF Snakes

  • Kontury mají problém opsat “tvar rohlíku” a nic je nenutí dostávat se do úzkých prostor

    width=600

  • Převedeme externí síly na nové pole sil, kterému říkáme Gradient Vector Flow

  • Chceme v podstatě přidat síly tam, kde nejsou hrany, aby celý prostor k hranám směřoval

    width=600

Hloupé na hadech (to, co jsem popisoval doteď) je, že závisí na parametrizaci (rozložení bodů na křivce). GAC naopak na parametrizaci nezávisí.

Pro křivku minimalizujeme funkcionál (Chceme co nejmenší délku křivky && co nejnižší celkový gradient): +

EGAC(C)=01g(GσI(C(q))Vyhlazenyˊ obraˊzek)C(q)Deˊlka krˇivkydqE_{GAC}(C) = \int_0^1 g(\underbrace{\lvert \nabla G_{\sigma} * I(C(q)) \rvert}_{\text{Vyhlazený obrázek}}) \underbrace{\lvert C'(q) \rvert}_{\text{Délka křivky}} dq

Funkce gg snižuje vliv délky na energii, pokud leží v místech s vysokým gradientem, protože je nepřímo úměrná rozmazanému gradientu: +

g(GσI(x,y))=11+GσI(x,y)g(\lvert \nabla G_{\sigma} * I(x, y) \rvert) = \frac{1}{1 + \lvert \nabla G_{\sigma} * I(x, y) \rvert}

Do teď jsme neřešili, co je uvnitř kontur (homogenita), jen co je přímo pod nimi. Chan-Vese model se snaží najít konturu, která rozdělí obrázek na dvě části s různými intenzitami.

ECV(CKrˇivka,c1Odhadpoprˇedıˊ,c2Odhadpozadıˊ)=μL(C)Deˊlka krˇivky+λ1Ω1(f(x)c1)2dxOdlisˇnost poprˇedıˊ od odhadu+λ2Ω2(f(x)c2)2dxOdlisˇnost pozadıˊ od odhaduE_{CV}(\underbrace{C}_{\text{Křivka}}, \underbrace{c_1}_{\substack{\text{Odhad} \\ \text{popředí}}}, \underbrace{c_2}_{\substack{\text{Odhad} \\ \text{pozadí}}}) = \underbrace{\mu L(C)}_{\text{Délka křivky}} + \underbrace{\lambda_1 \int_{\Omega_1} (f(x) - c_1)^2 dx}_{\text{Odlišnost popředí od odhadu}} + \underbrace{\lambda_2 \int_{\Omega_2} (f(x) - c_2)^2 dx}_{\text{Odlišnost pozadí od odhadu}}
  • Regularity term - délka dělící křivky
  • Data term (2x) - rozdíl částí (fg,bg; Ω1,Ω2\Omega_1, \Omega_2) od konstantního odhadu c1,c2c_1, c_2
  • Parametry μ,λ1,λ2\mu, \lambda_1, \lambda_2 určují vyhlazení, variabilitu pozadí a popředí respectively
  • Pro danou křivku C jsou optimální c1,c2c_1, c_2 s průměrnou hodnotou

Abychom se zbavili integrálů přes část oblasti, můžeme použít Harveside funkci, která vrací 0 pro nulu a 1 jinak. Tím dostaneme:

ECV(C,c1,c2)=μΩH(u(x))dx Deˊlka krˇivky+λ1Ω(f(x)c1)2H(u(x))dx Odlisˇnost poprˇedıˊ od odhadu+λ2Ω(f(x)c2)2(1H(u(x)))dx Odlisˇnost pozadıˊ od odhadu\begin{aligned} E_{CV}(C, c_1, c_2) &= \mu \int_{\Omega} \lvert \nabla H (u(x)) \rvert dx &\footnotesize{\text{ Délka křivky}}\\ &+ \lambda_1 \int_{\Omega} (f(x) - c_1)^2 H(u(x)) dx &\footnotesize{\text{ Odlišnost popředí od odhadu}}\\ &+ \lambda_2 \int_{\Omega} (f(x) - c_2)^2 (1 - H(u(x))) dx &\footnotesize{\text{ Odlišnost pozadí od odhadu}}\\ \end{aligned}

Diskretizujeme a řešíme.

Rozdělení metodou minimálního řezu využívá algoritmu pro hledání minimálního řezu v grafech. Uvažme graf, kde každý pixel je jeden vrchol a pro všechny pixely vytvoříme hrany ze zdroje (s) a spotřebiče (t) a hrany mezi sousedy.

width=600

  • Riemannovy metrické systémy
    Normálně, když počítáme délku mezi dvěma body, počítáme ji, jako d=uTud = \sqrt{u^T \cdot u}. Můžeme ale přidat tzv. metrický tensor MM, který může tuto délku upravit v závislosti na směru d=uTMud = \sqrt{u^T \cdot M \cdot u}. Díky tomu potom můžeme měřit nejen euklidovskou vzdálenost, ale třeba:

    • Euklidovskou vzdálenost - M=IM = I (identita)
    • Geodetickou vzdálenost - M=diag(g(I()))M = \text{diag}(g(\lvert \nabla I(\cdot) \rvert))
    • Anizotropní vzdálenost - M=g(I())E+(1g(I()))uuTM = g(\lvert \nabla I(\cdot) \rvert) \cdot E + (1 - g(\lvert \nabla I(\cdot) \rvert)) \cdot u \cdot u^T (kde u=IIu = \frac{\nabla I}{\lvert \nabla I \rvert} jednotkový vektor ve směru gradientu a EE je identita)

Pro aproximaxi vah hran v našem grafu můžeme použít Riemannovu metriku:

wkR=wkϵdetM(ukTMuk)pw_k^R = w_k^{\epsilon} \cdot \frac{\text{det} M}{(u_k^T \cdot M \cdot u_k)^{p}}

Můžeme rozdělit obrázek na dvě části pomocí minimálního řezu na dvě části podle geodetické vzdálenosti.

  • t-linky (hrany do s nebo do t) - nastavíme váhy na 0 vyjma pixelů s omezením (víme, že jsou popředí, nebo pozadí), které nastavíme na nekonečno (nebo velkou konstantu), čímž zabráníme výběr do minimálního řezu
  • n-linky (mezi sousedními pixely) - podle Riemannovy metriky

Můžeme rozdělit obrázek na dvě části, přičemž o některých pixelech víme, jestli patří do popředí, nebo do pozadí. Podobně, jako u Chan-Vese křivkami máme odhad c1,c2c_1, c_2 pro popředí a pozadí. Nastavme tedy váhy:

  • t-linky - wsi=λ2(f(i)c2)2wit=λ1(f(i)c1)2w_{si} = \lambda_2(f(i) - c_2)^2 \text{, } w_{it} = \lambda_1(f(i) - c_1)^2
  • n-linky - podle Riemannovy metriky

Bohužel wsi,witw_{si}, w_{it} jsou závislé na c1,c2c_1, c_2. Typicky řešíme tak, že si tipneme, provedeme iteraci, aktualizujeme a jedem znovu

Pro sekvenci obrázků f(x,y,t)f(x,y,t) chceme získat vektor posunutí (u(x,y,t),v(x,y,t))T(u(x,y,t),v(x,y,t))^T každého místa obrázku.

Máme 2 předpoklady, které platí pro všechny metody:

  • Zachování jasu (mezi snímky se nám nemění jas obrázku) = BCA - f(x+u,y+v,t+1)=f(x,y,t)f(x+u, y+v, t+1)=f(x,y,t) (posunutý pixel se rovná neposunutému)
  • Malé změny (objekty se posouvají jen o několik málo (třeba 5) pixelů) - Abychom mohli provést linearizaci Taylorovým rozvojem

Z rovnice pro BCA můžeme linearizací dostat

0=f(x+u,y+v,t+1)f(x,y,t)fx(x,y,t)u+fy(x,y,t)v+ft(x,y,t)0 = f(x+u, y+v, t+1) - f(x,y,t) \\ \approx f_x(x,y,t)u + f_y(x,y,t)v + f_t(x,y,t)

Tím dostáváme nejdůležitější rovnici linearised optical flow constraint (OFC) fxu+fyv+ft=0f_x u+f_y v + f_t = 0 (f(u,v)+ft=0\nabla f \cdot (u,v) + f_t = 0)

Problém Apertury = OFC ukazuje jen tok rovnoběžný s gradientem (protože jinak to není vidět)

width=500

Taky protože fxu+fyv+ft=0f(u,v)+ft=0f_x u+f_y v + f_t = 0 \sim \nabla f \cdot (u,v) + f_t = 0 , tedy máme projekci reálného směru na gradient
-> Nonuniqueness = Můžeme libovolně přidat část toku kolmou na gradient, aniž bychom změnili výsledek

Přidává 3. předpoklad:

  • Zachování jasu
  • Malé změny
  • Konstantnost toku v okolí BρB_\rho každého bodu

Chceme minimalizovat následující lokální energii (integrál přes okolí)

E(u,v)=12Bρ(x0,y0)(fxu+fyv+ft)2dxdyE(u,v) = \frac{1}{2} \int_{B_\rho(x_0, y_0)} (f_x u + f_y v + f_t)^2 dx dy

Parciálně zderivujeme a převedeme na matice

0=Eu=Bρfx(fxu+fyv+ft)dxdy=00=Ev=Bρfy(fxu+fyv+ft)dxdy=0(Bρfx2dxdyBρfxfydxdyBρfxfydxdyBρfy2dxdy)(uv)=(BρfxftdxdyBρfyftdxdy)\begin{aligned} 0 = \frac{\partial E}{\partial u} &= \int_{B_\rho} f_x(f_x u + f_y v + f_t) dx dy = 0 \\ 0 = \frac{\partial E}{\partial v} &= \int_{B_\rho} f_y(f_x u + f_y v + f_t) dx dy = 0 \end{aligned} \\ \begin{pmatrix} \int_{B_\rho} f_x^2 dxdy & \int_{B_\rho} f_x f_y dxdy \\ \int_{B_\rho} f_x f_y dxdy & \int_{B_\rho} f_y^2 dxdy \end{pmatrix} \begin{pmatrix} u \\ v \end{pmatrix} = \begin{pmatrix} -\int_{B_\rho} f_x f_t dxdy \\ -\int_{B_\rho} f_y f_t dxdy \end{pmatrix}

Můžeme nahradit integrály za Gaussovské rozmazání, čímž dostaneme “plynulejší okolí”

(Gσ(fx2)Gσ(fxfy)Gσ(fxfy)Gσ(fy2))(uv)=(Gσ(fxft)Gσ(fyft))\begin{pmatrix} G_{\sigma} * (f_x^2) & G_{\sigma} * (f_x f_y) \\ G_{\sigma} * (f_x f_y) & G_{\sigma} * (f_y^2) \end{pmatrix} \begin{pmatrix} u \\ v \end{pmatrix} = \begin{pmatrix} -G_{\sigma} * (f_x f_t) \\ -G_{\sigma} * (f_y f_t) \end{pmatrix}

Snadno vyřešíme Cramerovým pravidlem

Výhody

  • Jednoduché, rychlé
  • Stačí nám 2 snímky

Nevýhody

  • Problémové situace, kde neplatí podmínka konstantnosti - rotace, protijedoucí auta, …
  • Nepočítá husté pole toků (jen tam, kde víme gradient)

Pouze dva přepoklady:

  • Zachování jasu - fxu+fyv+ft=0f_x u+f_y v + f_t = 0
  • Malé změny (Gradient vektorového pole je malý) - Ω(u2+v2)dxdy\int_{\Omega} ( \lvert \nabla u \rvert^2 + \lvert \nabla v \rvert^2) dx dy je malé

Kombinací dostáváme funkcionál energie optického toku

EHS(u,v)=Ω((fxu+fyv+ft)2+α(u2+v2))dxdyE_{HS}(u, v) = \int_{\Omega} ((f_x u + f_y v + f_t)^2 + \alpha (\lvert \nabla u \rvert^2 + \lvert \nabla v \rvert^2)) dx dy

Tento funkcionál minimalizujeme na jediné globální řešení

Výhody

  • Jsme schopni řídit smoothness pomocí α\alpha
  • Produkuje dense flow fields

Nevýhody

  • Stále musíme splňovat předpoklad malých změn
  • Celkem pomalé a drahé