{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Laparoscopic Liver Surgery Simulations\n", "\n", "In this notebook, various simulation is conducted, to illustrate the potential accuracy of components of a laparoscopic liver surgery system.\n", "\n", "Aim: To simulate\n", "\n", "![System Diagram](LaparoscopicLiverSimulation/lap_liver_sim_diagram.jpg \"System Diagram\")\n", "\n", "and see the accuracy." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Jupyter notebook sets the cwd to the folder containing the notebook.\n", "# So, you want to add the root of the project to the sys path, so modules load correctly.\n", "import sys\n", "sys.path.append(\"../../\")" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# All imports.\n", "import random\n", "import copy\n", "import cv2\n", "import numpy as np\n", "from scipy.spatial.transform import Rotation as R\n", "import sksurgeryopencvpython as cvpy\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "def extract_rigid_body_parameters(matrix):\n", " t = matrix[0:3, 3]\n", " r = matrix[0:3, 0:3]\n", " r = R.from_matrix(r)\n", " euler = r.as_euler('zyx', degrees=True)\n", " return [euler[0], euler[1], euler[2], t[0], t[1], t[2]]\n", "\n", "def rigid_body_parameters_to_matrix(params):\n", " matrix = np.eye(4)\n", " r = (R.from_euler('zyx', [params[0], params[1], params[2]], degrees=True)).as_matrix()\n", " matrix[0:3, 0:3] = r\n", " matrix[0][3] = params[3]\n", " matrix[1][3] = params[4]\n", " matrix[2][3] = params[5]\n", " return matrix\n", "\n", "def convert_camera_point_to_world(point, camera_to_marker, marker_to_world):\n", " return marker_to_world @ camera_to_marker @ point\n", "\n", "def convert_4x1_to_1x1x3(p_41):\n", " p_113 = np.zeros((1,1,3))\n", " p_113[0][0][0] = p_41[0][0]\n", " p_113[0][0][1] = p_41[1][0]\n", " p_113[0][0][2] = p_41[2][0]\n", " return p_113\n", "\n", "def convert_1x2_to_1x1x2(p_12):\n", " p_112 = np.zeros((1,1,2))\n", " p_112[0][0][0] = p_12[0]\n", " p_112[0][0][1] = p_12[1] \n", " return p_112\n", "\n", "def project_camera_point_to_image(point, intrinsics, distortion):\n", " rvec = np.zeros((1,3))\n", " tvec = np.zeros((1,3))\n", " image_points, jacobian = cv2.projectPoints(convert_4x1_to_1x1x3(point), rvec, tvec, intrinsics, distortion)\n", " return image_points[0][0] # returns a list\n", " \n", "def convert_left_camera_to_right_camera(point, left_to_right):\n", " return left_to_right @ point\n", "\n", "def triangulate_points_to_3d(left_point, left_intrinsic, left_distortion, right_point, right_intrinsic, right_distortion, left_to_right):\n", " left_point_undistorted = cv2.undistortPoints(convert_1x2_to_1x1x2(left_point), left_intrinsic, left_distortion, None, left_intrinsic)\n", " right_point_undistorted = cv2.undistortPoints(convert_1x2_to_1x1x2(right_point), right_intrinsic, right_distortion, None, right_intrinsic)\n", " image_points = np.zeros((1,4))\n", " image_points[0][0] = left_point_undistorted[0][0][0]\n", " image_points[0][1] = left_point_undistorted[0][0][1]\n", " image_points[0][2] = right_point_undistorted[0][0][0]\n", " image_points[0][3] = right_point_undistorted[0][0][1]\n", " reconstructed = cvpy.triangulate_points_using_hartley(image_points,\n", " left_intrinsic,\n", " right_intrinsic,\n", " left_to_right[0:3, 0:3],\n", " left_to_right[0:3, 3]\n", " )\n", " result = np.ones((4,1))\n", " result[0][0] = reconstructed[0][0]\n", " result[1][0] = reconstructed[0][1]\n", " result[2][0] = reconstructed[0][2]\n", " return result \n", "\n", "def initialise_gold_standard(cam_x, cam_y, cam_z):\n", " point_in_left_camera_space = np.ones((4,1))\n", " point_in_left_camera_space[0][0] = cam_x\n", " point_in_left_camera_space[1][0] = cam_y\n", " point_in_left_camera_space[2][0] = cam_z\n", " point_in_right_camera_space = convert_left_camera_to_right_camera(point_in_left_camera_space, left_to_right)\n", "\n", " left_image_point = project_camera_point_to_image(point_in_left_camera_space, left_intrinsics, left_distortion)\n", " right_image_point = project_camera_point_to_image(point_in_right_camera_space, right_intrinsics, right_distortion)\n", " reconstructed = triangulate_points_to_3d(left_image_point, left_intrinsics, left_distortion, right_image_point, right_intrinsics, right_distortion, left_to_right)\n", "\n", " assert np.allclose(point_in_left_camera_space, reconstructed)\n", "\n", " gold_standard_world_point = convert_camera_point_to_world(point_in_left_camera_space, camera_to_marker, marker_to_world)\n", " return point_in_left_camera_space, gold_standard_world_point, left_image_point, right_image_point\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Load reference data.\n", "\n", "\n", "This data comes from the SmartLiver system." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "# Camera intrinsic calibration\n", "left_intrinsics = np.loadtxt('LaparoscopicLiverSimulation/lap_liver_left_intrinsics.txt')\n", "left_distortion = np.loadtxt('LaparoscopicLiverSimulation/lap_liver_left_distortion.txt')\n", "right_intrinsics = np.loadtxt('LaparoscopicLiverSimulation/lap_liver_right_intrinsics.txt')\n", "right_distortion = np.loadtxt('LaparoscopicLiverSimulation/lap_liver_right_distortion.txt')\n", "\n", "# Hand-eye or, marker-to-camera\n", "marker_to_camera = np.loadtxt('LaparoscopicLiverSimulation/lap_liver_sim_marker_to_camera.txt')\n", "camera_to_marker = np.linalg.inv(marker_to_camera)\n", "\n", "# Its a stereo laparoscope, so separation between right and left camera.\n", "left_to_right = np.loadtxt('LaparoscopicLiverSimulation/lap_liver_sim_l2r.txt')\n", "\n", "# And it's tracked, so, the marker to world transform.\n", "marker_to_world = np.loadtxt('LaparoscopicLiverSimulation/lap_liver_sim_marker_to_world.txt')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Create Simulation Parameters\n", "\n", "These parameters are the ones that define the simulation. As we vary them within certain ranges, we will see how much they change the accuracy of the result. \n" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "# First parameterise the main transformations. \n", "# The order is: rx, ry, rx, tx, ty, tz, where rotations are in degrees and translations in millimetres.\n", "marker2camera = extract_rigid_body_parameters(marker_to_camera)\n", "marker2world = extract_rigid_body_parameters(marker_to_world)\n", "left2right = extract_rigid_body_parameters(left_to_right)\n", "\n", "# Create some fixed parameters\n", "image_size = (1920, 1080) # width, height, in pixels\n", "distance_from_camera = 75 # millimetres" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Create Gold Standard\n", "\n", "We need to define the thing we are measuring or evaluating. Imagine we are trying to assess the accuracy of triangulating a point that is in front of the camera, and correctly positioning it in world space. So, first, generate a point in camera coordinates, project to 2D left and right images, triangulate back to 3D, to check we get the same point, and convert to world space." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Starting gold standard:\n", "Left camera point:[[ 5.]\n", " [10.]\n", " [75.]\n", " [ 1.]]\n", "World point:[[ 1.02579661e+02]\n", " [ 2.55747093e+00]\n", " [-1.76820496e+03]\n", " [ 1.00000000e+00]]\n", "Left image point:[904.48111684 790.14748833]\n", "Right image point:[1028.96090457 791.28608348]\n" ] } ], "source": [ "gold_left_camera_point, gold_world_point, gold_left_image_point, gold_right_image_point = initialise_gold_standard(5, 10, distance_from_camera)\n", "\n", "print(\"Starting gold standard:\")\n", "print(\"Left camera point:\" + str(gold_left_camera_point))\n", "print(\"World point:\" + str(gold_world_point))\n", "print(\"Left image point:\" + str(gold_left_image_point))\n", "print(\"Right image point:\" + str(gold_right_image_point))\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Experiment 1: Effect of noise on 2D feature detector\n", "\n", "So, given the gold standard we have established, we can start to investigate characteristics of the system. For example, if we are triangulating points, then how accurate will it be. Or put another way, what is the effect of errors in our feature detector?" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEGCAYAAACAd+UpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAjX0lEQVR4nO3de7yWY9r/8c/RTlI2bVAyajQRlWIxvEIlJEK2oSZKNZkhmjxjNzbzPGTTQ4SJqCFMQhSVFGpST9Fqx6oUP1NPEWWRNtO+4/fHedfTbrXuVeu6r3vzfb9e67XWvbuu446+neu8z+s4zd0REZHsVCbuAkREJDoKeRGRLKaQFxHJYgp5EZEsppAXEcli5eIuYEfVq1f3OnXqxF2GiEjGmDFjxo/uXqOox9Mq5OvUqUN+fn7cZYiIZAwzW7y3xzVdIyKSxRTyIiJZTCEvIpLF0mpOfk82bdrE0qVLWb9+fdylZJWKFStSu3ZtypcvH3cpIhKhtA/5pUuXUqVKFerUqYOZxV1OVnB3CgsLWbp0KXXr1o27HBGJUNpP16xfv55q1aop4EuRmVGtWjX9diSSA9I+5AEFfAT0ZyqSGzIi5EVEstbUqdC3b2SHV8iLiMTBHfr3h7PPhuefhzVrIjmNQj4JZcuWpUmTJjRs2JCLL76YlStXArBo0SLMjL/85S/bn/vjjz9Svnx5br75ZgAWLFhAixYtaNKkCQ0aNKB79+57PMeyZcto27Ztqde+ceNGzj77bDZv3lzqxxaRfbRmDVx3Hdx6K7RpA/n5ULlyJKeKPOTNrKyZzTKzUVGfKyoHHnggs2fPpqCggKpVq/Lss89uf6xu3bqMHj16++0333yTE088cfvtnj170qtXL2bPns38+fO55ZZb9niOJ554gm7dupV67RUqVKBVq1YMGzas1I8tIvvgyy/htNPgjTegTx8YMQIOPTSy06ViCeWtwHzg4P0+0m23wezZ+32YnTRpAk8+mfTTzzjjDD7//PPttytVqkSDBg3Iz88nLy+PYcOGcfXVV/Pdd98BYYReu3bt7c9v1KjRHo87fPhwHnzwQQBeeuklRowYwdq1a/nqq6+4/fbb2bhxI6+88goHHHAAY8aMoWrVqrRo0YKmTZvyySefsHbtWoYMGcLDDz/MF198Qfv27bcfr127dtx111106NChhH84IlKq3ngDbrwRDjwQxo2DVq0iP2WkI3kzqw1cBLwY5XlSZcuWLXz00UdccsklO91/zTXX8Prrr7NkyRLKli1LrVq1tj/Wq1cvzjnnHNq0aUO/fv22T/Xs6F//+heHHXYYBxxwwPb7CgoKePvtt5k+fTr33HMPlSpVYtasWZxxxhkMGTJk+/MqVKhAfn4+PXr04NJLL+XZZ5+loKCAl156icLCQgAaNmzI9OnTS/lPQ0SStmkT9OoF7dtDo0Ywc2ZKAh6iH8k/CfwZqFLUE8ysO9Ad4Fe/+lUxR3uy1AoriXXr1tGkSRO+/fZbGjRowHnnnbfT4xdccAH33nsvRxxxBO3bt9/psc6dO9O6dWvGjh3LyJEjef7555kzZ85Ogb5s2TJq1Ni5U2jLli2pUqUKVapU4ZBDDuHiiy8Gwm8CO/4mse0fnEaNGnHiiSdSs2ZNAH7961+zZMkSqlWrRtmyZalQoQKrV6+mSpUi/1OISBS++w6uvhqmTIGePcNKmgoVUnb6yEbyZtYWWO7uM/b2PHcf6O557p63a9Cli21z8osXL8bdd5qThzCaPuWUU3j88ce58sord3t9rVq16NKlCyNHjqRcuXIUFBTsdvxdL0za8R+BMmXKbL9dpkyZnT5E3fH+XV+z4/M2bNhAxYoVS/rWRWR/TJwITZuGaeahQ+Gpp1Ia8BDtdE0z4BIzWwS8DpxjZq9GeL7IVapUif79+/P444/vtlqld+/ePProo1StWnWn+8eOHcumTZsA+P777yksLOSoo47a6Tn169dn0aJFkdVdWFhI9erV1adGJFXc4bHHwpRM1arw2WdwzTWxlBJZyLv7Xe5e293rANcAH7t7x6jOlypNmzalcePGDB06dKf7TzzxRK6//vrdnj9u3DgaNmzISSedROvWrenbty9HHnnkTs856KCDOPbYY/n6668jqXnChAlcdNFFkRxbRHbxyy9w+eVwxx1wxRUh4E84Ib563D3yL6AFMKq4551yyim+q3nz5u12XzZ6++23/Z577onk2JdddpkvWLBgt/tz5c9WJGW++ca9Xj33cuXc+/Vz37o18lMC+b6XXE1JF0p3nwhMTMW5MtVll122fTVMadq4cSPt2rWjfv36pX5sEdnBihXQujUUFsKECXDmmXFXBGRAq2EIv23kQkOtrl27lvoxK1SoQKdOnXa7PwwARKRUrF0LbdvCkiXw4YfQrFncFW2X9m0NKlasSGFhoUKpFHmin7xW24iUgs2bw4eq+flhBU0aBTxkwEi+du3aLF26lBUrVsRdSlbZtjOUiOwHd7jpJhg1CgYMgHbt4q5oN2kf8uXLl9fuRSKSnv76V3jxRbjnHujRI+5q9ijtp2tERNLSCy+EkL/hBviv/4q7miIp5EVESmrUqDByb9MGBg6ENF4YopAXESmJadNCL5qTTw5dJdP8SnKFvIhIshYuDEsla9WC0aMj2+ijNCnkRUSS8f334WKnMmVg7Fg4/PC4K0pK2q+uERGJ3erVcOGFsHx56CxZr17cFSVNIS8isjcbN4ZGY59/Du+9B6eeGndFJaKQFxEpinvYrm/8ePj738NqmgyjOXkRkT35+eewk9Orr8KDD4b18BlII3kRkW3cYerUsPZ92DBYvx5uvhnuvjvuyvaZQl5E5Oefw4h94EAoKIAqVcLIvVu3sB4+gynkRSQ3ucP//E8I9jfeCKP2vLzQruCaazJiDXwyFPIiklt+/hleeSWE+9y5WTVq3xOFvIjkhrlzw+baWTxq3xOFvIhkN3d4/nm47TaoUCGrR+17opAXkey1cmUI9LfegvPPD9M0GdKOoLRonbyIZKfPPoOmTeGdd+CRR+D993Mu4EEhLyLZZutWePzxsNfq1q0waRLccUdoLJaDNF0jItnjxx/DnPvo0WG/1UGDoGrVuKuKVW7+0yYi2WfSJGjSJPSZ6d8f3n475wMeFPIikum2bAl7rLZsCQceGNoS3HJLWm/Jl0qarhGRzLVsGXTsCB9/DNddB889Fy5uku0U8iKSmcaNg9/9LmzoMWgQdO6s0fseaLpGRDLL+vXQu3fYiq96dZg+Hbp0UcAXQSN5Eckcs2aF0fvcudCjR1gqWalS3FWlNY3kRST9bd4MffrAb38LhYUwZgwMGKCAT4JG8iKS3r7+Gjp1CqtmrroqhHu1anFXlTE0kheR9LStsViTJjBvXtjUY9gwBXwJaSQvIuln2TLo2jVMy7RqFTbRPvrouKvKSBrJi0h6GT4cGjUKa9+feioslVTA7zOFvIikh19+CXPvV14JderAzJnQs2fONhYrLfrTE5F4/fQTvPtuGL3/4x9w333hQ9YGDeKuLCtoTl5EoucOS5fC/Pnw5Zfh+7av5cvDc+rXhylTwjJJKTUKeREpXStXwsSJOwf5l1/CmjX/95zDDgsj9bZtw/cGDeCcc0KDMSlVCnkRKT0LF4Zt9hYvDrdr14bjjw99ZbaFeYMGYYcmtSFIichC3swqApOAAxLnecvd74/qfCISs5kz4YILws/jxsHpp6sjZBqIciS/ATjH3deYWXlgspm97+7TIjyniMRh4kS45JKwScf48fCb38RdkSREtrrGg22TcOUTXx7V+UQkJiNGhBH8r34VPjhVwKeVSJdQmllZM5sNLAfGu/une3hOdzPLN7P8FStWRFmOiJS2wYPhiiugadOw/d5RR8Vdkewi0pB39y3u3gSoDZxmZg338JyB7p7n7nk1atSIshwRKU2PPQY33gjnnQcffqj9VNNUSi6GcveVwATgglScT0Qi5A5//jPccQe0bx8uZDrooLirkiJEFvJmVsPMDk38fCBwHvBlVOcTkRTYvDk0DuvbF/7wB3jtNahQIe6qZC+iXF1TE3jZzMoS/jF5w91HRXg+EYnS+vVw7bXhg9b77oMHHtBa9wwQWci7++dA06iOLyIptGoVXHppWCrZvz/cckvcFUmSdMWriOzd8uXQpg18/nnYuKNDh7grkhJQF0oR2TP3sCzyzDND/5mRIxXwGUghLyI7W7cOBg0K2+41bx4ajo0fDxdeGHdlsg8U8iIS/O//wl13haZiXbuGkfwLL8CiRdCsWdzVyT7SnLxILnOHTz4JH6aOGBFut2sXPlht3lyrZ7KAQl4kF61bB0OHhnCfMyf0d+/dO6x9P+aYuKuTUqSQF8klS5bA3/4WpmEKC6Fhw/DzdddBpUpxVycRUMiL5Ir33w/NxDZsCGvee/bUlEwOUMiL5ILhw8PVqg0bhp/r1o27IkkRra4RyXZDhsDVV8Opp8LHHyvgc4xCXiSbDRgA118PLVvCBx/AoYfGXZGkmEJeJFtt6xR58cUwahRUrhx3RRIDhbxItnGH++8PPd/btw9z8BUrxl2VxEQfvIpkE/ew3r1fP+jSBQYOhLJl465KYqSRvEi22LIFevQIAd+zZ1j/roDPeUmN5M3scKAZUAtYBxQA+e6+NcLaRCRZmzbBDTfAP/4Bd98NDz6o9e8CFBPyZtYSuBOoCswClgMVgXbAsWb2FvC4u6+KuE4RKcqGDWHufeRIePhhuPPOuCuSNFLcSP5CoJu7/++uD5hZOaAtYe/W4RHUJiLF+fe/4bLLYNw4ePppuPnmuCuSNLPXkHf3/9jLY5uBEaVdkIgkadUqaNsWpkyBwYOhc+e4K5I0lOyc/KFAJ6DOjq9x956RVCUiezd9erjI6auvQjfJq6+OuyJJU8kuoRwDTAO+APRhq0hcNmyAv/4VHnsMataEsWOhVau4q5I0lmzIV3T3P0VaiYjsXX5+WEEzd25YA//EE3DIIXFXJWku2XXyr5hZNzOraWZVt31FWpmIBBs2wF/+AqefDj//DKNHhz1YFfCShGRH8huBvsA9gCfuc+DXURQlIgkzZ4a594KCMIrv109NxqREkg353kA9d/8xymJEJGHjxnBBU58+cPjhocHYRRfFXZVkoGRD/mvg31EWIiIJs2aFUfvnn0OnTvDkk2EPVpF9kGzIrwVmm9kEYMO2O7WEUqQUbdwYRu4PPQTVq8O774Y2wSL7IdmQH4EufBKJTkEBdOwIc+aE7089BVW1tkH2X1Ih7+4vR12ISE5yD+2Ab7sNDj4YRowIm2yLlJKkllCaWVszm2VmP5nZKjNbbWZqSiayP1auDI3FevSAs88Oc/AKeCllyU7XPAlcDnzh7l7Mc0WkONOmwbXXwtKl8OijcPvtUEbbO0jpS/b/qiVAgQJeZD9t3RpaEpx1Vrj9ySdhmz4FvEQk2ZH8n4ExZvZPdl5d80QkVYlkox9+CEsix42DK68MOzfpwiaJWLIh/xCwhrBhSIXoyhHJUh9+GFbN/PILPPccdO+unZskJZIN+Vru3jDSSkSy0aZNcP/98MgjcPzxMH48NGoUd1WSQ5KdCBxjZudHWolItlm8GJo3D1vy3Xhj6AGvgJcUS3YkfxNwu5ltADYBBri7HxxZZSKZ7O23Q7Bv2RI29bjmmrgrkhyV7MVQVaIuRCQrrF0Lf/pTuMApLw9efx2OPTbuqiSH7XW6xszqFPO4mVntUq1IJFPNmhWCfeDAsCxyyhQFvMSuuJF8XzMrA4wEZgArCCts6gEtgVbA/cDSXV9oZkcDQ4AjCL3nB7r7U6VXukia2Lo19Jq5806oVi18uHruuXFXJQIUE/LufpWZnQB0ALoANQkth+cT9n19yN3XF/HyzUBvd59pZlWAGWY23t3nlV75IjH7/vuwqce4caElwYsvhg6SImmi2Dn5RCjfU9IDu/syYFni59VmNh84ClDIS3YYPRo6d4Y1a2DAAPj977X2XdJOSq6lTsztNwU+3cNj3c0s38zyV6xYkYpyRPbPunVwyy3Qti3UqgUzZoQmYwp4SUORh7yZVQaGA7e5+26dK919oLvnuXtejRo1oi5HZP8UFMBpp8Ezz0CvXvDpp9CgQdxViRSp2JBPrKA5el8ObmblCQH/mru/vS/HEEkL7iHY8/JgxQoYOxaeeAIOOCDuykT2qtiQT3SeHFPSA5uZAYOA+WpkJhlt+XK45JIwRXPuuaHve+vWcVclkpRkp2tmmtmpJTx2M+B3wDlmNjvxdWEJjyESD3eYOBF+9zs45piwLPLpp+G99+Dww+OuTiRpybY1+C3QwcwWEzb13tbWoHFRL3D3yYnniWSO776Dl16CwYPh//0/OOSQsIKmZ8/QYEwkwyQb8vrdVLLXpk1hOeSgQTBmTLi4qUULeOABuPxyqFQp7gpF9lmyvWsWm9lJQGI7Gz5x9znRlSWSAgsWhGAfMiRs6FGzJtxxB3TpAvXqxV2dSKlIKuTN7FagG7BthcyrZjbQ3Z+OrDKRKKxfH5qGvfhi6C1TtmxY737jjdCmDZRL9pdbkcyQ7P/RNwK/dfe1AGb2KDAVUMhLZnAPH5r26gXffAP164cNtDt1giOPjLs6kcgkG/IGbNnh9hb0oapkioUL4dZbw9r2Bg3C9/PP1xWqkhOSDfm/A5+a2TuJ2+0Ia+BF0tfq1fDgg9CvHxx4YLh46eaboXz5uCsTSZliQz7RangaMBE4M3F3Z3efFWFdIvvOHf7xD/iP/4Bly+CGG8Ieq0ccEXdlIimXTBfKrWb2rLs3BWamoCaRfTdrVrgydcoUOPVUeOcd+O1v465KJDbJXvH6kZldkWhVIJJ+CgvhpptCb5mFC8PqmWnTFPCS85IN+d8DbwIbzGyVma02s906Soqk3JYtoZd7/frwwgthzn3hwrAkskxKOmmLpLVkulCWAS5w9zLuXsHdD3b3Ku5+cArqE9kzdxg5Ek45Bf7wB2jcOEzVPPUUHHpo3NWJpI1kulBuBZ5JQS0ixdu6Fd56C5o2hXbtwgqaYcPg44+hUaO4qxNJO5qTl8ywZQsMHRpG7FddFXZnevnl0Jrg6qu15l2kCCWdk9+oOXlJqc2bQ2+ZE06A664L0zRDh8K8eeFqVbUhENmrZBuUVYm6EJGdbNwIr7wCffqENgQnnQRvvhm6QuoDVZGkJfW3JbEFYEczuzdx+2gzOy3a0iQnbdgAzz0XVst07QqHHRY+YJ01C668UgEvUkLJ/o35G3AGcF3i9hrg2UgqktzkHta2H3tsWO9es2bo7T59eth6T3PuIvsk6Z2h3P1kM5sF4O4/m1mFCOuSXFJYGHq4v/suNGsWdmZq1UrBLlIKkg35TWZWFnAAM6sBbI2sKskdEydCx45hs+x+/cI2e5qSESk1yf5t6g+8AxxuZg8Bk4E+kVUl2W/zZrj3XjjnHDjooNCC4LbbFPAipSzZ1TWvmdkMoBWhj3w7d58faWWSvRYtCsshp04Nm2T37w+VK8ddlUhWSnqRsbt/CXwZYS2SC954A7p3/7/17tdcE3dFIllNvxtLaqxdC926Qfv2YXem2bMV8CIpoJCX6M2ZE1oADxoEd98NkyZB3bpxVyWSExTyEh33MN9+2mmwahV8+CE89JC23xNJITX+kGj8+GP4UHXUKLj4Yhg8GKpXj7sqkZyjkJfSl58PV1wBP/wATz8Nf/yjLmwSiYmma6R0DR4MZ54ZQn3KlLBTkwJeJDYKeSkdGzbA738ftt0766wwmj/llLirEsl5CnnZf0uXwtlnw8CBcOedMHas5t9F0oTm5GX/TJwYdmZatw6GDw/93kUkbWgkL/vGHZ54As49F6pVg88+U8CLpCGFvJTcmjVw7bXQu3fo9f7pp+EqVhFJOwp5KZmvvoLTTw9b8T3ySJiiOfjguKsSkSJoTl6S9957ofd7uXLhw9Xzzou7IhEphkbyUry1a0Pv90sugXr1YMYMBbxIhtBIXopWUADPPw9DhoTeM9dfDwMGwIEHxl2ZiCRJIS87W78e3norhPvkyVChAlx1VbjQaduVrCKSMRTyEixcGC5m+vvf4aefwrRM375www26sEkkg0UW8mY2GGgLLHf3hlGdR/bDxo0wcmQYtX/0UfhAtV076NEDWrbUfqsiWSDKkfxLwDPAkAjPIftiyRJ47rmwiccPP8Axx4Q+7126wJFHxl2diJSiyELe3SeZWZ2oji/7wD2M2nv3DnPvF10URu2tW0PZsnFXJyIR0Jx8rvj++zBSf/99OP/8EPZ16sRdlYhELPZJVzPrbmb5Zpa/YsWKuMvJTu+8Aw0bwoQJ8Mwz4UImBbxITog95N19oLvnuXtejRo14i4nu6xaFUbvl18eQn3WLO3SJJJjYg95icgnn8BJJ8HLL4erVadOheOPj7sqEUmxyELezIYCU4HjzGypmd0Y1blkBxs3wl13QfPm4cPUyZPhP/8TypePuzIRiUGUq2uujerYUoS5c6FDB5gzB7p1C/3eK1eOuyoRiZGma7LB1q3w5JNhT9XvvgsXOA0cqIAXES2hzHiLF4fNsz/6KHSJfOEFOPzwuKsSkTShkXymWrUK7r4bjjsOpk0L4T5ihAJeRHaikM80mzeHqZjf/AYefjhsoj1/PnTtqqWRIrIbTddkkg8+CC0J5s6Fs86C0aMhLy/uqkQkjWkknwnmzoU2beCCC0LPmeHD4Z//VMCLSLEU8uls+XK46SZo3DjMuz/xBMybF65g1dSMiCRB0zXpaP36sCSyTx9Ytw5uvhnuuw+qVYu7MhHJMAr5dOIOw4bBnXeGpZGXXgqPPQb168ddmYhkKE3XpIuVK8NeqtdeC4cdBh9/HJZEKuBFZD9oJJ8OPvsM2reHpUvDvqq9emkTDxEpFRrJx2nrVnj8cWjWLEzVTJ4Mt9+ugBeRUqORfFx+/BGuvx7GjAmrZV58MUzTiIiUIo3k4zBpEjRpAh9+CM8+C2+9pYAXkUgo5FNpyxZ48EFo2RIqVQpr3//wB615F5HIaLomVZYtg44dw6qZDh1gwACoUiXuqkQkyynkU2H8+BDwq1fD4MFwww0avYtISmi6JkqbN4d2wK1bQ40akJ8PnTsr4EUkZRTyUZk3D84+O7QD7to1rIU/4YS4qxKRHKOQL20bNsD994fVMwsWwNChof97pUpxVyYiOUhz8qVp0iTo3j2Ee4cO0K9fmKYREYmJRvKlYeXKEO7Nm4eR/Nix8OqrCngRiZ1Cfn+4w5tvQoMGMGhQaElQUBA+aBURSQOartlXS5bAH/8I770HJ58ctuI7+eS4qxIR2YlG8iW1ZQs8/XRYKfPRR6HB2KefKuBFJC1pJF8SX3wB3bqFUG/dOly1Wrdu3FWJiBRJI/lkrF8P99wTRuvffAOvvQbvv6+AF5G0p5F8cSZPDhczLVgQWgM//rj2WhWRjKGRfFFWrw4baJ91VlgW+cEH8NJLCngRySgK+T15/3048UT429/g1lvDXPz558ddlYhIiSnkd1RYCJ06wYUXQuXKMGUKPPlk+FlEJAMp5CFc1DRsWLioaehQuPdemDULzjgj7spERPaLPnj99tuwO9O770JeXtiSr3HjuKsSESkVuTuSd4cXXggXNY0fD//93zB1qgJeRLJKbo7kFy6EHj1gwgRo0SKEfb16cVclIlLqcmskv25dmG9v1Ahmzgx93j/+WAEvIlkrd0byY8aEde//+lfYb7VvXzjyyLirEhGJVPaP5Jcsgcsvh4suggMOCCP3V15RwItITsjekN+0KXyY2qBB2MSjTx+YMwdatoy7MhGRlMnO6ZrJk+Gmm8IGHhdfDP37Q506cVclIpJykY7kzewCM1tgZl+b2Z1RnguAFSugc+fQb2bVKhgxIqx/V8CLSI6KLOTNrCzwLNAGOAG41sxOiORkW7eGlTLHHRf2Vr3zTpg3Dy69NJLTiYhkiiina04Dvnb3bwDM7HXgUmBeqZ7l55+hTZuwkUfz5qGp2AnR/FsiIpJpopyuOQpYssPtpYn7dmJm3c0s38zyV6xYUfKzHHooHHssDBkSLm5SwIuIbBf7B6/uPhAYCJCXl+clPoBZ2KlJRER2E+VI/lvg6B1u107cJyIiKRJlyE8HfmNmdc2sAnAN8G6E5xMRkV1ENl3j7pvN7GbgA6AsMNjd50Z1PhER2V2kc/LuPgYYE+U5RESkaNnb1kBERBTyIiLZTCEvIpLFFPIiIlnM3Et+/VFUzGwFsHgfX14d+LEUy8kkufzeIbffv9577tr2/o9x9xpFPSmtQn5/mFm+u+fFXUcccvm9Q26/f7333HzvkPz713SNiEgWU8iLiGSxbAr5gXEXEKNcfu+Q2+9f7z13JfX+s2ZOXkREdpdNI3kREdmFQl5EJItlfMinfLPwNGJmg81suZkVxF1LqpnZ0WY2wczmmdlcM7s17ppSycwqmtlnZjYn8f7/GndNqWZmZc1slpmNiruWVDKzRWb2hZnNNrP8Yp+fyXPyic3CFwLnEbYXnA5c6+6lu49smjKzs4E1wBB3bxh3PalkZjWBmu4+08yqADOAdjn0396Ag9x9jZmVByYDt7r7tJhLSxkz+xOQBxzs7m3jridVzGwRkOfuSV0Ilukj+e2bhbv7RmDbZuE5wd0nAT/FXUcc3H2Zu89M/LwamM8e9hDOVh6sSdwsn/jK3BFbCZlZbeAi4MW4a0l3mR7ySW0WLtnNzOoATYFPYy4lpRLTFbOB5cB4d8+l9/8k8Gdga8x1xMGBcWY2w8y6F/fkTA95yXFmVhkYDtzm7qvirieV3H2Luzch7J98mpnlxJSdmbUFlrv7jLhricmZ7n4y0Ab4Y2LatkiZHvLaLDyHJeaihwOvufvbcdcTF3dfCUwALoi5lFRpBlySmJt+HTjHzF6Nt6TUcfdvE9+XA+8Qpq2LlOkhr83Cc1Tig8dBwHx3fyLuelLNzGqY2aGJnw8kLD74MtaiUsTd73L32u5eh/B3/mN37xhzWSlhZgclFhpgZgcB5wN7XV2X0SHv7puBbZuFzwfeyKXNws1sKDAVOM7MlprZjXHXlELNgN8RRnGzE18Xxl1UCtUEJpjZ54TBznh3z6mlhDnqCGCymc0BPgNGu/vYvb0go5dQiojI3mX0SF5ERPZOIS8iksUU8iIiWUwhLyKSxRTyIiJZTCEvGcfMXjSzE1J8zieLu7LQzHqYWad9PP4iM6u+l8dfN7Pf7MuxJbdpCaVIMcysGmE98ukRnmMRe+ksaGbNgY7u3i2qGiQ7aSQvaStxdd/oRM/0AjNrn7h/opnlJX6+0cwWJnqrv2BmzyTuf8nMBpjZNDP7xsxaJPrvzzezl3Y4xwAzyy+mJ/sVwNgdXrPIzB5L9PT+zMzqJe5/wMxuN7NyZjbdzFok7n/YzB5K/Nwx8ZrZZvZ8ol12se8Z+AQ418zK7fcfrOQUhbykswuA79z9pES//J2u7DOzWsC9wOmEK2CP3+X1hwFnAL0I7S76AScCjcysSeI597h7HtAYaG5mjfdQRzNCv/od/eLujYBnCB0Rt0tciX0DMMDMzk28j7+aWQOgPdAs0VhsC9Ahmffs7luBr4GT9lCfSJEU8pLOvgDOM7NHzewsd/9ll8dPA/7p7j+5+ybgzV0ef8/DfOQXwA/u/kUiLOcCdRLPudrMZgKzCP8A7GmuvyawYpf7hu7w/YxdX5Bor/EKMAroktjvoBVwCjA90SK4FfDrErzn5UCtPdQnUiSFvKQtd18InEwIvgfN7L4SHmJD4vvWHX7edrucmdUFbgdauXtjYDRQcQ/HWbeH+72In3fUCFgJHJ64bcDL7t4k8XWcuz+w00H3/p4rJmoRSZpCXtJWYjrm3+7+KtCXEH47mk6YYjksMVd9RQlPcTCwFvjFzI4g9Ofek/lAvV3ua7/D96l7qP1yoCpwNvB0omPkR8CVZnZ44jlVzeyYXV63t/dcn2I6DorsSh/iSDprBPQ1s63AJuCmHR9092/NrA+hG99PhFa7u07pFMnd55jZrMTrlgBTinjqaOD37LzV3GGJDpAbgGt3fHJiKeQjhN8QliQ+DH7K3a83s78QdvUpk3hPfwQWF/eeE/8IrXP375N9fyKgJZSS4cyscmIz63KEDRQGu/s7EZxnMtDW3VcWt9wxCmbWC1jl7oNSdU7JDpqukUz3QOJDzALgX8CIiM7TG/hVRMdOxkrg5RjPLxlKI3kRkSymkbyISBZTyIuIZDGFvIhIFlPIi4hkMYW8iEgW+/+hqT24NBzjWwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "left_image_point = np.zeros((2))\n", "right_image_point = np.zeros((2))\n", "number_samples = 1000\n", "x_values = []\n", "y_values = []\n", "for sigma_counter in range(0, 25, 1):\n", " sigma = sigma_counter / 5\n", " rms = 0\n", " for i in range(number_samples):\n", " left_image_point[0] = gold_left_image_point[0] + random.normalvariate(0, sigma)\n", " left_image_point[1] = gold_left_image_point[1] + random.normalvariate(0, sigma)\n", " right_image_point[0] = gold_right_image_point[0] + random.normalvariate(0, sigma)\n", " right_image_point[1] = gold_right_image_point[1] + random.normalvariate(0, sigma)\n", " left_camera_point = triangulate_points_to_3d(left_image_point, \n", " left_intrinsics, \n", " left_distortion, \n", " right_image_point, \n", " right_intrinsics, \n", " right_distortion, \n", " left_to_right)\n", " left_world_point = convert_camera_point_to_world(left_camera_point,\n", " camera_to_marker,\n", " marker_to_world\n", " )\n", " diff = np.linalg.norm(gold_world_point - left_world_point)\n", " diff = diff * diff\n", " rms = rms + diff\n", " rms = rms / number_samples\n", " rms = np.sqrt(rms)\n", " x_values.append(sigma)\n", " y_values.append(rms)\n", " \n", "plt.plot(x_values, y_values, 'r', label='RMS (mm)') \n", "plt.legend(loc='upper left')\n", "plt.xlabel('sigma (pixels)')\n", "plt.ylabel('error (mm)')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Experiment 1: Discussion, Conclusion\n", "\n", "So as noise increases, so does the triangulation error in millimetres. So, we know that our feature detector for triangulation purposes must be sub-pixel accuracy to get triangulation accuracy around 1mm.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Experiment 2: Accuracy of stereo calibration\n", "\n", "What is the effect of inaccuracies in our stereo calibration?\n" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[-0.18663586637738608, -1.022194275896184, -0.31552983713490357, -6.05146, 0.058148, -0.143609]\n" ] } ], "source": [ "print(left2right)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "So, rotation about y axis, or translation in x-axis are of interest.\n" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEGCAYAAACJnEVTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAysElEQVR4nO3dd3hU1dbH8e+iht6rlNAsSCcgxUJHFBEpomABQZBLEURAinQRFFRAVHpEURABQQSkN+lI6EXahdB7SSgp+/1jD77IRRhgZs5MZn2eh8fJZGbOOkZ/Oeyz99pijEEppVRwSOR0AUoppXxHQ18ppYKIhr5SSgURDX2llAoiGvpKKRVEkjhdwM0yZ85sQkNDnS5DKaUCxsaNG08bY7K4+3q/Cv3Q0FA2bNjgdBlKKRUwROS/9/J6Hd5RSqkgoqGvlFJBRENfKaWCiF+N6d9OTEwMkZGRXL161elS7llISAi5cuUiadKkTpeilFJAAIR+ZGQkadKkITQ0FBFxuhy3GWM4c+YMkZGR5MuXz+lylFIKCIDhnatXr5IpU6aACnwAESFTpkwB+TcUpVTC5fehDwRc4N8QqHUrpRKugAh9pZRKsFavhk8/9dnhNPTdkDhxYkqUKEGRIkV44YUXOH/+vNMlKaUCnTHw1VfwzDMwahRcvuyTw2rouyFFihRERESwbds2MmbMyMiRI50uSSkVyKKjoWlTaNMGatSA9eshdWqfHNqroS8iHUVku4hsE5EfRSTEm8fzhfLly3PkyBH27dtHqVKl/n7+r7/++sfXSil1W/v3Q4UK8N130LcvzJoFGTL47PBem7IpIg8B7YHCxpgrIvIT8AoQft8f2qEDRER4orz/V6IEfPGFWy+Ni4tj0aJFNG/enAIFCpAuXToiIiIoUaIEEyZMoFmzZp6tTSmVsMyZA02a2MezZ8Nzz/m8BG8P7yQBUohIEiAlcNTLx/OKK1euUKJECbJnz86JEyeoXr06AC1atGDChAnExcUxZcoUGjdu7HClSim/FB9vr+pr14bQUNi40ZHABy9e6RtjjojIEOAQcAWYb4yZf+vrRKQl0BIgT548d/5QN6/IPe3GmH50dDQ1a9Zk5MiRtG/fnvr169O3b1+qVKlC6dKlyZQpkyP1KaX82Llz8Npr9ir/jTfg668hZUrHyvHalb6IZABeBPIBOYFUIvLara8zxow2xoQZY8KyZHG7JbQjUqZMyfDhwxk6dCixsbGEhIRQs2ZNWrdurUM7Sqn/FREBYWGwYIGdqRMe7mjgg3eHd6oBB4wxp4wxMcB0oIIXj+cTJUuWpFixYvz4448ANGnShESJElGjRg2HK1NK+ZXvvoPy5eHqVVi2DFq3Bj9YsOnN3juHgHIikhI7vFMVCMgdUi7fMn/2119//fvxypUradasGYkTJ/Z1WUopf3T9OnTsaK/sK1WCyZMhWzanq/qbN8f014rIz8CfQCywCRjtreM54aWXXmLfvn0sXrzY6VKUUv7gxAl46SW7yrZTJxg0CJL4V19Lr1ZjjOkN9PbmMZw0Y8YMp0tQSvmLgwehenU4cgSmTIGXX3a6otvyr19BSikViHbssCtro6Jg4UK7+MpPaRsGpZR6EOvWwVNPQVycvWHrx4EPGvpKKXX/Fi2CKlUgXTpYuRKKFXO6orvS0FdKqfsxY4ZdVZsvnw38AgWcrsgtOqZ/B2fOnKFq1aoAHD9+nMSJE3NjAdm6detIliyZk+UppZwyfjy8/TaULQu//QYZMzpdkds09O8gU6ZMRLgavPXp04fUqVPz/vvvO1uUUspZn31mp2NWr26v9lOlcrqie6LDO/dh/fr1FCtWjKtXrxIVFcXjjz/Otm3bnC5LKeVNxkDPnjbwGzSAX38NuMCHALvS7zCvAxHHIzz6mSWyl+CLZ7+4p/eUKVOGOnXq0LNnT65cucJrr71GkSJFPFqXUsqPxMVB27bwzTfQooX9Z4Cuwg+o0PcnvXr1okyZMoSEhDB8+HCny1FKecv167Y75pQp0LUrfPyxX/TQuV8BFfr3ekXuTWfOnOHy5cvExMRw9epVUgXgX/OUUncRHQ3168O8eTB4MHTp4nRFD0zH9O9Tq1at6N+/P02aNKFr165Ol6OU8rQ1a+ym5fPnw5gxCSLwIcCu9P3FxIkTSZo0KY0bNyYuLo4KFSqwePFiqlSp4nRpSqkHtW2bvWE7c6btjvnzz7aJWgIhxhina/hbWFiY2bDhn92Xd+7cyWOPPeZQRQ8u0OtXKmgcOAC9e8P330OaNPbK/t13IXVqpyu7IxHZaIwJc/f1eqWvlApuJ07ARx/9/4yczp3tDdsAWnB1LzT0lVLB6cIFGDIEPv/c7m7VvDn06gUPPeR0ZV4VEKFvjEECcIqUPw2dKaVcrlyBkSPt1MuzZ6FRI+jXDx5+2OnKfMLvZ++EhIRw5syZgAtQYwxnzpwhJCTE6VKUUgAxMXYWTqFCdginbFnYuNFuZxgkgQ8BcKWfK1cuIiMjOXXqlNOl3LOQkBBy5crldBlKBa+4ONsBc8oUmDYNTp60m5VPmmSnYwYhvw/9pEmTki9fPqfLUEoFivh4O8d+yhSYOhWOHYMUKaB2bXjzTdsOOQCHiz3F70NfKaXuyhjYsMEG/U8/weHDkDw51Kplx+xr1/b7qZe+oqGvlApMxsDmzf8f9Pv3Q9Kkdq/ajz6CF1+EtGmdrtLvaOgrpQLPunV2qGbXLju3vmpV6NHDrpzNkMHp6vyahr5SKrBs2gQ1a0L69HZBVb164NrRTt2dhr5SKnBs3253rEqbFpYuhbx5na4o4Pj9PH2llAJgzx47jJMsGSxapIF/n/RKXynl/w4csIEfHw9LlkDBgk5XFLA09JVS/i0y0gZ+VJQNfO1a+0A09JVS/uv4cRv4Z87AwoVQvLjTFQU8DX2llH86fdretI2MtLtXlSnjdEUJgoa+Usr/nD9vF1nt3Qu//QYVKzpdUYKhoa+U8i+XLtn2Cdu22S0LdRtSj9LQV0r5j+ho2ydn/XrbLK1WLacrSnA09JVS/uHqVahb17ZCnjQpQW1G7k809JVSzrt+HRo2hAULYMIEeOUVpytKsDT0lVLeEx8PFy/a/WjPn7f/vPnPjefWr7dtFb76Cpo2dbbmBE5DXynlOdeuQadOMGuWDfOLF+/+nuTJbWfMESOgdWvv1xjkNPSVUp5x5AjUrw9r19p/5soF6dLZP+nT///jW79OntzpyoOKhr5S6sH98Qc0aACXL9u9aOvVc7oi9S+0y6ZS6sGMHg2VK9vtCNes0cD3c14NfRFJLyI/i8guEdkpIuW9eTyllA9dvw7vvAOtWtn+OOvWweOPO12VugtvD+8MA+YZYxqISDIgpZePp5TyhWPH7HDOqlXwwQcwYIDdtlD5Pa9d6YtIOuBpYByAMea6Mea8t46nlPKRtWshLAwiIuym5B9/rIH/ALac2MLEzRN9djxvDu/kA04BE0Rkk4iMFZFUXjyeUsrbJkyAp5+2M25Wr4aXX3a6ooAUGx/Lzzt+5pnwZyj+TXHe+/09rsVe88mxvRn6SYBSwNfGmJJAFPDBrS8SkZYiskFENpw6dcqL5Sil7ltMDLRtC2+9BU89ZRdTFSvmdFUB51TUKQauGEi+YfloOLUhhy4c4tPqn7Kn3R6SJ/HN1FVvjulHApHGmLWur3/mNqFvjBkNjAYICwszXqxHKXU/Tp60LRKWL7cLrwYNgiQ62/tebDi6gS/XfcnkbZO5FneNavmrMfK5kTxf6HkSJ/Lt0JjXfnLGmOMiclhEHjHG7AaqAju8dTyllBesXAmNG8OpU/D999CkidMVBYzrcdf5ecfPjFg3gjWRa0iVNBXNSzanbdm2PJbFuS0fvf3ruh0wyTVzZz/QzMvHU0p5QnQ09OgBw4ZBaKhdfFWqlNNVBYRjl44xauMoRm0cxfHLxymYsSBf1PyCpiWaki4kndPleTf0jTERQJg3j6GU8rAVK+zY/d690KaNHc5JndrpqvzaqahTzNw9k593/MyiA4uIjY+lVsFatCvbjpoFa5JI/GcdrA7MKaWs6Gjo3h2GD7dX94sX25W26raOXTrGjF0zmLZzGksPLiXexJM/Q37eK/ceLUq1oFCmQk6XeFsa+kope3XfrBns22dn6Xz8sV7d38ahC4eYvnM603ZO449Df2AwPJLpEbo92Y0GhRtQPFtxRMTpMu9IQ1+pYBYVZa/uR4ywV/dLlkClSk5X5Vf2nd3HtJ3TmLZzGuuOrAOgaNai9KnUh/qP1adwlsJ+H/Q309BXKlgtX27H7vXq/n+cv3qeH7f+yLhN49h4bCMApXOU5uOqH1P/sfp+O3TjDg19pYJNVBR062av7vPntztWPfOM01U5zhjDykMrGbtpLFO3T+VK7BWKZyvOkOpDqPdYPfJlyOd0iR6hoa9UMFm+3I7d798P7drZq/tUwd0d5WTUSSZunsjYP8ey+8xu0iRLwxvF3+DtUm9TKkepgBq6cYeGvlLBYsoUu9AqNDTor+7j4uNYsH8BY/8cy8zdM4mNj6Vi7op88OQHNCzckFTJEu4vQg19pYLB1Kl2NW3FivDbb5AmjdMVOeLQhUNM2DSB8RHjOXThEJlSZKJ92fY0L9WcwlkKO12eT2joK5XQTZ8Or74K5crBnDlBebP26KWjvD//fSZvm4zBUD1/dYZUH0KdR+r4rNGZv9DQVyoh++UXaNQIypaFuXODLvBj42MZsXYEvZf25nrcdbpU7EKr0q0SzE3Z+6Ghr1RCNWuW7XcfFgbz5gXdkM7KQyv5z2//YevJrdQqWIsRtUZQIGMBp8tynIa+UgnR7Nl2O8OSJW3gp03rdEU+czLqJF0XdiU8IpzcaXMz/eXp1H20boKbhXO/NPSVSmjmzIH69aF4cfj9d0jnfGdHX4iLj2P0xtF0X9ydy9cv07ViVz58+sMEPRPnfmjoK5WQ/P471KsHRYrA/PmQPr3TFfnEhqMbaP1bazYc3UDl0MqMfG6koz3r/ZmGvlIJxfz58OKL8NhjsGABZMjgdEVed+7KObov6s6ojaPIljobk+pN4tUir+pQzh24FfoikhWoCOQErgDbgA3GmHgv1qaUctfChTbwH3nEPs6Y0emKvCouPo7vtnxHlwVdOHPlDO2faE/fSn39YpMSf3fH0BeRyth9bTMCm4CTQAhQFyggIj8DQ40xF71cp1Lq3yxeDHXqQKFCsGgRZMrkdEVec/zyccZvGs/ojaP574X/Uj5XeeY/P58S2Us4XVrAuNuV/nPA28aYQ7d+Q0SSALWB6sA0L9SmlLqbZcugdm3bOG3RIsic2emKPM4Yw9KDS/lm4zdM3zmd2PhYquSrwpAathGaP+1KFQjuGPrGmM53+F4s8IunC1JKuSEmBn79FV5/3fbSWbQIsmRxuiqPOnvlLBM3T+SbDd+w+8xuMoRkoF3ZdrQq3YpHMj/idHkBy90x/fTAG0Doze8xxrT3SlVKqf915Yq9QTttmg38c+fsTdvFiyFbNqer8whjDGuPrOWbDd8wZfsUrsZepVyucnxb91saFm5IiqQpnC4x4Lk7e2cOsAbYCujNW6V85dIlO+9+2jT7z6goOw3zhRfsXPwaNSBF4AfhpWuX+GHrD3yz8RsijkeQOllqmhZvSquwVjpe72Huhn6IMeY9r1ailLLOnLFX8tOn22mY165B1qzw2mt2Dn6lSpAsmdNVeoQxhh+2/kD7ee05e+UsxbIV4+vnv6ZJ0SakSR5cbSN8xd3Q/05E3gZmA9duPGmMOeuVqpQKNtHR8O23NuiXLIG4OMiTB1q3tkFfoQIkTux0lR51/PJx3pn9DjN3z6RcrnIMrTGU8rnK6xx7L3M39K8DnwI9AON6zgD5vVGUUkElMtJOudy0CR5+GLp0sUFfujQkwAC8cXXfbm47omOi+bT6p3Qs15HEiRLWLzV/5W7odwIKGmNOe7MYpYLO+vV2UdXly3ZI5/nnE2TQ33Dr1f2EFyfwaOZHnS4rqLgb+nuBaG8WolTQmToV3ngDsme3Y/dFijhdkdfo1b3/cDf0o4AIEVnCP8f0dcqmUvfKGBgwAHr1stsXTp9ub9QmUHp171/cDf1f0IVYSj24q1eheXP44Qe7sGrMGEieMLfr06t7/+RW6BtjvvV2IUoleMePQ926sHYtDBwIH3yQYMfv9eref7m7Irc20B/I63qPAMYYEzzb8Sj1ILZssT1yzpyxC63q1XO6Iq+4FnuN8Ihwui3qplf3fsrd4Z0vgHrAVmOMuctrlVI3+/VXePVVu5J2xQooVcrpijzu0rVLjN44ms/WfMbRS0epkLsC4+qM06t7P+Ru6B8GtmngK3UPjIGhQ+28+9KlYeZMyJnT6ao86nT0aYavHc6X677k3NVzVMlXhfAXw6mWv5ousvJT7oZ+F2COiCzjn7N3PvNKVUoFuuvX7Wra8eOhYUMID4eUKZ2uymMOXzjM0NVDGfPnGKJjoqn7aF0+qPgBT+R6wunS1F24G/ofAZexG6gkjKYfSnnL+fP2hu2yZfDhh9CnDyRKGD3fd53exeA/BvP9lu8BaFy0MV0rdqVwlsIOV6bc5W7o5zTGJNyVI0p5yokTULMm7NgB338PTZo4XZFHrD+ynkF/DGLGzhmEJAmhdVhrOpXvRN70eZ0uTd0jt1sri0gNY8x8r1ajVCA7eBCqV4ejR2H2bNv2OIAZY1i4fyGD/xjMogOLSB+Snh5P9aD9E+3JkiphbdgSTNwN/dbA+yJyDYhBp2wq9U87dtjAj462G5OXL+90Rfftetx1pmybwpDVQ9hyYgvZU2fnk2qf0CqsFWmT6//ygc7dxVna2Fqpf7NuHdSqZXvcL18ORYs6XdF9uXD1AqM3jmbY2mEcuXSEwlkKM77OeBoXbUzyJAlz1XAwumPoi0ioMebgHb4vwEPGmEhPF6ZUQFi82HbJzJLFXuHnD7xu44cvHGbY2mGM3jiaS9cvUSVfFca8MIZnCz6r0y4ToLtd6X8qIomAmcBG4BR2Bk9BoDJQFegNaOir4PPLL9Coke2B//vvATcHf9OxTQxdPZQp26dgjKFRkUZ0Kt+JUjkS3uIx9f/uGPrGmIYiUhhoArwF5MC2WN6J3Tf3I2PM1Tt9hogkBjYAR4wxtT1StVJOCw+3jdPKloXffoOMGZ2uyC3GGH7f9ztDVg1h0YFFpE6WmnZl29GhXAfypMvjdHnKB+46pm+M2YHdMet+vYv9JaF3gFTC8MUX0LGjvXE7fTqkTu10RW7ZfnI7r894nU3HN5EzTU4GVxtMy9ItSR+S3unSlA+5O3vnvohILuB57OIu3VhdBTZjbA/8AQOgQQM7Dz9A2iKHR4Tzn9/+Q9rkaQl/MZxXi75KssS6zjIYeTX0sY3augD/OvtHRFoCLQHy5NG/Xio/FR8P7dvDyJF2WGfUqIDYqDzqehRt5rTh283fUiVfFSbVm0T21NmdLks56K5rw8XKfa8f7GrHfNIYs/FOrzPGjDbGhBljwrJk0QUfyg/FxNgNT0aOhM6d7cYnARD4209up8yYMkzcPJE+z/Rh/mvzNfDV3UPf1Vlzzn18dkWgjogcBCYDVUTk+/v4HKWcs2wZVK5sd7r6+GP45JOA2PgkPCKcMmPKcPbKWRa+sZDelXprT3sFuBH6Ln+KSJl7+WBjTDdjTC5jTCjwCrDYGPPavRaolCOWLoVKleyfvXvh22/tTld+Lup6FE1/aUqzmc0on7s8Ee9EUCVfFafLUn7E3TH9J4AmIvJf7CbpN9owFPNaZUr5mjGwZAn07WtX1mbPDp9/Di1bBkRb5O0nt9NwakN2nd5Fn2f60PPpnnp1r/6Hu6Ff80EOYoxZCix9kM9QymuMgUWLbNivXGkXWQ0bBm+/DSlSOF2dWyZsmkCbOW1ImzwtC15fQNX8VZ0uSfkpd3vv/FdEigNPuZ5aYYzZ7L2ylPIBY2DBAhv2q1bBQw/BiBHQogWEhDhdnVt0do66V26N6YvIu8AkIKvrz/ci0s6bhSnlNcbAvHlQoYLtfX/okJ2Zs28ftG0bEIF/NfYqSw4s0dk56p65O7zTHHjCGBMFICKDgdXACG8VppTHGWObovXsaTtj5s4NX38NzZr5/SKro5eOsvrwalYdXsWqyFVsPLqRmPgYsqXKpsM56p64G/oCxN30dZzrOaUCw/r1dvbN4sWQJ49dXNW0qW2H7Gdi42PZemLr3wG/6vAqDp4/CEDyxMkp81AZOpbrSIXcFagUWol0IemcLVgFFHdDfwKwVkRmuL6uC4zzSkVKedLu3dCjB0ybBpkz274577zjd1f2m45tYvrO6ayKXMXayLVExUQBkDNNTirkrkD7su2pkLsCJXOU1PYJ6oHcNfRdrZXXYGffPOl6upkxZpMX61LqwURG2hu0EybYGTi9e0OnTpDGf/YDOnflHD9s/YFxm8ax6fgmEktiSmQvwVsl36JC7gpUyF2B3Glza0975VHudNmMF5GRxpiSwJ8+qEmp+3f2LAwaZGfhxMVBmzb2Sj9rVqcrAyDexLPs4DLGbRrHtJ3TuBp7lRLZSzCi1ggaF21MxhSB0aJZBS53h3cWiUh9YLqrLYNS/iU62s6tHzwYLl6E116Dfv0gNNTpygA4cvEI4RHhjI8Yz/5z+0mXPB1vlXiL5qWa66YlyqfcDf1W2NbIsSJyFd0YXfmLmBgYN84G/LFjULs2DBzoF/vUxsTFMHvPbMZuGsu8vfOIN/FUDq1Mv0r9qPdYPVIkDYyFXyphcXdM/1ljzB8+qEcp91y/DpMm2YDfuxcqVoSffoInn7z7e70s6noUA5YPYHzEeE5GnSRnmpx8UPED3ir5FgUyFnC6PBXk3B3T/xIo6YN6lLqzK1fslf0nn8Dhw1CiBMyaZa/w/eCG59YTW2n0cyN2nd5F3Ufr0rxkc2oWrEmSRN7eukIp9+iYvgoMFy/ahVSffQYnT9or+1Gj4Nln/SLsjTGM3jiaDr93IH1Ieha+sVC7Wyq/dK9j+nEicgUd01e+cvq0vUH75Zdw/rxtm9C9Ozz9tNOV/e381fO0/LUlU3dMpWaBmkx8aSJZU/nHbCGlbuVuwzX/mdysgsORIzB0qL2aj46GevWgWzcIC3O6sn9YG7mWV6a9QuTFSAZXG8z7Fd4nkbi7TYVSvudW6ItdHdIEyGeM6e/aPjGHMWadV6tTwWf/fjvtMjzczrNv3Ni2Tyhc2OnK/iHexDN01VC6L+7OQ2keYkWzFZTLVc7pspS6K3eHd74C4oEqQH/gMjASuKfdtJT6V9euQYcOMHo0JEkCb70FXbpAvnxOV/Y/TkWd4s1f3mTu3rnUf6w+Y+uMJX1IeqfLUsotbu+cZYwpJSKbAIwx50REG4Aozzh92g7frFgB7drZK/ucOZ2u6raWHFhCk+lNOHvlLF899xXvhL2jbRJUQHE39GNEJDFgAEQkC/bKX6kHs3s3PP+87ZXz44/wyitOV3RbsfGx9F/Wn/7L+/NwpoeZ22QuxbMXd7ospe6Zu6E/HJgBZBWRj4AGQE+vVaWCw5Il9go/aVL7uHx5pyu6rQPnDtB0ZlOW/3c5TUs05ctaX5IqWSqny1Lqvrg7e2eSiGwEqmKna9Y1xuz0amUqYRs3zrY4fvhhmD3br8buo2OiWfHfFczfN5/5++ez7eQ2UiVNxXcvfcdrxV5zujylHojbywSNMbuAXV6sRQWD+Hg79fKTT6BGDds6IZ2zm4DEm3i2nNhiQ37ffFYeWsm1uGskT5ycp/I+xRvF3qDh4w0JTR/qaJ1KeYKuDVe+Ex0Nr78O06fbq/wRI+xMHQccu3SMBfsXMH/ffBbsX8DJqJMAFMlahDZl2lCjQA2eyvsUKZOmdKQ+pbxFQ1/5xtGjUKcO/PknfP45vPuuz9snGGOYuHkiQ1YPYdvJbQBkTZWV6vmrU6NADarlr0bONP45a0gpT9HQV963ebNtiHbuHMycCS+84PMSDp4/SKvZrZi/bz6lc5RmcLXB1ChQg2LZiukKWhVUNPSVd82ebadhpk8PK1farpg+FBcfx8j1I+m+qDsiwpe1vqR1mdYa9Cpoaegr7zAGhg+H996zQf/rrz5fcLXj1A5azGrB6sjVPFvwWUbVHkWedHl8WoNS/kZDX3ne0aPQsaOdmVO3Lnz/PaTy3bz263HXGbxyMANWDCB1stR899J3NCnaRFfOKoWGvvKk2FgYORI+/NDubNW/v22DnMh3Qynrj6yn+azmbD25lVeKvMKwZ4dpm2OlbqKhrzxjzRpo3RoiImzP+y+/hIIFfXb46Jhoei3pxedrPid76uzMfGUmdR6p47PjKxUoNPTVgzl71i62GjMGcuSAqVOhfn2fTsdccmAJLX5twf5z+2lVuhWDqw0mXYizC76U8lca+ur+GAPffgudO9upmB07Qp8+kMY3++3ExMWw6vAqwjeHEx4RTsGMBVny5hIqhVbyyfGVClQa+urebd9uh3JWrLBN0r7+Gop7v+Pk4QuHmbd3HnP3zmXh/oVcun6JpImS8n759+lbua+unlXKDRr6yn1RUdCvn92cPG1aGDsWmjXz2o3aa7HXWHloJXP3zmXe3nlsP7UdgNxpc/NqkVepVagWVfJVIW1y3apZKXdp6Cv3/PILtG8Phw/bXa0GD4bMmT1+mAPnDvwd8osPLCYqJopkiZPxdN6naVaiGbUK1eKxzI/p9Eul7pOGvrqzo0ehbVuYMQOKFrUbnVSs6PHDrIlcQ4d5HVh7ZC0A+dLn483ib1KrUC0qhVYidbLUHj+mUsFIQ1/dXny87XnfubPdv3bQILu6NmlSjx7mxOUTdFvUjQkRE8iZJief1fiM5x9+nkIZC+nVvFJeoKGv/tdff0HLlrB0KVSubDcr9/Cc+9j4WEauG0mvpb24EnOFLhW60PPpnqRJ7pvZP0oFKw199f9iYuxN2j59IHlyO/e+eXOPz7lfenAp7ea2Y9vJbdQoUIPhzw7nkcyPePQYSqnb09BX1p9/2oCPiLCLq0aMsIutPCjyYiSdF3Rm8rbJ5E2XlxmNZvDiIy/qMI5SPuS10BeR3MBEIBtggNHGmGHeOp66T9HR9sr+s88ga1a7q9VLL3n0ENdir/H5ms8ZsHwAcSaO3s/0pmvFrqRImsKjx1FK3Z03r/RjgU7GmD9FJA2wUUQWGGN2ePGY6l4sXmzH7vftg7fftvvWpk/v0UPM2zuP9nPb89fZv6j7aF0+q/EZ+TL4zyboSgUbr7U/NMYcM8b86Xp8CdgJPOSt46l7cO4ctGgBVava8frFi+3NWg8G/p4ze6g7uS61JtUCYG6TucxoNEMDXymH+WRMX0RCgZLA2tt8ryXQEiBPHt3gwutmz7ZX9ydPQteu0Ls3pPDcMMuu07sYsHwAP277kRRJUjCo6iA6lOtA8iTJPXYMpdT983roi0hqYBrQwRhz8dbvG2NGA6MBwsLCjLfrCVrnz9umaOHhdpHV7NlQqpTHPn7HqR30X96fKdumkCJpCt4r9x7vV3ifbKmzeewYSqkH59XQF5Gk2MCfZIyZ7s1jqTv4/Xc7nHPsGPToAb16QbJkHvnobSe30X95f6Zun0rKpCnpUrELncp3IkuqLB75fKWUZ3lz9o4A44CdxpjPvHUcdQeXLkGnTna+/WOP2Zk5Zcp45KO3nNhCv2X9mLZzGmmSpaHbk93oWL4jmVN6vh+PUspzvHmlXxF4HdgqIhGu57obY+Z48ZjqhsWLbWO0w4ehSxfo2xdCQh74Yzcd20T/5f2ZsWsGaZOn5cOnP6RDuQ5kTJHRA0UrpbzNa6FvjFkJ6KobX7t8GT74wO5VW6gQrFxpe94/oI1HN9JveT9m7Z5FuuTp6P1Mb9594l0ypMjggaKVUr6iK3ITkhUroGlTOHAAOnSAjz6ClPe/schfZ/5iyvYp/LT9J7ae3EqGkAz0q9SP9k+01+0IlQpQGvoJQXS0vUE7bBjky2cbpT399H191IFzB/4O+k3HNwHwZJ4nGVFrBG8Uf0M3LFEqwGnoB7qlS6FVK9izB/7zH7u5Sep76z1/6MIhpm6fypTtU1h/dD0A5XKV4/Oan9OgcANypc3lhcKVUk7Q0A9Uu3fbG7SzZkHevLBwoV1h66ajl47+HfSrI1cDEJYzjE+qfcLLj79M3vR5vVW5UspBGvqB5tQpOxPnm2/seP3AgXb83s1VteuOrOPDJR+yYN8CDIbi2YozsMpAXn78ZQpkLODd2pVSjtPQDxRXrsDw4Tbko6LskE7v3rYzpht2ntpJj8U9mLFrBllSZqH3M715pcgr2sdeqSCjoe/v4uPtvrTdu8OhQ/DCC7Yb5qOPuvX2QxcO0XdpX8I3h5MqaSr6VepHh3IddIcqpYKUhr4/W7bMrqjduNH2yQkPt9sXuuF09Gk+XvExI9ePxGB494l36f5Ud10xq1SQ09D3R7t32w6YM2dC7tzw3XfQuDEkunsn7MvXL/P56s/5dNWnRMVE8WbxN+lTqQ950mkHU6WUhr5/OXcOPvzwvm7SXou9xuiNoxmwYgAno07y0qMvMaDKAApnKez9upVSAUND3x8YA9OmQdu2cPr0Pd2kjYuP44etP9BraS8Onj9IpdBKzHxlJuVylfNB4UqpQKOh77QjR6BNGzuUU6oUzJ0LJUve9W3GGOb8NYdui7qx9eRWSmYvyajXRlE9f3XdaFwp9a809J0SHw+jRtmx+9hYGDIE3n0Xktz9R7L68Gq6LuzKikMrKJChAD/W/5GXH3+ZROK13S+VUgmEhr4Tdu60G5H/8QdUq2bDP3/+u75tx6kddF/UnZm7Z5ItVTZGPjeSt0u9TdLESX1QtFIqIdDQ96Xr121vnAEDbH+c8HB44w27OfkdHL5wmD5L+/w9175/5f50KNeB1MnurceOUkpp6PvK6tX26n77dnj1Vfjii7veqD175Swfr/iYEetG6Fx7pZRHaOh726VLtu3xl19Crlx2Q/Lnn7/jW6Jjohm+djiDVg7i4rWLvF78dfpV6qdN0JRSD0xD35tmz7btjiMj7XTMjz6CNP/e/uBKzBW+2/IdfZf15eilo9R+uDYDqwykaLaiPixaKZWQaeh7w5490LEjzJkDjz8Oq1ZBudvPmzfGsPbIWsIjwpm8bTIXrl2gfK7yTK4/mafyPuXjwpVSCZ2GviddumRv0n7+ud2EfMgQaNcOkiX7n5ceu3SM77Z8R3hEODtP7yRFkhQ0KNyApiWaUjm0ss61V0p5hYa+J8THw/ff2zn3x49Ds2a2hUL27P942fW46/y6+1cmRExg3t55xJk4KuSuwJgXxvDy4y/rVoRKKa/T0H9Q69fbq/m1a+GJJ+zK2rJl//GSiOMRTNg0gUlbJ3HmyhlypslJ5wqdaVqiqfazV0r5lIb+/TpxArp1gwkTIFs2O+f+9df/7oR5LfYaEyImMGrjKCKOR5AscTJefORFmpVoRo0CNUicKLGz9SulgpKG/r26fh1GjIB+/exuVp07Q8+ekNYOzcTGxzJx80T6LuvLoQuHKJWjFCNqjaBx0cZkTJHR4eKVUsFOQ/9ezJtnWx3v3g3PPWdv2D78MADxJp6p26fSa2kv9pzZQ5mcZRj7wliq5a+mN2WVUn5DQ98dx4/befbTpkGhQv9YYGWMYfae2Xy45EM2n9hMkaxF+KXRL9R5pI6GvVLK72jo34kxMHGinXMfHW0XV3XqBMmTA7D4wGJ6LO7Bmsg1FMxYkB/q/UCjIo2026VSym9p6P+bQ4fsZibz5kHFijBuHDxiZ9qsiVxDj8U9WHxgMbnT5mbMC2N4s/ib2u1SKeX3NPRvdaPPfZcu9vHw4XaTk0SJ2Hx8Mz2X9GT2ntlkTZWVL2p+QauwVoQkCXG6aqWUcouG/s327oUWLWDZMqhaFcaMgXz52HRsEwNXDuTnHT+TPiQ9A6sMpN0T7bS1sVIq4GjoA8TF2VbHH34ISZPC2LHw1lusjlzDgB/aMuevOaRNnpaeT/WkU4VOpA9J73TFSil1XzT0t2+H5s3titoXXsB89RVLYvYwYGJVlhxcQqYUmRhQeQBtyrbRsFdKBbzgDf2YGBg0CPr3h7RpMZMmMad0Wgb83pA1kWvIkToHQ2sMpVXpVqRKlsrpapVSyiOCM/TXrYOWLWHzZuIbvcz0d2vw0ZZPiZgcQd50efn6+a9pWqKp3qBVSiU4wRX658/bXay+/prYHNn4cdy7DLzyO7vm/8TDmR5mwosTaFK0iU69VEolWMER+sbA5MnQsSNXz54k/P1KfJJjPwcOD6NYtmJMrj+ZBoUbaBM0pVSCl/BD/6+/oE0bLixfwNd1H+KLEhk5cW0JZdOUZdhzI6j9cG1tl6CUChoJN/SvXYPBgzk+7CO+KAdfdw/hojlCzVw1+eDJD3gm7zMa9kqpoJMwQ3/RIvZ1acGnOQ8S3jYRMYmhYeGGdK3YlZI5SjpdnVJKOcaroS8izwLDgMTAWGPMIG8ejxMniOjWjMEX5/JTbUiSOClNSzajc8XOFMxY0KuHVkqpQOC10BeRxMBIoDoQCawXkVnGmB2ePpaJi2P5yM4M+nME8/LFksYk4/0n2tLhyffJkSaHpw+nlFIBy5tX+mWBvcaY/QAiMhl4EfBo6F88/l9qDirCmgyXyfpQUgYW7Ujr53rp6lmllLoNb4b+Q8Dhm76OBJ649UUi0hJoCZAnT557PkjabHkokDgzb2SuT9O3vyJFspT3Wa5SSiV8jt/INcaMBkYDhIWFmXv+ABG+H3rA02UppVSC5M0tno4AuW/6OpfrOaWUUg7xZuivBwqJSD4RSQa8Aszy4vGUUkrdhdeGd4wxsSLSFvgdO2VzvDFmu7eOp5RS6u68OqZvjJkDzPHmMZRSSrnPm8M7Siml/IyGvlJKBRENfaWUCiIa+kopFUTEmHtfD+UtInIK+O99vj0zcNqD5QSSYD53CO7z13MPXjfOP68xJou7b/Kr0H8QIrLBGBPmdB1OCOZzh+A+fz334Dx3uP/z1+EdpZQKIhr6SikVRBJS6I92ugAHBfO5Q3Cfv5578Lqv808wY/pKKaXuLiFd6SullLoLDX2llAoiARf6IvKsiOwWkb0i8sFtvp9cRKa4vr9WREIdKNMr3Dj3p0XkTxGJFZEGTtToLW6c+3siskNEtojIIhHJ60Sd3uLG+b8jIltFJEJEVopIYSfq9Ia7nftNr6svIkZEEtQ0Tjd+9k1F5JTrZx8hIi3u+IHGmID5g23RvA/IDyQDNgOFb3nNf4BvXI9fAaY4XbcPzz0UKAZMBBo4XbOPz70ykNL1uHVC+bnfw/mnvelxHWCe03X76txdr0sDLAfWAGFO1+3jn31T4Et3PzPQrvT/3mzdGHMduLHZ+s1eBL51Pf4ZqCoi4sMaveWu526MOWiM2QLEO1GgF7lz7kuMMdGuL9dgd2pLKNw5/4s3fZkKSCgzNNz5fx6gPzAYuOrL4nzA3fN3W6CF/u02W3/o315jjIkFLgCZfFKdd7lz7gnVvZ57c2CuVyvyLbfOX0TaiMg+4BOgvY9q87a7nruIlAJyG2N+82VhPuLuf/v1XUObP4tI7tt8/2+BFvpK3ZGIvAaEAZ86XYuvGWNGGmMKAF2Bnk7X4wsikgj4DOjkdC0O+hUINcYUAxbw/yMdtxVooe/OZut/v0ZEkgDpgDM+qc67gnmjebfOXUSqAT2AOsaYaz6qzRfu9Wc/GajrzYJ86G7nngYoAiwVkYNAOWBWArqZe9efvTHmzE3/vY8FSt/pAwMt9N3ZbH0W8KbrcQNgsXHd7QhwwbzR/F3PXURKAqOwgX/SgRq9yZ3zL3TTl88Df/mwPm+647kbYy4YYzIbY0KNMaHY+zl1jDEbnCnX49z52ee46cs6wM47fqLTd6fv4272c8Ae7B3tHq7n+mF/0AAhwFRgL7AOyO90zT489zLYMb8o7N9utjtdsw/PfSFwAohw/ZnldM0+Pv9hwHbXuS8BHne6Zl+d+y2vXUoCmr3j5s/+Y9fPfrPrZ//onT5P2zAopVQQCbThHaWUUg9AQ18ppYKIhr5SSgURDX2llAoiGvpKKRVENPSV3xKRsb7uFikiX4jI07d5vpKIzPZlLXcjIm1F5C2n61CBRUNf+S1jTAtjzA5fHU9EMgHljDHLvXycxB76qPFAOw99lgoSGvrKcSKSSkR+E5HNIrJNRBq5nl96Yzm9iDQXkT0isk5ExojIl67nw0XkaxFZIyL7XVfk40Vkp4iE33SMr0Vkg4hsF5G+/1JKfWDeTe95VkR2icifQL1b6h3vqmWTiLzoej6liPzk6us/w7Wfw436L4vIUBHZDJQXkddc748QkVE3fhGISA0RWS12X4SpIpLa9fygm/YLGAJgbFfRgyJS1jM/CRUMNPSVP3gWOGqMKW6MKcJNwQsgIjmBD7F9VSoCj97y/gxAeaAjdon658DjQFERKeF6TQ9jTBh2v4FnRKTYbeqoCGx0HTMEGAO8gO1lkv2m1/XAtvcoi+3j/6mIpMLu5XDOGFPYVe/NPVBSAWuNMcWxq6UbARWNMSWAOKCJiGTGNkqrZowpBWwA3nP9DeQl7CrbYsCAmz53A/DUbc5FqdvS0Ff+YCtQXUQGi8hTxpgLt3y/LLDMGHPWGBODbbNxs1+NXVq+FThhjNlqjInHLk0Pdb3mZdcV+ybsL4Tb3SvIAZxyPX4UOGCM+cv12d/f9LoawAciEoFd9h8C5AGexDY7wxizDdhy03vigGmux1WxvxDWuz6jKnaTjHKuuv5wPf8mkBfbHvwqME5E6gHRN33uSSDnbc5FqdtK4nQBShlj9rh6oj8HDBCRRcaYfvfwETc6DMbf9PjG10lEJB/wPlDGGHPONewTcpvPufIvz99KgPrGmN3/ePLOe/VcNcbE3fT+b40x3W55/wvAAmPMq/9zQDuEUxXbRLAtUMX1rRBX3Uq5Ra/0leNcwzfRxpjvsX3wS93ykvXYIZkMrnbZ9e/xEGmxTeguiEg2oNa/vG4nUND1eBcQKiIFXF/fHMS/A+3ElfKuDp8AfwAvu54rDBT9l+MsAhqISFbXazOK3dN3DVBRRAq6nk8lIg+7xvXTGWPmYIewit/0WQ8D2+72L0CpG/RKX/mDothx8XggBrvH7d+MMUdEZCC2a+pZbCDfOgT0r4wxm0Vkk+t9h7HhfDu/Aa2AscaYqyLSEvhNRKKBFdje7WC35vsC2CJ2E48DQG3gK+BbEdnhOtb229VpjNkhIj2B+a73xwBtjDFrRKQp8KOIJHe9vCdwCZjpus8gwHs3fVxFoI+7/y6U0i6bKiCISGpjzGXXlf4MYLwxZoYXjrMSqG2MOX8f700MJHX9wiiAbff8iLF7m3qc628Y7xljXvfG56uESa/0VaDoI3ZnrBBgPvCLl47TCXtT9vx9vDclsEREkmKvyP/jrcB3yYydJaSU2/RKXymlgojeyFVKqSCioa+UUkFEQ18ppYKIhr5SSgURDX2llAoi/wdowFZN/gvmQwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "number_samples = 1000\n", "x_values = []\n", "y_values = []\n", "for sigma_counter in range(0, 25, 1):\n", " sigma = sigma_counter / 50\n", " rms = 0\n", " for i in range(number_samples):\n", " tmp_left2right = copy.deepcopy(left2right)\n", " tmp_left2right[1] = tmp_left2right[1] + random.normalvariate(0, sigma)\n", " tmp_stereo_matrix = rigid_body_parameters_to_matrix(tmp_left2right)\n", " left_camera_point = triangulate_points_to_3d(gold_left_image_point, \n", " left_intrinsics, \n", " left_distortion, \n", " gold_right_image_point, \n", " right_intrinsics, \n", " right_distortion, \n", " tmp_stereo_matrix)\n", " left_world_point = convert_camera_point_to_world(left_camera_point,\n", " camera_to_marker,\n", " marker_to_world\n", " )\n", " diff = np.linalg.norm(gold_world_point - left_world_point)\n", " diff = diff * diff\n", " rms = rms + diff\n", " rms = rms / number_samples\n", " rms = np.sqrt(rms)\n", " x_values.append(sigma)\n", " y_values.append(rms)\n", "\n", "rotation_values = copy.deepcopy(y_values)\n", "\n", "# Repeat of above, just doing translation.\n", "y_values = []\n", "for sigma_counter in range(0, 25, 1):\n", " sigma = sigma_counter / 50\n", " rms = 0\n", " for i in range(number_samples):\n", " tmp_left2right = copy.deepcopy(left2right)\n", " tmp_left2right[3] = tmp_left2right[3] + random.normalvariate(0, sigma)\n", " tmp_stereo_matrix = rigid_body_parameters_to_matrix(tmp_left2right)\n", " left_camera_point = triangulate_points_to_3d(gold_left_image_point, \n", " left_intrinsics, \n", " left_distortion, \n", " gold_right_image_point, \n", " right_intrinsics, \n", " right_distortion, \n", " tmp_stereo_matrix)\n", " left_world_point = convert_camera_point_to_world(left_camera_point,\n", " camera_to_marker,\n", " marker_to_world\n", " )\n", " diff = np.linalg.norm(gold_world_point - left_world_point)\n", " diff = diff * diff\n", " rms = rms + diff\n", " rms = rms / number_samples\n", " rms = np.sqrt(rms)\n", " y_values.append(rms)\n", "\n", "translation_values = copy.deepcopy(y_values)\n", "\n", "plt.plot(x_values, rotation_values, 'r', label='Ry') \n", "plt.plot(x_values, translation_values, 'g', label='Tx') \n", "plt.legend(loc='upper left')\n", "plt.xlabel('sigma (degrees)')\n", "plt.ylabel('error (mm)')\n", "plt.show() " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Experiment 2: Discussion and Conclusion\n", "\n", "These numbers are interesting. With only 4-6mm baseline, the system is very sensitive to changes in the stereo calibration parameters. The standard deviation on both rotation and translation, should be less than tenths of a millimetre or degree. This is exceptionally hard to achieve as each calibration gives numbers that vary by more than this. This suggests that we cannot be confident about the stereo calibration, and it should probably be done once by an engineer, under lab conditions, checked, and not done by a user in the field." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Experiment 3: Accuracy of hand-eye calibration\n", "\n", "Similarly for hand-eye. Harder to isolate calibration errors to individual parameters. So, we will try just adding noise to all parameters." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[88.04229429051159, -41.70296857104102, -149.50452705845592, 20.81503, 131.959603, -247.139331]\n" ] } ], "source": [ "print(marker2camera)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEGCAYAAACJnEVTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAgjklEQVR4nO3de5yWc/7H8ddHIqckZTeH3Ra7Sw6hqfjlHDY2sZvjalmn7IocclrsWpJFDjkryXFRSueEUiJKR51DsaseKCrpMJmpz++P7xUj08w1Ndd9zdzX+/l4zKOZe+77vj7XDO++fa/v9fmauyMiItmwRdoFiIhI7ij0RUQyRKEvIpIhCn0RkQxR6IuIZMiWaRdQUr169bxhw4ZplyEiUm1MmjTpK3evH/f5VSr0GzZsyMSJE9MuQ0Sk2jCz/1bk+ZreERHJEIW+iEiGKPRFRDJEoS8ikiEKfRGRDFHoi4hkiEJfRCRDFPoiImlZvRr69oW77srZIRX6IiK59N13MGQItGsHu+wCp58OjzwCRUU5OXyVuiNXRCQvFRfD6NHw0kvwyiuwdCnstBOcdVb4OOoo2DI3cazQFxFJwrp18O67IehffhkWLYLtt4dTTw1Bf/zxsNVWOS9LoS8iUpkmT4YXX4TeveGzz6BWLWjdOgT9SSfBNtukWp5CX0SkMqxeDR07Qs+eULMm/O538O9/Q5s2sMMOaVf3PYW+iMjmmjMnXJCdMQNuuAGuvRbq1k27qlIp9EVENsdzz8Hf/hambYYPDyP8KkxLNkVENsWqVXDhhXDuudCkCUydWuUDHxT6IiIVN2sWNGsGTz0FN98MI0fCbrulXVUsmt4REamIZ58N0znbbRemc044Ie2KKkQjfRGROFauhPPPh/POg6ZNw3RONQt8UOiLiJRv5swwnfPMM/CPf8CIEbDrrmlXtUk0vSMiUpann4ZLLw1r7V9/HY47Lu2KNotG+iIipSkuDqtzzj8fmjcP0znVPPBBoS8i8lNFRfCnP0GvXmF1zogR0KBB2lVVCk3viIiUVFQEZ58N/frBvffC1VenXVGlUuiLiKz33XehMVr//nD//XDllWlXVOkU+iIiEAL/jDNg4EB44IHQPC0PKfRFRNasCQ3TBg+Ghx6Cyy5Lu6LEKPRFJNvWrIHTTgtbGD78MHTokHZFiVLoi0h2FRZC27YwbBg8+mhor5DnFPoikk2FhfDHP8Krr8Ljj8Mll6RdUU4kvk7fzGqY2RQzG5L0sUREYiksDHvVvvoq9OiRmcCH3NycdQUwOwfHEREp3+rVcMopoaVCz55w8cVpV5RTiYa+me0O/B7omeRxRERiWbUq7Fn7xhsh8C+8MO2Kci7pkX434Dpg3caeYGbtzWyimU1cvHhxwuWISGatD/yRI0N7hQsuSLuiVCQW+mbWGljk7pPKep6793D3AncvqF+/flLliEiWTZsGhx8Ob74Zumb+5S9pV5SaJEf6LYA2ZvYp8BJwrJk9n+DxRER+rKgIOneGggJYuBAGDAh72mZYYqHv7n93993dvSFwFvCmu7dL6ngiIj8ybRoceij8859hLf7MmWF6J+PUWllE8kvJ0f2CBaFb5osvQr16aVdWJeTk5ix3Hw2MzsWxRCTDpk8P8/WTJ4dumQ89pLDfgEb6IlL9FRXB7bdDkyYa3ZdDbRhEpHrT6L5CNNIXkeqpqAi6dAmj+88+0+g+Jo30RaT6mT49bFg+aRKceWZoiaywj0UjfRGpPoqK4Lbbwuj+f/+Dvn3hpZcU+BWgkb6IVA9Tp4a5+w8+CBuXP/igwn4TaKQvIlXbd9+FG6yaNoUvvgiblr/wggJ/E2mkLyJV16RJYe5++nT485+hWzeoWzftqqo1jfRFpOopLIQbb4TmzeHrr8OG5c8+q8CvBBrpi0jVMn58GN3Pnh3+vO8+qFMn7aryhkb6IlI1rF4N110H//d/8O23YSvDXr0U+JVMI30RSd/YsWFTkw8/hPbtoWtXqF077arykkb6IpKeSZNCu+PDD4c1a8I2ht27K/ATpNAXkdybPDlsTl5QAO+8E5qlTZ8Oxx2XdmV5T9M7IpI7kyfDrbfCoEFhrr5zZ7j8cthxx7QrywyFvogkb8qUEPYDByrsU6bQF5HkTJ0awn7AgBD2t90GHTsq7FOk0BeRylcy7HfcMXx+xRUK+ypAoS8ilWftWrjsMnj88R/CvmNHrbWvQhT6IlI5iopCf5zeveHKK+GWWxT2VZBCX0Q2X2Fh2Mxk0CC4+2649tq0K5KNUOiLyOZZtQpOPTXcWPXII3DppWlXJGVQ6IvIplu+HFq3Dm0UnnoqbHIiVZpCX0Q2zZIl0KpVWIP/wgthekeqPIW+iFTcokVw/PEwZw706xf650i1oNAXkYpZuBBatgwbkw8ZEsJfqg2FvojE98knIfC/+gpeew2OOCLtiqSCFPoiEs/cuaEL5sqVMGIENGuWdkWyCRT6IlK+9W2P3WH0aDjwwLQrkk2kfvoiUraJE+Hoo2HLLWHMGAV+NaeRvoj81Ndfh81NxoyBnj2hbl0YORL23DPtymQzKfRFJKzIefvtEPJjxsDMmeHxrbeGI4+EJ5+EPfZIt0apFAp9kaxxh3nzfhzy8+eH722/PbRoAX/6U1iZ07Qp1KqVbr1SqRT6IlnxxRdw440wfDh8/nl4bOedQ7hfdlkY0TduHObuJW8l9ts1s1rAGGDr6Dh93f2WpI4nImXo2xf++tew3PIPfwgBf+SRsM8+sIXWc2RJkn+lrwGOdfcVZlYTeMfMXnX3cQkeU0RKWro0jOJfeCFM1Tz7bAh6yazE/or3YEX0Zc3ow5M6nohs4PXX4YADoE+fsDftu+8q8CXZdfpmVsPMpgKLgDfcfXwpz2lvZhPNbOLixYuTLEckG1auDD3tf/e7sGXh+PHwj39orl6AhEPf3de6+0HA7kAzM9u/lOf0cPcCdy+oX79+kuWI5L9334WDDgp71HbqBJMmwSGHpF2VVCE5uYLj7suAUUCrXBxPJHPWrAkrc444AoqLYdQouOceLbeUn0gs9M2svpnViT7fBjgemJPU8UQya9q00Pzs3/+GCy4IXx91VNpVSRWV5Ei/ATDKzKYBEwhz+kMSPJ5ItqxdC3fdBQUF8OWXMHgwPPEE7LBD2pVJFZbYlR13nwYcnNT7i2Ta+PFw+eUwYQKcdho89hjUq5d2VVINxAp9M9sFaAHsCqwGZgAT3X1dgrWJyIY+/xxuuCGst2/QIKy/P+ssMEu7Mqkmygx9MzsGuAGoC0whLL2sBZwK7GVmfYF73X15wnWKZNuaNdCtG9x+O3z3XQj+G2/UVI5UWHkj/ZOAi939fxt+w8y2BFoTLtD2S6A2EXEPc/VXXx2apLVpA/feC3vvnXZlUk2VGfrufm0Z3ysGBlR2QSISmT0brrwy3Fm7775hT9oTTki7Kqnm4s7p1wHOBRqWfI27d0ykKpEsW7YMbr0VHn4YttsO7r8fOnSAmjXTrkzyQNzVO8OAccB0QBdvRZKwdi306gU33QRffQUXXQRduoDuVJdKFDf0a7n71YlWIpJlc+aEjUumTIHDDw8979U+QRIQN/SfM7OLgSGElskAuPuSRKoSyZJ58+DYY0P7BC3BlITFDf3vgK7ATfzQHtkB7ZIssjkWLICWLcOSzLfegv1/0pNQpFLFDf1OwN7u/lWSxYhkypdfhsBfuhTefFOBLzkRN/Q/BlYlWYhIpixZEpZfLlgQlmI2aZJ2RZIRcUN/JTDVzEbx4zl9LdkUqajly+HEE8PF26FDw4VbkRyJG/oD0I1YIptv1So4+eSwuckrr8Bxx6VdkWRMrNB392eSLkQk761ZA23bwttvh1U6bdqkXZFkUKx++mbW2symmNkSM1tuZt+amZqsicRVXAxnnx3W3/fsGZZliqQg7vRON+CPwHR393KeKyIlrVsH558P/fvDAw+E3a1EUhJ356zPgBkKfJEKcodLL4Xnnw8tFTpq7YOkK+5I/zpgmJm9xY9X79yXSFUi+cAdrrkGuneHv/899L8XSVnc0O8CrCBsoLJVcuWI5JFbb4X77gvbGnbpknY1IkD80N/V3XW7oEhJxcWwenX4WLXqx3++8QbcdluYv+/WTb10pMqI3VrZzE5w99cTrUakqikshE6dQl+cDcO9qKjs1555JvToAVvEvXQmkry4of834BozWwMUAQa4u9dOrDKRtH3+OfzhDzB+PLRuDXXqwLbbwjbblP9n7dpQUKDAlyon7s1Z2n1ZsmXyZDjllNAj55VXQviL5IEyhyFm1rCc75uZ7V6pFYmkrW/f0A/HDMaOVeBLXinv355dzayfmZ1rZvuZ2S5m9gszO9bMOgNjgX1zUKdI8tzDxdfTT4eDDoIJE8KfInmkzOkddz/dzBoB5wAXAA0ILZZnE/bN7eLuhYlXKZK0VavCXbN9+sC554YLsFtvnXZVIpWu3Dl9d59F2DFLJD8tWBDm76dMgbvvDjdUaYml5Km4q3dE8tP778Opp8K338KgQWGVjkge03oyya4XXoAjj4RateC99xT4kgnlhn60QmePXBQjkhPr1sFNN8E550Dz5mG0r/1pJSPKDf2os+awHNQikrwVK8JGJnfcARddFNol1KuXdlUiORN3emeymTVNtBKRpC1YAEccEebuu3ULK3S2Uv9AyZa4F3KbA+eY2X8Jm6Svb8NwYGKViVSmKVPCnP3y5TBkSNiYXCSD4ob+7xKtQiRJgweH7Ql33jncYXugxiqSXbGmd9z9v0Ad4OToo0702EaZ2R5mNsrMZpnZTDO7YrOrFakI9zCNc8op0KhRaJymwJeMi7sx+hXAf4Bdoo/nzezycl5WDHRy90bAoUCH6O5ekeQVF8Nll8FVV4V1+KNHQ4MGaVclkrq40zsXAs3dfSWAmd0FvAc8tLEXuPvnwOfR59+a2WxgN2DWZlUsUp7ly0Mv++HD4dpr4c471eJYJBI39A1YW+LrtdFj8V4cunUeDIyPXZnIpvjf/8IF21mzwuqciy9OuyKRKiVu6D8FjDez/tHXpwJPxnmhmW0P9AOudPflpXy/PdAe4Be/+EXMckRKMWECnHxy2Nlq+HA47ri0KxKpcuLckbsFMA44H1gSfZzv7t1ivLYmIfD/4+6vlPYcd+/h7gXuXlC/fv2K1C7yg1degaOOCrtWvfeeAl9kI+J02VxnZo+4+8HA5LhvbGZG+NfAbHe/bzNqFNk4d7jnHrj++tBSYcAA+NnP0q5KpMqKe3VrpJm1jYI8rhbAn4FjzWxq9HFSxUsU2YiVK6F9e7juurDxyZtvKvBFyhF3Tv8S4Gqg2MwKibExuru/QwUu9orE5h46ZF5/PSxcCDfeCJ07a4WOSAxx5/RbufsW7r6Vu9d29x3KCnyRxEyYAC1aQLt28POfwzvvQJcuCnyRmOJ02VwHPJyDWkQ27vPPw3aGzZrB/PnQq1doidyiRdqViVQrSc7pi2y+NWvgrrvgN7/5YUrnww/DXwAa3YtUWEXn9Nea2WpizOmLbBb30AK5UyeYNy/0z7nnHth777QrE6nWYoW+u++QdCEi35s5E668EkaMCI3SXn8djj8+7apE8kLchmtmZu3M7B/R13uYWbNkS5PMWbIELr8cGjeGiRPhwQdh6lQFvkglijsp+ihwGPCn6OsVwCOJVCTZ9MYbYVT/6KNwySXw0UfhL4CaNdOuTCSvxA395u7eASgEcPelgPaZk81XVBQuzp5wQtjkZNIkeOQR7VsrkpC4F3KLzKwG4ABmVh9Yl1hVkg3z58PZZ4ell+3bw/33w7bbpl2VSF6LO9J/EOgP7GJmXYB3gDsSq0ryX+/ecPDBMHcu9OkD3bsr8EVyIO7qnf+Y2SSgJWG55qnuPjvRyiQ/rVwJV1wBTz4Jhx0W1t43bJh2VSKZEXd6B3efA8xJsBbJd9OmhR2t5s4N/XL+9S9dqBXJMd3SKMlzDxdnmzWDb74JK3W6dFHgi6Qg9khfZJMsWQIXXAADB8JJJ8HTT4M2yxFJjUb6kpwxY8KNVsOGwX33weDBCnyRlCn0pfItWhQu1h5zDNSqFbYvvOoqNUgTqQI0vSOVZ/lyuPfe8FFYCBdfDF27wg5q3SRSVSj0ZfMVFob2CXfcAV9/DWecEXay+s1v0q5MRDagf2/LpisuDuvtf/3r0AK5oCC0UejdW4EvUkUp9KXi3KFfP9h/f7joIth9dxg1CoYPh0MOSbs6ESmDQl8qZsSIsN7+tNOgRg0YMADefReOPjrtykQkBoW+xDNhAhx3XOhtv3hxWG8/bVrY0Uq7aIpUGwp9KVtREVx3XRjdT5sGDzwQ2iicd14Y6YtItaLVO7Jxn34KZ50F48fDX/8Kd9+t5Zci1ZxCX0rXrx9ceGG4aPvyy2EOX0SqPU3vyI8VFkKHDiHkf/vbsEetAl8kbyj05Qdz58Khh4Ybra65Bt5+G371q7SrEpFKpOkdCZ59Fi69FLbZBoYODR0xRSTvaKSfdStWhJU4550X7qidOlWBL5LHFPpZ9sEHIeiffx5uuQVGjoTddku7KhFJkKZ3ssgdHnsMrr4a6tYNYa87akUyQSP9rPnyS2jbNqzQOfbYMNpX4ItkhkI/K9zDNE6jRmEnq65dYcgQ7WQlkjGa3smCBQvCHbVDh8Jhh0GvXrDPPmlXJSIp0Eg/n7nDE0/AfvuF1sfduoW19wp8kcxKLPTNrJeZLTKzGUkdQ8rwySehI2b79tCkCUyfHvatVZM0kUxLcqT/NNAqwfeX0qxbBw8+GDY4ef996N49rM7Zc8+0KxORKiCxOX13H2NmDZN6fynF3LmhSdrYseEGq8cfhz32SLsqEalCUp/TN7P2ZjbRzCYuXrw47XKqp+JiuOsuaNwYZs0KLRWGDFHgi8hPpB767t7D3QvcvaC+lg9W3KRJoUnaDTdA69Yh9P/8Z+1mJSKlSj30ZRN98gmcc05oo/DZZ6Hnfd++8POfp12ZiFRhCv3q5uuvQ/uEffaB/v3hppvgo4/U815EYklyyeaLwHvAb81sgZldmNSxMmH1arjzTthrr7BP7bnnhrC//XaoXTvt6kSkmkhy9c7ZSb13pqxdGy7M/vOf4c7ak08O4d+oUdqViUg1pOmdqso99Mg56CC44ALYdVd46y0YNEiBLyKbTKFfFU2cCC1bwu9/H/as7dMHxo2DI49MuzIRqeYU+lXJF1/A2WdD06YwYwY89BDMnAmnn64lmCJSKdRls6ro2zd0wly5Em6+Ga69VhdoRaTSKfTTtmwZXH556HXftGm4aKsumCKSEE3vpGnkSDjgAHjxRfjXv0LPHAW+iCRII/00rF4d2iY8+GAI+XHjwp21IiIJ00g/1yZMgEMOCYHfsSNMnqzAF5GcUejnSlER3Hpr2K5wxQoYMSLcWbvNNmlXJiIZoumdXJgzJ7RNmDAB2rULSzHr1Em7KhHJII30k7R+F6uDD4b580MnzOeeU+CLSGo00k/K/Plhf9qRI8MuVj17QoMGaVclIhmnkX5lKy6Ge+75YY/aHj3CLlYKfBGpAjTSr0xTp4Y9aidPhjZt4NFHYbfd0q5KROR7GulXhvXr7gsKYOHCMHc/YIACX0SqHI30N9eoUWHu/uOPwyi/a1fYaae0qxIRKZVG+ptq6VK4+GI49tjQ+37kyHCxVoEvIlWYQr+i3ENHzH33haeeguuug2nTQviLiFRxmt6piIULoUMHGDgwrL0fNiy0VBARqSY00o/DHZ54ImxT+NprcPfdYTmmAl9EqhmN9MuzYAFcdFEI+2OOCeG/115pVyUiskk00t8Y97Chyf77w9tvw8MPhyZpCnwRqcY00i/NF1+ErQsHDoQWLeDpp2HvvdOuSkRks2mkv6E+fcLofvjw0E7hrbcU+CKSNxT66331FZx5ZvjYc0+YMgU6dYIaNdKuTESk0ij0IUzj7L8/9O8PXbrAu++GdfgiInkm23P6y5bBFVeEC7aNG4cVOo0bp12ViEhisjvSf+21MLr/z3/g5pvDunsFvojkueyF/rJloTFaq1ZQuza89x507gxbbZV2ZSIiictW6A8cGO6qfeYZuP760Pe+adO0qxIRyZlszOkvWgQdO0Lv3nDggTB4MDRpknZVIiI5l98jffcwZ9+oUViZ07kzTJigwBeRzMrfkf6CBeGu2qFDoXlz6NUrhL+ISIbl30jfPWxGvt9+8OabcN99MHasAl9EhIRD38xamdlcM/vYzG5I8lgAzJsHLVvCJZeE/WpnzICrrtJdtSIikcRC38xqAI8AJwKNgLPNLJnh9tq1YUR/wAEwaVIY6Y8YEdopiIjI95Kc028GfOzu8wHM7CXgFGBWpR5l6VI48UQYPx5at4bHHoPdd6/UQ4iI5Iskp3d2Az4r8fWC6LEfMbP2ZjbRzCYuXry44kepUyf0uH/hBRg0SIEvIlKG1FfvuHsPoAdAQUGBV/gNzMKyTBERKVeSI/2FwB4lvt49ekxERFKSZOhPAH5tZr8ys62As4BBCR5PRETKkdj0jrsXm9llwGtADaCXu89M6ngiIlK+ROf03X0YMCzJY4iISHz5d0euiIhslEJfRCRDFPoiIhmi0BcRyRBzr/j9UEkxs8XAfzfx5fWAryqxnOoky+cO2T5/nXt2rT//X7p7/bgvqlKhvznMbKK7F6RdRxqyfO6Q7fPXuWfz3GHTz1/TOyIiGaLQFxHJkHwK/R5pF5CiLJ87ZPv8de7ZtUnnnzdz+iIiUr58GumLiEg5FPoiIhlS7UK/vM3WzWxrM+sdfX+8mTVMocxExDj3I81sspkVm9lpadSYlBjnfrWZzTKzaWY20sx+mUadSYlx/n81s+lmNtXM3klsP+oUlHfuJZ7X1szczPJqGWeM3/1fzGxx9LufamYXlfmG7l5tPggtmucBewJbAR8AjTZ4zqXA49HnZwG90647h+feEDgQeBY4Le2ac3zuxwDbRp//LV9+7xU4/9olPm8DDE+77lyde/S8HYAxwDigIO26c/y7/wvwcNz3rG4j/e83W3f374D1m62XdArwTPR5X6ClmVkOa0xKuefu7p+6+zRgXRoFJijOuY9y91XRl+MIO7Xlizjnv7zEl9sB+bJCI87/8wCdgbuAwlwWlwNxzz+26hb6cTZb//457l4MfAPsnJPqkhVro/k8VdFzvxB4NdGKcivW+ZtZBzObB9wNdMxRbUkr99zN7BBgD3cfmsvCciTuf/tto6nNvma2Rynf/151C32RMplZO6AA6Jp2Lbnm7o+4+17A9cDNadeTC2a2BXAf0CntWlI0GGjo7gcCb/DDTEepqlvox9ls/fvnmNmWwI7A1zmpLllZ3mg+1rmb2XHATUAbd1+To9pyoaK/+5eAU5MsKIfKO/cdgP2B0Wb2KXAoMCiPLuaW+7t3969L/PfeE2hS1htWt9CPs9n6IOC86PPTgDc9utpRzWV5o/lyz93MDga6EwJ/UQo1JinO+f+6xJe/Bz7KYX1JKvPc3f0bd6/n7g3dvSHhek4bd5+YTrmVLs7vvkGJL9sAs8t8x7SvTm/C1eyTgA8JV7Rvih67jfCLBqgFvAx8DLwP7Jl2zTk896aEOb+VhH/dzEy75hye+wjgS2Bq9DEo7ZpzfP4PADOjcx8F7Jd2zbk69w2eO5o8Wr0T83f/7+h3/0H0u9+nrPdTGwYRkQypbtM7IiKyGRT6IiIZotAXEckQhb6ISIYo9EVEMkShLzllZj1z3QHSzLqZ2ZGlPH60mQ3JZS1xmFkDM3s9gfcdYWY7Vfb7SvWi0JeccveL3H1Wro5nZjsDh7r7mISPU6MS364V8Folvt96zxG60EqGKfQlEWa2nZkNNbMPzGyGmZ0ZPT56/S3yZnahmX1oZu+b2RNm9nD0+NNm9piZjTOz+dGIvJeZzTazp0sc4zEzm2hmM83s1o2U0hYYXuI1rcxsjplNBv64Qb29olqmmNkp0ePbmlmfqFd//2iPhvX1rzCze83sA+AwM2sXvX6qmXVf/xeBmZ1gZu9Z2OvgZTPbPnr8TvthD4B7StTcCng1Ou+3zGxg9HO408zOiY4x3cz2qsjPi3An59kV/21KXkn7bjN95OcHIWyfKPH1jtGfowkN0XYFPgXqAjWBt4l6ggNPE/rHGKGN7HLgAMIgZRJwUPS8utGfNaL3PbCUOp4BTo4+r0XoWPjr6L37AEOi790BtIs+r0O4A3I74Bqge/T4/kAx0R2fhPbFZ0Sf70tofFUz+vpR4FygHqHP+3bR49cD/yR0fp3LD/tU1ylxLlOjz48GlgENgK0JPVdujb53BdCtIj+v6LkfATun/d+HPtL70EhfkjIdON7M7jKzI9z9mw2+3wx4y92XuHsRoXVGSYPd3aP3+dLdp7v7OsLt5g2j55wRjdinAPsBpV0raAAsjj7fB/jE3T+K3vv5Es87AbjBzKYS/gKpBfwCOJwQqLj7DGBaidesBfpFn7ckNLqaEL1HS8LGF4dGdY2NHj8P+CWh5Xch8KSZ/RFYvxdAc2B8iWNMcPfPPTTUmgesn+ufXuLnEPfnBbCI8BeuZNSWaRcg+cndP4z6nJ8E3G5mI939tgq8xfqugetKfL7+6y3N7FeEUXhTd18aTWPUKuV9Vm/k8Q0Z0Nbd5/7owbL33yl097UlXv+Mu/99g9efDLzh7j+ZVjGzZoS/HE4DLgOOBU6kxHQUPz33kj+XLUt5Xqk/rxJf1yL8TCSjNNKXRJjZrsAqd3+e0Nv+kA2eMgE4ysx2stACu20FD1Gb0FjuGzP7GSEsSzMb2Dv6fA7QcP1cOD+e334NuNyilI+6dgKMBc6IHmtEmDYpzUjgNDPbJXpuXQv79I4DWpjZ3tHj25nZb6J5/R3dfRhwFdA4ep+WhOZxlS46t58TptUkozTSl6QcAHQ1s3VAEWHf2u+5+0Izu4PQCXUJIZA3nALaKHf/wMymRK/7jBDOpRkKXAL0dPdCM2sPDDWzVYTrCDtEz+sMdAOmWdiY4xOgNWFu/hkzmxUda2Zpdbr7LDO7GXg9en0R0MHdx5nZX4AXzWzr6Ok3A98CA82sFuFfCVebWX3Cvx6+jftzqKAmwDgPO8pJRqnLpqTGzLZ39xXRSL8/0Mvd+ydwnHeA1u6+bBNeW4NwcbYw+hfCCOC3HvYrrVQWdv3a3d3vrOz3jt7/AULL6ZFJvL9UDxrpS5r+ZWG3q1qEC5QDEjpOJ8JF2WWb8NptgVFmVpMwIr80icAHiKbCkjRDgS8a6YuIZIgu5IqIZIhCX0QkQxT6IiIZotAXEckQhb6ISIb8P41o/vGuto7FAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "number_samples = 1000\n", "x_values = []\n", "y_values = []\n", "for sigma_counter in range(0, 25, 1):\n", " sigma = sigma_counter / 50\n", " rms = 0\n", " for i in range(number_samples):\n", " tmp_marker2camera = copy.deepcopy(marker2camera)\n", " for j in range(6):\n", " tmp_marker2camera[j] = tmp_marker2camera[j] + random.normalvariate(0, sigma)\n", " tmp_marker2camera_matrix = rigid_body_parameters_to_matrix(tmp_marker2camera)\n", " tmp_camera2marker = np.linalg.inv(tmp_marker2camera_matrix)\n", " left_world_point = convert_camera_point_to_world(gold_left_camera_point,\n", " tmp_camera2marker,\n", " marker_to_world\n", " )\n", " diff = np.linalg.norm(gold_world_point - left_world_point)\n", " diff = diff * diff\n", " rms = rms + diff\n", " rms = rms / number_samples\n", " rms = np.sqrt(rms)\n", " x_values.append(sigma)\n", " y_values.append(rms)\n", "\n", "plt.plot(x_values, y_values, 'r') \n", "plt.xlabel('sigma (degrees/mm)')\n", "plt.ylabel('error (mm)')\n", "plt.show() " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Experiment 3: Brief Discussion and Conclusion\n", "\n", "So, it can be seen that any noise on the calibration parameters of hand-eye calibration can have a huge affect on the final accuracy of the ability to measure the location of a world/tracker coordinate. The x-axis shows sigma, which is is used, assuming a Gaussian noise model, to add noise to the 3 rotation and 3 translation parameters of the hand-eye matrix. You can see that these parameters must be correct to within a sigma of 0.2 millimetres or degrees in order to keep the accuracy of measurement within 2mm. So, is this realistic, in the field? You should do a user study and look at the variation in calibration parameters when many users do many calibrations. If the variability in the calibration parameters that you see in the field, is larger than this, could we ever assume that system accuracy is good enough in practice? If not, then how would you design your tracking marker and laparoscope? Is calibration something that the user should be doing? If nott, what should the vendor do?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "mphy0026", "language": "python", "name": "mphy0026" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.8" } }, "nbformat": 4, "nbformat_minor": 2 }