Search
\[ \def\_#1{\mathbf{#1}} \def\R{\mathbb{R}} \def\dist{\mathop{\text{dist}}} \def\iverson#1{[\![#1]\!]} \]
(V.Voráček)
Hledejme kružnici, která 'nejlépe' prokládá dané body $(x_1,y_1),\ldots,(x_n,y_n)\in\R^2$. Všimněte si, že tato úloha není přesně specifikovaná, protože slovem 'nejlépe' může každý myslet něco trochu jiného. Budeme řešit dvě podúlohy, z nichž každá bude odpovídat jiné formalizaci tohoto slova.
V jedné možné formalizaci úlohy hledáme kružnici, která minimalizuje součet čtverců vzdálenosti daných bodů od této kružnice. Kružnice se středem $(x_0,y_0)$ a poloměrem $r$ je popsána rovnicí \begin{equation}\label{circ} (x-x_0)^2+(y-y_0)^2=r^2 . \end{equation} To odpovídá úloze \begin{equation}\label{fit} \min_{x_0,y_0,r} \; \sum_{i=1}^n \dist(x_i,y_i,x_0,y_0,r)^2 \end{equation} kde $\dist(x,y,x_0,y_0,r)$ označuje vzdálenost bodu $(x,y)$ od kružnice (vymyslete, jak se tato vzdálenost spočítá!). Funkce \eqref{fit} má ale obvykle mnoho lokálních minim a budeme ji (lokálně) minimalizovat v příští domácí úloze iteračními metodami. Teď budeme hledat přibližné minimum funkce \eqref{fit}.
Za tím účelem popíšeme kružnici trochu jinou rovnicí než \eqref{circ}. Kružnice je speciální případ kuželosečky, která je popsána rovnicí \begin{equation}\label{conic} ax^2 + by^2 + cxy + dx + ey + f = 0. \end{equation} Pro kružnici máme $a=b$ a můžeme položit $c=0$, tedy \eqref{conic} se zjednoduší na \begin{equation}\label{circ-hom} ax^2 + ay^2 + dx + ey + f = 0. \end{equation} Protože $a\neq0$, můžeme rovnici \eqref{circ-hom} vydělit číslem $a$, čímž dostaneme \begin{equation}\label{circ-nhom} x^2 + y^2 + dx + ey + f = 0 \end{equation} (kde tedy $d,e,f$ značí jiná čísla než v \eqref{circ-hom}). Uvědomte si, že \eqref{circ} se liší od \eqref{circ-hom} a \eqref{circ-nhom} doplněním na čtverec (viz kapitola o kvadratických funkcích ve skriptech).
Levá strana rovnice \eqref{circ-hom} příp. \eqref{circ-nhom} je nulová právě když bod $(x,y)$ leží na kružnici. Proto můžeme doufat, že neuděláme velkou chybu, když úlohu \eqref{fit} nahradíme úlohou \begin{equation}\label{fit-nhom} \min_{d,e,f} \; \sum_{i = 1}^n (x_i^2 + y_i^2 + dx_i + ey_i +f)^2 \end{equation} nebo úlohou \begin{equation}\label{fit-hom} \min_{a,d,e,f} \; \sum_{i = 1}^n (ax_i^2 + ay_i^2 + dx_i + ey_i + f)^2 \quad\text{za podmínky}\quad a^2+d^2+e^2+f^2=1 . \end{equation} Úloha \eqref{fit-nhom} resp. \eqref{fit-hom} je přibližné řešení přeurčené nehomogenní příp. homogenní lineární soustavy, viz skripta. Všimněte si, že podmínku v úloze \eqref{fit-hom} nemůžeme vynechat.
Uvědomte si, že výraz $|x^2 + y^2 + dx + ey +f|$ příp. $|ax^2 + ay^2 + dx + ey +f|$ (v počítačovém vidění nazývaný algebraická vzdálenost) není přesně roven (či přímo úměrný) vzdálenosti bodu $(x,y)$ od kružnice, je to jen aproximace této vzdálenosti.
Jako úkoly napište tyto funkce:
d = dist(X, x0, y0, r)
[x0 y0 r] = quad_to_center(d,e,f)
[d e f] = fit_circle_nhom(X)
[d e f] = fit_circle_hom(X)
Následně experimentálně porovnejte obě metody. To znamená, že pro optimální řešení úloh \eqref{fit-nhom} a \eqref{fit-hom} spočítáte hodnotu kritéria původní neaproximované úlohy \eqref{fit}. Rozdíl mezi metodami by měl být patrný zejména když body jsou podél optimální kružnice rozmístěny nerovnoměrně.
Formulace \eqref{fit} předpokládá, že pozorované body $(x_1,y_1),\dots,(x_n,y_n)$ všechny leží na nějaké kružnici, až na (i.i.d. gaussovské) chyby v měření. Může se ale stát, že jen část bodů (tzv. inliers) patří kružnici a zbylá část bodů s hledanou kružnicí vůbec nesouvisí (tzv. outliers, česky vychýlené body). Kdybychom z takových dat chtěli odhadnout kružnici metodami popsanými výše, odhadnutá kružnice by se mohla libovolně lišit od skutečné kružnice. Pro takové případy je třeba použít metody robustního odhadování. V této části se zaměříme na jednoduchou a široce užívanou robustní metodu, známou jako RANSAC (RAndom SAmple Consensus).
Úloha \eqref{fit} měla tvar \begin{equation}\label{fit-rho} \min_{x_0,y_0,r} \; \sum_{i=1}^n \rho(\dist(x_i,y_i,x_0,y_0,r)) \end{equation} kde funkce $\rho{:}\ \R\to\R$ byla $\rho(t)=t^2$. Toto je obecný tvar tzv. M-estimátoru (estimace = odhad parametrů modelu, zde kružnice). Robustní M-estimátory používají jiné funkce $\rho$. My použijeme funkci \begin{equation}\label{rho} \rho(t) = \begin{cases} -1 & \text{pro } |t|\le\theta \\ 0 & \text{pro } |t|>\theta \end{cases} \end{equation} kde $\theta>0$ je zvolený práh. Vyřešit přesně úlohu \eqref{fit-rho} pro tuto volbu $\rho$ je ovšem obtížné (účelová funkce je nekonvexní a dokonce nespojitá), proto se uchylujeme k přibližnému řešení algoritmem RANSAC.
Algoritmus předpokládá, že optimální kružnice prochází nějakou trojicí z bodů $(x_1,y_1),\dots,(x_n,y_n)$ (uvědomte si, že kružnice je jednoznačně definovaná třemi nekolineárními body). Ovšem nezkouší všech $n\choose3$ trojic, ale jen $k$ náhodně vybraných trojic, kde $k\ll{n\choose3}$ je zvolené podle relativního množství inlierů. Přesně, algoritmus opakuje $k$-krát tuto operaci:
Přibližné optimum úlohy \eqref{fit-rho} je ta kružnice, která v průběhu algoritmu měla nejvíce inlierů (této množině inlierů říkáme concensus set).
Úkol: Implementujte tento algoritmus ve funkci [x0 y0 r] = fit_circle_ransac(X,num_iter,threshold), kde $\_X$ je matice s body, $(x_0,y_0,r)$ jsou parametry kružnice s nejvíce inliery. proměnné num_iter a threshold odpovídají v popisu $k$ a $\theta$ respektive.
[x0 y0 r] = fit_circle_ransac(X,num_iter,threshold)
num_iter
threshold
Poznámka: Nalezenou nejlepší množinu inlierů můžeme později použít pro přesnější odhad kružnice nerobustními metodami z prvního podúkolu nebo z příští domácí úlohy.
Templaty, včetně skriptu pro testování si stáhněte zde: template pro matlab, template pro python. Pro python spouštějte skript circlefit.py, pro matlab circlefit.m.
circlefit.py
circlefit.m
Obrázky ukazují výsledky prokládání kružnice správně implementovanými funkcemi fit_circle_nhom, fit_circle_hom a fit_circle_ransac:
fit_circle_nhom
fit_circle_hom
fit_circle_ransac