(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), kde řádky matice $\_X\in\R^{n\times 2}$ jsou body $(x_i, y_i)$ a $\_d\in\R^{1\times n}$ je řádkový vektor orientovaných (uvnitř kružnice záporné, vně kladné) vzdáleností bodů od kružnice se středem $(x_0, y_0)$ a poloměrem $r$.
• [x0 y0 r] = quad_to_center(d,e,f), která přepočítá reprezentaci kružnice rovnicí \eqref{circ-nhom} na reprezentaci rovnicí \eqref{circ}.
• [d e f] = fit_circle_nhom(X), kde $\_X$ je matice s body a $(d,e,f)$ je optimální řešení úlohy \eqref{fit-nhom}.
• [d e f] = fit_circle_hom(X), kde $\_X$ je matice s body a $(d,e,f)$ je optimální řešení úlohy \eqref{fit-hom}, ve kterém čísla $d,e,f$ jsou vydělená číslem $a$ ve shodě s rovnicí \eqref{circ-nhom}

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:

• Vyberte náhodně 3 body z bodů $(x_1,y_1),\dots,(x_n,y_n)$ a najděte kružnici jimi procházející.
• Spočítejte vzdálenosti všech bodů $(x_1,y_1),\dots,(x_n,y_n)$ od této kružnice. Pokud je vzdálenost bodu od kružnice menší než $\theta$, prohlašte ho za inlier.

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.

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.

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: