Zagadnienia

7.1 Dopasowanie prostej
- 7.1.1 Uproszczenia i analiza problemu
- 7.1.2 Implementacja
7.2 Dopasowanie okręgu
7.3 Uwagi i uzupełnienia

7. Linie i okręgi

Podczas kursu matematyki obliczeniowej zetknęliśmy się z zadaniem: mając dany zestaw punktów na płaszczyźnie:

$P=\{ p_{i}=(x_{i},y_{i})\in\mathbb{R}^{2}:i=0,\ldots,N\},$

znaleźć prostą $L(x)=ax+b$ taką, że

$\sum _{{i=0}}^{N}|L(x_{i})-y_{i}|^{2}=\min!$

Było to po prostu inne sformułowanie liniowego zadania najmniejszych kwadratów na współczynniki $a,b$ . Zwróćmy jednak uwagę, że takie sformułowanie zadania nierówno traktuje zmienne $x_{i}$ oraz $y_{i}$ .

W wielu zastosowaniach, na przykład: w analizie obrazów, w energetyce, czy w komunikacji, mamy do czynienia z zadaniem minimalizacji średniej odległości prostej od zestawu $P$ punktów na płaszczyźnie. Rzeczywiście, jest to matematyczne sformułowanie praktycznego zadania w rodzaju: Jak poprowadzić (prostą) autostradę tak, by dojazd do niej z pobliskich miejscowości był możliwie szybki? Albo: Jak, na podstawie nieco rozmytego konturu budynku na zdjęciu, zdefiniować jego brzeg?

Przez analogię do zadania najmniejszych kwadratów określimy zadanie geometrycznego dopasowania prostej do punktów, polegające na wskazaniu takiej prostej $L$ , dla której

$\sum _{{i=1}}^{N}\operatorname{dist}(p_{i},L)^{2}=\min!$

przy czym odległość będziemy mierzyć w normie euklidesowej.

Powyższe zadanie jest najprostszą wersją całej klasy zadań geometrycznego dopasowania, w których krzywą $\gamma$ określonego typu chcemy dopasować do zadanego zestawu punktów $P$ tak, by

$\sum _{{i=1}}^{N}\operatorname{dist}(p_{i},\gamma)^{2}=\min!$

W przypadku budynków na przykład, szukalibyśmy prostokąta spełniającego powyższy warunek (por. uwagę na końcu wykładu). W niniejszym rozdziale zajmiemy się na początek, w oparciu o [6], geometrycznym dopasowaniem prostej, a w dalszej części wykładu, korzystając z materiałów zawartych w [1] — dopasowaniem okręgu do zadanych punktów.

Liczba punktów pomiarowych $N$ w praktyce nie jest bardzo duża, dlatego będziemy mieli do czynienia z zadaniami niezbyt intensywnymi obliczeniowo: głównym problemem będzie tutaj właściwe sformułowanie zadania w taki sposób, by można było je skutecznie zaatakować dostępnymi metodami.

7.1. Dopasowanie prostej

Nietrudno zauważyć, że reprezentowanie prostej $L$ w postaci $L(x)=ax+b$ nie jest zbyt zręczne, bo na przykład czasem najlepsze dopasowanie może zachodzić dla prostej $x=0$ , której nie jesteśmy w stanie tak przedstawić. Dlatego korzystniej będzie rozpatrzyć ogólne równanie prostej $L$ postaci

$Ax+By+C=0,$

oczywiście uzupełnione jakimś warunkiem normującym (współczynniki $A,B,C$ są dane z dokładnością do mnożnika). Tutaj najrozsądniej będzie przyjąć, że $A^{2}+B^{2}=1$ : wektor $[A,B]$ jest wektorem normalnym do prostej $L$ . Wtedy odległość punktu $p_{i}=(x_{i},y_{i})$ od $L$ jest dana równością

$\operatorname{dist}(p_{i},L)=\frac{|Ax_{i}+By_{i}+C|}{\sqrt{A^{2}+B^{2}}}=|Ax_{i}+By_{i}+C|,$

zatem mamy znaleźć $A,B,C$ takie, że $A^{2}+B^{2}=1$ oraz wyrażenie

$f(A,B,C)=\sum _{{i=1}}^{N}(Ax_{i}+By_{i}+C)^{2}$

przyjmuje najmniejszą wartość.

7.1.1. Uproszczenia i analiza problemu

Ponieważ na parametr $C$ nie ma ograniczeń, możemy łatwo go wyeliminować z rozważań przyrównując do zera pochodną

$\dfrac{\partial f}{\partial C}=2\sum _{{i=1}}^{N}(Ax_{i}+By_{i}+C)=0\qquad\implies\qquad C=-(A\bar{x}+B\bar{y}),$

gdzie $\bar{x}=\dfrac{1}{N}\sum _{{i=1}}^{N}x_{i}$ i analogicznie $\bar{y}=\dfrac{1}{N}\sum _{{i=1}}^{N}y_{i}$ .

W ten sposób zadanie zredukowało się do znalezienia wektora $u=(A,B)^{T}$ takiego, że

$\| Mu\| _{2}^{2}=\min!,\text{ przy czym }\| u\| _{2}^{2}=1,$

gdzie

$M=\begin{pmatrix}x_{1}-\bar{x}&y_{1}-\bar{y}\\ x_{2}-\bar{x}&y_{2}-\bar{y}\\ \vdots&\vdots\\ x_{N}-\bar{x}&y_{N}-\bar{y}\\ \end{pmatrix}.$

Jest to więc nic innego, jak zadanie wyznaczenia wektora szczególnego¹⁵Więcej na temat wektorów i wartości szczególnych można dowiedzieć się np. w wykładzie z Matematyki Obliczeniowej II, w rozdziale dotyczącym rozkładu SVD macierzy. macierzy $M$ , odpowiadającego najmniejszej wartości szczególnej. Rzeczywiście, jeśli $M=U\Sigma V^{T}$ i macierze $U,V$ są ortogonalne, a

$\Sigma=\begin{pmatrix}\sigma _{1}&0\\ 0&\sigma _{2}\\ 0&0\\ \vdots&\vdots\\ 0&0\end{pmatrix},$

przy czym dla ustalenia uwagi $\sigma _{1}\geq\sigma _{2}\geq 0$ , to

$\| Mx\| _{2}^{2}=\| U\Sigma V^{T}x\| _{2}^{2}=\|\Sigma V^{T}x\| _{2}^{2}=x^{T}V\begin{pmatrix}\sigma _{1}^{2}&0\\ 0&\sigma _{2}^{2}\end{pmatrix}V^{T}x\geq\sigma _{2}^{2}\| V^{T}x\| _{2}^{2}=\sigma _{2}^{2}\| x\| _{2}^{2}=\sigma _{2}^{2},$

a równość zachodzi gdy $V^{T}x=[0,1]^{T}$ , czyli dla $x=v_{2}$ .

7.1.2. Implementacja

Ponieważ w Octave znajduje się gotowa funkcja wyznaczająca pełny rozkład SVD zadanej macierzy, możemy pokusić się o implementację funkcji wyznaczenia $A,B,C$ na podstawie punktów, których współrzędne na osiach $OX$ i $OY$ podane są, odpowiednio, w talicach $x$ i $y$ :

function ABC = linefit(x,y)
x = x(:); y = y(:);
ABC = NaN(3,1);
xm = mean(x); ym = mean(y);
[U, S, V] = svd([x-xm, y-ym], 0); % ekonomiczny rozkład SVD
ABC(1:2) = V(:,end); % [A, B]
ABC(3) =  - ABC(1)*xm - ABC(2)*ym;
end

Aby przetestować działanie naszego kodu, zaburzymy losowo punkty prostej $y=5x+6$ i sprawdzimy, co wyjdzie. Porównamy także nasz (geometryczny) fit z dopasowaniem na podstawie standardowej metody najmniejszych kwadratów.

W celu wygodnego rysowania prostej zadanej w postaci uwikłanej, naprędce opracujemy funkcję lineplot, której listing zamieszczamy poniżej.

function lineplot(ABC, X, Y, color)
if nargin < 4
	color = '';
end
xmin = X(1); xmax = X(2);
ymin = Y(1); ymax = Y(2);
A = ABC(1); B = ABC(2); C = ABC(3);
if abs(A/B) <= 1
	plot(X, (-C-A*X)/B, ['-',color]);
else
	plot((-C-B*Y)/A, Y, ['-',color]);
end
end

Teraz możemy rozpocząć testy, korzystając na przykład z kodu poniższej postaci:

a = 5; b = 6;
x = linspace(1,3,300);
y = a*x + b;
x = x+0.6*(2*rand(size(x))-1);
y = y+0.6*(2*rand(size(y))-1);
% fit geometryczny
ABC = linefit(x,y);
% fit LZNK
ba = [ones(size(x)); x]'\y';
plot(x,y,'or');
hold on;
lineplot(ABC,[0,4],[min(y)-1,max(y)+1],'r');
lineplot([ba(2),-1,ba(1)],[0,4],[min(y)-1,max(y)+1],'k');
lineplot([a,-1,b],[0,4],[min(y)-1,max(y)+1],'g');
hold off;
legend('zaburzone dane','fit geometryczny','fit LZNK','oryginalna prosta');
grid on;

Dopasowanie w sensie geometrycznym jest optycznie znacznie lepsze!

Rys. 7.1. Przykładowy wynik dopasowania prostej w sensie LZNK i w sensie geometrycznym.

7.2. Dopasowanie okręgu

Analogicznie, można rozważyć zadanie geometrycznego dopasowania okręgu (lub, w większej ogólności, elipsy) $S$ do danych punktów na płaszczyźnie. Zastosowania mogą dotyczyć [1] medycyny (zwłaszcza — okulistyki), kontroli jakości wykonania obrotowych części mechanicznych, wyznaczania torów cząstek elementarnych, robotyki, a nawet — archeologii (wyznaczenie średnicy stadionu na podstawie zachowanych fragmentów). W wymienionych przypadkach zazwyczaj mamy do czynienia z zestawem punktów znajdujących się na łuku stanowiącym jedynie niewielki fragment okręgu. Z powodu zaburzeń danych, spowodowanych na przykład błędami pomiaru, zadane punkty nie układają się idealnie na poszukiwanej krzywej.

Jeśli okrąg $S$ reprezentować przez środek $O=(a,b)$ oraz promień $r$ :

$S=\{(x,y)\in\mathbb{R}^{2}:(x-a)^{2}+(y-b)^{2}=r^{2}\},$

to odległość $d_{i}$ zadanego punktu $p_{i}=(x_{i},y_{i})$ od $S$ jest równa, jak łatwo sprawdzić,

$d_{i}=|\| p_{i}-O\| _{2}-r|.$

W takim razie, musimy zminimalizować wartość funkcjonału

$f(a,b,r)=\sum _{{i=1}}^{N}\left(\sqrt{(x_{i}-a)^{2}+(y_{i}-b)^{2}}-r\right)^{2}.$

7.2.1. Atak na wprost

Ponieważ zadanie minimalizacji tym razem dotyczy bardziej skomplikowanego funkcjonału nieliniowego, użyjemy standardowej procedury optymalizacji nieliniowej w Octave, sqp.

function d = residgeom(X,x,y)
	% a = X(1); b = X(2); r = X(3);
	d = sumsq( sqrt((x-X(1)).^2 + (y-X(2)).^2) - X(3));
end

function [S,dist,info] = fitcircle1(S0,x,y)
% S = (a,b,r)
[S,dist,info] = sqp(S0,@(X)residgeom(X,x,y));
end

Funkcja residgeom wyznacza wartość funkcjonału $f(a,b,r)$ , przy czym przekazujemy jej wszystkie parametry funkcjonału w pierwszym argumencie, w postaci wektora X. Musimy tak zrobić, bo sqp spodziewa się funkcji jednego argumentu — taką tworzymy w postaci funkcji anonimowej

@(X)residgeom(X,x,y)

przekazywanej do sqp. Ponieważ w zadaniu nieliniowym wymagany jest sensowny punkt startowy $S_{0}$ — taki więc musimy przekazać, prócz zadanych współrzędnych punktów $(x_{i},y_{i})$ , do funkcji fitcircle1. Zadaniem użytkownika jest właściwy dobór $S_{0}=(a_{0},b_{0},r_{0})$ tak, by metoda iteracyjna była zbieżna do właściwego minimum.

7.2.2. Zmiana sformułowania zadania

Jakkolwiek dokonany przez nas wybór parametrów zadania: $a,b,r$ , jest niewątpliwie kuszący, to jednak, gdy punkty $(x_{i},y_{i})$ są prawie współliniowe, wyznaczany promień $r$ może być bardzo duży w porównaniu do współrzędnych środka $(a,b)$ . Ponadto małe zaburzenia położenia punktów $(x_{i},y_{i})$ mogą prowadzić do dużych skoków wielkości $r$ (a także — położenia środka okręgu), dlatego warto zadanie sformułować inaczej.

Okrąg $S$ możemy opisać za pomocą standardowego równania krzywej drugiego stopnia,

$A(x^{2}+y^{2})+Bx+Cy+D=0.$

Ponieważ wynika stąd, że

$\left(x+\dfrac{B}{2A}\right)^{2}+\left(y+\dfrac{C}{2A}\right)^{2}=\dfrac{B^{2}+C^{2}-4AD}{4A^{2}},$

Jak podaje [1], Pratt zasugerował warunek normalizacyjny

$B^{2}+C^{2}-4AD=1.$

Ponadto wówczas odległość $d_{i}$ punktu $(x_{i},y_{i})$ od okręgu jest równa

$d_{i}=2\dfrac{r_{i}}{1+\sqrt{1+4Ar_{i}}},\text{ gdzie }r_{i}=A(x_{i}^{2}+y_{i}^{2})+Bx_{i}+Cy_{i}+D,$

zatem musimy minimalizować funkcjonał

$f(A,B,C,D)=\sum _{i}d_{i}^{2}$

z ograniczeniem

$g(A,B,C,D)\equiv B^{2}+C^{2}-4AD-1=0.$

To zadanie po raz kolejny zrealizujemy za pomocą funkcji sqp:

function d = residgeom2(X,x,y)
	A = X(1); B = X(2); C = X(3); D = X(4);
	P = A*(x.^2 +y.^2) + B*x + C*y + D;
	d = sumsq(2*(P./(1+sqrt(1+4*A*P))));
end

function d = constr(X)
d = X(2)^2 + X(3)^2 - 4*X(1)*X(4) - 1;
end

function [S,dist,info] = fitcircle2(S0,x,y)
% S = (A,B,C,D)
c = zeros(4,1);
c(1) = 1/(2*S0(3)); c(2:3) = -2*c(1)*S0(1:2); c(4) = (c(2)^2 + c(3)^2 - 1)/(4*c(1));

[c,dist,info] = sqp(c, @(X)residgeom2(X,x(:),y(:)), @constr);
S = [-c(2:3)/(2*c(1)); sqrt(c(2)^2 + c(3)^2 - 4*c(1)*c(4))/(2*abs(c(1)))];
end

Jakkolwiek ostatnią linijkę kodu moglibyśmy teoretycznie zastąpić tańszym S = [-c(2:3); 1]/(2*c(1)), to jednak bezpieczniej będzie zrobić tak, jak powyżej — bo warunek $B^{2}+C^{2}-4AD=1$ może nie być idealnie spełniony.

7.2.3. Wybór przybliżenia początkowego jako nowe zadanie

Pozostaje jeszcze pytanie, w jaki sposób będziemy mogli sensownie wyznaczyć początkowe przybliżenie $S_{0}$ . Nadspodziewanie dobrym kandydatem na $S_{0}$ jest (patrz [1]) okrąg spełniający tzw. algebraiczny warunek dopasowania, tzn. minimalizujący, przy ograniczeniu $B^{2}+C^{2}-4AD=1$ , sumę kwadratów residuum równania okręgu:

$\mathcal{F}(A,B,C,D)=\sum _{i}r_{i}^{2},$

gdzie $r_{i}$ jest zadane przez $r_{i}=A(x_{i}^{2}+y_{i}^{2})+Bx_{i}+Cy_{i}+D$ .

Jest to więc — bardzo podobnie jak widzieliśmy to w przypadku fitowania prostej — zadanie najmniejszych kwadratów dla $u=(A,B,C,D)$

$\| Mu\| _{2}^{2}=\min!,\text{ z ograniczeniem }u^{T}Ku=1,$

gdzie

$M=\begin{pmatrix}x_{1}^{2}+y_{1}^{2}&x_{1}&y_{1}&1\\ \vdots&\vdots&\vdots&\vdots\\ x_{n}^{2}+y_{n}^{2}&x_{n}&y_{n}&1\end{pmatrix},\qquad K=\begin{pmatrix}&&&-2\\ &1&&\\ &&1&\\ -2&&&\end{pmatrix}_{{4\times 4}}.$

Niestety, macierz $K$ nie jest dodatnio określona, co komplikuje algorytm rozwiązywania tego zadania.

Funkcja Lagrange'a jest postaci

$L(u,\lambda)=u^{T}M^{T}Mu-\lambda(u^{T}Ku-1),$

zatem koniecznym warunkiem dla ekstremum jest

$L_{u}(u,\lambda)=0\text{ oraz }u^{T}Ku=1.$

Ponieważ $L_{u}(u,\lambda)=M^{T}Mu-\lambda Ku$ , to pierwszy warunek oznacza po prostu, że para $(\lambda,u)$ jest parą własną powyższego uogólnionego zadania własnego,

$\displaystyle M^{T}Mu-\lambda Ku=0.$

(7.1)

Wprowadzając (ekonomiczny) rozkład QR macierzy $M$ ,

$M=QR,$

mamy równoważnie $R^{T}Ru-\lambda Ku$ i stąd $(Ru,\lambda)$ musi być parą własną macierzy $RK^{{-1}}R^{T}$ . Ponieważ macierz $K^{{-1}}$ ma wartości własne $\{ 1,1,1/2,-1/2\}$ , to macierz $RK^{{-1}}R^{T}$ (na mocy tw. Sylvestera) także ma dokładnie jedną wartość własną ujemną, a pozostałe trzy — dodatnie.

Mnożąc stronami (7.1) przez $u^{T}$ dostajemy

$\lambda=\dfrac{u^{T}M^{T}Mu}{u^{T}Ku}=\dfrac{\| Mu\| _{2}^{2}}{u^{T}Ku}=\| Mu\| _{2}^{2},$

na mocy warunku $u^{T}Ku=1$ . Ponieważ $\| Mu\| _{2}^{2}\geq 0$ , interesująca nas $\lambda$ nie może być ujemna. Ze względu na to, że poszukiwane przez nas $u$ musi minimalizować $\| Mu\| _{2}^{2}$ , znaczy to, że minimum zostanie przyjęte równe najmniejszej dodatniej wartości własnej $B$ . W takim razie szukana przez nas para własna to para odpowiadająca najmniejszej dodatniej wartości własnej macierzy $RB^{{-1}}R^{T}$ .

Stąd jednym ze sposobów implementacji powyższego algorytmu wyznaczenia $u$ — a tym samym okręgu spełniającego wymieniony na początku tej sekcji algebraiczny warunek dopasowania — może być funkcja:

function [S,dist,info] = fitcircle3(x,y)
x = x(:); y = y(:);
invB = [0 0 0 -0.5; 0 1 0 0; 0 0 1 0; -0.5 0 0 0]; % $B^{-1}$ explicite 
X = [x.^2+y.^2, x, y, ones(size(x))];
[Q R] = qr(X,0);
[V L] = eig(R*invB*R'); 
[L i] = sort(diag(L));
c = R \ V(:,i(2)); % interesuje nas druga najmniejsza w.wl, a raczej: wektor
S = [-c(2:3)/(2*c(1)); sqrt(c(2)^2 + c(3)^2 - 4*c(1)*c(4))/(2*abs(c(1)))];
dist = residgeom(S,x,y);
info = 1;
end

Ćwiczenie 7.1

Przeredaguj funkcje fitcircle1 i fitcircle2 tak, by jeśli użytkownik nie poda przybliżenia startowego, zostało wybrane przybliżenie wyznaczanie funkcją fitcircle3.

Rozwiązanie:

Musimy zmienić kolejność argumentów tak, by przybliżenie początkowe było ostatnim. Wtedy początek funkcji, np. fitcircle1, należałoby zapisać np. w takiej formie:

function [S,dist,info] = fitcircle1(x,y,S0)
if nargin < 3
	S0 = fitcircle3(x,y);
end
% ... tu dalsze instrukcje ...
end

Poniżej przytaczamy przykładowy skrypt dopasowujący okrąg do zadanych punktów trzema opisanymi wcześniej metodami.

t = linspace(0,0.2*pi,100)';
X = 3*sin(t)+0.2*rand(size(t));
Y = 3.35*cos(t)+0.1*rand(size(t));

XY = [X,Y];
e = fitcircle3(X,Y);
c = fitcircle1(e,X,Y);
d = fitcircle2(e,X,Y);
[e c d]

t = linspace(0,2*pi,300);
plot(X,Y,'ro');
axis('square');
hold on;
plot(c(1) + c(3)*sin(t),c(2)+c(3)*cos(t),'g-');
plot(d(1) + d(3)*sin(t),d(2)+d(3)*cos(t),'r-');
plot(e(1) + e(3)*sin(t),e(2)+e(3)*cos(t),'k-');
plot(X,Y,'ro');
hold off;
legend('dane', 'fc1', 'fc2', 'fc3');
[residgeom(c,X,Y), residgeom(d,X,Y), residgeom(e,X,Y)]

Wszystkie trzy okręgi praktycznie są nie do odróżnienia

Rys. 7.2. Przykładowy wynik dopasowania okręgu trzema opisywanymi metodami, gdy zadane punkty prawie leżą na pewnym okręgu. W tym (typowym) przypadku, każda z metod daje praktycznie identyczny wynik!.

7.3. Uwagi i uzupełnienia

W przypadku, gdy wiadomo, które punkty $P$ należą do jakich boków prostokąta, zadanie geometrycznego dopasowania prostokąta jest w miarę prostym uogólnieniem zadania dopasowania linii, zob [6, Chapter 6, str. 88]. [1] jest prawdziwą kopalnią wiedzy na temat dopasowywania prostych, okręgów i elips.

Treść automatycznie generowana z plików źródłowych LaTeXa za pomocą oprogramowania wykorzystującego LaTeXML.

Obliczenia naukowe