\[
\def\_#1{\mathbf{#1}}
\def\R{\mathbb{R}}
\def\dist{\mathop{\text{dist}}}
\def\iverson#1{[\![#1]\!]}
\]

===== Prokládání bodů kružnicí =====

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.


==== Minimalizace algebraické vzdálenosti ==== 

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ě.

/*
Obecněji lze ukázat, že když souřadnice $x_i,y_i$ bodů mají velmi různé absolutní hodnoty a jejich těžiště je daleko od počátku, tak úlohy \eqref{fit-nhom} a \eqref{fit-hom} mohou být velmi špatnou aproximací úlohy \eqref{fit}. Proto je běžnou praxí, že se před řešením úloh body posunou a oškálují, tj.
\begin{align*}
x_i &:= \alpha x_i+u,\\
y_i &:= \alpha y_i+v,
\end{align*}
tak, aby $\sum_ix_i=\sum_iy_i=0$ a $\max_i\max\{x_i,y_i\}=-\min_i\min\{x_i,y_i\}=1$.
*/


==== Robustní prokládání metodou RANSAC ====

/*
Lze ukázat, že úloha \eqref{fit} je statisticky optimální odhad kružnice (přesněji, odhad ve smyslu maximální věrohodnosti), jestliže pozorované body $(x_1,y_1),\dots,(x_n,y_n)$ vznikly jako
\begin{align*}
x_i &= \bar x_i + \epsilon_x \\
y_i &= \bar y_i + \epsilon_y
\end{align*}
kde 
$(\bar x_1,\bar y_1),\dots,(\bar x_n,\bar y_n)$ jsou (nepozorovatelné) body ležící na nějaké (neznámé) kružnici a $\epsilon_x,\epsilon_y$ jsou náhodné proměnné (chyby měření, šum) s i.i.d. normálním (neboli gaussovským) rozdělením.
*/

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 //[[https://en.wikipedia.org/wiki/Random_sample_consensus|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. [[https://en.wikipedia.org/wiki/M-estimator|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.

/*
Opakujte ''num_iter''-krát:
  * Vyberte $3$ náhodné body z $X$ a proložte jimi kružnici.
  * Spočítejte čtverce vzdáleností všech bodů od této kružnice. Pokud je čtverec vzdáleností bodu od kružnice menší než ''threshold'', prohlašte ho za inlier.
  * Vezměte všechny inliery a proložte jimi kružnici funkcí ''fit_circle_nhom''. Zjistěte kolik bodů má čtverec vzdálenosti od této nové kružnice menší než ''threshold''.
Výstup funkce je ta kružnice, kterou jste nalezli v průběhu algoritmu ve třetím bodu a měla nejvíce inlierů.
*/


==== Templaty ====

Templaty, včetně skriptu pro testování si stáhněte zde: {{./matlab_template_circlefit.zip|template pro matlab}}, {{./python_template_circlefit.zip|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'':

{{circle-fit-nhom.svg?300}}
{{circle-fit-hom.svg?300}}
{{circle-fit-ransac.svg?300}}