Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 63 Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 65 Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 67 Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 69 Notice: Undefined variable: base in /home/misc/mst/public_html/lecture.php on line 36 Obliczenia naukowe – 7. Linie i okręgi – MIM UW

Zagadnienia

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ólnego15Wię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.

Projekt współfinansowany przez Unię Europejską w ramach Europejskiego Funduszu Społecznego.

Projekt współfinansowany przez Ministerstwo Nauki i Szkolnictwa Wyższego i przez Uniwersytet Warszawski.