neděle 30. prosince 2012

Převod quaternionu na Eulerovy úhly

V tomto příspěvku odvodím vztahy pro převod quaternionu na Elerovy úhly používané v letectví, tedy úhly, které postupně reprezentují rotaci okolo osy $z$ , $y$ a $x$. Tyto úhly značím postupně jako $\alpha$, $\beta$ a $\gamma$. Prvek matice $\mathbf{R}$ na řádku $i$ a ve sloupci $j$ značím jako $\mathbf{R}_{ij}$.

Nejprve převedeme quaternion na rotační matici a z té posléze získáme Eulerovy úhly. Rotační matici z quaternionu $\mathbf{q}=(q_w, q_x, q_y, q_z)$ získáme pomocí následujícího vztahu [1]:
$$
\mathbf{R}_q =

\begin{bmatrix}
    1 - 2 q_y^2 - 2 q_z^2 & 2 q_x q_y - 2 q_z q_w & 2 q_x q_z + 2 q_y q_w \\
    2 q_x q_y + 2 q_z q_w & 1 - 2 q_x^2 - 2 q_z^2 & 2 q_y q_z - 2 q_x q_w \\
    2 q_x q_z - 2 q_y q_w & 2 q_y q_z + 2 q_x q_w & 1 - 2 q_x^2 - 2 q_y^2
\end{bmatrix}
$$ Rotační matice pro Eulerovy úhly definované v úvodním odstavci při značení $\sin\!\left(\alpha\right)=s_\alpha$ (pro kosinus a další úhly analogicky) má následující podobu
$$ \mathbf{R}_E =
\begin{bmatrix} c_\alpha\, c_\beta & c_\alpha\, s_\beta\, s_\gamma - c_\gamma\, s_\alpha & s_\alpha\, s_\gamma + c_\alpha\, c_\gamma\, s_\beta\\ c_\beta\, s_\alpha & c_\alpha\, c_\gamma + s_\alpha\, s_\beta\, s_\gamma & c_\gamma\, s_\alpha\, s_\beta - c_\alpha\, s_\gamma\\ - s_\beta & c_\beta\, s_\gamma & c_\beta\, c_\gamma \end{bmatrix}
$$ Z matice $\mathbf{R}_E$ určíme nejprve úhel $\beta$ jako
$$ \beta = \mathrm{asin}\,\left(-\mathbf{R}_{31} \right).
$$ Dále se postup v závislosti na hodnotě prvku $\mathbf{R}_{31}$ dělí na tři případy.
  1. Pokud platí, že $\mathbf{R}_{31} = 1$, pak má rotační matice nutně následující podobu
$$
\mathbf{R}_E' =
\begin{bmatrix} 0 & -c_\alpha\, s_\gamma - c_\gamma\, s_\alpha & s_\alpha\, s_\gamma - c_\alpha\, c_\gamma\\
0 & c_\alpha\, c_\gamma - s_\alpha\,  s_\gamma & -c_\gamma\, s_\alpha\, - c_\alpha\, s_\gamma\\
1 & 0 & 0 \end{bmatrix}
=
\begin{bmatrix} 0 & -s_{\gamma + \alpha} & c_{\gamma + \alpha}\\
0 & c_{\gamma + \alpha} & -s_{\gamma + \alpha}\\
1 & 0 & 0 \end{bmatrix}.
$$ Řešení  pro úhly $\alpha$ a $\gamma$ je nekonečně mnoho. Je možné získat pouze součet těchto úhlů. Z nekonečně mnoha řešení tedy vybereme jedno. \begin{align*}
\alpha &= 0\\
\gamma &= \mathrm{atan2}\!\left( \mathrm{-\mathbf{R}_{12}}\mathrm{\mathbf{R}_{22}} \right)
\end{align*}
  1. Pokud platí, že $\mathbf{R}_{31} = -1$, pak má rotační matice nutně následující podobu
$$
\mathbf{R}_E' =
\begin{bmatrix} 0 & c_\alpha\, s_\gamma - c_\gamma\, s_\alpha & s_\alpha\, s_\gamma + c_\alpha\, c_\gamma\\
0 & c_\alpha\, c_\gamma + s_\alpha\,  s_\gamma & c_\gamma\, s_\alpha\, - c_\alpha\, s_\gamma\\
-1 & 0 & 0 \end{bmatrix}
=
\begin{bmatrix} 0 & s_{\gamma - \alpha} & c_{\gamma - \alpha}\\
0 & c_{\gamma - \alpha} & -s_{\gamma - \alpha}\\
-1 & 0 & 0 \end{bmatrix}.
$$ Obdobně jako v předchozím případě je řešení  pro úhly $\alpha$ a $\gamma$ nekonečně mnoho. Opět vybereme jedno řešení, které má například následující podobu \begin{align*}
\alpha &= 0,\\
\gamma &= \mathrm{atan2}\!\left( \mathrm{\mathbf{R}_{12}}, \mathrm{\mathbf{R}_{22}} \right).
\end{align*}
  1. Pro jiné hodnoty $\mathbf{R}_{31}$ určíme zbývající úhly podle následujících vztahů
\begin{align*} \alpha &= \mathrm{atan2}\!\left(\mathrm{\mathbf{R}_{21}}, \mathrm{\mathbf{R}_{11}} \right),\\
\gamma &= \mathrm{atan2}\!\left( \mathrm{\mathbf{R}_{32}}, \mathrm{\mathbf{R}_{33}} \right).
\end{align*} Rotační matici s konkrétními hodnotami jsme získali převodem z quaternionu a známe vztahy pro výpočet jednotlivých Eulerových úhlů z této matice, které jsme odvodili z obecného tvaru rotační matice pro tyto úhly. Již tedy zbývá jen do těchto vztahů dosadit. Pro jednoduchost zde uvádím jen dosazené vztahy pro poslední případ.\begin{align*} \alpha &= \mathrm{atan2}\!\left( 2 q_x q_y + 2 q_z q_w, 1 - 2 q_x^2 - 2 q_y^2 \right),\\
\beta &= \mathrm{asin}\,\left(2 q_y q_w - 2 q_x q_z\right),\\
\gamma &= \mathrm{atan2}\!\left( 2 q_y q_z + 2 q_x q_w, 1 - 2 q_y^2 - 2 q_z^2 \right).
\end{align*} Převod quaternionů na Eulerovy úhly implementovaný v programovacím jazyce C# vypadá následovně:
        public void ToEulerAngles(out double yaw, out double pitch, out double roll)
        {
            double test = q[1] * q[3] - q[0] * q[2];

            if (test > 0.49999)
            {
                roll = Math.Atan2(2 * q[1] * q[0] - 2 * q[2] * q[3], 2 * q[1] * q[3] + 2 * q[2] * q[0]);
                pitch = -Math.PI / 2;
                yaw = 0;
                return;
            }

            if (test < -0.49999)
            {
                roll = Math.Atan2(2 * q[2] * q[3] - 2 * q[1] * q[0], 2 * q[1] * q[3] + 2 * q[2] * q[0]);
                pitch = Math.PI / 2;
                yaw = 0;
                return;
            }

            double sqx = q[1] * q[1];
            double sqy = q[2] * q[2];
            double sqz = q[3] * q[3];

            roll = Math.Atan2(2 * q[2] * q[3] + 2 * q[0] * q[1], 1 - 2 * sqx - 2 * sqy);
            pitch = Math.Asin(-2 * test);
            yaw = Math.Atan2(2 * q[1] * q[2] + 2 * q[3] * q[0], 1 - 2 * sqy - 2 * sqz);
        }

Zdroje

Individuální projekt (2)

Vnější souřadnicová soustava

Vnější souřadnicová soustava nosiče vůči Zemi má počátek v hmotném středu nosiče, osa $x$ směřuje na sever, osa $z$ do středu Země a osa $y$ je orientována tak, aby byl souřadnicový systém pravotočivý.

Vnitřní souřadnicová soustava nosiče

Souřadnicová soustava nosiče je znázorněna na Obrázku 1. Jedná se o v letectví běžně používaný souřadnicový systém, kde počátek je ve středu nosiče (je tedy souhlasný s počátkem vnější souřad. soustavy), osa $x$ je orientována ve směru letu, osa $y$ prochází osou křídel a osa $z$ je orientována tak, aby směřovala směrem dolů a systém splňoval pravidlo pravotočivé souřadnicové soustavy.

Obrázek 1: Souřadnicová soustava nosiče


Poloha a orientace nosiče

Polohu nosiče získáme přímo z polohy určené GPS přijímačem v inerciální jednotce. Tato poloha sice neodpovídá poloze nosiče, jelikož je dána polohou antény, ale vzhledem k tomu, že rozdíl v těchto polohách je minimální vzhledem k ostatním uvažovaným vzdálenostem ve výpočtu přímé kinematické úlohy, lze tento rozdíl zanedbat, případně by byl kompenzován přímo v inerciální jednotce.

Orientaci nosiče získáme opět přímo z inerciální jednotky. Inerciální jednotka může být umístěna v rámci nosného objektu tak, že její orientace se neshoduje s orientací nosného objektu. Tato situace je však kompenzována přímo v inerciální jednotce. Orientace je v inerciální jednotce uchovávána v podobě quaternionů.

Popis samotné inerciálně stabilizované platformy

Samotná inerciálně stabilizovaná platforma má čtyři klouby. Dva zajišťují rotaci platformy v azimutální ose a zbylé dva v elevační ose. Pro rotaci okolo jedné osy jsou tedy vždy použity dva klouby. Jeden je schopen rotace o $n\cdot360^\circ$ a je vůči druhému relativně pomalý, druhý je schopen rotace jen ve velmi omezeném rozsahu, ale je velmi rychlý. Pro samotnou stabilizaci se tedy s výhodou využívají převážně ony rychlé klouby.

Pro naše účely je ale výhodnější reprezentovat kamerovou platformu jako otevřený kinematický řetězec s pouze dvěma klouby. Pro přímou kinematickou úlohu pouze přičteme k pomalému kloubu pro danou osu otočení rychlého kloubu. Pro inverzní kinematickou úlohu budeme uvažovat pouze pomalé klouby. Takovýto kinematický řetězec můžeme popsat například Denavit-Hartenbergovou notací. Pro náš případ je také výhodné, že v této notaci se bude shodovat osa kamery a dálkoměru s osou $x$ souřadnicové soustavy posledního kloubu.

Parametry Denavitovy-Hartenbergovy notace
$n$ $\theta_n$ $d_n$ $a_n$ $\alpha_n$
1 $q_1$ $d_1$ $0$ $\pi/2$
2 $q_2$ $0$ $0$ $0$
Parametry $q_1$ a $q_2$ jsou kloubové souřadnice robota, které odpovídají natočení azimutálního a elevačního kloubu. Tato natočení jsou vyčítána řídící jednotkou kamerové platformy z inkrementálních rotačních čítačů. Ty však bohužel nemají shodnou orientaci s kloubovými souřadnicemi. Na Obrázku 2 jsou znázorněny orientace jednotlivých čítačů. Z nich snadno získáme převodní vztah do kloubových souřadnic.

\begin{align*}
q_1 &= -(\beta_0 + \beta_1) + \pi \\
q_2 &= \alpha_0 + \alpha_1 - \pi/2
\end{align*}
Obrázek 2: Orientace kloubových úhlů


Individuální projekt (1)

Tato série příspěvků se bude věnovat úloze, kterou jsem řešil v rámci předmětu Individuální projekt (A3B35IND). Úloha se týká projektu inerciálně stabilizované kamerové platformy, tedy projektu, jehož se okrajově týkaly i předchozí příspěvky o inerciální jednotce. Podrobnější popis projektu inerciálně stabilizované kamerové platformy je na stránkách skupiny AA4CC.

Cílem práce je vyřešit přímou a inverzní kinematickou úlohu pro kamerovou hlavici na nějakém nosném objektu. Přímá část práce spočívá v určení místa na mapě, jenž je ve středu zorného pole kamery. Motivací pro tuto úlohu je fakt, že operátor kamerové hlavice má aktuálně k dispozici aplikaci, v níž je vykreslována pozice nosného objektu a obraz z kamery, bohužel však již nemá k dispozici informaci o tom, na jaké místo na mapě je kamerová hlavice zaměřena. Inverzní úloha spočívá v určení kloubových souřadnic kamerové platformy takových, aby zadané místo na mapě bylo v zorném poli kamery. Zde je motivace zřejmá. Místo nepohodlného ručního zaměřování určitého místa na mapě přenechat tuto úlohu řídícímu systému a celý tento proces tak zautomatizovat.

sobota 29. prosince 2012

C# třída pro práci s (transformačními) maticemi

Jelikož programovací jazyk C# nenabízí ve standardních třídách vhodnou třídu pro práci s maticemi, natož pak s transformačními maticemi, musel jsem si vytvořit vlastní třídu. Třída obsahuje mimo základní operace s maticemi, jako je sčítání, maticové násobení, násobení skalárem a transpozici i operace spojené s transformačními maticemi. Jmenovitě se jedná o generování rotačních matic okolo osy $x$, $y$ a $z$, rotační matice podle Eulerových úhlů v notaci $ZYX$, translační matice a transformační matice podle parametrů z Denavit-Hartnebergovy notace. Rotační matice jsou dimenze $4\times 4$, tedy v homogenních souřadnicích. Třída přetěžuje operátory sčítání, násobení a indexového přístupu pro maximálně intuitivní práci s jednotlivými instancemi třídy. Její aktuální verze s podrobnější dokumentací je přístupná na mém repositáři na githubu.

Seznam funkcí implementovaných ve třídě:
        public Matrix(int rows, int cols)
        public Matrix Copy()
        public static Matrix ZeroMatrix(int rows, int cols)
        public static Matrix IdentMatrix(int n)
        public override string ToString()
        public static Matrix Transpose(Matrix m)
        public static Matrix Multiply(Matrix m1, Matrix m2)
        private static Matrix ScalarMultiply(double n, Matrix m)
        private static Matrix Add(Matrix m1, Matrix m2)
        public static Matrix EulerRotaionMatrix(double yaw, double pitch, double roll)
        public static Matrix RotaionMatrixZ4(double angle)
        public static Matrix RotaionMatrixX4(double angle)
        public static Matrix RotaionMatrixY4(double angle)
        public static Matrix TranslationMatrix(double xTrans, double yTrans, double zTrans)
        public static Matrix DHTransMatrix(double[][] DHparams)