diff --git a/README.md b/README.md
index a13a1d60a4c257e5b442fb6a3c851927a2435386..a0a57e543031943db63ccd69f0427d33996a5121 100644
--- a/README.md
+++ b/README.md
@@ -66,6 +66,18 @@ python setup.py install
```
## Blender-python version
+Navipy can be interfaced with blender. It is highly recommended to use the same version of packages of blender when doing so, in order to reduce problem of compatibility.
+To determine the packages that you will need, you can run the script: `navipy/script/check_blender_versions_pip.py` in blender or via commandline:
+
+```
+blender -b -P check_blender_versions_pip.py
+```
+
+It will create a textfile containing all packages used by blender. They can be installed in your virtualenvironment (prior to navipy) by doing:
+```
+pip install -r requirement.txt
+```
+
| Blender version | Python version |
| --------------- | -------------- |
| 2.79b | 3.5.3 |
diff --git a/doc/source/tutorials/02-recording-animal-trajectory.ipynb b/doc/source/tutorials/02-recording-animal-trajectory.ipynb
index 557203069aa8a3301d4584f28f772fcc2fa59692..2b74d09a62c3f7ca2eac4a721a1308c764078877 100644
--- a/doc/source/tutorials/02-recording-animal-trajectory.ipynb
+++ b/doc/source/tutorials/02-recording-animal-trajectory.ipynb
@@ -6,23 +6,125 @@
"source": [
"# Recording animal trajectory\n",
"\n",
- "## Camera calibration\n",
+ "## Camera recording\n",
"\n",
- "A camera is an optical intrument for aquiring images. Cameras are composed of a sensors (converting photon to electrical signal) and a lens (foccussing light rays on the sensors). \\\n",
- "Many experimental paradigm require the recording of animal motion with one or more camera. \\\n",
- "The experimenter is more interested by the position of the animal in his or her arena (world \\\n",
+ "A camera is an optical intrument for aquiring images. Cameras are composed of a sensors (converting photon to electrical signal) and a lens (foccussing light rays on the sensors). \n",
+ "Many experimental paradigm require the recording of animal motion with one or more camera. \n",
+ "The experimenter is more interested by the position of the animal in his or her arena (world \n",
"coordinate system) than by the position of the animal within the image (camera-coordinate system). Therefore, we need to transform the position of the animal in the image to its position in the world. To do so, we need to calibrate our camera. This calibration require two steps:\n",
"\n",
"* remove distortion from the lens (intrinsic calibration)\n",
"* determine the position and orientation (i.e. pose) of the camera in the environment (extrinsic calibration)\n",
"\n",
- "### Intrinsic calibration\n",
+ "### Direct Linear Transformation (DLT)\n",
+ "\n",
+ "When we observe an animal through a camera, we only observe the projection of the animal. For example the center of mass of the animal in three dimension is projected in two dimension on the camera, i.e. in the camera plane, and only these two dimensions are visible from the camera. Mathematically, the point in 3D space is transformed into the camera space. This transformation can be described as follow:\n",
+ "\n",
+ "$$\n",
+ "\\begin{align}\n",
+ "\\begin{bmatrix} u-u_0 \\\\ v-v_0 \\\\ -d \\end{bmatrix} &= T \\begin{bmatrix} x-x_0 \\\\ y-y_0 \\\\ z-z_0 \\end{bmatrix} \\\\\n",
+ "&=c \\begin{bmatrix} r_{11} & r_{12} & r_{13} \\\\ r_{21} & r_{22} & r_{23} \\\\ r_{31} & r_{32} & r_{33} \\end{bmatrix}\\begin{bmatrix} x-x_0 \\\\ y-y_0 \\\\ z-z_0 \\end{bmatrix}\n",
+ "\\end{align}\n",
+ "$$\n",
+ "\n",
+ "with:\n",
+ "* $u$, $v$ the position of the animal on the camera\n",
+ "* $x$, $y$, $z$ the position of the animal in the environment\n",
+ "* $u_0$, $v_0$, the center of the camera\n",
+ "* $x_0$, $y_0$, $z_0$ the center of the environment (origin)\n",
+ "* $d$, the distance between the point on the camera and the point in the environment\n",
+ "* $c$, a constant of colinarity\n",
+ "\n",
+ "**Note** The unit of $u$, $v$ , $u_0$, and $v_0$ have the same units than $x,y,z$. However, $u$ and $v$ will be usually measured in pixels, and thus we need to further transform $u$ and $v$ by a constant of proportionality. $u-u_0 \\Rightarrow \\lambda_u(u-u_0)$, and $v-v_0 \\Rightarrow \\lambda_v(v-v_0)$.\n",
+ "\n",
+ "#### Calibration\n",
+ "\n",
+ "We ultimatly want to have the position of the animal within its environment from two or more cameras. To be able to reconstruct (triangulate) the animal position, we first need to have a method to describe the transformation from camera to world and vice versa, i.e. obtained the parameters of the transformation from known points. This process is called camera calibration.\n",
+ "\n",
+ "We need to express $u$ and $v$ as a function of $x$, $y$, and $z$. \n",
+ "\n",
+ "$$\n",
+ "u = \\frac{L_1x+L_2y+L_3z+L_4}{L_9x+L_{10}y+L_{11}z+L_1}\\text{, and }\n",
+ "v = \\frac{L_5x+L_6y+L_7z+L_8}{L_9x+L_{10}y+L_{11}z+L_1}\n",
+ "$$\n",
+ "\n",
+ "the coefficients $L_1$ to $L_{11}$ are the DLT parameters that reflect the relationships between the environment frame of our animal and the image frame. \n",
+ "\n",
+ "The system of equation above can rewritten as a product of matrices:\n",
+ "$$\n",
+ "\\begin{bmatrix} u \\\\ v \\end{bmatrix} = \\begin{bmatrix} x & y & z & 1 & 0 & 0 & 0 & 0 & -ux & -uy & -uz \\\\ 0 & 0 & 0 & 0 & x & y & z & 1 & -vx & -vy & -vz \\end{bmatrix}\\begin{bmatrix} L_1 \\\\ L_2 \\\\ . \\\\ L_{10} \\\\ L_{11} \\end{bmatrix}\n",
+ "$$\n",
+ "\n",
+ "We can expand this equations for $n$ points visible on the camera:\n",
+ "\n",
+ "$$\n",
+ "\\begin{bmatrix} \n",
+ "u_1 \\\\ v_1 \\\\\n",
+ "u_2 \\\\ v_2 \\\\\n",
+ "... \\\\\n",
+ "u_{n-1} \\\\ v_{n-1} \\\\\n",
+ "u_n \\\\ v_n \\\\\n",
+ "\\end{bmatrix} = \n",
+ "\\begin{bmatrix} \n",
+ "x_1 & y_1 & z_1 & 1 & 0 & 0 & 0 & 0 & -u_1x_1 & -u_1y_1 & -u_1z_1 \\\\ \n",
+ "0 & 0 & 0 & 0 & x_1 & y_1 & z_1 & 1 & -v_1x_1 & -v_1y_1 & -v_1z_1 \\\\\n",
+ "x_2 & y_2 & z_2 & 1 & 0 & 0 & 0 & 0 & -u_2x_2 & -u_2y_2 & -u_2z_2 \\\\ \n",
+ "0 & 0 & 0 & 0 & x_2 & y_2 & z_2 & 1 & -v_2x_2 & -v_2y_2 & -v_2z_2 \\\\\n",
+ ". \\\\\n",
+ "x_{n-1} & y_{n-1} & z_{n-1} & 1 & 0 & 0 & 0 & 0 & -u_{n-1}x_{n-1} & -u_{n-1}y_{n-1} & -u_{n-1}z_{n-1} \\\\ \n",
+ "0 & 0 & 0 & 0 & x_{n-1} & y_{n-1} & z_{n-1} & 1 & -v_{n-1}x_{n-1} & -v_{n-1}y_{n-1} & -v_{n-1}z_{n-1} \\\\\n",
+ "x_n & y_n & z_n & 1 & 0 & 0 & 0 & 0 & -u_nx_n & -u_ny_n & -u_nz_n \\\\ \n",
+ "0 & 0 & 0 & 0 & x_n & y_n & z_n & 1 & -v_nx_n & -v_ny_n & -v_nz_n \\end{bmatrix}\n",
+ "\\begin{bmatrix} L_1 \\\\ L_2 \\\\ . \\\\ L_{10} \\\\ L_{11} \\end{bmatrix}\n",
+ "$$\n",
+ "\n",
+ "\n",
+ "If we have 11 equations or more we can could deterime the DLT parameters from the system of equations. Since each points in the environment will lead to two variables $u$ and $v$, we only need 6 points to be able to obtain the DLT parameters\n",
+ "\n",
+ "#### Reconstruction\n",
+ "\n",
+ "\n",
+ "#### Exercises\n",
+ "* Derive from the $r_11$ to $r_33$, $d$ and $c$ what are $L_1$ to $L_{11}$\n",
+ "\n",
+ "#### 16 DLT parameters\n",
+ "The DLT method exposed above does not correct for optical distortion nor de-centering distortion. \n",
+ "\n",
+ "$$\n",
+ "\\frac{1}{R}\\begin{bmatrix} u \\\\ v \\end{bmatrix} = \n",
+ "\\frac{1}{R}\\begin{bmatrix} \n",
+ "x & y & z & 1 & 0 & 0 & 0 & 0 & -ux & -uy & -uz & \\zeta r^2R & \\zeta r^4R & \\zeta r^6R & (r^2+2\\zeta^2)R & \\zeta\\eta R \\\\ \n",
+ "0 & 0 & 0 & 0 & x & y & z & 1 & -vx & -vy & -vz & \\eta r^2R & \\eta r^4R & \\eta r^6R & \\zeta\\eta R & (r^2+2\\eta^2)R\\\\ \\end{bmatrix}\\begin{bmatrix} L_1 \\\\ L_2 \\\\ . \\\\ L_{15} \\\\ L_{16} \\end{bmatrix}\n",
+ "$$\n",
+ "\n",
+ "with \n",
+ "* $R=L_9x + L_{10}y + L_{11}z +1$, \n",
+ "* $r^2 = \\eta^2 + \\zeta^2$, \n",
+ "* $\\zeta = u-u_0$, and \n",
+ "* $\\eta = v-v_0$ "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Alternative method\n",
+ "\n",
+ "A camera is an optical intrument for aquiring images. Cameras are composed of a sensors (converting photon to electrical signal) and a lens (foccussing light rays on the sensors). \n",
+ "Many experimental paradigm require the recording of animal motion with one or more camera. \n",
+ "The experimenter is more interested by the position of the animal in his or her arena (world \n",
+ "coordinate system) than by the position of the animal within the image (camera-coordinate system). Therefore, we need to transform the position of the animal in the image to its position in the world. To do so, we need to calibrate our camera. This calibration require two steps:\n",
+ "\n",
+ "* remove distortion from the lens (intrinsic calibration)\n",
+ "* determine the position and orientation (i.e. pose) of the camera in the environment (extrinsic calibration)\n",
+ "\n",
+ "#### Intrinsic calibration\n",
"\n",
"`With Matlab <https://de.mathworks.com/help/vision/ug/single-camera-calibrator-app.html>`_\n",
"\n",
"`With Opencv <https://docs.opencv.org/3.3.1/dc/dbb/tutorial_py_calibration.html>`_\n",
"\n",
- "### Extrinsic calibration\n",
+ "#### Extrinsic calibration\n",
"\n",
"To obtain the pose of the camera (i.e. the extrinsic parameters), we need \n",
"to find a transformation such that the projection of physical points matche\n",
@@ -33,7 +135,7 @@
"Using Ransac decrease the bad effect that outliers may have on the pose\n",
"estimation.\n",
"\n",
- "#### Building your own manhattan\n",
+ "## Building your own manhattan\n",
"\n",
"The manhattan should contains at least 7 towers, and not all towers should\n",
"share the same plane or line. For best results it is recommended that the\n",
diff --git a/doc/source/tutorials/02c-orientation-2markers.ipynb b/doc/source/tutorials/02c-orientation-2markers.ipynb
index 9f7177ef167ec8a8a3acda0656837762fa5e2a9a..b0ceeba4a707ef68c42ccfbd1c62f426f3d42c70 100644
--- a/doc/source/tutorials/02c-orientation-2markers.ipynb
+++ b/doc/source/tutorials/02c-orientation-2markers.ipynb
@@ -13,61 +13,12 @@
" R.v^{ref}=v^{bee}\n",
" \\end{equation}\n",
"$$\n",
- " \n",
+ "\n",
"here $v^{bee}$ is the vector given by the 3D coordinates of the two markers, $v^{ref}$ is the vector given by the 3D coordinates of the two markers when the orientation of the agent is null, i.e. the orientation of the matrix is equal to the identity matrix. \n",
"\n",
"The system has 3 equations and 9 unknown variables: the elements of the orientation matrix.\n",
"\n",
- "**Note** Why 9 unknown variables if we have three rotation angles? The rotation angles appear in the orientation matrix in cosine and sine function, both nonlinear function. Therefore the sine and the cosine of the rotation angle have to be treated separatly. Moreover the multiplication of two cosines, two sine, or one cosine and one sine has to be treated has variables, because multiplication is... nonlinear.\n",
- "\n",
- "$$\\begin{align}\n",
- " v_x^{bee} & = +v_x^{ref} \\cos\\alpha \\cos\\beta + v_y^{ref}(\\cos\\alpha \\sin\\beta \\sin\\gamma - \\sin\\alpha \\cos\\gamma) + v_z^{ref} (\\cos\\alpha \\sin\\beta \\cos\\gamma + \\sin\\alpha \\sin\\gamma) \\\\\n",
- " v_y^{bee} & = +v_x^{ref} \\sin\\alpha \\cos\\beta + v_y^{ref}(\\sin\\alpha \\sin\\beta \\sin\\gamma + \\cos\\alpha \\cos\\gamma) + v_z^{ref} ( \\sin\\alpha \\sin\\beta \\cos\\gamma - \\cos\\alpha \\sin\\gamma )\\\\\n",
- " v_z^{bee} & = -v_x^{ref}\\sin\\beta + v_y^{ref} \\cos\\beta \\sin\\gamma + v_z^{ref} \\cos\\beta \\cos\\gamma\n",
- "\\end{align}$$\n",
- "\n",
- "or equivalently we can look at $v^{ref}=R^Tv^{bee}$, because $R^T=R^{-1}$\n",
- "\n",
- "$$\\begin{align}\n",
- " v_x^{ref} & = +v_x^{bee} \\cos\\alpha \\cos\\beta +v_y^{bee} \\sin\\alpha \\cos\\beta -v_z^{bee}\\sin\\beta \\\\\n",
- " v_y^{ref} & = +v_x^{bee}(\\cos\\alpha \\sin\\beta \\sin\\gamma - \\sin\\alpha \\cos\\gamma) +v_y^{bee}(\\sin\\alpha \\sin\\beta \\sin\\gamma + \\cos\\alpha \\cos\\gamma)+ v_z^{bee} \\cos\\beta \\sin\\gamma \\\\\n",
- " v_z^{ref} & = +v_x^{bee}(\\cos\\alpha \\sin\\beta \\cos\\gamma + \\sin\\alpha \\sin\\gamma) + v_y^{bee} ( \\sin\\alpha \\sin\\beta \\cos\\gamma - \\cos\\alpha \\sin\\gamma )+v_z^{bee} \\cos\\beta \\cos\\gamma\n",
- " \\end{align}$$\n",
- "\n",
- "To simplify the problem, we need to remove terms in the system of equations. Removing terms is easily done by letting certain variables to be zero. We want to determine $\\alpha$, $\\beta$, and $\\gamma$. Thus we can only set to zeros $v_x^{ref}$, $v_y^{ref}$, or $v_z^{ref}$.\n",
- "\n",
- "### A simple case\n",
- "\n",
- "If we assume that the two markers are aligned with the roll axis,\n",
- "i.e. $v^{ref}=(1,0,0)^T$, the problem can easily be solve as follow:\n",
- "\n",
- "\n",
- "$$\n",
- " \\begin{align}\n",
- " v_x^{bee} & = +\\cos\\alpha \\cos\\beta \\\\\n",
- " v_y^{bee} & = +\\sin\\alpha \\cos\\beta \\\\\n",
- " v_z^{bee} & = -\\sin\\beta\n",
- " \\end{align}\n",
- "$$\n",
- "\n",
- "\n",
- "$$\n",
- " \\begin{align}\n",
- " \\tan\\alpha & = \\frac{\\pm v_y^{bee}}{\\pm v_x^{bee}} &\\quad \\text{from L1 and L2} \\\\\n",
- " \\tan\\beta & = \\frac{-v_z^{bee}}{\\pm\\sqrt{v_y^{bee}+v_x^{bee} }} &\\quad\\text{from all}\n",
- " \\end{align}\n",
- "$$\n",
- "\n",
- "We remark that the solution does not depend of the angle $\\gamma$, \n",
- "the roll angle. Thus, when we do not care about the full orientation of \n",
- "a body but simply care about pitch and yaw, two markers are sufficients.\n",
- "\n",
- "### Other cases\n",
- "\n",
- "The two markers may be aligned with the pitch (or yaw axis) of the body. \n",
- "In such situation, the problem is slightly more complex. We can, still \n",
- "solve the system of equation, but can not find a solution independent of \n",
- "the pitch (or yaw angle), we need to guess the pitch (or yaw angle)."
+ "**Note** Why 9 unknown variables if we have three rotation angles? The rotation angles appear in the orientation matrix in cosine and sine function, both nonlinear function. Therefore the sine and the cosine of the rotation angle have to be treated separatly. Moreover the multiplication of two cosines, two sine, or one cosine and one sine has to be treated has variables, because multiplication is... nonlinear."
]
},
{
@@ -76,6 +27,10 @@
"source": [
"### Angles from two markers\n",
"\n",
+ "#### Simples cases\n",
+ "\n",
+ "In navipy, only few cases have been implemented (see below to derive the equations for a custom cases). The rotation $R$ can be written as a series of rotation along x, y, and z-axis (see: Euler angles). Simple case to extract the euler angles with only two markers arise when the first applied rotation is aligned with our two markers. For example for $R=R_zR_yR_x$, and x-aligned marker, the decomposition is simple. Indeed our measure markers can be obtained by rotating the x-axis, here the first rotation along the x-axis does not change the orientation of the x-axis. Thus, the euler angle around the x-axis will not affect the euler angles around the two other axis. \n",
+ "\n",
"Here we illustrate, how the solution of the system in three cases: \n",
"roll-aligned, pitch-aligned, and yaw-aligned, varies as a function of \n",
"the a priori known angle.\n",
@@ -85,241 +40,170 @@
"\n",
"* The rotation around x, with convention $R_zR_yR_x$, for x-aligned markers\n",
"* The rotation around y, with convention $R_zR_xR_y$, for y-aligned markers\n",
- "* The rotation around z, with convention $R_yR_xR_z$, for x-aligned markers\n",
- "\n",
- "Thus the none of the yaw pitch roll angles may match of the a priori \n",
- "known angle (except for x-aligned markers, because the internal convention is the yaw pitch roll convention)"
+ "* The rotation around z, with convention $R_yR_xR_z$, for x-aligned markers"
]
},
{
"cell_type": "code",
- "execution_count": 1,
+ "execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"import pandas as pd\n",
"import numpy as np\n",
- "import matplotlib.pyplot as plt\n",
- "from navipy.maths.homogeneous_transformations import compose_matrix\n",
- "import navipy.trajectories.transformations as mtf\n",
- "from navipy.trajectories.triangle import Triangle"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Similar to the notebook about the background of the orientation, we place a triangle at a known position orientation"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 2,
- "metadata": {},
- "outputs": [],
- "source": [
- "yaw_t = +3*np.pi/4\n",
- "pitch_t = -1*np.pi/6\n",
- "roll_t = -1*np.pi/12\n",
- "angles = [yaw_t, pitch_t, roll_t]\n",
- "euler_axes = 'rzyx'\n",
- "# Create a triangle in a given orientation\n",
- "# and get the x,y,z axis used as our two markers\n",
- "triangle_mode = 'x-axis=median-from-0'\n",
- "transform = compose_matrix(angles=angles,\n",
- " axes=euler_axes)\n",
- "markers = pd.Series(data=0,\n",
- " index=pd.MultiIndex.from_product(\n",
- " [[0, 1, 2], ['x', 'y', 'z']]))\n",
- "markers.loc[(0, 'x')] = -1\n",
- "markers.loc[(2, 'y')] = np.sin(np.pi / 3)\n",
- "markers.loc[(1, 'y')] = -np.sin(np.pi / 3)\n",
- "markers.loc[(1, 'x')] = np.cos(np.pi / 3)\n",
- "markers.loc[(2, 'x')] = np.cos(np.pi / 3)\n",
- "equilateral = Triangle(markers.loc[0],\n",
- " markers.loc[1],\n",
- " markers.loc[2])\n",
- "equilateral.transform(transform)\n",
- "_, x_axis, y_axis, z_axis = mtf.triangle2bodyaxis(\n",
- " equilateral, triangle_mode)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "If can only can get the position of two markers out of three, we have only access to a vector in 3D space and not a triangle. \n",
+ "from navipy.trajectories.transformations import twomarkers2euler\n",
"\n",
- "We need to assume that this vector is align to one axis, and that one of the three Euler angles is known. When we do not know the angle, we can look at all possible angles, for our known angles."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 3,
- "metadata": {},
- "outputs": [],
- "source": [
- "known_angles = np.linspace(-np.pi, np.pi, 180)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Let's assume that the two markers are aligned with one the axes\n",
- "1. along the x-axis,\n",
- "2. along the y-axis, or\n",
- "3. along the z-axis."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 4,
- "metadata": {},
- "outputs": [],
- "source": [
+ "mark0 = pd.Series(index=['x','y','z'], data=np.random.rand(3))\n",
+ "mark1 = pd.Series(index=['x','y','z'], data=np.random.rand(3))\n",
+ "known_angle = 2*np.pi*np.random.rand()\n",
+ "# x-aligned\n",
"axis_alignement = 'x-axis'\n",
- "mark0 = pd.Series(data=0, index=['x', 'y', 'z'])\n",
- "mark1 = pd.Series(x_axis, index=['x', 'y', 'z'])\n",
- "solution_x_axis = pd.DataFrame(index=known_angles, columns=['yaw',\n",
- " 'pitch',\n",
- " 'roll'])\n",
- "for known_angle in known_angles:\n",
- " angles_estimate = mtf.twomarkers2euler(\n",
- " mark0, mark1, axis_alignement, known_angle, euler_axes)\n",
- " solution_x_axis.loc[known_angle, :] = angles_estimate"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 5,
- "metadata": {},
- "outputs": [],
- "source": [
+ "euler_axes = 'zyx'\n",
+ "alpha, beta, gamma = twomarkers2euler(mark0, mark1, axis_alignement,\n",
+ " known_angle, euler_axes)\n",
+ "# y-aligned\n",
"axis_alignement = 'y-axis'\n",
- "mark0 = pd.Series(data=0, index=['x', 'y', 'z'])\n",
- "mark1 = pd.Series(y_axis, index=['x', 'y', 'z'])\n",
- "solution_y_axis = pd.DataFrame(index=known_angles, columns=['yaw',\n",
- " 'pitch',\n",
- " 'roll'])\n",
- "for known_angle in known_angles:\n",
- " angles_estimate = mtf.twomarkers2euler(\n",
- " mark0, mark1, axis_alignement, known_angle, euler_axes)\n",
- " solution_y_axis.loc[known_angle, :] = angles_estimate"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 6,
- "metadata": {},
- "outputs": [],
- "source": [
+ "euler_axes = 'zxy'\n",
+ "alpha, beta, gamma = twomarkers2euler(mark0, mark1, axis_alignement,\n",
+ " known_angle, euler_axes)\n",
+ "\n",
+ "# z-aligned\n",
"axis_alignement = 'z-axis'\n",
- "mark0 = pd.Series(data=0, index=['x', 'y', 'z'])\n",
- "mark1 = pd.Series(z_axis, index=['x', 'y', 'z'])\n",
- "solution_z_axis = pd.DataFrame(index=known_angles, columns=['yaw',\n",
- " 'pitch',\n",
- " 'roll'])\n",
- "for known_angle in known_angles:\n",
- " angles_estimate = mtf.twomarkers2euler(\n",
- " mark0, mark1, axis_alignement, known_angle, euler_axes)\n",
- " solution_z_axis.loc[known_angle, :] = angles_estimate"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "The possible two other angles are shown below, as a function of the known angle"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 7,
- "metadata": {},
- "outputs": [
- {
- "data": {
- "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3oAAAEWCAYAAAAw+uVwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XdcleX/x/HXxZ6CDBFFRJyAe+DOnaZm2rDhHllmmdavsr2s7JulmavcKzXTsnKUe6SmmBv3QFAcgKjI5ly/Pw6OyoF64D4cPs/Hg0fCuc99v493XpzPuZbSWiOEEEIIIYQQwnbYGR1ACCGEEEIIIYRlSaEnhBBCCCGEEDZGCj0hhBBCCCGEsDFS6AkhhBBCCCGEjZFCTwghhBBCCCFsjBR6QgghhBBCCGFjpNATBU4pdUIp1Tr3z28ppSYbkKG5UiquoK8rhLBu0j4JIayRUmqtUqp/7p+7KaX+MCBDiFJKK6UcCvra4t7IjRKG0lp/anQGIYS4GWmfhBDWSGs9B5hjdA5h/aRHTwghhBBCCCFsjBR6Ik+UUuWVUklKqdq535dSSp1XSjW/xbGrlVKJSqkEpdQcpZT3Lc77gVJq9g3f91RKxeQ+991/DaP6QCn1g1JqplLqslJqn1Kq7g3PLaWUWpib67hSavANj7kqpaYrpS4opaKBepb72xFCGEkp9ZpSauG/fjZGKfX1TY6V9kkIUWCUUk8qpVJu+MpQSq29yXHFlVK/5bYRF3L/HHSLc/ZWSm284fsHlVIHlVIXlVLjlVLrbhjm2VsptVEpNTL3vMeVUg/d8FwvpdQUpVS8UuqUUmq4Uso+9zH73OclKKWOAR0s/fcj8pcUeiJPtNZHgTeA2UopN2AaMENrvfYmhyvgM6AUEAaUAT640zWUUuHAeKAbEAh4AaX/dVgnYB7gDfwCjM19rh3wK7Ar9zmtgCFKqba5z3sfKJ/71RbodedXLYQoJGYD7a4WbLnzR54CZt7kWGmfhBAFRms9X2vtobX2wNzuHAPm3uRQO8zvrcoCwUAauW3I7Sil/IAfgTcBX+Ag0Ohfh9XP/bkf8D9gilJK5T42HcgGKgC1gAeB/rmPPQt0zP15XeDxO75gYVWk0BN5prWeBBwB/sL8RuftWxx3RGu9QmudobU+D3wFNMvDJR4HftVab9RaZwLvAfpfx2zUWi/VWucAs4AauT+vB/hrrT/SWmdqrY8BkzC/2QPoCnyitU7SWscCY/L6uoUQ1k1rHQ+sB57I/VE7IEFrvf0mx0r7JIQocLkf+HwPrNVaf/vvx7XWiVrrhVrrVK31ZeAT8tY2tQf2aa0Xaa2zMbcfZ/51TIzWelJu2zQD83u4AKVUQO7zh2itr2itzwGj+GfbNFprHau1TsL8IZkoRGQxFnG3JmH+pHqA1jpDKdUUWJb7WIzWOiK34fgaaAp4Yv5A4UIezl0KiL36jdY6VSmV+K9jbmy8UgGX3E/vywKllFLJNzxuD2y42bmBmDzkEUIUHjOAgZjbqO7ALGmfhBBW5BPMbc5gpVQwEH31Aa21R+5oqVGYP6gqnvuQp1LKPrdAu5V/t01a/XfV3jM3PJ6a25nnAfgAjkD89Q4+7G44n7RNhZwUeiLPlFIewGhgCvCBUmqh1noD5sbiRp9i/qS7mtY6SSnVmTwMPwDigco3XM8V8zCEvIgFjmutK97m3GWAfbnfB+fxvEKIwuFnYIJSqirmoUava61PIu2TEMJgSqmngKeBelrrLOBmbdOrmNuY+lrrM0qpmsAOzMPNbyceuDaXL3dI5k3n9t1ELJAB+OX2Bt7s3GVu+F7apkJGhm6Ku/E1EKW17g8sASbe4jhPIAW4qJQqDbyWx/P/CDyslGqklHLCPG/mTg3cVVuBy0qpN3IXNrBXSlVVSl1d1OAH4M3cyc5BwEt5PK8QohDQWqdjbkO+B7bmFnk3I+2TEKLAKKVqAd8AnXOHi9+KJ+Z5eclKKR/Mc3fzYglQTSnVOXcEwSCgZF6emDvs/Q/gS6VUMaWUnTIvWHV1yOgPmHsgg5RSxYFhecwkrIQUeiJPlFKPYB5OMDD3R68AtZVS3W5y+IdAbeAi5gZoUV6uobXeh/kNzjzMnyKlAOcwf9p0p+fmYP4UvyZwHEgAJmNeMOFqppjcx/7APH9GCGFbZgDVuP2/b2mfhBAF6RHMQzE33rDy5rKbHDcacMXcPmwBlufl5FrrBMzzk/8HJALhQBR5aJty9QScMA8lvYD5Q63A3McmAb9jXkjqb/LYXgrrobT+91xyIaxD7lDRZKCi1vq40XmEENYtd97LAaCk1vpSPl9L2ichhNXJXfQlDuimtV5jdB5hLOnRE1ZFKfWwUspNKeUOjAT2ACeMTSWEsHa5b25eAeblV5En7ZMQwhoppdoqpbyVUs7AW5iHlW8xOJawAlLoCWvzCHA696si8JSWbmchxG3kFl6XgDbkfV7LvZD2SQhhjRoCRzEP+3wY83zANGMjCWsgQzeFEEIIIYQQwsZIj54QQgghhBBC2JhCtY+en5+fDgkJMTqGEMKCtm/fnqC19jc6x/2QtkkI2yTtkxDCGuW1bSpUhV5ISAhRUVFGxxBCWJBSKsboDPdL2iYhbJO0T0IIa5TXtkmGbgohhBBCCCGEjZFCTwghhBBCCCFsjBR6QgghhBBCCGFjCtUcPSGEEELkv6ysLOLi4khPTzc6SoFwcXEhKCgIR0dHo6MUCLm/QhQNUugJIYQQ4h/i4uLw9PQkJCQEpZTRcfKV1prExETi4uIoV66c0XEKhNxfIYoGGbophBBCiH9IT0/H19fX5osAAKUUvr6+RaZ3C+T+ClFUSKEnhBBCiP8oCkXAVUXptV5VlF5zUXqtQtxICj0hhEX9dPgn1sWuMzqGEEIIIYTVM5k0o1Yc4uCZyxY/txR6QgiL0FozYdcE3tv0HosOLzI6jrByB89c5mJqltExhBBCCMNkZOcw9IedfL3qMEt2n7b4+WUxFiHEfcsyZTF8y3AWHV5Ep/Kd+KDRB0ZHElau7ej1FHNxYPcHbY2OIoQQQhS4i2lZPD9rO5uPJfJ6u8oMbFbe4tcwrEdPKeWilNqqlNqllNqnlPrQqCxCiHuXmpXK4NWDWXR4Ec9Vf47hjYfjaCdLWIs7u5SeTXaOyegYwgq99957jB49+tr3b7/9Nl9//TWtWrWidu3aVKtWjcWLFwPwxRdfMGbMGACGDh1Ky5YtAVi9ejXdunUr+PDijuT+iqLudHIaXSduJiomiVFP1uCF5hXyZS6pkT16GUBLrXWKUsoR2KiUWqa13mJgJiHEXUhIS2DQqkEcSDrAew3f44lKTxgdSRQyUTEXaBDqa3QMcRsf/rqP6NOXLHrO8FLFeP/hiFs+3rdvXx599FGGDBmCyWRi3rx5bNq0iT59+lCsWDESEhJo0KABnTp1omnTpnz55ZcMHjyYqKgoMjIyyMrKYsOGDTzwwAMWzW2L5P4KUbAOnLlE76nbSMnIZnqfSBpX8Mu3axlW6GmtNZCS+61j7pc2Ko8Q4u4cv3icgSsHkpSexJgWY2hWppnRkUQhtDL6rBR64j9CQkLw9fVlx44dnD17llq1auHj48PQoUNZv349dnZ2nDp1irNnz1KnTh22b9/OpUuXcHZ2pnbt2kRFRbFhw4ZrPUHCusj9FUXVpiMJPDdrO27O9vzwXEPCSxXL1+sZOkdPKWUPbAcqAOO01n/d5JgBwGuAt7+/fwEnFELczM5zO3lx9YvYK3umtp1KVb+qRkcqcNI23R9nBzsysk0s3RPPm+3DsLeT5c+t1e16ZvJT//79mT59OmfOnKFv377MmTOH8+fPs337dhwdHQkJCSE9PR1HR0fKlSvH9OnTadSoEdWrV2fNmjUcOXKEsLAwQ7Ib7W7aJ7m/QhSMn3ec4rUfd1HOz53pfSIp5e2a79c0dNVNrXWO1romEAREKqX+825Ra/2d1rqi1to/ODi44EMKIf5hZcxK+v/RH29nb2Y/NLtIFnkgbdP90kA5P3dOX0xn1f6zRscRVqhLly4sX76cbdu20bZtWy5evEiJEiVwdHRkzZo1xMTEXDu2adOmjBw5kgceeICmTZsyceJEatWqVWT3TysM7ZPcX1FUaK0Zv/YIQ+bvpE7Z4ix4vlGBFHlgJdsraK2TgTVAO6OzCCFubc7+Obyy9hUq+1Rm1kOzKFOsjNGRRCGltebB8AACvVyYuTnmzk8QRY6TkxMtWrSga9eu2Nvb061bN6KioqhWrRozZ86kSpUq145t2rQp8fHxNGzYkICAAFxcXGjatKmB6cWdyP0VRUGOSfPu4r38b/lBHq5Rihl9I/FyLbgF6wwbuqmU8geytNbJSilXoA3wuVF5hBC3ZtImRm0fxfR902lZpiUjHhiBq0PBfBolbJNJg4O9onuDsnzx+0F2xyVTPcjb6FjCiphMJrZs2cKCBQsA8PPzY/PmzTc9tlWrVmRlXd+X8dChQwWSUdw7ub/C1qVl5vDS3B2s3H+W55qF8kbbKtgV8DQFI3v0AoE1SqndwDZghdb6NwPzCCFuIjMnk2HrhzF933SeqvwUXzX/Soo8cd+01tgpRc+GZfFxd2LEsgOY1+gSAqKjo6lQoQKtWrWiYsWKRscRFib3V9i6xJQMnp60hVUHzvJhpwjefCiswIs8MHbVzd1ALaOuL4S4s0uZl3h59ctEnY1iaJ2h9InoI3MihEWYNCil8HRx5KWWFfjw12hWRJ/lwYiSRkcTViA8PJxjx44ZHUPkE7m/wpbFJF6h19StxF9MZ0K3OrSratzvNauYoyeEsD7xKfH0WtaLned3MqLpCPpW7StFnrCIqz13V/9v6la/LFVKevL2z3tJTs00LpgQQghxH3bGJvPo+E0kp2Xx/bP1DS3yQAo9IcRNHEw6SPel3Tlz5Qzftv6WDqEdjI4kbMjVEZp2uR8cODnYMfKJGly4kslrP+7GZJIhnEIIIQqXldFneeq7zbg527NwYCPqlPUxOpIUekKIf9p8ejO9lvdCKcWMh2YQGRhpdCRhY0y5ld6N0xWqlvbizfZhrIg+y+hVhw1KJoQQQty9OX/FMGBWFBVLeLJoYGPK+3sYHQkweMN0IYR1+fXor7z353uU8y7H+FbjKeku86WE5V3tsPv3SOC+jUM4EH+JMasO4+PmSO/G5Qo+nBBCCJFHWmtG/nGQcWuO0qKyP2OfqY27s/WUV9KjJ4RAa82k3ZN4a+Nb1Amow4x2M6TIE/lGkztH71+VnlKKTx+txoPhAXzwazSjVx6SlTjFf/Tv35/o6GgAPv300zse37t3b3788cf8jiUsRO6vKCwys028umAX49Yc5al6ZZjUs65VFXkghZ4QRV62KZvhW4YzZscYOoR2YELrCXg6eRodS9iwf8/Ru5GjvR1jn6nN43WCGL3yMEPm7yQ9K6eAEwprNnnyZMLDw4G8FQKicJH7KwqDy+lZ9J2+jUV/n+KVNpX47NFqONhbX1llfYmEEAUmNSuVoWuG8sOhH+hXtR+fNvkUR3tHo2MJG3d1jt6tFnF1crDji8er81rbyizeeZrO4/5k76mLBZhQWIMTJ05QpUoVunXrRlhYGI8//jipqak0b96cqKgohg0bRlpaGjVr1qRbt24AzJw5k+rVq1OjRg169Ohx7Vzr16+nUaNGhIaGSu+PlZD7KwqrMxfTeWLiZrYcS+SLx6szuFVFq12V3Lr6F4UQBSYxLZGXVr/EvsR9vF3/bZ6q8pTRkUQRcb1H79bHKKUY1KICYYGeDFu4h0fG/ckLzcszsHl53JzkV1eBWjYMzuyx7DlLVoOHRtzxsIMHDzJlyhQaN25M3759GT9+/LXHRowYwdixY9m5cycA+/btY/jw4WzatAk/Pz+SkpKuHRsfH8/GjRs5cOAAnTp14vHHH7fs6ynM5P4KkWeHzl6m99StXEzLYkrvejSr5G90pNuSHj0hiqCYSzH0WNaDwxcOM6r5KCnyRIG6vurmnT8BbVklgBVDm/FIjVJ8s/oIzb9Yy9ytJ8nOMeV3TGEFypQpQ+PGjQHo3r07GzduvOWxq1ev5oknnsDPzw8AH5/rS5t37twZOzs7wsPDOXv2bP6GFnkm91cUJpuPJvLYhE1kmTTzn2to9UUeSI+eEEXOrvO7eGnVSwBMbjuZGv41DE4kipq73SbPy82Rr56syTP1g/ls2QHeXLSH8WuP0L9JKE/UDZIevvyWh56Z/HKzBXvuhbOz87U/ywI//yL3V4g7+mXXaf7vh12U8XFlRt9Igoq7GR0pT6RHT4giZM3JNfT/vT8eTh7Maj9LijxhjNssxnI7dUN8+PH5hkzuWZcSni68/8s+Gny6ind/3sueuIvyBs8GnTx5ks2bNwPw/fff06RJk3887ujoSFZWFgAtW7ZkwYIFJCYmAvxjaJ+wTnJ/hbXTWjNp/TEGz91BjTJeLBzYqNAUeSCFnhBFxvwD8xmydggVvCsw66FZlC1W1uhIooi62YbpeaWUonV4AAsHNmLhwIY0r1yC+VGxPDx2Iw99vYFv1x0lJvGKhRMLo1SuXJlx48YRFhbGhQsXGDhw4D8eHzBgANWrV6dbt25ERETw9ttv06xZM2rUqMErr7xiUGqRV3J/hTXLMWk+/DWaT5bup321kszqVx9vNyejY90VVZg+Aa1bt66OiooyOoYQhYpJmxjz9xim7J1Cs6Bm/O+B/+HmaD2fRimltmut6xqd435I23R3ElMyqDN8JR92iqBXo5D7Pt/FtCx+3XWaBdvj2BWbDECVkp60jShJyyolqFraC/t7qSqLsP379xMWFmZohhMnTtCxY0f27t1bINe72Wu21fZJ7q8Qt5eelcOQeTtZvu8M/ZqU4+32YdhZ0e+RvLZNNjmx4cNf9xF9+pLRMYQwnCab044zuGj/F8WzH+Dckafoc2RXvl83vFQx3n84It+vIwqnqx8vWup3pperI90blKV7g7LEJqXy+74z/L7vDGNWH+brVYfxdnOkcXk/Glfwo2lFP8r4WM8HHUIIIazLhSuZ9J8Zxd8nL/BOhzD6Nw01OtI9s8lCTwgBOaQS6ziRVPsDlMjqgm9OOxTW82mUKLqu76Nn+f8fy/i40b9pKP2bhpKQksGfRxLYcDiBjYcTWLInHoDS3q7UDSlO3RAf6pYtTqUAT+nxs0IhISEF1tsjCp7cX2GNYpNS6TV1K3HJaYx9ujYdqgcaHem+2GShJz0Joqg7c+UML6x6gczk43za+FMeLv+w0ZGEuObqjIH83l/Wz8OZR2qW5pGapdFac/R8ChsOJ7DtRBKbjiayeOdpADydHahdtji1g4tTvYwXNYK88XEvXPMwhBBC3J89cRfpM30rWTma2f3qE1nO585PsnI2WegJUZQdvnCYgSsHkpKVwrjW42hUqpHRkYT4B32Pq27eD6UUFUp4UqGEJ30al0NrTWxSGlExSWw7cYGoE0msP3z+Wrag4q7UKONNjSAvqgd5U7W0Fx7O8itTCCFs0ZqD5xg052+Kuzkxb0A9KpTwNDqSRchvLSFsyNb4rQxZMwRXB1dmtJtBZZ/KRkcS4j/uZ9VNS1FKEezrRrCvG4/WDgLgcnoWe05dZHfcRXbHJbPzZDJLdsfnHg8V/D2oHuRNjdxevyqBnjg72Bv3IoQQQty3+dtO8tZPe6lS0pNpvetRopiL0ZEsRgo9IWzEkmNLeOfPdyjrWZYJrScQ6FG4x5UL23Vtjp6VzRn1dHGkUXk/GpX3u/azhJQM9sRdZFdcMrvjLrLu0DkW/h0HgJO9HWGlilGrjDc1y3hTo4w3Ib5u+TL3UAghhGVprRm18jBjVh3mgUr+jO9W2+ZGbtjWqxGiCNJaM23fNEZtH0XdgLqMbjEaL2cvo2MJcUsFNUfPEvw8nGlRpQQtqpQAzP/eTl9MZ1dsMrtik9kRm8z8bbFM33QCMK8AWiO38KuZ2/Pn6+Fs4CuwfSEhIURFReHn54eHhwcpKSlGRxIWJPdX5IesHBNvLdrDgu1xPFEniE8frYajve1tLy6FnhCFWI4phxFbRzDv4DweCnmI4U2G42Qvi0gI62bEHD1LUUpR2tuV0t6utK9m7jXPzjFx5HwKO08msysumR0nkxm7+jCm3NdZoYQH9cv50CDUl/qhPpTwtJ1hQQVFa43WGjs723sjJuT+ioKVkpHNC3P+Zv2h87zcqiJDWle02ZEYUugJUUilZafxxvo3WBO7ht4RvRlaZyh2Sn5JCut3fXsFg4NYiIO9HVVKFqNKyWI8FRkMQGpmNnviLrL95AW2Hk9i8c7TzPnrJAChfu7UD/WhWSV/mlT0t7mhQpZy4sQJ2rZtS/369dm+fTuvv/46I0eORGtNhw4d+Pzzz42OKO6D3F9hhHOX0ukzfRsHzlxmxKPVrrXZtsqw3y5KqTLATCAA8/6532mtvzYqjxCFyYX0C7y4+kX2nN/DsMhhdAvrZnQkIfLs+obpNlLp3YSbkwP1Q32pH+rLC83NvX7R8Zf461gSfx1P5Lfd8czdGoujvSKynA8tKpegVVgA5fzcjY7+H59v/ZwDSQcses4qPlV4I/KNOx53+PBhZsyYQXBwMA0aNGD79u0UL16cBx98kJ9//pnOnTtbNFdRJPdXFBVHzl2m19RtJF3JZHLPuteG5NsyIz9GzAZe1Vr/rZTyBLYrpVZoraMNzCSE1Yu9HMvAlQM5c+UMXzX/itZlWxsdSYi7Yms9ennhYG9H9SBvqgd58+wDoWTnmNgec4HVB86x+sA5hi/Zz/Al+wkPLEbnWqXoVKM0Jb1kiGfZsmVp0KABixcvpnnz5vj7+wPQrVs31q9fL4VAISf3VxSUbSeS6D8jCkd7xfznGlA9yNvoSAXCsEJPax0PxOf++bJSaj9QGpBCT4hb2Juwl0GrBpGjc5j04CRqlahldCQh7pq+VugVoUrvXxzs7a71+L3ZPozYpFT+iD7LLztP8enSA3y27ACNy/vRvUEwrcMCcDBwkYC89MzkF3d36+vhtDVyf4WtW7onniHzdxLk7cr0PpEE+7oZHanAWMWEHqVUCFAL+Osmjw1QSh1WSp0/efJkQUcTwmqsj1tP39/74urgyqyHZkmRZzBpm+7d9cVYjM1hTcr4uNGvSTkWv9iE1a82Y3DLihxPuMLzs/+m6f/WMGn9MVIzs42OaZjIyEjWrVtHQkICOTk5zJ07l2bNmhkdy2oVtvZJ7q/IL1M2HmfQ939TrbQXCwc2KlJFHlhBoaeU8gAWAkO01pf+/bjW+jutdUWttX9wsG1PmBTiVhYcWsBLq1+inFc5ZrefTTmvckZHKvKkbbp3pkK86mZBCPX3YGibSqx/vQWTetalnJ87nyzdT9PP1/Dd+qNkZOcYHbHABQYGMmLECFq0aEGNGjWoU6cOjzzyiNGxrFZha5/k/gpLM5k0H/8Wzce/RfNgeABz+tenuHvRW5Xc0KW+lFKOmIu8OVrrRUZmEcIaaa0Zu3Ms3+3+jialm/Blsy9xcyxan0YJ23N9w3RxO/Z2ijbhAbQJD2B7TBKjVx7m06UHmLs1lg86RdCskr/REfNVSEgIe/fuvfb9008/zdNPP/2f406cOHHtz7LHWuEh91fkl/SsHF5dsIslu+Pp3SiEdzuGY19Eh5AY1qOnzJMzpgD7tdZfGZVDCGuVZcrinT/f4bvd3/FoxUcZ03KMFHnCJlzfML1o/uK9F3XK+jCrX31m9o1EAb2mbuWV+Tu5klF0h3MKIcS/Jadm0nPKVpbsjuet9lV4/+GiW+SBsT16jYEewB6l1M7cn72ltV5qYCYhrEJKZgqvrH2FzfGbeaHmCzxf/Xl5UyxsxtUevSL8u/eePVDJn2VDmjJ+zVG+WX2YnXHJjO9WmyolixkdTQghDBV3IZXe07ZxMjGVMU/XolONUkZHMpyRq25uREbuCPEf51LPMWjVIA5fOMxHjT6iS8UuRkcSwqKkR+/+ODvYM7RNJRqE+jJ43g46j/uTsU/XpnV4gEWvo7UuMvfo6kqwRYncX2FL9p66SJ/p20jPymFG30galvc1OpJVMHwxFiHEdUeTj9J9aXdOXjrJuFbjpMgTNkkjPXqW0LC8L0sHN6VygCfPzd7Or7tOW+zcLi4uJCYmFok3yFprEhMTcXEpOvsWyv0VtmT9ofM8+e1mHO0UCwc2kiLvBoYuxiKEuC7qTBSD1wzG2d6Zae2mEe4bbnQkIfKFrLppOf6eznz/bAP6TNvG0Pk7cXe2p2WV++/ZCwoKIi4ujvPnz1sgpfVzcXEhKCjI6BgFRu6vsBULomJ5c9EeKpTwYHqfSEp6SUF/Iyn0hLACy08s560NbxHkGcSE1hMo7VHa6EhC5BvTtbGbxuawFe7ODkzpXZenJ21h8Nyd/DyoERVKeN7XOR0dHSlXTrZxsVVyf0Vhp7Vm7OojfLniEI0r+DKxex08XRyNjmV1ZOimEAabuW8mr617jap+VZn10Cwp8oTN09KjZ3GeLo5816MuLo52PD/7b9Kzit5ee0KIoiE7x8RbP+3hyxWH6FKrNNN6R0qRdwtS6AlhkBxTDp9v/Zwvor6gTdk2THpwEl7OXkbHEiLfaVl1M1+U8nZl1JM1OXIuhVErDhkdRwghLO5KRjbPzoxi7tZYBrUoz1dda+DkIOXMrcjfjBAGSM9O57X1rzF7/2y6h3VnZLORONs7Gx1LiAJhujZyUyo9S2ta0Z9n6gfz3YZj7I+/ZHQcIYSwmPOXM3h60hbWHTrP8M5Vea1tlSKzcuy9kkJPiAKWnJ7MgBUDWBmzktfqvsYbkW9gp+Sfoig6pEcvf73etjKezg58vvyA0VGEEMIijp1P4dEJf3Lo7GW+61GX7g3KGh2pUJB3l0IUoLjLcfRY1oN9Cfv4otkX9IzoaXQkIQqcSfbRy1febk4MalGBtQfPs/V4ktFxhBDivmyPucBjEzaRmpHDvAENLb5nqC2TQk+IAhKdGE33pd1JTE/kuwe/o21IW6MjCWGIqz16Uufln16NQiju5si0P48bHUUIIe7Z7/vO8MykLXi5OrJwYCNqlvE2OlKhIoWeEAVg46mN9F7eGyd7J2Y9NIs6AXWMjiSEYa5u0SyrbuYfF0fZ9/i1AAAgAElEQVR7utYrwx/RZzmdnGZ0HCGEuGszN5/g+dnbqRJYjIUDGxHi5250pEJHCj0h8tlPh3/ixVUvUrZYWWa3n0157/JGRxLCUCaZo1cgutcvi0lr5m09aXQUIYTIM5NJ89my/by3eB+tqpRg3rMN8PWQBevuhRR6QuQTrTUTdk7gvU3vUT+wPtPbTaeEWwmjYwlhuOtz9IzNYevK+LjRoJwvy/edMTqKEELkSUZ2DkPm7+Tbdcfo3iCYid3r4Opkb3SsQksKPSHyQZYpi/c3vc/4XeN5pPwjjG01FndHGXIgBFzv0ZPFWPJfm/AADp1NISbxitFRhBDiti6mZdFr6lZ+2XWa19tV5uNHquJgL6XK/ZC/PSEsLDUrlZdWv8RPR37i+RrP83Hjj3G0czQ6lhDWI7dHT+bo5b82uavTrYg+a3ASIYS4tdPJaTwxcRPbYy4w6skavNC8gnwYaAEORgcQwpYkpCXwwsoXOHThEO83fJ/HKz1udCQhrI7M0Ss4ZXzcqFLSk5XRZ+jvsRlWvAvpRXwj9aF7wbOk0SmEELn2x1+iz7RtXMnIZnqfSBpX8DM6ks2QQk8ICzl+8TgDVw4kKT2JMS3H8EDQA0ZHEsIqXZujh1R6BaFRqBeVoz6ExasguBEENzA6krEc3YxOIITI9eeRBJ6ftR13Zwd+eL4hYYHFjI5kU6TQE8ICdpzbwUurX8Je2TOt7TQi/CKMjiSE1ZJ99ApQ5hWejXuXQLt1JNV6EZ+HPwI7WdhACGG8n3bE8fqPuynn5870PpGU8nY1OpLNuWWhp5R6JQ/Pv6K1/taCeYQodFbErGDY+mEEegQyofUEyniWMTqSEFbNJHP0CkbKefj+CUqe38VbWf2oFzyQLlLkCSEMprVm/NqjfPH7QRqE+vBtj7p4ucpaBvnhdouxvAZ4AJ63+Xo1vwMKYc3m7J/Dq2tfJcw3jFkPzZIiT4g8kB69ApB4FKa0hnMHMHWdzSK7Nuw9VcTn5gkhDJdj0ry7eC9f/H6QTjVKMaNvpBR5+eh2Qzdnaa0/ut2TlVKyXrwokkzaxFdRXzEjegatglsxoukIXBxcjI4lRKGQ26EnPXr5JXYbfN8VlB30XoJ9UB3CAv9k76mLRicTQhRhaZk5vDR3Byv3n+X5ZuV5vW1l7GRVrnx1y0JPa/36nZ6cl2OEsDWZOZm8vfFtlp9YztNVnuaNem9gL8OhhMgzWXUzHx1YAj/2hWKloNuP4FsegKqlvPh5xylMJi1vrIQQBS4xJYN+M6LYFZfMR49E0LNhiNGRioR7nqOntf7K8nGEsG4XMy4yZM0Qos5G8UqdV+gd0Vv2eRHiLl1bdVP+6VjW1kmw7HUoVRuemQ/u15cor1TSk8sZ2Zy7nEFJLxl9IIQoOCcSrtB72lbiL6YzoVsd2lWV7U0Kyu2Gbnrm/rcyUA/4Jff7h4Gtlri4Umoq0BE4p7WuaolzCpFf4lPiGbhyIDGXY/i86ee0D21vdCQhCqXrc/Sk0rMIkwlWfQh/jobK7eGxKeD0zy0E/D2cAEi8IoWeEKLg7Dh5gX4zotBa8/2zDahTtrjRkYqU2w3d/BBAKbUeqK21vpz7/QfAEgtdfzowFphpofMJkS8OJB3ghZUvkJ6dzretvyUyMNLoSEIUWlpW3bSc7AxYPAj2LIC6/aD9FzfdPsHH3RmApCuZBZ1QCFFErYg+y0tz/8bf05kZfSIJ9fcwOlKRk5d99AKAG38zZOb+7L5prdcrpUIscS4h8sum05t4Ze0reDh6MOOhGVQsXtHoSEIUalfn6EmZd5/SkmF+dzixAVq9D02G3nI8rI+7eVU7KfSEEAVh9pYY3lu8l6qlvZjSqx7+ns5GRyqS8lLozQS2KqV+yv2+MzAj/yIJYT1+OfoL7//5PqHeoYxvNZ4Ad4t8xiFEkSY9ehZwMQ7mPAEJh6HLd1DjydsefrVHLzFFCj0hRP7RWjPyj4OMW3OUllVKMPaZWrg55aXcEPnhjn/zWutPlFLLgSa5P+qjtd6Rv7GuU0oNwLynn7e/v39BXVYUcVprJu2ZxDc7vqF+YH1GNR+Fp5PnnZ8oigxpm+6dSfbRuz9n98HsxyEzBbr/CKHN7/gUb1dH7JT06BUV0j4JI2Rmmxi2cDeLdpzi6cgyfPxIVRzsb7dlt8hvefrb11pvB+YCPwGJSqngfE31z2t/p7WuqLX2Dw4usMuKIizblM3HWz7mmx3f0CG0AxNaTZAiT/yHtE33Tsuqm/fu2DqY2s785z7L8lTkAdjZKYq7OZEohV6RIO2TKGiX0rPoM30ri3ac4tU2lfi0SzUp8qzAHXv0lFKdgC+BUsA5IBg4AETkbzQhCl5qViqvr3+ddXHr6F+tP4NrDZaVAYWwMM3VffTk39Zd2f0D/PwC+FU075HnVfqunu7j7kTSlYx8CieEKKrOXEyn97StHDmXwhePV+eJumWMjiRy5aXU/hhoABzSWpcDWgNbLHFxpdRcYDNQWSkVp5TqZ4nzCnEvEtMS6fd7Pzac2sA79d/h5dovS5EnRD4wyRy9u6M1bPgKFj0LwQ3MPXl3WeTB1UJPevSEEJZz8Mxluoz/k9ikVKb2ridFnpXJy+zILK11olLKTillp7Veo5QabYmLa62ftsR5hLhfMZdieH7F8ySkJTC6+WhaBLcwOpIQNkvm6N0FUw4sfQ2ipkDVx6HzeHC4t9XrfD2cOHDmsoUDCiGKqs1HExkwKwpXR3t+eL4hEaW8jI4k/iUvhV6yUsoDWA/MUUqdA67kbywhCs6u87t4cdWLKBRT2k6hun91oyMJYdNkjl4eZabCwv5wcAk0HmLeQsHu3ue8SI+eEMJSftl1mv/7YRfBvm5M71OPoOJuRkcSN5GXQu8RIA0YCnQDvICP8jOUEAVl9cnVvL7+dUq4lWBi64kEF5NJ60LkN61ljt4dXUmA75+EU9uh/UiIfPa+T+nj7kxyahbZOSZZJEEIcU+01kzacIxPlx4gMsSHST3r4uXmaHQscQu3LfSUUvbAb1rrFoAJ2T9P2JB5B+bx2dbPiPCN4JuW3+Dr6mt0JCGKhKtz9KTMu4WkYzD7Mbh0Gp6cDWEdLXJaX3cnAC6kZsnmxUKIu5Zj0nz8WzTTN52gQ7VAvuxaAxdHe6Njidu4baGntc5RSpmUUl5a64sFFUqI/GTSJr7++2um7p1K86DmfP7A57g5ypADIQqK9OjdRtx2+L4raBP0+hXKRFrs1D7XCr1MKfSEEHclPSuHl+ft4Pd9Z+nfpBxvtQ/Dzk7acGuXl6GbKcAepdQKbpibp7UenG+phMgnmTmZvPvnuyw9vpSulbryZv03cbDLyz8DIYSlyKqbt3BwGSzoAx4loPsi8Ktg0dNf7dHbcfICWTkmi567MKkU4ImjDF0VIs8uXMmk34xt7IhN5t2O4fRrUs7oSCKP8vIOd1HulxCF2qXMSwxdM5StZ7bycu2X6Ve1n2yfIIQBTFrGbv5H1FRY8ioE1oBnfjAXexZW0ssFgDcW7rH4uQuTrW+1okQxF6NjCFEoxCal0mvqVuKS0xj3TG3aVws0OpK4C3cs9LTWMi9PFHpnrpxh4MqBnLh4gk+bfMrD5R82OpIQRZ6M+sG8BOnq4bBhJFRsC09MAyf3fLlUqL8H8wc0IDktK1/OX1gUc5WFI4TIi91xyfSdvo2sHM2c/vWpF+JjdCRxl25Z6CmlvtNaD7jdk/NyjBBGO3ThEANXDuRK1hXGtx5Pw1INjY4kRJFmkjl6ZtmZ8Otg2DUXaveCDl+Bff4OJa8fKotOiSJGa/MqtklHIfkkXDkPKefMP7tyHtIuQHY6ZGdATob5v6YcsHcCByewdzb/19Ed3H3BzQ/c/cDdHzwCwCfU/OXsYejLzDHlcPrKaU5eOknMpRjOpp4lMS2RpPQkktKTSMtOIzMnk0xTJtmmbJzsnXC2d8bZ3hk3Bzf8XP3wdfXFz9WPku4lSbrgxcglSRR3Kc68AfWoUMLT0Ncn7s3tfqN0Vkql3+ZxBciu0sKq/RX/F0PWDMHNwY0Z7WZQ2aey0ZGEKPJMso8epF+CH3rAsbXQ8h1o+n9F/C9ECAtITYLTOyB+J5yNNhd3iUch49I/j7NzNA+PdvcD1+Lg5msu5hxcwMEZlB3kZOUWf5nmr8wrkHAErmyGtCTzgkk38ggAn/LmubWBNSCwJgREgKOrxV9mlimLQ0mH2Je4j+jEaPYl7uNo8lGyTNd76x3sHPBx8cHXxRcfVx9KOZTCyd4JJzsnHOwczEVfTibpOelcybrC8YvHiTobRXJG8rVz2JcDnLwYsaMKEX4RRPhGUM2vGoEeMnyzsLhdofdaHp6/wVJBhLC0JceW8M6f7xBSLIQJrSdQ0r2k0ZGEEFzfML3I9uhdOg1znoDzB6DzBKj5jNGJhCh8tDYXcTF/mr9Objb32F3lFQy+5aH6k+b/+lYA77LmAs/F6/4+WDHlmHsCL502b4eSdNT838RjsP83+Hum+ThlDwHhENwIyuZ+3cP8W5M2EZ0YzV/xf7Ht7DZ2nN1Banaq+WU6exHuE073sO6ULVb22pefq99dr0OgtWbkimjGb4iierlMOtZx4NSVGPYn7Wdm9EyyTdkAlHIvRb2S9YgMjKR+yfoEuAfc9WsSBeOWhZ7MzROFldaaqXunMvrv0dQrWY/RLUZTzKmY0bGEELmuDt0sknXeuf0w+3FITzYvulKhldGJhCg8Mi6be8EP/Q5HVsLlePPP3f3NRVTdflCqprlHzbV4/uWws88dvukHgdX/+ZjWcOkUxO+C0zsh9i9z4bf1W/PjAVWhYhuo+CAERd5yuHZadhp/xf/F2ti1rItbR0JaAgDlvcrzcPmHqRtQl6p+VSntUdoiC8tl5Zh4c9EeftwexxN1avLpo9X+sTptZk4mhy4cYtf5XUSdiWJt3FoWH10MQJhPGC2CW9CiTAsqF68sC91ZEVlXXtiUHFMOn239jPkH5/NQyEMMbzIcJ3sno2MJIW5wdR89VdSW3Ty+AeZ1Mw/l6rPsv28QhRD/lZYMB36DfT/B8fXmYZTOxaB8SwhtDmUbg19F6/nkSCnwCjJ/Velg/ll2prnwO7EBjqyCP8fAxlHmnsUqHSGiC4Q2J01nsy52HUuPL2Xz6c2k56Tj4ehB49KNaRbUjIalGuLn6mfxyCkZ2bww52/WHzrPy60qMqR1xf8Ua072TlT1q0pVv6p0C+uGSZs4dOEQf576k7Wxa5mwcwLjd44n0D2QlsEtaV+uPdX8qknRZzAp9ITNSMtO4431b7Amdg19IvowpM4Q7JTslSSEtbk+dNPYHAVqz4/w80AoXg66/wjewUYnEsJ6ZWfCwaWwax4cXWUu7ryDIXIAVGoHwQ3AvhCtnurgBGXqmb+avgLpF+HoGji4jOz9v7Ll4EKWeBVnlasLaeRQwq0Ej1V6jGZBzagbUBfHfHyt5y6l02f6Ng6cucznj1XjyXp5a5vslB1VfKpQxacK/ar1IyEtgQ1xG1gdu5oFBxcwZ/8cyniWoUNoBzqU60CIV0i+vQZxa3ku9JRSblrr1PwMI8S9upB+gRdXv8ie83sYFjmMbmHdjI4khLiFIrVhutaweSz88Y655+GpOfk7pEyIwizhCPw9A3Z+D6kJ4FnKXNxFPAqla1tPr939cvEitkwdFqaf4OesgySmJ+KJPe1TrtDh8kXqFHPHrow3eFfO14L2yLnL9Jq6jQupmUzuVZcWle99/04/Vz+6VOxCl4pduJx5mZUxK1lybAnf7vqWibsmUt2vOo9Veox2Ie1wc3Sz4KsQt3PHQk8p1QiYDHgAwUqpGsBzWusX8jucEHkReymWgasGcubKGUY1H0WrsjLnRQhrVmTm6Jly4Pe34K+J5qFZnSeCo2zULcQ/mHLMQzO3TjIPbVT2UPkhqNPbPDzTzt7ohBaTlZPFmtg1LDi0gC3xW7BX9jQNakrnCp1pWropTlnpsHehudhd/gaseA/CH4GGg8xzDy1o6/Eknp0ZhaO9HfMHNKRakJfFzu3p5Hmt6Dt75SzLji/jpyM/8f6m9/li2xd0DO3IE5WfoFLxSha7pri5vPTojQLaAr8AaK13KaUeyNdUQuTR3oS9DFo1iBydw+QHJ1OzhGUbQiGE5eV26Nn23I2sNFj0LOz/FRq+CG0+BjsZSi7ENZmpsHMObB4HF46bV8Rs9R7U7AaetrVKdnxKPPMPzuenIz+RlJ5EoHsgg2oOokuFLv9csdLeCer2MX/F7zYv4rJrHuz5Aco9AI0GQ4XW9/0p2ZLd8Qz9YSdBxV2Z0SeSMj7518MW4B5A76q96RXRi7/P/c2CQwtYdHgR8w7Oo4Z/DbpW7krbkLY42zvnW4aiLE9DN7XWsf/6hZyTP3GEyLt1set4bf1r+Lj4MKH1BMp5lTM6khAiD7TWtj0/LzUJ5j4FsVuh7WfQUAbACHFNahL89S1smwSpiVC6LrT50LwoiQ313gHsOr+LWdGzWBmzEo2mWVAzHq/0OI1LNcb+Tq81sDp0GAmt3oXtM2DLBJjzOPiHQaOXoHrXexrWOWXjcYYviaZ2cHEm96xLcfeCWbBOKUWdgDrUCajDsHrDWHx0MT8e+pG3N77Nl1Ff8lTlp+hauSu+rr4FkqeoyEuhF5s7fFMrpRyBl4H9+RtLiNv74eAPfPLXJ1TxqcK4VuPyZRUqIUT+MGltu715F07A7McgORaemA4RnY1OJIR1SLtg7r3bMgEyU6DSQ9B4MAQ3tKlx3NmmbFadXMXM6JnsPr8bT0dPeoT34Jkqz9zbRuMuXua/p/rPm1ce3fQNLH4B1n8BzV6Hal1vuUXDjUwmzSdL9zNl43HaRZRk9FM1cXE0prD2dvGmV0Qveob35K8zfzErehbjd41n8p7JdCzfkR5hPahQvIIh2WxNXgq954GvgdLAKeAPYFB+hhLiVrTWfLPjGybtmUTT0k0Z2WykTOoVopDR2kZX3Dy9w7wRek4W9FwMZRsanUgI46VfNBd3m8dBxiUI7wzNh0GJMKOTWdTlzMssOryI7/d/z+krpynjWYZhkcPoXKEz7o7u938BByeo8aS5J+/Q77DmE/NKvutHQrM3oNrjt+wRTc/K4dUfdrFkTzy9G4Xwbsdw7K2gEVZK0SCwAQ0CG3Ds4jHmRM/hl6O/sOjwIhqVakTP8J40KtXIdj8YLAB3LPS01gmALGEoDJeVk8UHmz/gl6O/8FjFx3inwTs42MkOIUIUNiZtg/PzDq+AH3qBmy/0Xgr+ssiAKOKyM8wLEW34CtKTzUMzm78JJasancyiEtMSmb1/NvMOzCMlK4W6AXV5I/INmgU1u/PwzHuhFFRuB5XawoElsHYE/DQANoyE1h9A5fb/6CFNTs3k2ZlRbDtxgbfbh9G/aTmrbH9DvUJ5t+G7vFTrJRYcWsDcA3N5fuXzVPCuQL9q/WgX0k7e892DW/6NKaW+4fqc+f/QWg/Ol0RC3ERKZgqvrH2FzfGbGVRzEM9Vf84qGyohxJ1prW1rq/S/Z8KvQyAgArotsLmFJIS4K1pD9GLzipHJMVChDbR8x+KrRhotPiWe6fums/DwQjJzMnkw5EH6VO1DhG9EwQRQCsI6mgu7A7/Cqo9h3jPmbVwe/BhK1yE2KZXe07YSm5TGmKdr0alGqYLJdh+8Xbx5tvqz9I7ozbITy5i2dxpvbniTsTvG0ieiD50rdpaFW+7C7UrjqAJLIcRtnEs9xwsrX+Bo8lE+bvwxnSvInBchCjONjeyhpzWs/QzWfQ7lW0HXGeDsaXQqIYxzegcsfwtOboIS4dB9EVSwrS2PTlw8wZS9U/jt6G8APFz+YfpU7WPcgnB2duYtGCq3N2/LsHYETGpJcvlOvBjTjvPZfszsF0mD0MK1yImjvSOdyneiY2hH1sWuY/KeyQz/azgTd0+kR3gPulbqioeTh9Exrd4tCz2t9YyCDCLEzRy5cISBqwZyKeMSY1uNpXHpxkZHEkLcJ5PJBlbdzMky9+LtnA01u8PDo/N1Y2MhrNqleFj1Eeyaax6+3HEU1OqZp0VCCosDSQeYvGcyf5z4A2d7Z56s8iS9wnvd2wIr+cHeEer1h+pPEvPLZ5TYO4kFaikpNfriExhpdLp7ZqfsaBHcguZlmrPtzDYm75nMqO2jmLxnMk9Vforu4d3xcfExOqbVysuG6b/y3yGcFzH3+H2rtU6/14srpdphXujFHpistR5xr+cStmfbmW28vOZlnO2dmd5uOmG+tjVxW4iiqtDP0cu4bJ6Pd3QVNBtmXliiML8eIe5VZqp5Fcg/R4Mp27w6ZNNXzStF2ogd53YwafckNpzagIejB/2r9adbWDer3QZgwZ5k3tzRhAb+9ZgU9Ds+uybBoQXm+ZF1+xbaD6SUUkQGRhIZGMm+hH1M3jOZyXsmMyt6Fo9VeozeEb0p6S7D5v8tLx+1HAP8gbm53z8JXAYqAZOAHvdyYaWUPTAOaAPEAduUUr9oraPv5XzCtiw/vpy3Nr5FkGcQE1tPpJSH9Y8rF0LkjUYX3rro8hnzyppn90Gnb6B2T6MTCVHwTCbYswBWfQiXTpmHDrb+EHxsYz9brTWbTm9i0p5JbD+7neLOxRlcazBPVnmSYk7FjI53U1prvll9hK9WHKJJBT8mdK+Nq0sXaDoIfn8Llr0O2ybDg8Oh4oOF+sOpCL8IRrUYxbGLx5i6ZyrzD8xn/oH5PFz+YfpW7UuIV4jREa1GXgq9Rlrrejd8/6tSapvWup5Sat99XDsSOKK1PgaglJoHPAJIoVeEaa2ZGT2TkVEjqV2iNmNajsHL2XY+GRRCXN1eoRC+yTh/EGY/bt7k+Zn5ULGN0YmEKHgnt8DyN+H03xBYEx6bDGUbGZ3KIkzaxKqTq5i0exL7k/YT4BbAsMhhPFrxUVwdXI2Od0vZOSbe+Xkv87bF8mit0ox4rDpODnbmB0tWg56/wKHl8Mc78H1XCG0BbT8xLyBViIV6hTK8yXAG1Rx0bWGcn4/8TJuybehfrb+MBCNvhZ6HUipYa30SQCkVDFyd/Zh5H9cuDcTe8H0cUP/fBymlBgCvAd7+/v73cTlh7XJMOYyMGsns/bNpU7YNnzX9TFZWElZL2qZ7Z94w3egUdylmM8x9CuydoM8SKFXL6ERC3FK+tE8XYmDl++ZNuz0DofNEqP6keTGQQi7LlMWy48uYvGcyxy8ep2yxsnzU6CM6hnbE0cqHOl7JyObF7/9mzcHzvNiiAq8+WOm/Q+OVgsoPmReNippiXrBlYhPziIQWb4NHCWPCW0igRyBv1n+TAdUHXNvq4o+YP2hSugnPVnuW2gG1jY5omLwUeq8CG5VSRwEFlANeUEq5A/m+YIvW+jvgO4C6devecrsHUbilZ6fz1sa3WBGzgh7hPfi/uv+HnSr8vzyE7ZK26d4Vuh69fT/DogHgHQzdf4TiIUYnEuK2LNo+pV+CjV/B5vGg7MzzUhsPBicLbAJusPTsdH468hPT907n9JXTVC5emS+afUGb4Db5sweehZ2/nEHf6dvYd/oin3SpSrf6ZW//BAcnaDDQXKCv+9w8lHPPQmj6CjR4ARxdCiZ4PvF19eXl2i/Tp2of5h+Yz6zoWfRa3ovaJWrzbPVnaVyqceGeH34P8rJh+lKlVEWgSu6PDt6wAMvo+7j2KaDMDd8H5f5MFDHJ6ckMXjOYned28lrd1+gZIXNehLBlJl2IVt3cPN48v6VMJDw9D9xkdTdRRORk5y7X/xlcOQ81noaW74JXaaOT3bfLmZeZf9BcCCSlJ1HTvyZvN3ibpqWbFppC4Nj5FHpN20rC5Uwm9axLq7CAvD/ZzQce+ty8Sucf75rnWm6fZp5nGdGlUM/fAyjmVIxnqz9L9/DuLDq8iGl7pzFw5UDCfMLoX60/rYJbFYpC3hLyuu5tHSAk9/gaSim01jPv89rbgIpKqXKYC7yngGfu85yikIm7HMfAlQM5nXKaL5p9QduQtkZHEkLkM5MGrH3LdJPJPJ9lyzgIexgenQSO1jtHRwiL0RoO/Q4r3oWEQ+YNuJ/5AUoX/uFviWmJzNk/h7kH5pKSlULjUo3pX60/dQLqFJoCD2B7TBL9Z0RhpxRzBzSgZhnvezuRX0V4Zh4cWwu/vw0/9oG/JkLbzyCojkUzG8HVwZVuYd3oWqkrvx37jSl7p/DqulcJKRZCv2r96BDaAUc76x6ae7/ysr3CLKA8sBPIyf2xBu6r0NNaZyulXgR+x7y9wlSt9f0s7iIKmX2J+xi0chBZpiy+e/A76gQU/kZFCJEXVt6jl5UOPz0H0T9D/eeh7adQRD79FUXc6Z3mDzhObADfCvDUXPPcrkJUBN1MfEo80/dNZ9HhRWTkZNC6bGv6V+tPuG+40dHu2vK9Z3h53g4CvVyY0TeSsr4WGEIb2hyeWw87ZsPq4TC5JVTrCq3fB6+g+z+/wRztHelSsQudyndixckVTN49mXf/fJfxO8fTO6I3j1Z8FBcHKxi2qrXF/63lpUevLhCutbb4HBSt9VJgqaXPK6zfhrgNvLruVYo7F2dq26mEeocaHUkIUUBMJiueo5eaBPO6wclN5mXIG75Y6N/kCnFHWsPiQbBzjnnD8/YjoU7vQrvn2lWHLhxi5r6ZLDm2BICO5TvSt2pfynkVzm0gZmw6wQe/7qNGkDdTetXF18OCC9bZ2UOdXlD1Udg4CjaNhf2/QqOXzHMynT0tdy2D2NvZ0y6kHW3LtmXjqY1M2jOJz7Z+xre7v6VHeA+erPwknk4GvM6TW8wFdq3uUOMpi546L4XeXqAkEG/RK4sia9HhRXy0+SMqFa/EuFbj8HeTFQuFKEqsZdXNjJwMktKSSEpPIiUrhdTkGNLWjSAtNYHUps+S5uFM5s5xmLQJjb7236ufezraOeJk74SzvTNO9k7X/uzh6EExp2J4Onni5exFMadiuJTyVHEAACAASURBVDm6GfxqhbgNpcDJAxoPMS/MUYg3PL+6B96MfTPYHL8ZF3sXulbuSu+I3gR6BBod756YTPr/27vv8KiqtY3Dv5VGQqgp1AABpXcIiGgABRXsHLGiUkQUFRHbsXusx4L6WbHTBRF7wUITFBBCh1CkE0ogCSRACCmzvj8mcFCBBJLJnvLc15WLTDLZ+0kgL/PutfZavPjjGt6bvZEeTavz5vVtiQjz0CyDchWh+xPuRn/af2D2S+5FWzoPhY6DoVyFoo7g9YwxJMYlkhiXyKLURXyw/ANeX/w6H6/4mKsaXcX1Ta4vm/2bty+CGc/BhukQWQ1P3NJQnEYvBkg2xiwADh/5oLX28lJPI37NWsvIZSMZuWwk59Q6h1e6vUJkqO+v2iUip8bi+RE9ay3pOelsydpCyv4UUg6ksH3/drYf2M6eQ3vIyMngYN7Bf35hBBBRGVJ+cr8BwSYYYwwGQ5AJIsgEYa0l15WLy7qKlScyNJJq5atRvXz1o3/WrVSX+Erx1K9cX/uFivMufsnpBCVyuOAw32/8nnHJ41i/bz2xEbEMazeMqxtd7dO/X4fzC7j/s+V8u2wHN3aqy1OXtyC4LOa+V6kLfT6GTne6F+SZ/hTMews63w0db/WLVVcB2ldvT/sL2rMqfRWjVo5iXPI4xiWPo3vd7tzU7CZax7Yu/fs3dy53/0zX/gARUXDB0+6FcTzwMy1Oo/efUj+rBJw8Vx5Pz3uar9Z/xZVnXskTZz/h9zfAisjxeWJEL/NwJkm7klietpy1GWtZnbGajJyMo583GKpHVqd2hdq0iGlBdHg0UeFRRIVHUXXvNirOeY3yoRWIuOx1yldvQURIBOVDyhe5h1a+K5/cglz3myuXw/mH2Z+3n6zcLLIOZ5GVm0Xm4Uz2HNrD7uzdpB5M5Y+df5B2KI0CW3D0ONHh0TSLbkbLmJY0j2lO69jWPv3iVKSspB9KZ/K6yUxaM4mMnAwaV23M8+c+T8/4nl6/B15RMg/lcdu4JOZvzODfPZtwe9cGZb9oTFx797Yy2xa6m5NpT8LcN6HT7ZBwi9+sRNw8ujkjuo5g54GdTFw7kSnrpvDzlp9pEd2CG5vdyIXxF5bsdau17ntff38d1k9zj5qf/5j7PnAPTostzvYKv3rs7BIQDuYd5L5Z9/H7jt8Z0noIQ1oP8anVrUSkdJXG/eYFrgIW717Mr9t+ZcGuBazJWIPFEhIUQsMqDekS14XGVRtTv3J94irGUTOyJmHBYf880JIJ8NMzENsU+n4GlU5taldIUAghQSGnPDUzz5XH9v3b2Zy1mc2Zm/lz358kpyfz2/bfsFiCTBAtYlpwbu1zSaydSLPoZtpbVKSQtZak1CQmr53MtK3TyHfl0yWuCzc3u5mONTr6xWuMHfsO0X/UAjalHeT/rm3DlW0d3taiTge46QvYtsC9B9+MZ2HOq9DmBvcefNFnOJuvlNSsUJN729/L7a1u55sN3zBh9QQemvMQrya9yjWNr6F3w95UK38KG8wX5MPqb9wN3s6l7ima5z/uHsGLOM3VUk+BKWqNFWNMJ+BNoCkQhnuFzIPW2koeT/c3CQkJNikpqaxPKyWQdiiNO6bdwbq963i80+Nc1egqpyOJlzHGLLLWJjidoyRUm07NsElLWLptH78+cN4pf+2mzE1MXjuZHzf/SNqhNMKCwmhTrQ0danSgY42OtIxpWbyr+NbC7Jdh5nPuFeeuGQfhZf7f2j8czDtIcnoyC3ct5Lftv7EybSUWS63IWlx+5uVcfsbl1KlYp+gDSalQffIumYcz+XbDt0xeN5lNmZuoFFaJK868gqsbXe2zC6wcz+qdWfQftYDswwW8d1N7Op8Z43Skf0pNhnlvw4rJUJAHTS5xT+mM7wJB/nNRymVd/Lb9N8Ylj2P+zvkEm2AS4xLp07AP59Y+98T78e1PhaXjYdFo2LcVos5wL2rT6rpS2Zi+uLWpOI1eEu497j7DvQLnzUAja+3DJU55ivypWAWCjZkbuWPaHWTkZDCi6wi6xHVxOpJ4Ib2QCjx3T1zCiu2ZzLy/W7G/ZlHqIt5f/j5zd8wlJCiErnFd6Rnfky5xXU59oZOCfPj+Xvdm0K2ug8vfhJDjjPZ5gYycDH7f/jvfbfyOeTvmAXBBvQu4tdWtNIlq4nA6/6f65LwCVwELdi3gu43f8fPmn8kpyKFVbCuuaXQNF8Vf5B3L4pei39encdu4RVQoF8KoAR1oWtP5C1AntT8VFrwPSR/Bob1QpZ579cg2N/jF1gzH2pq1lc///Jyv139Nek461ctXp3fD3vQ+s7d78RaXCzbOdDd3a38AVz7EJ7oXsWlySalu01OqjZ61NsEYs9xa26rwY0ustW1LKWux+XqxCiSLUxczdMZQQoJCeKf7OzSPae50JPFSeiEVeO76ZDHJO7OYcV+3Ip+7KXMTry16jZnbZhITEcO1ja+lT6M+xESc5hXuwwfcmwL/+TMk3u++R8JHpnntOriLSWsmMWntJA7mHaRXfC8e7Pjg6f8spEiqT85Zm7GW7zZ+xw8bf2D3od1UDK1Iz/o9uabxNX57keOLxSk8OGU5Z8RWYPTADtSsHOF0pOLLy4E137kvoG2aDSYIzugOra6Fxj39YnuGI/Jcefy67Vem/DmFudvnYrG0C4uh1740LtyTQlR4VXej264/xJzpkQzFrU3FWYwl2xgTBiw1xryEe5sF/xmTlVL3y5ZfeGj2Q9SqUIt3eryjaUYi8hfWFr2ItLWWiWsm8krSK4QGh3J327u5sdmNRISU4IXPgd3wyTWwcxlc+hokDDz9YzmgRmQN7ml/DwNaDGBs8lhGrRzFbzt+4/6E++l9Zm+/uC9JApe1luT0ZKZvnc6MrTPYkLmBEBPCuXHn8u8G/6Zrna6UCy7FfeO8iLWWd2Zt4OWf1nJ2g2jevak9lSN8bCGZ0HBo2cf9lrHJvSfjkgnwxSAILgcNL4BmV8KZ3X1+AZdQgukREkWP4LpsP5DEd3mpTI3M47nyobxQrw6danbiovrt6RJZlWiHsxan0bsJd2N3FzAcqAPoRis5rvHJ43lp4Uu0jm3Nm+e/SZVwz99oKiK+xb3YyImbkszDmTz+++PM3DaTxNqJPH3O0yUftUpbD+P/BQf3wHUT3VeYfVTlcpUZ2nYolzS4hKfnPc2Tc5/kj51/8FTnp/xuGpv4t0P5h1iSuoRZKbOYsXUGqdmpBJtg2ldvz7VNrqVnfE+qhld1OqZH5Re4ePKbVUz4YytXtKnFS31aUS7EQ3vklZWo+u7ZEt0egW1/QPJXkPy1e8TPBEHtBHfjd0Z3qNkagovTjjhsfyps+Q02zHTPCDmQCiaI2nU7c1uLIQxuejnrcjOYumkqUzdN5Ym5T2AwtIxtSde4rkcXCCvrC3JFTt30Jr46/SAQuKyLV5JeYWzyWLrX7c4LiS/oBYcUi6ZGBZ7bxy1iU9pBfhr+z/t292TvYfAvg9mctZnh7YZzU7ObSv4f49Y/YOJ17hcYN0x2LxfuJ1zWxYcrPuStJW/Rtlpb3ur+FhXD/GeKlNNUn0pXdl42q9JXsTh1MfN3zmfZnmXkufIIDw6nc63OnF/3fLrGdQ2Yi8TZufncPXEJ01bvZki3M3jgwsYElcUeeU5wudwbhK//Bf78BXYsASyERrpX9KzbGeIS3I1fpMPT0QvyYc8ad8btSbD5d0j/0/25cpXdo5KNerqb1eOMTlprSc5IZnbKbGZvm83K9JUAVIuoRvsa7UmonkD76u1pUPn0t8sozambvmfqQ7BrhdMpAsZhLI+adH4y2dxgK/Dg5nUEb9agb0Cr0RJ6veB0CvFSJ9pHL/NwJoN+HsTOgzsZ2WMknWp2KvnJVn8Lnw+CSrWg7xS/WQL8iCATxOBWg6lbsS4Pz3mYO6ffyfsXvK8LbXLKFu5aSJAJok7FOsRGxJboAou1ll0Hd7ExcyMb9m1g/b71rEhbwcbMjbisC4OhSVQT+jbtS6eanWhXvV3JpmX7oLQDh7llTBLLU/bx9BXNufnseKcjeVZQkLuhq9MBznsEDqa57+XbOg+2zHPv0Ufh4FPFWlCzFcQ0co8ORjWAqvXdi7uU4oIm5B50TzNNX/+/t7Q/IXUV5B9yP6dcJajbCdrdBPHnQo2iRyCNMTSPbk7z6OYMaT2EtENpzEmZw9wdc0nalcTUTVPd32ZoRRpHNaZJVBOaRjelbbW2pX67k382elJmMnExzOxhkTnMfa4q9KMipsi7b0QkkFn4x4vIPFced8+4m237t/HeBe/RoUaHkp/oj/dh6oNQuz3c8KnzV4k9qGf9nhhjeODXB3jkt0d4pesrumdPTsnLC19mdcZqAMKDw4mrGEdMRAxVy1WlSngVKperTGhQKMEmmJAg98vHQ/mHyM7P5lDeIfYd3sfu7N3sObSH3dm7OVxw+Oixq5SrQvOY5lxQ7wJaxLSgVUyrgBm1O57NaQfpN2oBuzJzePfG9lzUvIbTkcpeZAy0+Jf7DSAn033/9M5lsHO5e8BmwwwoyP3f1wSFQoVq7lG08jFQPtr9FlYegsMgONT9Z1AouPLcC8TkF77lHYLsdPe92gd3u//MPfDXTBVruS8GJgyEWm3db1ENSrxdRExEjHt1zoa9sdaSsj+FpNQkVqatZE3GGqasm0JOQQ6DWw1maNuhJTrX35200TPGBAMvWmvvL9WzeppGEsrEjgM7GDJtCNv2u3jp3JfoVb+X05FExAdYa/n77KS3l7zN4t2LeSHxhZI3eS4XTP+Pe4PaxhfDVR+5Xwj4uYviL2LngZ28sugVJq6ZyA1Nb3A6kviQV7q+wpb9W9i2fxsp+1PYtn8bGTkZ7Dy4k705e8nKzTru14UGhRIREkHlcpWpVr4aLWJaUC2iGnUr1aV+5fqcUeUMosJ9e/GN0rRk615uGZOEtZZPbu1E+3r+fQ9isYVXhvpd3G9HuFywfwdkbHSPvO3d5G7QstPdI4J7N0N2hnv07diG8FhBIRAS4V4sJiLK3SjWbAMVqkOFWKgaD9Fnuve5K1fB49+mMYY6lepQp1IdejfsDUC+K58tWVsoH1L6/0+dtNGz1hYYY84t9bOKz1uTsYY7pt1BTn5O6V19F5GA4LJ/3dFg6e6lfLzyY65qeBWXNLikZAfPPwxf3QErp0DCLXDxy6U71cfL9Wvej4WpCxmRNIJzap9DvUr1nI4kPuLIi88TcVkXBa4C8m0+Ba4CLJaIkIijo3tStF+SUxk6cTHVKoYzekAHGsR6vrHwaUFB7umaleP+2gAej7XufesKct1vQaEQEu4TC72EBIVwRhXP3FZQnLHIJcaYb4wxNxlj/nXkzSNpxCfM3T6XflP7ERwUzNheY9XkicgpcY/ouTs9l3Xx3wX/JTYilgc7PFiyAx/aB+Ovcjd5Pf4Dl7wSUE0euK8WP9X5KcKCwxixcITTccSPBJkgQoPdo3cVwipQMayimrxTMH7+Fm4bl0Tj6hX54o7OavJKmzHuqZthkRBR1T065wNNnqcVp9ELB9KB84HLCt8u9WQo8V5fr/+aO6ffSVzFOMb3Gs+ZVT2zEaSI+C/3iJ670ft+4/ckpydzb8K9lA8twbSVzBT4uCdsnQ//+gDOHe4zG6GXtpiIGG5rdRuzUmYxf+d8p+OIBDRrLS/9uIbHvlrJeY2rMXFwJ2Iq+Od+gOJ9imx1rbUDyiKIeDdrLe8vf5+3lr7FWTXP4v+6/R8VwnQ1SkROnctaDO66MmrVKBpWbcjF9S8+/QPuWgkT+rhXULtxCjToVkpJfVffpn0ZlzyOUStHlc7qpSJyynLzXfz78+V8uWQ713esyzNXNCckuGQLe4iciiL/tRljGhljphtjVhY+bmWMeczz0cRb5LvyeXr+07y19C0ubXApI7uPVJMnIiUSZGDeznn8ufdP+jXrd/orRG6c5R7Jw8CAqWryCoUFh3F9k+uZu2Muf+790+k4IgEnKyePAaMX8OWS7dx/YSOe791CTZ6UueL8i/sAeBjIA7DWLgeu82Qo8R7ZedkMmzmMKeumcGvLW3n+3OcJDQ51OpaI+DBX4T16k9ZMIiYi5vRX7F32KYzvA1XqwKBpUKNF6Qb1cVc3uprw4HA+Xfup01FEAsquzByueXcef2zMYMTVrbnr/Iba7kQcUZxGr7y1dsHfPpbviTDiXdIOpTHwp4H8tv03Hu/0OHe3u1uFSkRKzOUCl8nh9+2/0zO+J2HBYad2AGthzivw5WD3RrYDpkLl2p4J68OqhFehS1wXpm2ZRoGrwOk4IgFh7a799H7nd1L2HmLUgA70aR/ndCQJYMVp9NKMMWdQuF29MaYPsNOjqcRxmzM3c9MPN7Fh3wZeP+91rml8jdORRMRPWCzZwSvJdeXSo16PU/vignz4/l6Y/jS06AM3fg4RgbvxclEuqHcB6TnpLNuzzOkoIn5v3oZ0+rw7lwKX5dPbOpHYMNbpSBLgirPu6J3A+0ATY8x2YBNwo0dTiaOW7l7K0BlDCTJBfHzRx7SMbel0JBHxIy4L+0OWEBUeRZvYNsX/wtyDMOUWWDcVzrkHuj/p3mdJTigxLpHQoFCmbZ1Gu+rtnI4j4re+XrqdBz5bTr3o8owe2JHaVSKcjiRS9IietXajtbYHEAs0sdaea63dXJKTGmOuNsasMsa4jDEJJTmWlK7pW6cz6OdBVAqrxLhe49TkiUipc9kCDgSt4rw65xFc3H3uDuyBMZfBuh/h4hFwwVNq8oohMjSSs2udzeyU2U5HEfFL1lre+3UDwyYtpU3dKky5vbOaPPEaJxzRM8bce4KPA2CtfbUE510J/At4rwTHkFI2cc1EXljwAs2jm/NW97eICo9yOpKI+KHDpOIyh2hbrW3xviB9g3sj9P074drx0FRbuZ6KttXaMjtlNlm5WVQKq+R0HBG/UeCyPP3tKsbM28IlrWryytWtCQ8t5sUrkTJwsqmbFT11UmvtakALe3gJl3Xx+uLX+Xjlx3SL68ZLXV8iIkRXo0TEMw4FbQageXTzop+ckgSfXONegKXft1Cno2fD+aFm0c0ASE5P1p56IqUkJ6+AYZOW8NOqVG5NrM/DvZoSFKTXteJdTtjoWWufKssg4ozcglwe//1xftj0A9c2vpaHOj5ESFBxbt0UETk9OWYLQTaM+pXrn/yJa36AKQOhQjW48QuIObNsAvqZIw31qrRVavRESkHGwVwGjVnIkm37eOLSZgw8t4haJuKQIl/RG2NGUbji5rGstQOL+LppQI3jfOpRa+3XxQ1ojBkMPABUiY3V6kWlKSs3i+Ezh7Ng1wKGtRvGLS1u0SirSDGpNp2+nOAtRFDv5PfnLfwIfrgfaraGGya7mz05LZXLVaZ2hdqsSl/ldBQpI6pPnrM1PZv+oxaQsu8Q79zQjl4tazodSeSEijN0890x74cDvYEdRX1R4QIuJWatfR/3qp8kJCT8o+GU07Pr4C6GTBvC5qzNPH/u81x2xmVORxLxKapNpyfflU9u0DYq0+34T7DWvXXCb69Cw4vg6lEQFlmmGf1R8+jmavQCiOqTZyxP2cfA0QvJK7B8MugsEuK1loF4tyIbPWvt58c+NsZMBH7zWCLxuHV71zFk2hCy87IZ2WOkpvKISJnZkrUFa/KItPX++cn8XPjmLlj+KbTrB5e8CsGaSl4amkU34+ctP2tBFpHTNHPNbu6YsJioyDAmDe7ImdUqOB1JpEinszZ1Q6BEc2iMMb2NMSnA2cD3xpifSnI8Kb4/dv5Bv6n9ABjdc7SaPBEpU7uzdwMQbqL/+omcTJjQx93knf8YXPa6mrxSVDPSPb0sLTvN4SQivmfigq0MGpvEGdUi+fLOzmryxGcU5x69/bjv0TOFf+4C/l2Sk1prvwS+LMkx5NR9t/E7Hv/9ceIrxTOyx0hqRB7vFkoREc/Zm7MXgDCOGVXK2gETroY9a+DKkdDmBofS+a+q4VUByMjJoAENHE4j4hustbz2yzremLGero1iebtvOyqU0wUo8R3FmbrpsW0WpGxYa/lo5Ue8vvh1OtTowP+d93+auiMijsjIyQAg9Eijl5rsHsnLyXQvunJmdwfT+a8j+6Ie+fmLyMnlFbh4+IsVTFmUwjUJcTzXuyWhwaczEU7EOcUZ0TNAX6C+tfYZY0xdoIa1doHH00mJFbgK+O+C//Lp2k/pVb8Xz57zLGHBYU7HEpEAlZGTATaIsKBI2DQHJvWF0AgYMBVqtnI6nt860ugdGVEVkRM7cDifIeMXMefPNO7p0ZBh3RtqVXLxScUZf34HcAHnA88A+4HPgQ4ezCWl4FD+IR6c/SCzts1iQIsB3NPuHoKMrkaJiHMycjIIspGclT0Txr8MVevDjVOgSl2no/m1KuFVAI3oiRRld1YO/UctZG3qfl66qhXXdKjjdCSR01acRu8sa207Y8wSAGvtXmOMhoS8XEZOBkOnD2VF2goe7vgwNzTVPS8i4ry9OXupUuBi4K7noN45cN0EiKjqdCy/FxoUSqWwSmr0RE5i/e799Pt4IXuzc/mwXwLnNdb+neLbitPo5RljgincNN0YE4t7hE+81Lasbdw+7XZSs1N5rdtrdK+ne15EvEJ2Bvz5s9MpHJWR8gcN89NZVKEb7W/8FELDnY4UMKLCo9ToiZzAgk0ZDBqzkLCQYD4dfDYt4yo7HUmkxIrT6L2Be4XMasaY54A+wGMeTSWnbcWeFdw14y5c1sWHF35Im2ptnI4kIkfs2wJf3uZ0CkftjatFmK3L2FpP0F5NXpmKCo9i72Hdoyfyd98v38nwyUuJqxrBmAEdqRNV3ulIIqWiOKtuTjDGLAK6495i4Upr7WqPJ5NTNmvbLB749QGiI6J5t8e7xFeOdzqSiByrWjO4e4nTKRyV8WNfcjObUCdI9wuXtajwKDZlbnI6hohX+XDORp77YTXt61blg5sTqBqpu5PEfxRrMxBr7RpgjYezSAlMXjuZ5/54jqZRTXmr+1vERMQ4HUlE/i6kHEQF7h5meQV57M87QLirAkFawa7MVQ2vyuLdi52OIeIVXC7Ls9+v5uPfN9GrRQ1eu7YN4aHBTscSKVXa9dHHWWt5c8mbfLDiAxJrJzKi6wjKh2rKgYh4n6P3hxVUQH1e2YsKj2Jvzl4KXAUEB+kFrQSunLwC7p28lB9W7KJ/53gev7QZwUEqSuJ/1Oj5sLyCPJ6c+yTfbvyWqxpexWOdHiMkSH+lIuKdjtwfZlwVMOhFVVmrGl4ViyUzN/PovnoigWZfdi63jk1i4ea9PHZJU245t772yBO/pa7ARx3IPcDwWcOZv3M+d7W5i8GtBqtQiYhXyzj0vxE9XTwve9Hh0YD770GNngSibRnZ9B+1gG0Zh3jz+rZc1rqW05FEPEqNng9KPZjKHdPvYOO+jTxzzjNceeaVTkcSkWI4mHeQ5PRkp2M45uj9YQWRukfPAVXD3fsVzt85P6BX32wd25qwYC24EWhWbs9kwOiFHM4rYNwtHTmrQbTTkUQ8To2ej1m/dz1Dpg8h63AWb3d/m861OzsdSUSKaXPWZgb+NNDpGI4KCQrB5lfSPXoOqBXpHr14ceGLDidx1oyrZxBbPtbpGFKGfl23hzvGL6JK+TA+GXQWDatXdDqSSJlQo+dDFu5ayLAZwwgPCWdMrzE0iWridCQROQXxleL5+KKPnY7hqOiIaK59a6OmmjugTqU6TLlsClm5WU5HcVTlctoIO5BMTtrGw1+soFH1iowe0IHqlbR/pwQONXo+YuqmqTz626PUqViHkT1GUquC5pWL+JrI0Eg61OjgdAzHWbtB9+g5pHFUY6cjiJQJay1vTF/Pa9PWkdgwhnf6tqNieKjTsUTKlBo9L2etZWzyWEYkjaBdtXa8cf4buhopIj7NZdHUTRHxmLwCF499uZJPk7bxr3a1eeFfrQgLCXI6lkiZU6PnxQpcBbyc9DITVk/gwnoX8nzi85QLLud0LBGRErHWajEWEfGIg4fzufOTxcxau4eh55/JvRc00lRxCVhq9LxUTn4OD895mGlbp3Fzs5u5L+E+goyuRomI73NZ1OiJSKnbs/8wA0cvZNWOTJ7v3ZIbzqrrdCQRR6nR80L7cvYxdMZQlu1ZxoMdHuSmZjc5HUlEpNS4rHU6goj4mQ17DtB/1ALS9ufywc0JdG9a3elIIo5To+dlUvanMGTaEHYc2MGIriO4MP5CpyOJh+Tl5ZGSkkJOTo7TUcpEeHg4cXFxhIbqZviApxE9ESlFi7ZkcMuYJIKNYdLgTrSuU8XpSCJeQY2eF1mVvoo7pt1BviufDy78gHbV2zkdSTwoJSWFihUrEh8f7/f3D1hrSU9PJyUlhfr16zsdRxzmslarbopIqfhx5U6GTVpKrSoRjB7QgXrRkU5HEvEauunLS8xJmcOAHwcQHhzOuF7j1OQFgJycHKKjo/2+yQMwxhAdHR0wo5dyclp1U0RKw+jfNzFkwmKa1arElNvPVpMn8jca0fMCX/z5BU/Pe5pGVRvxdve3iS0f63QkKSOB0OQdEUjfq5ycRatuisjpc7ksL/y4hvdnb+SCZtV547q2RIQFOx1LxOs40ugZY14GLgNygQ3AAGvtPieyOMlayzvL3uHdZe9yTq1zeKXbK0SG6mqUiPg394ieGj0ROXWH8wu4/7PlfLtsBzd1qsd/Lm9OsOaCixyXU1M3fwFaWGtbAeuAhx3K4Zg8Vx5PzH2Cd5e9y5VnXsmb3d9UkyciAcFaq6mbInLKMg/lcfNHC/h22Q4e6tWEp69QkydyMo6M6Flrfz7m4XygjxM5nHIw7yD3zbqP33f8zpDWQxjSeoiubotIwLAWLcYiIqdk+75DDBi1gE1pB3n9ujZc0aa205FEvJ433KM3EPj0RJ80xgwGHgCqxMb6/r1re7L3cOf0O1m3dx1PdX6KfzX8l9ORxAs89e0qkndkHS7ZagAAF65JREFUleoxm9WqxJOXNT/pc5544gmioqK45557AHj00UepVq0a33zzDXv37iUvL49nn32WK664gpdffply5cpx9913M3z4cJYtW8aMGTOYMWMGH330ERMmTCjV/N7O32pTWXKvuqlOT8RT/K0+Je/IYsDoBWQfLmDMwI50PiPG6UgiPsFjUzeNMdOMMSuP83bFMc95FMgHTvgK0Vr7vrW2obU2tm7dup6KWyY27tvIjT/cyOaszbx5/ptq8sRxAwcOZOzYsQC4XC4mTZrEddddx5dffsnixYuZOXMm9913H9ZaEhMTmTNnDgBJSUkcOHCAvLw85syZQ5cuXZz8NhzhT7WprLksqM0T8Rx/qk+//ZnGNe/Nw2D4bMjZavJEToHHRvSstT1O9nljTH/gUqC7tdZ6Koe3WJy6mKEzhhIaFMqonqNoHn3ykRYJLEWNvHlKfHw80dHRLFmyhNTUVNq2bUtUVBTDhw9n9uzZBAUFsX37dlJTU2nfvj2LFi0iKyuLcuXK0a5dO5KSkpgzZw5vvPGGI/nF9xwp95quLiJF+WJxCg9OWc4ZsRUYPbADNStHOB1JxKc4tepmT+BBoKu1NtuJDGXp580/8/Cch6lVoRYje4wkrmKc05FEjho0aBCjR49m165dDBw4kAkTJrBnzx4WLVpEaGgo8fHx5OTkEBoaSv369Rk9ejSdO3emVatWzJw5k/Xr19O0aVOnvw3xEUcu62nqpoiciLWWd2Zt4OWf1nJ2g2jeu7k9lcJDnY4l4nOcWnXzLaAi8IsxZqkx5l2HcnjcuORx3P/r/TSLbsa4XuPU5InX6d27Nz/++CMLFy7koosuIjMzk2rVqhEaGsrMmTPZsmXL0ecmJiYyYsQIunTpQmJiIu+++y5t27bV6IwUm6uw09NiLCJyPPkFLh77aiUv/7SWK9vUYszAjmryRE6TU6tununEecuSy7oYkTSCccnj6FG3B/9N/C/hIeFOxxL5h7CwMM477zyqVKlCcHAwffv25bLLLqNly5YkJCTQpEmTo89NTEzkueee4+yzzyYyMpLw8HASExMdTC++xlU4oqdrAyLyd9m5+dw9cQnTVu9mSLczeODCxgTpqpDIafOGVTf9zuGCwzwy5xF+3vIzfZv25YGEBwgOCnY6lshxuVwu5s+fz2effQZATEwM8+bNO+5zu3fvTl5e3tHH69atK5OM4j8sukdPRP4p7cBhbhmTxIqUfTxzRXNuOjve6UgiPk+NXinLPJzJ3TPuZvHuxdyfcD83N7tZL2jEayUnJ3PppZfSu3dvGjZs6HQcCQC6R09E/m5z2kH6jVpAalYO797Yngub13A6kohfUKNXinYc2MGQaUPYtn8bL3V5iV71ezkdSeSkmjVrxsaNG52OIQHEdXTVTYeDiIhXWLJ1L7eMSQLgk1s70a5uVYcTifgPNXqlZHX6au6YfgeHCw7z3gXv0aFGB6cjiYh4nf+N6DmbQ0Sc90tyKkMnLqZ6pXBGD+hI/ZhIpyOJ+BWnVt30K3O3z6X/j/0JCQphbM+xavJERE7gf6tuqtMTCWTj5m/htnFJNK5ekc+HdFaTJ+IBGtEroa/Wf8VTc5/ijCpn8E6Pd6hWvprTkUREvNaRVTdFJDC5XJaXf17LyFkb6N6kGm/e0JbyYXo5KuIJ+s06TdZa3lv+Hm8vfZtONTvxWrfXqBBWwelYIiLeTYuxiASs3HwXD05ZxldLd3DDWXV5+vLmhARrcpmIp+i36zTku/J5at5TvL30bS5rcBnvdH9HTZ74lUGDBpGcnAzA888/X+Tz+/fvz5QpUzwdS/yANkwXCUxZOXn0H7WAr5bu4IGLGvPclS3U5Il4mEb0TlF2Xjb3/3o/c7bP4daWtzK07VBtnyB+58MPPzz6/vPPP88jjzziYBr/sS87l1/X7XE6hqNCC1/YqW6KBI6dmYcYMGoh63cf4JWrW3NV+zinI4kEBDV6pyDtUBp3Tb+L1RmrebzT41zT+BqnI4m/mPoQ7FpRuses0RJ6vXDSp2zevJmePXvSvn17Fi9eTPPmzRk7diwXX3wxI0aMYMqUKRw6dIg2bdrQvHlzJkyYwNixYxkxYgTGGFq1asW4ceMAmD17Nq+++iq7du3ipZdeok+fPqX7/fiBbRmHGDZpqdMxvIJG9EQCw9pd++k/agH7c/IZNaADiQ1jnY4kEjDU6BXT5szN3D7tdjJyMnjjvDfoWqer05FESsXatWv56KOPOOeccxg4cCDvvPPO0c+98MILvPXWWyxd6m5OVq1axbPPPsvcuXOJiYkhIyPj6HN37tzJb7/9xpo1a7j88svV6B1Hw+oVmHFf4NaOGWt28+z3qwGN6IkEgrkb0rht3CLKhwUz+bazaVarktORRAKKGr1iWLp7KUNnDCXIBPHRhR/RMral05HE3xQx8uZJderU4ZxzzgHgxhtv5I033jjhc2fMmMHVV19NTEwMAFFRUUc/d+WVVxIUFESzZs1ITU31bGgfFR4aTIPYwL2fd/u+Q0ffV58n4t++Xrqd+z9bRnx0JKMHdqR2lQinI4kEHN0FW4TpW6Yz6OdBVAqrxPhe49Xkid/5+8jK6Y60lCtX7uj71moNffmnqMiwo+9r1U0R/2St5d1fNzBs0lLa1a3KlNs7q8kTcYgavZOYuGYiw2cNp3HVxoy7eBx1KtVxOpJIqdu6dSvz5s0D4JNPPuHcc8/9y+dDQ0PJy8sD4Pzzz+ezzz4jPT0d4C9TN0WKEh35v4sBukdPxP8UuCz/+WYVL0xdw6WtajL2lo5ULh/qdCyRgKVG7zhc1sWri17l+T+ep2udrnx40YdEhUcV/YUiPqhx48a8/fbbNG3alL179zJkyJC/fH7w4MG0atWKvn370rx5cx599FG6du1K69atuffeex1KLb6oauT/XvAZ1OmJ+JOcvAKGjF/EmHlbGNylAW9c15ZyIcFOxxIJaLpH729yC3J57PfHmLppKtc2vpaHOz5McJAKlfivkJAQxo8f/5ePzZo16+j7L774Ii+++OLRx/369aNfv35/ef7o0aP/8vjAgQOlnlN837Ev+jRzU8R/ZBzMZdCYhSzZto8nL2vGgHPqOx1JRFCj9xdZuVncM/MeFu5ayD3t7mFgi4FaGU5EpBQZA9bqHj0Rf7E1PZt+oxawY98hRvZtR88WNZ2OJCKF1OgV2nVwF0OmDWFz1mb+m/hfLm1wqdORRDwuPj6elStXOh1DAkhU+TDSD+ZqRE/EDyzbto9bxiwk32WZMOgsEuJ1m4uIN1GjB6zNWMsd0+4gOz+bkT1G0qlmJ6cjiYj4pSrlQ0k/mKsRPREfN2NNKndOWEJ0hTBGD+jImdUCd+sYEW8V8IuxzN85n/4/9gcDo3uOVpMnIuJBVcq7t1hQnyfiuyYu2MqtYxdxRrVIvrijs5o8ES8V0CN63274lifmPkF8pXhG9hhJjcgaTkcSEfFrVY82eur0RHyNtZbXflnHGzPW07VRLO/0bUdkuYB+KSni1QLyt9Nay0crP+L1xa/TsUZHXjvvNSqFVXI6loiI36tauKdW9uF8h5OIyKnIK3Dx0Ocr+HxxCtcm1OHZ3i0IDQ74iWEiXs2R31BjzDPGmOXGmKXGmJ+NMbXK6tz5rnyenf8sry9+nYvrX8zIHiPV5IkUIT4+nrS0NAAqVNAUHTl9VSPdI3p7s/McTiIixbU/J4+Boxfy+eIUhvdoxAtXtVSTJ+IDnBrRe9la+ziAMeZu4Angdk+fNDsvm3/P/jezUmYxsMVAhrUbRpBRoRIB90i3tZagIP1OiOccmbq5LzvX4SQiUhypWTkMGLWQtan7ealPK65JqON0JBEpJkcaPWtt1jEPIwHr6XNm5GRw1/S7WJm2kkfOeoTrm1zv6VOKFNuLC15kTcaaUj1mk6gm/Lvjv0/6nM2bN3PRRRdx1llnsWjRIh588EFGjBiBtZZLLrnkLxuli5SGBrGRAFSKCHU4iYgU5c/U/fQftZC92bl81C+Bbo2rOR1JRE6BY/foGWOeA24GMoHzTvK8wcADQJXY2NjTOtfWrK0MmTaE1OxUXjvvNbrX7X5axxHxR3/++Sdjxoyhbt26dOrUiUWLFlG1alUuvPBCvvrqK6688kqnI3ql0qhNgejCZtV576b2nN9ELxhFPKU06tOCTRkMGrOQsJBgJt92Ni1qVy7VjCLieR5r9Iwx04DjLWP5qLX2a2vto8CjxpiHgbuAJ493HGvt+8D7AAkJCac88rd8z3Lumn4XFsuHF35Im2ptTvUQIh5X1MibJ9WrV49OnTrx9ddf061bN468KOjbty+zZ89Wo3cCJa1NgcoYw0XNtcKxiCeVtD59v3wnwz9dSlxUBGMGdKROVPlSzyginuexRs9a26OYT50A/MAJGr2SmLl1Jg/OfpCYiBhG9hhJfOX40j6FiM+LjIx0OoKIiHiJD+ds5NnvV5NQryof9ks4uveliPgep1bdbHjMwyuA0r05CZi8djL3zLqHM6ucyfiLx6vJEylCx44d+fXXX0lLS6OgoICJEyfStWtXp2OJiEgZKHBZnvp2Fc9+v5peLWowftBZavJEfJxT9+i9YIxpDLiALZTyiptvL32bd5e9S5e4Lrzc5WXKh2rKgUhRatasyQsvvMB55513dDGWK664wulYIiLiYQUuy9CJi/lhxS4GnBPP45c0IyjIOB1LRErIqVU3r/Lk8eMrxXN1o6t55KxHCAkKyD3hRYolPj6elStXHn18/fXXc/31/1yRdvPmzUffP3DgQFlEExGRMhIcZKgfE8ljlzRlUGIDp+OISCnxyy7okgaXcEmDS5yOISIiIuITHrioidMRRKSUaWdkERERERERP6NGT8RB1gbOqvyB9L2KiIiIOE2NnohDwsPDSU9PD4gGyFpLeno64eHhTkcRERERCQh+eY+eiC+Ii4sjJSWFPXv2OB2lTISHhxMXF+d0DBEREZGAoEZPxCGhoaHUr1/f6RgiIiIi4oc0dVNERERERMTPqNETERERERHxM2r0RERERERE/IzxpRX/jDF7gC0lOEQMkFZKcUqDN+XxpiygPCfjTVmg5HnqWWtjSyuME0qhNoF3/b16UxbwrjzelAWU52RKI4vqk3f9nYJ35fGmLOBdebwpC/hfnmLVJp9q9ErKGJNkrU1wOscR3pTHm7KA8pyMN2UB78vjq7zp5+hNWcC78nhTFlCek/GmLL7M236O3pTHm7KAd+XxpiwQuHk0dVNERERERMTPqNETERERERHxM4HW6L3vdIC/8aY83pQFlOdkvCkLeF8eX+VNP0dvygLelcebsoDynIw3ZfFl3vZz9KY83pQFvCuPN2WBAM0TUPfoiYiIiIiIBIJAG9ETERERERHxe2r0RERERERE/EzANXrGmGeMMcuNMUuNMT8bY2o5nOdlY8yawkxfGmOqOJjlamPMKmOMyxjj2BK0xpiexpi1xpj1xpiHnMpRmOVjY8xuY8xKJ3MUZqljjJlpjEku/Hsa5nCecGPMAmPMssI8TzmZx9epNhWZx/H6pNp0wiyqTX7Om+qTatNxM6g2nYA31ScnalPA3aNnjKlkrc0qfP9uoJm19nYH81wIzLDW5htjXgSw1v7boSxNARfwHnC/tTbJgQzBwDrgAiAFWAhcb61NLusshXm6AAeAsdbaFk5kOCZLTaCmtXaxMaYisAi40sGfjQEirbUHjDGhwG/AMGvtfCfy+DrVpiLzOFqfVJtOmkW1yc95U31SbfrH+VWbTp7Ha+qTE7Up4Eb0jhSqQpGAo52utfZna21+4cP5QJyDWVZba9c6df5CHYH11tqN1tpcYBJwhVNhrLWzgQynzn8sa+1Oa+3iwvf3A6uB2g7msdbaA4UPQwvfAuvKUSlSbSoyj9P1SbXpBFSb/J831SfVpn9QbToJb6pPTtSmgGv0AIwxzxljtgF9gSecznOMgcBUp0M4rDaw7ZjHKTj4gsFbGWPigbbAHw7nCDbGLAV2A79Yax3N4+tUm7yaalMxqDb5Ly+tT6pNqk3F5g31qaxrk182esaYacaYlcd5uwLAWvuotbYOMAG4y+k8hc95FMgvzORoFvFuxpgKwOfAPX+7ylrmrLUF1to2uK+odjTGOD5Nw5upNpU8j3gv1Sbf5k31SbVJSpu31Keyrk0hnjy4U6y1PYr51AnAD8CTHoxTZB5jTH/gUqC79fBNk6fws3HKdqDOMY/jCj8mQOGc7s+BCdbaL5zOc4S1dp8xZibQE/CKG7C9kWpTyfI4TLXpJFSbfJ831SfVplOi2lQEb6xPZVWb/HJE72SMMQ2PeXgFsMapLOBeKQl4ELjcWpvtZBYvsRBoaIypb4wJA64DvnE4k1covIn3I2C1tfZVL8gTawpXOzPGROC+EdzR3ydfptrk9VSbTkC1yf95U31SbfoH1aaT8Kb65ERtCsRVNz8HGuNeIWkLcLu11rErH8aY9UA5IL3wQ/MdXMmqN/AmEAvsA5Zaay9yIMfFwP8BwcDH1trnyjrDMVkmAt2AGCAVeNJa+5FDWc4F5gArcP/7BXjEWvuDQ3laAWNw/z0FAZOttU87kcUfqDYVmcfx+qTadMIsqk1+zpvqk2rTcTOoNp04j9fUJydqU8A1eiIiIiIiIv4u4KZuioiIiIiI+Ds1eiIiIiIiIn5GjZ6IiIiIiIifUaMnIiIiIiLiZ9ToiYiIiIiI+Bk1egHIGDPXA8eMN8bcUNrHPc55uhljvjvF52caY0q0jK4xpr8x5q3C94cbY7YeeSwipUO16bTOq9okUgZUn07rvKpPDlOjF4CstZ09cNh4wOPF6jTNsdZe/PcPGmNCTudg1trXgCdKnEpE/kK1yU21ScT7qD65qT75FjV6AcgYc6Dwz27GmFnGmCnGmDXGmAnGGFP4uc3GmJeMMSuMMQuMMWcWfny0MabP348FvAAkGmOWGmOG/+18FYwx040xiwuPd0Xhx+ONMauNMR8YY1YZY342xkQUfq6DMWZ54fFeNsasPM73EWmM+bgw35Ijxy3ie+9mjJljjPkGSC782FfGmEWFGQYf89wBxph1xpgFwDmn8jMWkVOn2qTaJOKtVJ9Un3yRGj1pC9wDNAMa8NdfykxrbUvgLeD/ijjOQ7iv/rQpvGpzrBygt7W2HXAe8MqRogg0BN621jYH9gFXFX58FHCbtbYNUHCCcz4KzLDWdiw87svGmMgicgK0A4ZZaxsVPh5orW0PJAB3G2OijTE1gadw/zzOxf3zEZGyo9qk2iTirVSfVJ98gho9WWCtTbHWuoCluKcRHDHxmD/PLsE5DPC8MWY5MA2oDVQv/Nwma+3SwvcXAfHGmCpARWvtvMKPf3KC414IPGSMWQrMAsKBusXIs8Bau+mYx3cbY5YB84E6uAvoWcAsa+0ea20u8GkxjisipUe1SbVJxFupPqk++YTTmmcrfuXwMe8X8Nd/E/Y47+dTeIHAGBMEhBXjHH2BWKC9tTbPGLMZd2E53vkjip3cXQSvstauPYWvATh49ADGdAN6AGdba7ONMbOOySYizlFtUm0S8VaqT6pPPkEjenIy1x7z55ErRJuB9oXvXw6EFr6/H6h4guNUBnYXFqrzgHonO6m1dh+w3xhzVuGHrjvBU38Chh4zN77tyY57kmx7CwtVE6BT4cf/ALoWTkUIBa4+jWOLiGeoNqk2iXgr1SfVJ6+hRk9OpmrhlIFhwJGbhD/A/Uu8DPeUhCNXeJYDBcaYZX+/oRiYACQYY1YANwNrinHuW4APCqcWRAKZx3nOM7iL5XJjzKrCx6fqRyDEGLMa903R8wGstTuB/+Au0r8Dq0/j2CLiGapNqk0i3kr1SfXJaxhrbdHPkoBTOEUgwVqb5tD5K1hrj6xw9RBQ01o77DSO0w2431p7aSnn64/753NXaR5XRE5OtanI4/ZHtUnEEapPRR63P6pPZUojeuKtLilcHnglkAg8e5rHyQVamBJu+nmswqtuDwNZpXVMEfEZqk0i4q1Un+QvNKInIiIiIiLiZzSiJyIiIiIi4mfU6ImIiIiIiPgZNXoiIiIiIiJ+Ro2eiIiIiIiIn1GjJyIiIiIi4mf+H1q817FJG35VAAAAAElFTkSuQmCC\n",
- "text/plain": [
- "<Figure size 1080x288 with 3 Axes>"
- ]
- },
- "metadata": {},
- "output_type": "display_data"
- }
- ],
- "source": [
- "fig, axarr = plt.subplots(1, 3, figsize=(15, 4),\n",
- " sharey=True)\n",
- "ax = axarr[0]\n",
- "solution_x_axis.plot(ax=ax)\n",
- "ax.set_title('x-aligned')\n",
- "ax = axarr[1]\n",
- "solution_y_axis.plot(ax=ax)\n",
- "ax.set_title('y-aligned')\n",
- "ax = axarr[2]\n",
- "solution_z_axis.plot(ax=ax)\n",
- "ax.set_title('z-aligned')\n",
- "\n",
- "for ax in axarr:\n",
- " ax.set_xlabel('input angle [rad]')\n",
- " ax.set_ylabel('euler angle [rad]')\n"
+ "euler_axes = 'yxz'\n",
+ "alpha, beta, gamma = twomarkers2euler(mark0, mark1, axis_alignement,\n",
+ " known_angle, euler_axes)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "We observe that when the markers are aligned with the x-axis, the yaw and pitch can be determined (without having to guess the roll), because yaw and pitch are constant. However we can not obtain the roll angle. \n",
+ "#### Advance cases\n",
+ "\n",
+ "The few cases above may be a bit limiting. We will illustrate, how the solution of the system can be done for other axis. Note that, although the method will not expose a general solution, it will hopefully serve as basis for deriving euler angles in custom cases. \n",
+ "\n",
+ "We focus on the following euler angle sequence: $R=R_zR_yR_x$.\n",
+ "\n",
+ "We start with the following equation:\n",
+ "$v^{bee}=R.v^{ref}$\n",
+ "\n",
+ "here:\n",
+ "* $v^{bee}$ is a vector formed by our two markers\n",
+ "* $v^{ref}$ is the vector formed by our two markers when the bee will have a null orientation (yaw=pitch=roll=0)\n",
+ "* $R$ the rotation matrix\n",
+ "\n",
+ "We write the transformation with a matrix of cosine and sine of euler angles\n",
+ "\n",
+ "$$\\begin{align}\n",
+ "v^{bee}&=Rv^{ref}\\\\\n",
+ " &=R_zR_yR_xv^{ref}\\\\\n",
+ " &=\\begin{pmatrix}\\cos\\gamma & \\sin\\gamma & 0 \\\\\n",
+ " -\\sin\\gamma & \\cos\\gamma & 0 \\\\\n",
+ " 0 & 0 & 1 \\end{pmatrix}\n",
+ " \\begin{pmatrix}\\cos\\beta & 0 & -\\sin\\beta \\\\\n",
+ " 0 & 1 & 0 \\\\\n",
+ " \\sin\\beta & 0 & \\cos\\beta \\end{pmatrix}\n",
+ " \\begin{pmatrix}1 & 0 & 0 \\\\\n",
+ " 0 & \\cos\\alpha & \\sin\\alpha \\\\\n",
+ " 0 & -\\sin\\alpha & \\cos\\alpha \\end{pmatrix}v^{ref}\\\\\n",
+ " &=\\begin{pmatrix}\n",
+ "\\cos\\gamma \\cos\\beta & \\cos\\gamma \\sin\\beta \\sin\\alpha + \\sin\\gamma \\cos\\alpha & -\\cos\\gamma \\sin\\beta\\cos\\alpha + \\sin\\gamma \\sin\\alpha \\\\\n",
+ "-\\sin\\gamma \\cos\\beta & -\\sin\\gamma \\sin\\beta \\sin\\alpha + \\cos\\gamma \\cos\\alpha & \\sin\\gamma \\sin\\beta \\cos\\alpha + \\cos\\gamma \\sin\\alpha \\\\\n",
+ "\\sin\\beta & -\\cos\\beta \\sin\\alpha & \\cos\\beta \\cos\\alpha\\end{pmatrix}v^{ref}\\\\\n",
+ "\\end{align}$$\n",
"\n",
- "When the two markers are aligned with the y-axis or z-axis, the pitch or the yaw needs to be apriori known to get the two other angles, repectively. "
+ "To simplify the problem, we need to remove terms in the system of equations. Removing terms is easily done by letting certain variables to be zero. We want to determine $\\alpha$, $\\beta$, and $\\gamma$. Thus we can only set to zeros $v_x^{ref}$, $v_y^{ref}$, or $v_z^{ref}$.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "### Solving pitch aligned markers with yaw, pitch, roll convention\n",
- "\n",
- "$v^{ref}=(0,-1,0)^T$\n",
+ "##### Case 1: $v^{ref}$ is y-aligned\n",
"\n",
- "$$\\begin{align}\n",
- " v_x^{bee} & = -\\cos\\alpha \\sin\\beta \\sin\\gamma + \\sin\\alpha \\cos\\gamma &\\quad \\text{from L1}\\\\\n",
- " v_y^{bee} & = -\\sin\\alpha \\sin\\beta \\sin\\gamma - \\cos\\alpha \\cos\\gamma &\\quad \\text{from L2}\\\\\n",
- " v_z^{bee} & = -\\cos\\beta \\sin\\gamma &\\quad \\text{from L3}\n",
- " \\end{align}$$\n",
+ "We know that the two markers are aligned with the y-axis and thus: $v^{ref}=[0, 1, 0]^T$. The system of equation is then simplified as follow:\n",
"\n",
- "#### for $v_z^{bee}\\neq0 \\Rightarrow \\sin\\gamma\\neq0$\n",
+ "$$\n",
+ "v^{bee} = \\begin{pmatrix} x \\\\ y \\\\ z \\end{pmatrix}^{bee}= \\begin{pmatrix} \\cos\\gamma\\sin\\beta\\sin\\alpha + \\sin\\gamma\\cos\\alpha \\\\ \n",
+ "-\\sin\\gamma\\sin\\beta\\sin\\alpha + \\cos\\gamma\\cos\\alpha \\\\\n",
+ "-\\cos\\beta\\sin\\alpha \\end{pmatrix}\n",
+ "$$\n",
"\n",
- "$$\\begin{align}\n",
- " \\cos\\gamma & = v_x^{bee}\\sin\\alpha -v_y^{bee}\\cos\\alpha &\\quad \\text{from L1 and L2} \\\\\n",
- " \\cos\\beta & = -\\frac{v_z^{bee}}{\\sin\\gamma } &\\quad\\text{from L3}\n",
- " \\end{align}$$\n",
+ "The simplest equation is $z=-\\cos\\beta\\sin\\alpha$. If we know $\\beta$ and that $\\beta$ is different of $\\pm\\pi/2$, we have the following relation:\n",
"\n",
- "From the two last equation $\\gamma$ and $\\beta$ can be found as a function of $\\alpha$ and $v^{bee}$\n",
+ "$$ \\alpha = \\arcsin\\frac{z}{-\\cos\\beta} + k\\pi ,\\quad k\\in\\mathbf{Z}$$ \n",
"\n",
- "#### for $v_z^{bee}=0$\n",
- "$$\\begin{align}\n",
- " \\tan\\alpha & = -\\frac{v_x^{bee}}{v_y^{bee}} &\\quad \\text{from L1 and L2}\\\\\n",
- " \\cos\\gamma & =\\pm\\sqrt{ \\left(v_x^{bee}\\right)^2 + \\left(v_y^{bee}\\right)^2 } &\\quad \\text{from L1 and L2}\\\\\n",
- " \\end{align}$$\n",
+ "We found $\\alpha$, know we can look for $\\gamma$ by using the two other equations. We recognise that the two equations share many term and that we can probably extract $\\sin\\gamma$ and $\\cos\\gamma$ out of the two equations. \n",
"\n",
- "for $\\beta$ known, and $\\beta\\neq\\pm\\pi/2$:\n",
+ "$$\n",
+ "\\begin{align}\n",
+ "x\\sin\\beta\\sin\\alpha + y\\cos\\alpha &= \\cos\\gamma\\sin^2\\beta\\sin^2\\alpha + \\sin\\gamma\\sin\\alpha\\sin\\beta\\sin\\alpha - \\sin\\gamma\\sin\\alpha\\sin\\beta\\cos\\alpha + \\cos\\gamma\\cos^2\\alpha \\\\\n",
+ " &= \\cos\\gamma(\\sin^2\\beta\\sin^2\\alpha + \\cos^2\\alpha) \\\\\n",
+ "x\\cos\\alpha - y\\sin\\beta\\sin\\alpha &= \\cos\\alpha\\cos\\gamma\\sin\\beta\\sin\\alpha + \\sin\\gamma\\cos^2\\alpha + \\sin\\gamma\\sin^2\\beta\\sin^2\\alpha - \\cos\\gamma\\cos\\alpha\\beta\\sin\\alpha \\\\\n",
+ " &= \\sin\\gamma(\\cos^2\\alpha + \\sin^2\\alpha\\sin^2\\beta)\n",
+ "\\end{align}$$\n",
"\n",
- "$$\\begin{align}\n",
- " \\sin\\gamma & = \\frac{v_z^{bee}}{-\\cos\\beta} \\\\\n",
- " \\tan\\left(\\alpha+\\theta\\right) & = \\frac{v_x^{bee}}{-v_y^{bee}} \\\\\n",
- " \\tan\\theta&=\\frac{v_z^{bee}\\tan\\beta}{\\pm\\sqrt{1-\\sin^2\\gamma}}\n",
- " \\end{align}$$\n",
+ "The two equations above yield to:\n",
"\n",
- "for $\\gamma$ known, and $\\gamma\\neq 0 + k2pi$:\n",
+ "$$\\tan\\gamma = \\frac{\\sin\\gamma}{\\cos\\gamma}$$\n",
"\n",
- "$$\\begin{align}\n",
- " \\sin\\left(\\alpha+\\theta\\right) & = \\frac{\\cos\\gamma}{\\sqrt{ \\left(v_x^{bee}\\right)^2+\\left(v_y^{bee}\\right)^2}} \\\\\n",
- " \\cos\\beta & = -\\frac{v_z^{bee}}{\\sin\\gamma } \\\\\n",
- " \\tan\\theta&=\\frac{v_y^{bee}}{v_x^{bee}}\n",
- " \\end{align}$$"
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 120,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "-0.06706145865876341 -0.06706145865876341\n",
+ "0.3759553148789634 0.3759553148789634\n",
+ "-0.5011653731009601 -0.5011653731009601\n"
+ ]
+ }
+ ],
+ "source": [
+ "import pandas as pd\n",
+ "import numpy as np\n",
+ "from navipy.trajectories.transformations import twomarkers2euler\n",
+ "from navipy.maths.euler import matrix as euler_rot\n",
+ "\n",
+ "def twomarker_yaligned_zyx(vbee, beta_th):\n",
+ " # Get alpha knowing beta\n",
+ " beta = beta_th\n",
+ " alpha = np.arcsin(vbee.z/(-np.cos(beta)))\n",
+ " # Get gamma knowing the two others\n",
+ " nominator = vbee.x*np.cos(alpha) - vbee.y*np.sin(beta)*np.sin(alpha)\n",
+ " nominator /= np.cos(alpha)**2 + (np.sin(alpha)**2) * (np.sin(beta)**2)\n",
+ " denominator = vbee.x*np.sin(beta)*np.sin(alpha)+vbee.y*np.cos(alpha)\n",
+ " denominator /= np.cos(alpha)**2+(np.sin(alpha)**2)*(np.sin(beta)**2)\n",
+ " gamma = np.arctan2(nominator, denominator)\n",
+ " return alpha, beta, gamma\n",
+ "\n",
+ "# Get the position of two markers\n",
+ "# knowing the euler angles\n",
+ "# in order to test our function above\n",
+ "euler_axes = 'zyx'\n",
+ "alpha_th = np.pi*(np.random.rand()-0.5)\n",
+ "beta_th = np.pi*(np.random.rand()-0.5)\n",
+ "gamma_th = np.pi*(np.random.rand()-0.5) #2*np.pi*np.random.rand()\n",
+ "rotmat = euler_rot(gamma_th, beta_th, alpha_th, axes=euler_axes)[:3, :3]\n",
+ "axis_alignment = np.array([0, 1, 0]).transpose()\n",
+ "\n",
+ "mark0 = pd.Series(index=['x','y','z'], data=0)\n",
+ "mark1 = pd.Series(index=['x','y','z'], \n",
+ " data=rotmat.dot(axis_alignment))\n",
+ "# Call the function\n",
+ "vbee = mark1-mark0\n",
+ "alpha, beta, gamma = twomarker_yaligned_zyx(vbee, beta_th)\n",
+ "\n",
+ "# Check the results\n",
+ "print(alpha_th, alpha)\n",
+ "print(beta_th, beta)\n",
+ "print(gamma_th, gamma)"
]
}
],
diff --git a/navipy/arenatools/cam_dlt.py b/navipy/arenatools/cam_dlt.py
index 30877498b1ec8b42a94950c5aead39825ef600fd..baa8bd3283a171f0f61ae2fe29370343a9eb0cbc 100644
--- a/navipy/arenatools/cam_dlt.py
+++ b/navipy/arenatools/cam_dlt.py
@@ -101,6 +101,53 @@ units are [u,v] i.e. camera coordinates or pixels
return xyz, rmse
+def dlt_principal_point(coeff):
+ normalisation = np.sum(coeff[8:11]**2)
+ u_0 = np.sum(coeff[0:3]*coeff[8:11])/normalisation
+ v_0 = np.sum(coeff[4:7]*coeff[8:11])/normalisation
+ return u_0, v_0
+
+
+def dlt_principal_distance(coeff):
+ normalisation = np.sum(coeff[8:11]**2)
+ return 1/np.sqrt(normalisation)
+
+
+def dlt_scale_factors(coeff):
+ u_0, v_0 = dlt_principal_point(coeff)
+
+ normalisation = np.sum(coeff[8:11]**2)
+ du = (u_0*coeff[8] - coeff[0])**2
+ du += (u_0*coeff[9] - coeff[1])**2
+ du += (u_0*coeff[10] - coeff[2])**2
+ du /= normalisation
+
+ dv = (v_0*coeff[8] - coeff[4])**2
+ dv += (v_0*coeff[9] - coeff[5])**2
+ dv += (v_0*coeff[10] - coeff[6])**2
+ dv /= normalisation
+ return du, dv
+
+
+def dlt_transformation_matrix(coeff):
+ u_0, v_0 = dlt_principal_point(coeff)
+ d = dlt_principal_distance(coeff)
+ du, dv = dlt_scale_factors(coeff)
+ transform = np.zeros((3, 3))
+ transform[0, 0] = (u_0*coeff[8]-coeff[0])/du
+ transform[0, 1] = (u_0*coeff[9]-coeff[1])/du
+ transform[0, 2] = (u_0*coeff[10]-coeff[2])/du
+
+ transform[1, 0] = (v_0*coeff[8]-coeff[4])/dv
+ transform[1, 1] = (v_0*coeff[9]-coeff[5])/dv
+ transform[1, 2] = (v_0*coeff[10]-coeff[6])/dv
+
+ transform[2, 0] = coeff[8]
+ transform[2, 1] = coeff[9]
+ transform[2, 2] = coeff[10]
+ return d*transform
+
+
def dlt_inverse(coeff, frames):
"""
This function reconstructs the pixel coordinates of a 3D coordinate as
@@ -123,10 +170,68 @@ def dlt_inverse(coeff, frames):
coeff[5]+frames[:, 2]*coeff[6]+coeff[7]
uv[:, 0] /= normalisation
uv[:, 1] /= normalisation
+ # Apply distortion
+ delta_uv = np.zeros((frames.shape[0], 2))
+ u_0, v_0 = dlt_principal_point(coeff)
+ zeta = uv[:, 0] - u_0
+ eta = uv[:, 1] - v_0
+ rsq = zeta**2 + eta**2
+ if coeff.shape[0] > 11:
+ delta_uv[:, 0] += zeta*(coeff[11]*rsq)
+ delta_uv[:, 1] += eta*(coeff[11]*rsq)
+ if coeff.shape[0] > 13:
+ delta_uv[:, 0] += zeta*(coeff[12]*(rsq**2) + coeff[13]*(rsq**3))
+ delta_uv[:, 1] += eta*(coeff[12]*(rsq**2) + coeff[13]*(rsq**3))
+ if coeff.shape[0] > 15:
+ delta_uv[:, 0] += coeff[14]*(rsq + 2*(zeta**2)) + coeff[15]*eta*zeta
+ delta_uv[:, 1] += coeff[15]*(rsq + 2*(eta**2)) + coeff[14]*eta*zeta
+ # print(eta, zeta, rsq)
+ uv += delta_uv
return uv
-def dlt_compute_coeffs(frames, campts):
+def _dlt_matrices_calib(vframes, vcampts, nparams=11, l1to11=np.zeros(11)):
+ # re arange the frame matrix to facilitate the linear least
+ # sqaures solution
+ if nparams < 11:
+ nparams = 11
+ matrix = np.zeros((vframes.shape[0]*2, 16)) # 16 for the dlt
+ u_0, v_0 = dlt_principal_point(l1to11)
+ for num_i, index_i in enumerate(vframes.index):
+ # eta, zeta, Rsq depends on L9to11
+ zeta = vcampts.loc[index_i, 'u'] - u_0 # -u_0 = 0 ??
+ eta = vcampts.loc[index_i, 'v'] - v_0 # -v_0 = 0
+ R = np.sum(l1to11[-3:] * vframes.loc[index_i, ['x', 'y', 'z']])+1
+ rsq = eta**2 + zeta**2
+ # populate matrix
+ matrix[2*num_i, 0:3] = vframes.loc[index_i, ['x', 'y', 'z']]
+ matrix[2*num_i+1, 4:7] = vframes.loc[index_i, ['x', 'y', 'z']]
+ matrix[2*num_i, 3] = 1
+ matrix[2*num_i+1, 7] = 1
+ matrix[2*num_i, 8:11] = vframes.loc[index_i,
+ ['x', 'y', 'z']]*(-vcampts.loc[index_i, 'u'])
+ matrix[2*num_i+1, 8: 11] = vframes.loc[index_i,
+ ['x', 'y', 'z']]*(-vcampts.loc[index_i, 'v'])
+ # 12th parameter
+ matrix[2*num_i, 11] = zeta*rsq*R
+ matrix[2*num_i+1, 11] = eta*rsq*R
+ # 13th and 14th parameters
+ matrix[2*num_i, 12] = zeta*(rsq**2)*R
+ matrix[2*num_i+1, 12] = eta*(rsq**2)*R
+ matrix[2*num_i, 13] = zeta*(rsq**3)*R
+ matrix[2*num_i+1, 13] = eta*(rsq**3)*R
+ # 15th and 16th parameters
+ matrix[2*num_i, 12] = (rsq + 2*(zeta**2))*R
+ matrix[2*num_i+1, 12] = eta*zeta*R
+ matrix[2*num_i, 13] = eta*zeta*R
+ matrix[2*num_i+1, 13] = (rsq + 2*(eta**2))*R
+
+ matrix[2*num_i, :] /= R
+ matrix[2*num_i+1, :] /= R
+ return matrix[:, : nparams]
+
+
+def dlt_compute_coeffs(frames, campts, nparams=11, niter=100):
"""
A basic implementation of 11 parameters DLT
@@ -142,36 +247,44 @@ single plane.
# remove NaNs
valid_idx = frames.dropna(how='any').index
valid_idx = campts.loc[valid_idx, :].dropna(how='any').index
-
# valid df
vframes = frames.loc[valid_idx, :]
vcampts = campts.loc[valid_idx, :]
-
- # re arange the frame matrix to facilitate the linear least
- # sqaures solution
- matrix = np.zeros((vframes.shape[0]*2, 11)) # 11 for the dlt
- for num_i, index_i in enumerate(vframes.index):
- matrix[2*num_i, 0:3] = vframes.loc[index_i, ['x', 'y', 'z']]
- matrix[2*num_i+1, 4:7] = vframes.loc[index_i, ['x', 'y', 'z']]
- matrix[2*num_i, 3] = 1
- matrix[2*num_i+1, 7] = 1
- matrix[2*num_i, 8:11] = \
- vframes.loc[index_i, ['x', 'y', 'z']]*(-vcampts.loc[index_i, 'u'])
- matrix[2*num_i+1, 8:11] = \
- vframes.loc[index_i, ['x', 'y', 'z']]*(-vcampts.loc[index_i, 'v'])
-
- # re argen the campts array for the linear solution
- vcampts_f = np.reshape(np.flipud(np.rot90(vcampts)), vcampts.size, 1)
- print(vcampts_f.shape)
- print(matrix.shape)
+ # Get the matrices for calib
+ matrix = _dlt_matrices_calib(vframes, vcampts)
+ vcampts_f = vcampts.values.flatten() # [u_1, v_1, ... u_n, v_n]
# get the linear solution the 11 parameters
coeff = np.linalg.lstsq(matrix, vcampts_f, rcond=None)[0]
# compute the position of the frame in u,v coordinates given the linear
# solution from the previous line
-
matrix_uv = dlt_inverse(coeff, vframes)
-
# compute the rmse between the ideal frame u,v and the
# recorded frame u,v
rmse = np.sqrt(np.mean(np.sum((matrix_uv-vcampts)**2)))
- return coeff, rmse
+ if nparams == 11:
+ return coeff, rmse
+ # Now we can try to guess the other coefficients
+ if nparams in [12, 14, 16]:
+ for _ in range(niter):
+ # 9th to 11th parameters are index 8 to 10 (0 being the 1st param)
+ l1to11 = coeff[: 11]
+ # Get the matrices for calib
+ matrix = _dlt_matrices_calib(vframes, vcampts,
+ nparams=nparams, l1to11=l1to11)
+ vcampts_normed = vcampts_f.copy()
+ for num_i, index_i in enumerate(vframes.index):
+ normalisation = np.sum(
+ l1to11[-3:] * vframes.loc[index_i, ['x', 'y', 'z']])+1
+ vcampts_normed[2*num_i] /= normalisation
+ vcampts_normed[2*num_i + 1] /= normalisation
+ coeff = np.linalg.lstsq(matrix, vcampts_normed, rcond=None)[0]
+ print(coeff)
+ # compute the position of the frame in u,v coordinates given the linear
+ # solution from the previous line
+ matrix_uv = dlt_inverse(coeff, vframes)
+ # compute the rmse between the ideal frame u,v and the
+ # recorded frame u,v
+ rmse = np.sqrt(np.mean(np.sum((matrix_uv-vcampts)**2)))
+ return coeff, rmse
+ else:
+ raise ValueError('nparams can be either [11,12,14,16]')
diff --git a/navipy/scripts/check_blender_versions_pip.py b/navipy/scripts/check_blender_versions_pip.py
new file mode 100644
index 0000000000000000000000000000000000000000..105511ab312d6ff6183720e0bc16d769bec26b2f
--- /dev/null
+++ b/navipy/scripts/check_blender_versions_pip.py
@@ -0,0 +1,17 @@
+"""
+List all installed packages in a python installation
+by using pip.
+
+It can be used to install a virtual environment, so that it matches
+another installation (for example blender)
+"""
+import pip
+
+
+installed_packages = pip.get_installed_distributions()
+installed_packages = sorted(["%s==%s" % (i.key, i.version)
+ for i in installed_packages])
+with open('requirement.txt', 'w') as cfile:
+ for line in installed_packages:
+ print(line)
+ cfile.write(line+'\n')
diff --git a/navipy/scripts/dlt_calibrator.py b/navipy/scripts/dlt_calibrator.py
index ed9b12a6d91b30acb4496bd8b1b5258028e6cf06..409986bbbd029a33e11f09f34da8f6c6ce7d2485 100644
--- a/navipy/scripts/dlt_calibrator.py
+++ b/navipy/scripts/dlt_calibrator.py
@@ -5,6 +5,10 @@ import cv2
import os
from navipy.arenatools.cam_dlt import dlt_inverse
from navipy.arenatools.cam_dlt import dlt_compute_coeffs
+from navipy.arenatools.cam_dlt import dlt_principal_point
+from navipy.arenatools.cam_dlt import dlt_principal_distance
+from navipy.arenatools.cam_dlt import dlt_scale_factors
+from navipy.arenatools.cam_dlt import dlt_transformation_matrix
keybinding = dict()
@@ -42,6 +46,16 @@ def parser_dlt_calibrator():
parser.add_argument('-s', '--scale',
default=0.5,
help=arghelp)
+ arghelp = 'number of dlt parameters (coeff)'
+ parser.add_argument('-c', '--ndltcoeff',
+ type=int,
+ default=11,
+ help=arghelp)
+ arghelp = 'number of iteration for calibration'
+ parser.add_argument('-e', '--epoque',
+ type=int,
+ default=100,
+ help=arghelp)
return parser
@@ -71,6 +85,8 @@ def main():
# The image may need to be scale because screen may have
# less pixel than the image
scale = args['scale']
+ ndlt_coeff = args['ndltcoeff']
+ niter = args['epoque']
imageref = cv2.imread(args['image'])
# define some constants
showframe_ref = 50
@@ -144,9 +160,17 @@ def main():
break
if key == ord("c"):
print('calibrate')
- coeff, rmse = dlt_compute_coeffs(frames, campts)
+ coeff, rmse = dlt_compute_coeffs(
+ frames, campts, nparams=ndlt_coeff,
+ niter=niter)
print(rmse)
print(coeff)
+ print('principal points: {}'.format(dlt_principal_point(coeff)))
+ print('principal distance: {}'.format(
+ dlt_principal_distance(coeff)))
+ print('scale factor: {}'.format(dlt_scale_factors(coeff)))
+ print('transform:')
+ print(dlt_transformation_matrix(coeff))
coeff = pd.Series(data=coeff)
coeff.to_csv(os.path.splitext(args['points'])[0]+'_coeff.csv')