Przekształcenia wymierne – testy

Niestety , metoda Newtona nie znajduje przekształcenia dla rzeczywistego obrazuPerspektywa
Dla pięciu punktów nie jest zbieżna z powodu łuku pomiędzy dwoma górnymi lub dolnymi punktami. Dla czterech i odpowiedniego wyboru punktów, które mają leżeć na linii daje się znaleźć transformację, tylko w odwrotną stronę niż potrzeba.
Ciekawa własność dla 4 punktów + punkt na linii : mianowniki obu wyrażeń wymiernych są takie same. Można to wykorzystać dla uproszczenia – 4 punkty, żadnych dodatkowych warunków na to, by punkt leżał na linii, a jedynie równość mianowników wyrażeń. Przez co będziemy mieli układ równań rozmiaru 8 a nie 10.
$x’ = (ax+by+c)/(1+dx+ey)$
$y’ = (fx+gy+h)/(1+dx+ey)$
Nazywając zmienne podobnie jak będą wykorzystywane w programie:
$Bx = (x0*Ax+x1*Ay+x2)/(1+x3*Ax+x4*Ay)$
$By = (x5*Ax+x6*Ay+x7)/(1+x3*Ax+x4*Ay)$
Pochodne cząstkowe $(x0*Ax+x1*Ay+x2)/(1+x3*Ax+x4*Ay)-Bx$
Ax/(1+Ax x3+Ay x4)
Ay/(1+Ax x3+Ay x4)
1/(1+Ax x3+Ay x4)
-((Ax (Ax x0+Ay x1+x2))/(1+Ax x3+Ay x4)^2)
-((Ay (Ax x0+Ay x1+x2))/(1+Ax x3+Ay x4)^2)
0
0
0
Pochodne czątkowe $(x5*Ax+x6*Ay+x7)/(1+x3*Ax+x4*Ay)-By$
0
0
0
-((Ax (Ax x5+Ay x6+x7))/(1+Ax x3+Ay x4)^2)
-((Ay (Ax x5+Ay x6+x7))/(1+Ax x3+Ay x4)^2)
Ax/(1+Ax x3+Ay x4)
Ay/(1+Ax x3+Ay x4)
1/(1+Ax x3+Ay x4)

Wynik:
Transformacja
Perspektywa poprawiona, jednak zakrzywienie jeszcze bardziej widoczne – y powinno zależeć od x w kwadracie , dodamy czynniki x^2 i x*y dla y; natomiast x pozostawiamy bez zmian.

$(x0*Ax+x1*Ay+x2)/(1+x3*Ax+x4*Ay) – Bx$
$(x5*Ax+x6*Ay+x7 +Ax^2*x8 +Ax*Ay*x9)/(1+x3*Ax+x4*Ay) – By$
Pochodne cząstkowe pierwszej, tak jak było:
Ax/(1+Ax x3+Ay x4)
Ay/(1+Ax x3+Ay x4)
1/(1+Ax x3+Ay x4)
-((Ax (Ax x0+Ay x1+x2))/(1+Ax x3+Ay x4)^2)
-((Ay (Ax x0+Ay x1+x2))/(1+Ax x3+Ay x4)^2)
0
0
0
0
0
Pochodne czątkowe drugiej:
0
0
0
-((Ax (Ax x5+Ay x6+x7+Ax^2 x8+Ax Ay x9))/(1+Ax x3+Ay x4)^2)
-((Ay (Ax x5+Ay x6+x7+Ax^2 x8+Ax Ay x9))/(1+Ax x3+Ay x4)^2)
Ax/(1+Ax x3+Ay x4)
Ay/(1+Ax x3+Ay x4)
1/(1+Ax x3+Ay x4)
Ax^2/(1+Ax x3+Ay x4)
(Ax Ay)/(1+Ax x3+Ay x4)
Wynik:
Fragment

Widać dosyć pogięte linie, transformacja nie uwzględnia że wygięcie postępowało w stronę góry, tak ze po odtworzeniu oryginału na samej górze jeszcze widać to wygięcie a na dole jest w przeciwną stronę.

Wniosek – trudności występują z powodu zagięcia na walcowatej powierzchni, gdyby nie to, mogło by zadziałać przekształcenie 4-punktowe.

Znajdowanie przekształcenia wymiernego

Dwuwymiarowa macierz obrotu ma postać

$ \begin{bmatrix}
\cos \theta & -\sin \theta \\
\sin \theta & \cos \theta \\
\end{bmatrix}$

Jednak sam obrót wokół punktu (0,0) bez przesunięcia (translacji) na wiele się nie przyda. Do przekształceń $\mathbb{R}^2 \Rightarrow \mathbb{R}^2$ potrzebna jest macierz przekształceń afinicznych rozmiaru 3×3.
$ \begin{bmatrix}
\cos \theta & -\sin \theta & 0 \\
\sin \theta & \cos \theta & 0\\
0 & 0 & 1 \\
\end{bmatrix}$

i macierz translacji
$ \begin{bmatrix}
1 & 0 & dx \\
0 & 1 & dy\\
0 & 0 & 1 \\
\end{bmatrix}$

macierz skalowania
$ \begin{bmatrix}
S_x & 0 & 0 \\
0 & S_y & 0\\
0 & 0 & 1 \\
\end{bmatrix}$

I inne, np. skrzywienie (skew,shearing)
$ \begin{bmatrix}
1 & F_x & 0 \\
F_y & 1 & 0\\
0 & 0 & 1 \\
\end{bmatrix}$

Macierze mnożymy ze sobą i na koniec mnożymy przez wektor współrzędnych

$$\begin{bmatrix} x’ \\ y’ \\1 \end{bmatrix} \text{=} \begin{bmatrix} 1 & 0 & d_1x \\ 0 & 1 & d_1y\\ 0 & 0 & 1 \\ \end{bmatrix} \begin{bmatrix} \cos \theta & -\sin \theta & 0 \\ \sin \theta & \cos \theta & 0\\ 0 & 0 & 1 \\ \end{bmatrix} \begin{bmatrix} 1 & 0 & d_2x \\ 0 & 1 & d_2y\\ 0 & 0 & 1 \\ \end{bmatrix}\begin{bmatrix} x \\ y \\1 \end{bmatrix}$$

Ogólnie, dowolne przekształcenie afiniczne to 6 liczb rzeczywistych bo nie liczy się ostatni wiersz macierzy.

Ale czasem takie przekształcenie nie wystarczy. Mamy powierzchnię obróconą w przestrzeni i rzutowaną na inną powierzchnię. Na przykład zdjęcie strony książki czy zdjęcie lotnicze/ z drona.

Perspektywa – projekcja ortograficzna

Oprócz perspektywy należy uwzględnić wszelkie przekształcenia linowe, którym mogło ulec zdjęcie w przestrzeni i afiniczne na płaszczyźnie na którą rzutujemy. W ten sposób dla przekształconego punktu mamy 10 współczynników:
$x’ = \frac{ax+by+c}{dx+ey+1}$
$y’ = \frac{fx+gy+h}{ix+jy+1}$

W mianowniku dla wyrazu wolnego wystarczy jedynka, ponieważ wyrażenia wymierne można zawsze skrócić przy założeniu że wyraz wolne będzie różny od zera. W istocie dla faktycznych przekształceń liczby d,e,i oraz j są niewielkie w porównaniu z jedynką.

Aby rozwiązać potrzebujemy układu dziesięciu równań nieliniowych które można uzyskać z przekształcenia pięciu punktów.

Użyjemy do tego wielowymiarowej metody Newtona. Metoda ta nie zawsze jest zbieżna, ale dla sensownej transformacji i punktu startowego – braku transformacji znajduje rozwiązanie.

Metoda ta jest w Rosetta w C#

            do
            {
                for (int i = 0; i < size; i++)
                    F[i] = fun.F(i, X);
                for (int i = 0; i < size; i++)
                    for (int j = 0; j < size; j++)
                        J[i, j] = fun.df(i, j, X);
                J.ElimPartial(F);
                X -= F;
                //need weight vector because different coordinates can diffs by order of magnitudes
            } while (F.norm(fun.weights()) > 1e-12);

Do zakończenia używam wektota wag, bo niektóre współrzędne mają przyjmują o wiele większe wartości niż inne.
Nie trzeba odwracać macierzy, wystarczy eliminacja Gaussa, którą mamy w Rosetta eliminacja w C#

        //with partial pivot
        internal void ElimPartial(Vector B)
        {
            for (int diag = 0; diag < rows; diag++)
            {
                int max_row = diag;
                double max_val = Math.Abs(this[diag, diag]);
                double d;
                for (int row = diag + 1; row < rows; row++)
                    if ((d = Math.Abs(this[row, diag])) > max_val)
                    {
                        max_row = row;
                        max_val = d;
                    }
                SwapRows(diag, max_row);
                B.SwapRows(diag, max_row);
                double invd = 1 / this[diag, diag];
                for (int col = diag; col < cols; col++)
                    this[diag, col] *= invd;
                B[diag] *= invd;
                for (int row = 0; row < rows; row++)
                {
                    d = this[row, diag];
                    if (row != diag)
                    {
                        for (int col = diag; col < cols; col++)
                            this[row, col] -= d * this[diag, col];
                        B[row] -= d * B[diag];
                    }
                }
            }
        }

Skoro
$x’ = \frac{ax+by+c}{dx+ey+1}$
$y’ = \frac{fx+gy+h}{ix+jy+1}$
To (x,y) należą do punktów $A_i$ a (x’,y’) do $B_i$. Natomiast a,c,d,e,f,h,h,i,j są szukane i nazwijmy je $x_k$.

Szukamy zera funkcji F, dla współrzędnej x:
$(x_0 * A_x + x_1 * A_y + x_2) / (1 + x_3 * A_x + x_4 * A_y) – B_x = 0$
Dla współrzędnej y:
$(x_5 * A_x + x_6 * A_y + x_7) / (1 + x_8 * A_x + x_9 * A_y) – B_y = 0$
Należy teraz policzyć pochodne cząstkowe, w czym może pomóc Mathematica. Kod jest na GitHubie
Weżmy transformację
$x’ = \frac{1.2x + 0.2y – 50}{1 + 0.0001x + 0.0002y}$
$y’ = \frac{0.3x + 1.1y – 200}{1 – 0.0005x + 0.00025y}$
p.x = () / ();

Dla punktów (0, 0), (0, 200), (200, 200), (200, 0), (100, 100)
przekształcenie daje
(-50, -200),
(-9.61538461538462, 19.0476190476191),
(216.981132075472, 84.2105263157895),
(186.274509803922, -155.555555555556),
(87.378640776699, -61.5384615384615)
Metoda Newtona znajduje to przekształcenie.
Jednak czasem wskazanie pięciu punktów może być niewygodne, zamieńmy warunek na cztery punkty + 2 punkty leżące w środku dwóch ramion czworokąta początkowego muszą lezeć na linii. Z Wikipedii mamy warunek konieczny i wystarczający współliniowości:
$\frac{x_3 – x_1}{x_2-x_1}= \frac{y_3 – y_1}{y_2 – y_1}$
Czyli
fx = (x0*Cx + x1*Cy + x2)/(1 + x3*Cx + x4*Cy)
fy = (x5*Cx + x6*Cy + x7)/(1 + x8*Cx + x9*Cy)
i (fx – Bx1)/(Bx2 – Bx1) = (fy – By1)/(By2 – By1) co oznacza że ma być równa zeru funkcja
(fx – Bx1)/(Bx2 – Bx1) – (fy – By1)/(By2 – By1)

$\frac{\frac{Cx \cdot x0+Cy \cdot x1+x2}{Cx \cdot x3+Cy \cdot x4+1}-Bx1}{Bx2-Bx1}-\frac{\frac{Cx \cdot x5+Cy \cdot x6+x7}{Cx \cdot x8+Cy \cdot x9+1}-By1}{By2-By1}$

gdzie $C_x = (A_1x+A_2x)/2$ i analogicznie C_y.
Teraz policzymy pochodne czątkowe po $x_0$ .. $x_9$:
$\frac{Cx}{(Bx2-Bx1) (Cx \cdot x3+Cy \cdot x4+1)}$
$\frac{Cy}{(Bx2-Bx1) (Cx \cdot x3+Cy \cdot x4+1)}$
$\frac{1}{(Bx2-Bx1) (Cx \cdot x3+Cy \cdot x4+1)}$
$-\frac{Cx \cdot (Cx \cdot x0+Cy \cdot x1+x2)}{(Bx2-Bx1) (Cx \cdot x3+Cy \cdot x4+1)^2}$
$-\frac{Cy \cdot (Cx \cdot x0+Cy \cdot x1+x2)}{(Bx2-Bx1) (Cx \cdot x3+Cy \cdot x4+1)^2}$
$-\frac{Cx}{(By2-By1) (Cx \cdot x8+Cy \cdot x9+1)}$
$-\frac{Cy}{(By2-By1) (Cx \cdot x8+Cy \cdot x9+1)}$
$-\frac{1}{(By2-By1) (Cx \cdot x8+Cy \cdot x9+1)}$
$\frac{Cx \cdot (Cx \cdot x5+Cy \cdot x6+x7)}{(By2-By1) (Cx \cdot x8+Cy \cdot x9+1)^2}$
$\frac{Cy \cdot (Cx \cdot x5+Cy \cdot x6+x7)}{(By2-By1) (Cx \cdot x8+Cy \cdot x9+1)^2}$
Wyniki:

NazwaWartość
x0
1.08658974903739
x10.199206687385285
x2-50
x3-0.000508833978851917
x40.000282504511930326
x50.301374174488111
x61.10061913356058
x7-200
x8-0.000508833978852153
x90.000282504511930329

Zauważamy że (z dokładnościa do tolerancji metody Newtona) $x_3 = x_8$ a $x_4 = x_9$ czyli wielomiany w mianowniku są równe.