{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "done\n"
     ]
    }
   ],
   "source": [
    "# from opencv documentation (tutorial_py_calibration)\n",
    "# actual calibration images not committed to repo, since they're ~10 MB each\n",
    "\n",
    "import numpy as np\n",
    "import cv2\n",
    "import glob\n",
    "\n",
    "length = 11\n",
    "width = 8\n",
    "\n",
    "# termination criteria\n",
    "criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 30, 0.001)\n",
    "# prepare object points, like (0,0,0), (1,0,0), (2,0,0) ....,(6,5,0)\n",
    "objp = np.zeros((length*width,3), np.float32)\n",
    "objp[:,:2] = np.mgrid[0:length,0:width].T.reshape(-1,2)\n",
    "# Arrays to store object points and image points from all the images.\n",
    "objpoints = [] # 3d point in real world space\n",
    "imgpoints = [] # 2d points in image plane.\n",
    "images = glob.glob('*.jpg')\n",
    "for fname in images:\n",
    "    img = cv2.imread(fname)\n",
    "    gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)\n",
    "    # Find the chess board corners\n",
    "    ret, corners = cv2.findChessboardCornersSB(gray, (length,width), None)\n",
    "    # If found, add object points, image points (after refining them)\n",
    "    if ret == True:\n",
    "        objpoints.append(objp)\n",
    "        corners2 = cv2.cornerSubPix(gray,corners, (11,11), (-1,-1), criteria)\n",
    "        imgpoints.append(corners)\n",
    "        # Draw and display the corners\n",
    "        #cv2.drawChessboardCorners(img, (length,width), corners2, ret)\n",
    "        #cv2.imshow('img', img)\n",
    "        #cv2.waitKey(500)\n",
    "cv2.destroyAllWindows()\n",
    "\n",
    "ret, mtx, dist, rvecs, tvecs = cv2.calibrateCamera(objpoints, imgpoints, gray.shape[::-1], None, None)\n",
    "\n",
    "np.save('ret.npy', ret)\n",
    "np.save('mtx.npy', mtx)\n",
    "np.save('dist.npy', dist)\n",
    "np.save('rvecs.npy', rvecs)\n",
    "np.save('tvecs.npy', tvecs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# test calibration\n",
    "\n",
    "import numpy as np\n",
    "import cv2\n",
    "\n",
    "ret = np.load('cal_results/ret.npy')\n",
    "mtx = np.load('cal_results/mtx.npy')\n",
    "dist = np.load('cal_results/dist.npy')\n",
    "rvecs = np.load('cal_results/rvecs.npy')\n",
    "tvecs = np.load('cal_results/tvecs.npy')\n",
    "\n",
    "img = cv2.imread('aruco/setup1_img1.jpg')\n",
    "h,  w = img.shape[:2]\n",
    "newcameramtx, roi = cv2.getOptimalNewCameraMatrix(mtx, dist, (w,h), 1, (w,h))\n",
    "\n",
    "# undistort\n",
    "dst = cv2.undistort(img, mtx, dist, None, newcameramtx)\n",
    "# crop the image\n",
    "x, y, w, h = roi\n",
    "dst = dst[y:y+h, x:x+w]\n",
    "cv2.imwrite('aruco_test.jpg', dst)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# undistort aruco test images\n",
    "\n",
    "import numpy as np\n",
    "import cv2\n",
    "import glob\n",
    "\n",
    "ret = np.load('cal_results/ret.npy')\n",
    "mtx = np.load('cal_results/mtx.npy')\n",
    "dist = np.load('cal_results/dist.npy')\n",
    "rvecs = np.load('cal_results/rvecs.npy')\n",
    "tvecs = np.load('cal_results/tvecs.npy')\n",
    "\n",
    "images = glob.glob('aruco/raw/*.jpg')\n",
    "for fname in images:\n",
    "    img = cv2.imread(fname)\n",
    "    h, w = img.shape[:2]\n",
    "    newcameramtx, roi = cv2.getOptimalNewCameraMatrix(mtx, dist, (w,h), 1, (w,h))\n",
    "    # undistort\n",
    "    dst = cv2.undistort(img, mtx, dist, None, newcameramtx)\n",
    "    # crop the image\n",
    "    x, y, w, h = roi\n",
    "    dst = dst[y:y+h, x:x+w]\n",
    "    fname_short = fname.rsplit(\"/\")[2]\n",
    "    cv2.imwrite('aruco/undistorted/' + fname_short, dst)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'topLeft' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-3-665cb3fe1e5e>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     30\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     31\u001b[0m     \u001b[0;31m# draw the bounding box of the ArUCo detection\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 32\u001b[0;31m     \u001b[0mcv2\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mline\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mimg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtopLeft\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtopRight\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m255\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m6\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     33\u001b[0m     \u001b[0mcv2\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mline\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mimg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtopRight\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbottomRight\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m255\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m6\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     34\u001b[0m     \u001b[0mcv2\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mline\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mimg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbottomRight\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbottomLeft\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m255\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m6\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'topLeft' is not defined"
     ]
    }
   ],
   "source": [
    "# aruco marker detection test w/thresholding and marking\n",
    "# thx https://www.pyimagesearch.com/2020/12/21/detecting-aruco-markers-with-opencv-and-python/\n",
    "\n",
    "import numpy as np\n",
    "import cv2\n",
    "import glob\n",
    "\n",
    "aruco_dict = cv2.aruco.Dictionary_get(cv2.aruco.DICT_6X6_100)\n",
    "# marker created via cv2.aruco.drawMarker(aruco_dict, 1, 8)\n",
    "images = glob.glob('aruco/undistorted/test/*.jpg')\n",
    "for fname in images:\n",
    "    img = cv2.imread(fname)\n",
    "    gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)\n",
    "    ret,gray_thresh = cv2.threshold(gray,127,255,cv2.THRESH_BINARY)\n",
    "    fname_short = fname.rsplit(\"/\")[2]\n",
    "    aruco_params = cv2.aruco.DetectorParameters_create()\n",
    "    (corners, ids, rejected) = cv2.aruco.detectMarkers(gray_thresh, aruco_dict, parameters=aruco_params)\n",
    "\n",
    "    # loop over the detected ArUCo corners\n",
    "    for (markerCorner, markerID) in zip(corners, ids):\n",
    "    # extract the marker corners (which are always returned in\n",
    "    # top-left, top-right, bottom-right, and bottom-left order)\n",
    "        corners = markerCorner.reshape((4, 2))\n",
    "        (topLeft, topRight, bottomRight, bottomLeft) = corners\n",
    "        # convert each of the (x, y)-coordinate pairs to integers\n",
    "        topRight = (int(topRight[0]), int(topRight[1]))\n",
    "        bottomRight = (int(bottomRight[0]), int(bottomRight[1]))\n",
    "        bottomLeft = (int(bottomLeft[0]), int(bottomLeft[1]))\n",
    "        topLeft = (int(topLeft[0]), int(topLeft[1]))\n",
    "        \n",
    "    # draw the bounding box of the ArUCo detection\n",
    "    cv2.line(img, topLeft, topRight, (0, 255, 0), 6)\n",
    "    cv2.line(img, topRight, bottomRight, (0, 255, 0), 6)\n",
    "    cv2.line(img, bottomRight, bottomLeft, (0, 255, 0), 6)\n",
    "    cv2.line(img, bottomLeft, topLeft, (0, 255, 0), 6)\n",
    "    # compute and draw the center (x, y)-coordinates of the ArUco\n",
    "    # marker\n",
    "    cX = int((topLeft[0] + bottomRight[0]) / 2.0)\n",
    "    cY = int((topLeft[1] + bottomRight[1]) / 2.0)\n",
    "    cv2.circle(img, (cX, cY), 12, (0, 0, 255), -1)\n",
    "    # draw the ArUco marker ID on the image\n",
    "    cv2.putText(img, str(markerID),(topLeft[0], topLeft[1] - 15), cv2.FONT_HERSHEY_SIMPLEX,0.5, (0, 255, 0), 2)\n",
    "    print(\"[INFO] ArUco marker ID: {}\".format(markerID))\n",
    "    # show the output image\n",
    "    cv2.imshow(\"Image\", cv2.resize(img,None,fx=0.5, fy=0.5, interpolation = cv2.INTER_AREA))\n",
    "    cv2.waitKey(0)\n",
    "    cv2.destroyAllWindows()\n",
    "    cv2.imwrite('test_inaccurate.jpg', img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('setup4_img2.jpg', 3314.5, 1382.5)\n",
      "('setup5_img5.jpg', 1966.0, 1592.5)\n",
      "('setup5_img7.jpg', 1899.5, 1492.5)\n",
      "('setup2_img7.jpg', 3158.0, 2333.5)\n",
      "('setup3_img9.jpg', 842.0, 1569.5)\n",
      "('setup5_img2.jpg', 2153.5, 1879.5)\n",
      "('setup1_img6.jpg', 1118.5, 964.5)\n",
      "('setup3_img8.jpg', 839.5, 1637.5)\n",
      "('setup2_img5.jpg', 2966.0, 2331.5)\n",
      "('setup1_img4.jpg', 923.0, 962.5)\n",
      "('setup5_img8.jpg', 1854.5, 1421.5)\n",
      "('setup4_img1.jpg', 3313.5, 1486.5)\n",
      "('setup5_img9.jpg', 1790.5, 1326.5)\n",
      "('setup4_img3.jpg', 3316.5, 1267.0)\n",
      "('setup2_img6.jpg', 3060.5, 2332.5)\n",
      "('setup1_img9.jpg', 1426.0, 968.5)\n",
      "('setup5_img1.jpg', 2197.0, 1945.5)\n",
      "('setup2_img3.jpg', 2802.5, 2330.5)\n",
      "('setup4_img5.jpg', 3319.0, 1059.5)\n",
      "('setup5_img6.jpg', 1930.0, 1538.0)\n",
      "('setup3_img4.jpg', 833.5, 2008.0)\n",
      "('setup3_img5.jpg', 835.5, 1930.5)\n",
      "('setup1_img7.jpg', 1228.0, 966.0)\n",
      "('setup3_img1.jpg', 829.0, 2291.5)\n",
      "('setup3_img3.jpg', 831.5, 2112.5)\n",
      "('setup3_img7.jpg', 839.0, 1717.0)\n",
      "('setup1_img8.jpg', 1315.5, 966.5)\n",
      "('setup5_img3.jpg', 2100.0, 1796.5)\n",
      "('setup1_img3.jpg', 836.0, 961.5)\n",
      "('setup2_img9.jpg', 3345.5, 2335.5)\n",
      "('setup4_img6.jpg', 3321.0, 960.0)\n",
      "('setup3_img6.jpg', 836.0, 1836.5)\n",
      "('setup2_img4.jpg', 2873.5, 2331.0)\n",
      "('setup4_img4.jpg', 3317.5, 1166.5)\n",
      "('setup4_img8.jpg', 3320.5, 819.0)\n",
      "('setup2_img1.jpg', 2620.0, 2328.5)\n",
      "('setup4_img7.jpg', 3320.5, 906.5)\n",
      "('setup1_img1.jpg', 688.5, 961.0)\n",
      "('setup4_img9.jpg', 3323.5, 716.0)\n",
      "('setup2_img2.jpg', 2708.5, 2329.5)\n",
      "('setup3_img2.jpg', 831.0, 2180.0)\n",
      "('setup2_img8.jpg', 3256.0, 2334.0)\n",
      "('setup1_img5.jpg', 1015.5, 962.5)\n",
      "('setup5_img4.jpg', 2050.5, 1722.0)\n",
      "('setup1_img2.jpg', 770.5, 962.0)\n"
     ]
    }
   ],
   "source": [
    "# aruco marker center detection\n",
    "# thx https://www.pyimagesearch.com/2020/12/21/detecting-aruco-markers-with-opencv-and-python/\n",
    "\n",
    "import numpy as np\n",
    "import cv2\n",
    "import glob\n",
    "\n",
    "aruco_dict = cv2.aruco.Dictionary_get(cv2.aruco.DICT_6X6_100)\n",
    "# marker created via cv2.aruco.drawMarker(aruco_dict, 1, 8)\n",
    "images = glob.glob('aruco/undistorted/*.jpg')\n",
    "for fname in images:\n",
    "    img = cv2.imread(fname)\n",
    "    gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)\n",
    "    ret,gray_thresh = cv2.threshold(gray,50,255,cv2.THRESH_BINARY)\n",
    "    fname_short = fname.rsplit(\"/\")[2]\n",
    "    aruco_params = cv2.aruco.DetectorParameters_create()\n",
    "    corners = cv2.aruco.detectMarkers(gray_thresh, aruco_dict, parameters=aruco_params)\n",
    "    cX = (corners[0][0][0][0][0] + corners[0][0][0][2][0]) / 2.0\n",
    "    cY = (corners[0][0][0][0][1] + corners[0][0][0][2][1]) / 2.0\n",
    "    cv2.circle(img, (int(cX), int(cY)), 12, (0, 0, 255), -1)\n",
    "    print((fname_short,cX,cY))\n",
    "    #cv2.imshow(\"Image\", cv2.resize(img,None,fx=0.5, fy=0.5, interpolation = cv2.INTER_AREA))\n",
    "    #cv2.waitKey(0)\n",
    "    #cv2.destroyAllWindows()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0   -0.003616\n",
      "1   -0.003807\n",
      "2    0.003670\n",
      "3    0.002580\n",
      "4   -0.001951\n",
      "5    0.013098\n",
      "6   -0.008544\n",
      "7    0.007945\n",
      "8   -0.009375\n",
      "Name: error (mm), dtype: float64\n",
      "0.007475426524096852\n",
      "[0.01249389 0.46241624]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAFzCAYAAADmJtp4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAzhklEQVR4nO3deZzNZf/H8dc1B82MtWwlbkvLRGTsJDWhVUh3i10hrZaIFqWfdNddSWhRlrShRWSpLGEia4NRlqS63WHUSDF2s1y/P66TW1KNOWfme5b38/GYh3O+58w5n8ss77mu7/e6LmOtRURERApWjNcFiIiIRCMFsIiIiAcUwCIiIh5QAIuIiHhAASwiIuIBBbCIiIgHChXkm5UpU8ZWqVKlIN9SRETEM6tXr/7ZWlv2ZI8VaABXqVKFlJSUgnxLERERzxhj/vtnj/3tELQx5jVjTLoxZv1xx84wxsw3xmzx/3t6sIoVERGJBrk5B/w6cPUJxx4EFlhrzwMW+O+LiIhILv1tAFtrFwO/nHC4LfCG//YbwPXBLUtERCSy5fUccHlr7U7/7R+B8n/2RGNML6AXwD/+8Y8/PJ6Zmcn27ds5fPhwHkuRYIiNjaVixYoULlzY61JERKJCwBdhWWutMeZPd3Sw1o4FxgLUr1//D8/bvn07xYsXp0qVKhhjAi1H8sBay+7du9m+fTtVq1b1uhwRkaiQ13nAPxljzgLw/5ue1wIOHz5M6dKlFb4eMsZQunRpjUKIiBSgvAbwTKCb/3Y3YEYgRSh8vaevgYhIwcrNNKQpwHIgwRiz3RjTA/g3cIUxZgvQ0n9fREREcik3V0F3sNaeZa0tbK2taK2dYK3dba1tYa09z1rb0lp74lXS+SY7x7Jg00+MXrCFBZt+IjvnT08/h7Ubb7yR77//Puive//997Nw4cKgv66IiJyaAl0JK1DZOZYuE1aSum0Ph45mE1fER2KlUrzVoxG+mPAcQrXWYq0lJuZ/fwtt2LCB7OxsqlWrFvT36927N7fffjvNmzcP+muLiEjuhdVmDMmb00ndtoeDR7OxwMGj2aRu20Py5jxfA8aBAwdo1aoVtWvXpmbNmrz77rsArF69mssuu4x69epx1VVXsXOnm3WVlJTEAw88QMOGDTn//PNZsmQJ4EKzYcOGJCYmctFFF7FlyxYARowYQc2aNalZsyYjR44EYOvWrSQkJNC1a1dq1qzJtm3bflfTpEmTaNu27bH7xYoVY+DAgVx44YW0bNmSVatWkZSURLVq1Zg5cyYAr7/+Otdffz1XXHEFVapU4cUXX2TEiBHUqVOHxo0b88svbpCicuXK7N69mx9//DHP/2ciIpHGi9HVsArgDWkZHDqa/btjh45mszEtI8+vOWfOHCpUqMC6detYv349V199NZmZmfTu3ZupU6eyevVqunfvzuDBg499TlZWFqtWrWLkyJEMHToUgFdeeYW+ffuSmppKSkoKFStWZPXq1UycOJGVK1eyYsUKxo0bx9q1awHYsmULd999Nxs2bKBy5cq/q2np0qXUq1fv2P0DBw7QvHlzNmzYQPHixXnkkUeYP38+06dPZ8iQIceet379eqZNm8YXX3zB4MGDiY+PZ+3atTRp0oQ333zz2PPq1q3L0qVL8/x/JiISSX4bXX3x+Q/4ZNJcek9ZS5cJK/M9hMMqgC+sUIK4Ir7fHYsr4qNGhRJ5fs1atWoxf/58HnjgAZYsWULJkiXZvHkz69ev54orriAxMZEnnniC7du3H/ucG264AYB69eqxdetWAJo0acKTTz7J008/zX//+1/i4uL4/PPPadeuHUWLFqVYsWLccMMNx3rMlStXpnHjxietaefOnZQt+7/NM4oUKcLVV199rN7LLruMwoULU6tWrWPvD3D55ZdTvHhxypYtS8mSJWnduvWxzzn+eeXKlSMtLS3P/2ciIpEk+eufqDHtTd55rR+PLhgXlNHV3AirAE5KKEdipVLEF/FhgHj/OeCkhHJ5fs3zzz+fNWvWUKtWLR555BEef/xxrLVceOGFpKamkpqayldffcW8efOOfc5pp50GgM/nIysrC4COHTsyc+ZM4uLiuPbaa//2QqeiRYv+6WNxcXG/m5NbuHDhY9OEYmJijr1/TEzMsfc/vq6/e97hw4eJi4v76/8YEZFosHs3VW7vzCNzxrC0cm3ubfsAEPjoam6EVQD7Ygxv9WjECx3q0P+K83mhQ52AL8BKS0sjPj6ezp07M3DgQNasWUNCQgK7du1i+fLlgFsuc8OGDX/5Ot9//z3VqlWjT58+tG3bli+//JJmzZrx4YcfcvDgQQ4cOMD06dNp1qzZ39ZUvXp1vv322zy36e9888031KxZM99eX0QkLCxeDImJVF21mKeu7EX3Gx/jl/iSQOCjq7kRVldBgwvhFtXL06L6ny4/fUq++uorBg4cSExMDIULF2bMmDEUKVKEqVOn0qdPH/bu3UtWVhb9+vXjwgsv/NPXee+993jrrbcoXLgwZ555Jg8//DBnnHEGt956Kw0bNgSgZ8+e1KlT53fDwSfTqlUrkpOTadmyZVDaeLzMzEy+/fZb6tevH/TXFhEJC1lZ8MQTMGwYnHMOdvlyvlqbSfwJM2wCGV3NDWNtwc2jrV+/vk1JSfndsU2bNlG9evUCqyEcHDp0iMsvv5ylS5fi8/n+/hNOwfTp01mzZg3Dhg37w2P6WohIxNu2DTp1giVLoGtXePFFKF6c7BxL8uZ0NqZlUKNCCZISygVleqsxZrW19qQ9nrDrAUeDuLg4hg4dyo4dO066g1QgsrKyGDBgQFBfU0QkLMyYAd27w9Gj8NZb0LnzsYeCPbqaGwrgEHXVVVfly+vedNNN+fK6IiIh6/BhuP9+eOklqFcPpkyB887zuqrwughLRETklGzaBI0aufDt3x+WLQuJ8AX1gEVEJBJZC6+9Bn36QNGi8NFHcO21Xlf1O+oBi4hIZNm7Fzp0gJ49oUkTWLcu5MIXFMAiIhJJVq6EOnVg6lR48kmYOxfOOsvrqk5KAXwCay05OTlelyEiIqciJweefhouucTdXrIEHnoIgjyVM5gUwPxxd6Jhw4bRoEEDLrroIh577LFjzxs2bBgJCQlccskldOjQgeHDh3tYtYiIAPDjj3D11fDgg9CuHaSmuqHnEBdaF2H16+f+44IpMRH82wD+lS1btvDGG2+QkZHB1KlTWbVqFdZa2rRpw+LFi4mLi+ODDz5g3bp1ZGZmUrdu3d/tWCQiIh6YO9ctqLFvH4wd6877mvDYHz60AthDv+1OdP/99zNv3jzq1KkDwP79+9myZQv79u2jbdu2xMbGEhsbe2ynIRER8cDRozB4MAwfDjVrwqJFUKOG11WdktAK4Fz0VPPLb7sTWWt56KGHuOOOO373+EgPaxMRkeN89527yvmLL+Cuu+C55yAMd3jTOeATXHXVVbz22mvs378fgB07dpCenk7Tpk2ZNWsWhw8fZv/+/cyePdvjSkVEotDkye4q5y1b4IMP4OWXwzJ8IdR6wCHgyiuvZNOmTTTxn8AvVqwYb7/9Ng0aNKBNmzZcdNFFlC9fnlq1alGyZEmPqxURiRL790Pv3vD669C0qQviIK+VX9C0G9Ip2L9/P8WKFePgwYNceumljB07lrp163pdVtCE09dCRKJIairccovr9T7yCAwZAoXCo/+o3ZCCpFevXmzcuJHDhw/TrVu3iApfEZGQYy288AIMHAhlysDChZCU5HVVQaMAPgWTJ0/2ugQRkejw889u68BZs6B1a7euc5kyXlcVVLoIS0REQktyMtSu7eb4jhrl9vGNsPCFEAnggjwPLSenr4GIeC4ry53fbd4cihWDFSvcbkZhsrDGqfI8gGNjY9m9e7cCwEPWWnbv3k1sbKzXpYhItPrhB3d+d9gw6NYNVq92040imOfngCtWrMj27dvZtWuX16VEtdjYWCpWrOh1GSISjaZNgx49IDsbJk2Cjh29rqhAeB7AhQsXpmrVql6XISIiBe3QIRgwAMaMgfr14Z134JxzvK6qwHg+BC0iIlFowwZo2NCF78CBsHRpVIUvhEAPWEREooi1MG6c2/2ueHGYMweuusrrqjyhHrCIiBSMPXvcilZ33AGXXALr1kVt+IICWERECsLy5W5/9unT4emnXc/3zDO9rspTCmAREck/2dnw5JPQrBnExMDnn8OgQe52lNM5YBERyR9padC1KyxY4IaeX30VtIvcMQpgEREJvo8/dgtqHDwIEybAbbdF7IpWeaUxABERCZ4jR6B/f2jVCipUgJQUt6mCwvcP1AMWEZHg2LIF2reHNWvg3nvh2WdBS9z+KQWwiIgE7q234O67oUgR+PBDaNvW64pCnoagRUQk7/btcxdade0KdetCaqrCN5cUwCIikjerV0O9em4Dhf/7P1i4ECpV8rqqsKEAFhGRU2MtPP88NGniNlRYtAgeewx8Pq8rCys6BywiIrm3axfcequbZtS2rZtiVLq011WFJfWARUQkdxYuhNq13cIaL77olpVU+OaZAlhERP5aZiYMHgwtW7qVrFauhHvu0dzeAGkIWkRE/tzWrdCxo9tMoUcPGDUKihb1uqqIoAAWEREAsnMsyZvT2ZCWwYUVSpD01WJ8vW53F129845bz1mCRgEsIiJk51i6TFhJ6rY92AMHGJo8Ad+aT7CNGmEmT4Zq1bwuMeIogEVEhOTN6aRu20PFHd/xwsxnSPj5B8ZdfBPnjBlB82oVvS4vIimARUSEDTv20m7lLB5dOJ59p8XT+eZhLK1ah/67DtHc6+IilAJYRCTa/for7Z/tT7l5H/FZ1boMaHUfPxc9nfgiPmpUKOF1dRFLASwiEs2WLoWOHSmblsakG+/lyQuu4WCmJb6Ij8RKpUhKKOd1hRFLASwiEo2ys+Gpp9wazpUrY5Yto329+py5OZ2NaRnUqFCCpIRy+GI01ze/KIBFRKLNjh3QuTMkJ7s5vmPGQIkS+IAW1cvTonp5ryuMCgpgEZFoMnu2W8v50CGYOBG6ddOKVh4JaClKY8x9xpgNxpj1xpgpxpjYYBUmIiJBdOQI9OsHrVu7LQPXrHFBrPD1TJ4D2BhzNtAHqG+trQn4gPbBKkxERILkm2/c1oGjRkGfPrBiBSQkeF1V1At0CLoQEGeMyQTigbTASxIRkaCwFt58022cEBsLM2e6HrCEhDz3gK21O4DhwA/ATmCvtXbeic8zxvQyxqQYY1J27dqV90pFRCT3MjKgSxc3zFy/Pqxbp/ANMYEMQZ8OtAWqAhWAosaYzic+z1o71lpb31pbv2zZsnmvVEREciclBerWhSlT4PHH3f69Z5/tdVVygkAuwmoJ/Mdau8tamwlMAy4OTlkiInLKcnLguefg4ovh6FH47DN49FHw+byuTE4ikHPAPwCNjTHxwCGgBZASlKpEROTUpKe7KUVz5kC7djB+PJxxhtdVyV8I5BzwSmAqsAb4yv9aY4NUl4iI5Nann0Lt2rBokVtU44MPFL5hIKCroK21jwGPBakWERE5FZmZMGQIPP00VK8O8+ZBrVpeVyW5pJWwRETC0X/+Ax06wMqV0KsXPP88xMd7XZWcAgWwiEi4efddF7rGwHvvwU03eV2R5EFAS1GKiEgBOnAAevaE9u3hwgshNVXhG8YUwCIi4eDLL92CGq+9Bg8/7KYYVanidVUSAAWwiEgosxZeegkaNoQ9e2D+fPjXv6BwYa8rkwApgEVEQtUvv8ANN8C990KLFq4X3KKF11VJkCiARURC0ZIlbm7vRx/BiBEwaxZoOd+IogAWEQkl2dkwdCgkJbkdjJYvh/vugxj9uo40moYkIhIqtm+HTp1g8WK3k9FLL0Hx4l5XJflEASwiEgpmzIDu3eHIEbeHb5cuXlck+UxjGiIiXjp8GHr3huuvd9OK1q5V+EYJBbCIiFe+/hoaN4YXX3TneZctg/PO87oqKSAaghYRKWjWwsSJrucbH++udL72Wq+rkgKmHrCISEHauxc6doQePVzvd906hW+UUgCLiBSUVaugTh14/323mtW8eVChgtdViUcUwCIi+S0nB555Bpo2dbcXL3brOft8XlcmHtI5YBGRIMvOsSRvTmdDWgZ1Ch/mkmEDMJ/OhxtvhHHjoFQpr0uUEKAAFhEJouwcS5cJK0ndtocGX6+i48fPk3n0IIXGvELMHf49fEXQELSISFAlb05nw9Zd9Jk7jjfef4zdcSW4sftoFl12vcJXfkc9YBGRINqesp43Jw6g9s4tvJ14DcOa9+Ro4dPYmJZBi+rlvS5PQogCWEQkWCZPptNdd3AgG+68/iHmJDQFIL6IjxoVSnhcnIQaBbCISKD273eLarz+Or6mTXn0ugEsPhiHOZpNXBEfiZVKkZRQzusqJcQogEVEApGaCu3bwzffwCOPYB57jOdjfCRvTmdjWgY1KpQgKaEcvhid/5XfUwCLiOSFtfDCCzBwIJQpAwsWwOWXA+ADWlQvr3O+8pcUwCIip+rnn93WgbNmwXXXuXWdy5TxuioJM5qGJCJyKpKToXZtmDsXRo6EmTMVvpInCmARkdzIyoIhQ6B5cyhaFFasgL59NbdX8kxD0CIif+eHH6BTJ/j8c+jWze3fW6yY11VJmFMAi4j8lenT3daBmZnw9tsuiEWCQEPQIiInc+gQ3H033HADnHMOrF2r8JWgUgCLiJxo40Zo1AjGjIEBA2DpUjj3XK+rkgijIWgRkd9YC+PHu4urihWDjz+Ga67xuiqJUOoBi4gA7NkDt9wCvXpB06awbp3CV/KVAlhEZPlySEx0F1z9+99uju9ZZ3ldlUQ4BbCIRK+cHHjqKWjWzM3nXbIEHngAYvSrUfKfzgGLSHTauRO6dHFrON98M7z6KpQq5XVVEkUUwCISfT75xC2osX+/u+iqe3etaCUFTuMsIhKxsnMsCzb9xOgFW1iw6SeyDx9x04quvdad41292i2yofAVD6gHLCIRKTvH0mXCSlK37eHQ0Wwu2PcjVT8aTrX/fg333APDh0NsrNdlShRTAItIRErenE7qtj0cPJrN9RsW8cS8l8mO8bFu9ERq977V6/JENAQtIpFpQ1oGMfv28dxHIxg5+zk2lKvGNbe9wOIaTb0uTQRQD1hEIlTjPVtp8+Z9VPp1J6Mu7sDopu05LbYINSqU8Lo0EUABLCKRxloYNYoGgwbxS9FS3NblKZacdSFxRXwkVipFUkI5rysUARTAIhJJdu2CW2+Fjz/GtGlDqfET6PZzNg3SMqhRoQRJCeXwxeiKZwkNCmARiQwLF0LnzvDLL/DCC3DPPfiMoUVZaFG9vNfVifyBLsISkfCWmQmDB0PLllCiBKxcCffeq7m9EvLUAxaR8LV1K3Ts6DZT6N4dRo+GokW9rkokVxTAIhKepk6Fnj3dhgpTpkD79l5XJHJKNAQtIuHl4EG44w646SZISIDUVIWvhCUFsIiEj/XroWFDGDsWBg2Czz+HatW8rkokTxTAIhL6rIVXXoEGDeDnn2HuXHj6aShc2OvKRPJMASwioe3XX91w8113waWXwrp1cOWVXlclEjAFsIiErqVLITERZsyAZ55x+/iW15xeiQwKYBEJPdnZ8MQTcNllUKiQC+KBAyFGv7IkcmgakoiElh073IpWycnQoYM791tCGyhI5Anoz0ljTCljzFRjzNfGmE3GmCbBKkxEotDs2VC7NqxaBRMnwqRJCl+JWIGO54wC5lhrLwBqA5sCL0lEos6RI9CvH7RuDRUrwurVblMFLScpESzPQ9DGmJLApcCtANbao8DR4JQlIlHjm2/cQhpr10KfPm56UWys11WJ5LtAesBVgV3ARGPMWmPMeGPMHxZhNcb0MsakGGNSdu3aFcDbiUhEsRbeeAPq1oUffoCZM2HUKIWvRI1AArgQUBcYY62tAxwAHjzxSdbasdba+tba+mXLlg3g7UQkYmRkQJcubpi5fn03t7d1a6+rEilQgQTwdmC7tXal//5UXCCLiPy5lBTX650yBR5/HBYsgLPP9roqkQKX5wC21v4IbDPGJPgPtQA2BqUqEYk8OTnw3HNw8cVw9Ch89hk8+ij4fF5XJuKJQOcB9wYmGWOKAN8DtwVekohEnJ9+csPNc+ZAu3YwfjyccYbXVYl4KqAAttamAvWDU4qIRKT589353j174OWX4c47Nb1IBK2EJSJBlJ1jSd6czoa0DGqWjSNp0gvEPPssXHCBC+JatbwuUSRkKIBFJCiycyxdJqwkddseSqfv4LLZzxKzYzM5PXsSM2oUxMd7XaJISNHK5iISFMmb00ndtofm6xbx0cQ+VP15O/f982EW9X9C4StyEuoBi0hQbP7uRx6b8Ty3fDmPNRUS6NNmEDtKlqdaWgYtqmsLQZETKYBFJHBffkm3vjcR9/23vNT4Jp6/pBNZvkLEF/FRo4I2UxA5GQ1Bi0jeWQsvvQQNGxJ/cB9P3TeKl67oTrY/fBMrlSIpoZzXVYqEJPWARSRvfvkFevSADz+Ea67BvP46D5YpS+PN6WxMy6BGhRIkJZTDF6MpRyInowAWkVO3ZAl07OgW2HjuObeVYEwMPqBF9fI65yuSCxqCFpHcy86GoUMhKcntWrR8OfTvDzH6VSJyqtQDFpHc2b4dOnWCxYuhc2e3qlXx4l5XJRK2FMAi8vdmzIDu3eHIEbeHb9euXlckEvY0biQif+7wYejdG66/HipXhjVrFL4iQaIAFpGT+/praNwYXnzRXWS1fDmcf77XVYlEDA1Bi8jvWQsTJ7qeb3w8zJ4NrVp5XZVIxFEPWET+Z+9eN72oRw9o1AjWrVP4iuQTBbCIOKtWQZ068P778K9/ue0DK1TwuiqRiKUAFol2OTnwzDPQtKmb57t4MTz8MPh8XlcmEtF0Dlgkmv34o7uqef58+Oc/Ydw4OP10r6sSiQrqAYtEq7lzoXZtt6zkq6+6oWeFr0iBUQCLRJujR2HQILj6aihbFlJSoFcvMNo0QaQgaQhaJJp89x106ABffAF33gkjRkBcnNdViUQlBbBItJgyBe64w11cNXWqO+crIp7RELRIpNu/363j3LEj1KoFqakKX5EQoAAWiWSpqVC/Prz+OjzyCHz2mVvTWUQ8pwAWiUTWwujRbjWrfftgwQIYNgwK6ayTSKjQT6NIpPn5ZzfkPGsWXHedW9e5TBmvqxKRE6gHLBJJkpMhMdHN8R05EmbOVPiKhCgFsEgkyMqCIUOgeXO3g9Hy5dC3r+b2ioQwDUGLhLsffoBOneDzz6FbN7d/b7FiXlclIn9DASwSzqZPd1sHZmbC22+7IBaRsKAhaJFwdOgQ3H033HADVKsGa9cqfEXCjAJYJNxs3OimF40ZAwMGwLJlcO65XlclIqdIQ9Ai4cJaGD/eXVxVrBh8/DFcc43XVYlIHqkHLBIO9uyBW25xuxY1bQrr1il8RcKcAlgk1C1fDnXqwLRp8NRTbo7vWWd5XZWIBEgBLBKqcnJc4DZr5u5//jk8+CDE6MdWJBLoHLBIKNq5E7p0cWs433wzvPoqlCrldVUiEkQKYJFQ88knbkGN/fth3Dg3z1crWolEHI1liYSKo0fdtKJrr4Uzz4SUFOjZU+ErEqHUAxYJBd9+C+3bw+rVboGN4cMhLs7rqkQkHymARbz29ttw111QuLC70rldO68rEpECoCFoEa/s3+/O9Xbp4rYQTE1V+IpEEQWwiBfWrIG6dV3vd8gQWLQI/vEPr6sSkQKkABYpSNbCyJHQuDEcPAgLF8LQoVBIZ4NEoo1+6kUKyq5dcNtt8NFH0Lo1TJwIpUt7XZWIeEQ9YJGCsGgR1K4N8+fD6NEwY4bCVyTKKYBF8lNWFjzyCLRoASVKwMqV0Lu35vaKiIagRfLNf/8LHTu6/Xq7d3c936JFva5KREKEAlgkP3zwgVvFKjsbJk+GDh28rkhEQoyGoEWC6dAhuPNOuPFGOO88WLtW4SsiJ6UAFgmW9euhQQO3c9GgQW77wHPO8boqEQlRGoIWCZS1MHYs9OvnLrSaOxeuvNLrqkQkxKkHLBKIX3+Fm25yw86XXgpffqnwFZFcUQCL5NXSpW4N5xkz4Jln3D6+5ct7XZWIhAkFsMipys6Gf/0LLrvMLSG5dCkMHAgx+nESkdwL+BywMcYHpAA7rLXXBV6SSAhLS4POnd3KVh06wJgxULKk11WJSBgKxp/sfYFNQXgdkdA2ezZcdJFbzeq112DSJIWviORZQAFsjKkItALGB6cckRB05Ii7wrl1a6hYEVavdpsqaDlJEQlAoD3gkcAgICfwUkRC0DffQJMmMGqUW8N5xQq44AKvqxKRCJDnADbGXAekW2tX/83zehljUowxKbt27crr24kUvDffhLp13ZrOM2a4tZxjY72uSkQiRCA94KZAG2PMVuAdoLkx5u0Tn2StHWutrW+trV+2bNkA3k6kgOzbB126QLduUK8erFsHbdp4XZWIRJg8B7C19iFrbUVrbRWgPbDQWts5aJWJeCElBerUcRsoDB0KCxe6874iIkGmpSglqmXnWJI3p7Nh+x6umT+Fc0c+iSlfHpKToVkzr8sTkQgWlAC21iYDycF4LZGCkp1j6TJhJT98vZUnpg/nvP+s5ovES6k7fxq+MqW9Lk9EIpx6wBK1kjenE/fZIqZ9+CwlD+/nkSvuYlqj1rywK4sWZbyuTkQindbOk+iUmUnJx4cwbtJg9sYWo23XEbxdtxWHMnPYmJbhdXUiEgXUA5bo85//QMeO1F+xgvfqXM1jST05VMRNL4or4qNGhRIeFygi0UABLNHlvffg9tsByJnyDh/uq4zZtgdzNJu4Ij4SK5UiKaGcx0WKSDRQAEt0OHgQ+vaF8eOhcWOYPJmYqlV5y38V9Ma0DGpUKEFSQjl8MVpiUkTynwJYIt9XX8Ett8DXX8NDD7n5vYULA+CLMbSoXp4W1bWPr4gULF2EJZHLWnj5ZWjQAH79FebNgyefPBa+IiJeUgBLZPrlF/jnP+Gee6B5c7ecZMuWXlclInKMAlgiz5IlkJjo9u997jn3bzldWCUioUUBLJEjOxsefxySkqBIEVi2DPr3hxh9m4tI6NFFWBIZtm+Hzp3hs8+gUyd37reE5vOKSOhSAEv4mzkTbrsNjhyBN96Arl29rkhE5G9pbE7C1+HD0KcPtG0LlSvDmjUKXxEJGwpgCU9ff+0W1HjhBejXD5Yvh/PP97oqEZFc0xC0hBdr4fXX4d57IT4eZs2C667zuioRkVOmHrCEj7173QVW3btDo0Zubq/CV0TClAJYwsOqVVCnjttM4YknYP58qFDB66pERPJMASyhLScHnnkGmjZ183wXL4bBg8Hn87oyEZGA6BywhK6ffnJXNc+b55aVHDcOTj/d66pERIJCPWAJTfPmwUUXuR7vK6/A++8rfEUkoiiAJbQcPQqDBsFVV0HZsvDFF3DHHWC0R6+IRBYNQUvo+P576NDBXXB1xx0wYoSbaiQiEoEUwBIapkxxoevzwdSp7pyviEgE0xC0eOvAATevt2NHqFULUlMVviISFRTA4p3UVKhXz61sNXiw28mocmWvqxIRKRAKYCl41ro1nBs1gowM+PRTt7hGIZ0REZHood94UrB273ZDzjNnQqtWMHGiu9pZRCTKqAcsBeezz6B2bZgzB0aOdBspKHxFJEopgCX/ZWXBY49B8+ZuWtHy5dC3r+b2ikhU0xC05K9t29wORkuWQLdu7txv8eJeVyUi4jkFsOSfDz9053szM+Gtt6BzZ68rEhEJGRqCluA7dAjuuQfatYNq1WDtWoWviMgJFMASXJs2uelFL78MAwbAsmVw7rleVyUiEnI0BC3BYS1MmAB9+kCxYvDxx3DNNV5XJSISshTAkmfZOZbkzels+WYH7V4ZSvk5M6FlS3jzTTjrLK/LExEJaQpgyZPsHEuXCSvJWb6cZ6c9TemMXUxpdxc3v/cCvkI+r8sTEQl5CmDJk+RNP1L/nbH0Tn6TH4uX4aZOz7C56oWU2/IzLaqX97o8EZGQpwCWU7dzJ+d1vZEWa5YxO+ESHr76XjJii2GOZrMxLUMBLCKSCwpgOTVz5kDXrpy9bz+PXteXt2q0PLaiVVwRHzUqlPC4QBGR8KBpSJI7R4/C/fe7K5vPPBO++ILv2rQn/rRCGCC+iI/ESqVISijndaUiImFBPWD5e99+Cx06QEoK3H03DB+OLy6Ot2q4q6A3pmVQo0IJkhLK4YvR+s4iIrmhAJa/NmkS3HknFC4M06a51a38fDGGFtXL65yviEgeaAhaTm7/frj1VreEZGIipKb+LnxFRCQwCmD5o7VroV49t4HCkCGwaBH84x9eVyUiElEUwPI/1sKoUdC4MRw4AAsXwtChUEhnKkREgk2/WcXZtQtuuw0++ghat4aJE6F0aa+rEhGJWOoBixtirl0b5s+H0aNhxgyFr4hIPlMAR7OsLHj0UWjRAkqUgJUroXfvYwtriIhI/tEQdLT673+hY0e3X2/37q7nW7So11WJiEQNBXA0+uAD6NkTsrNh8mS3yIaIiBQoDUFHk0OH3KIaN94I553nphspfEVEPKEAjhYbNkCDBvDqqzBoEHz+OZxzjtdViYhELQ1BRzprYexY6NfPXWg1dy5ceaXXVYmIRD31gCPZr7/CzTe7YedLL4Uvv1T4ioiECAVwpFq2zK3h/OGH8Mwz8MknUF6bJoiIhIo8B7AxppIxZpExZqMxZoMxpm8wC5M8ys6Gf/3L9XgLFYKlS2HgQIjR31oiIqEkkHPAWcAAa+0aY0xxYLUxZr61dmOQapNTlZbmdi9atMhd3TxmDJQs6XVVIiJyEnnuFllrd1pr1/hv7wM2AWcHqzA5RR995JaTXLkSXnvN7eOr8BURCVlBGZc0xlQB6gArg/F6cgqOHIH77oPrroOzz4bVq92mClpOUkQkpAUcwMaYYsAHQD9rbcZJHu9ljEkxxqTs2rUr0LeT433zDVx8MYwc6dZwXrECLrjA66pERCQXAgpgY0xhXPhOstZOO9lzrLVjrbX1rbX1y5YtG8jbyfHefBPq1oWtW93uRaNHQ2ys11WJiEguBXIVtAEmAJustSOCV5L8pX37oEsX6NYN6tWDdeugTRuvqxIRkVMUSA+4KdAFaG6MSfV/XBukuuRkVq92vd7Jk2HoUFi4ECpW9LoqERHJgzxPQ7LWfg7oSp+CkJPjzvM++KBbTCM5GZo187oqEREJgNaCDnXp6XDrrW4lq+uvhwkT4IwzvK5KREQCpOWRQtmCBW5u78KF8NJLMG2awldEJEIogENRZiY8/DBccQWcfjqsWgV33625vSIiEURD0KFm61a3jOSKFdCzpzv3W7So11WJiEiQKYBDyfvvw+23uz1833kHbrnF64pERCSfaAg6FBw8CL16ub17L7gAUlMVviIiEU4B7LWvvoIGDWD8eDfNaMkSqFrV66pERCSfaQjaK9bCK69A//5QqhTMmwctW3pdlYiIFBD1gL3wyy/wz3+6K5uTktxykgpfEZGooh5wPsvOsSRvTmdDWgYXVihB0q5v8HXuBD/+CMOHu60EY/R3kIhItFEA56PsHEuXCStJ3baHI4eP0nfVVC5fPAlbrSpm2TKoX9/rEkVExCMK4HyUvDmd1G17KL77JybMfo4mP3zFrJqXU/y1sSTVP9fr8kRExEMK4Hy0IS2Dizcu45mPR3Fa1lEGXHsf02o2p3+GJcnr4kRExFMK4Pxy+DDtXn+GPh+MZ335c+jdZhD/OeNs4ov4qFGhhNfViYiIxxTA+WHzZmjfnkqpqXzS4mYebNiJjBwf8UV8JFYqRVJCOa8rFBERjymAg8laeOMNuPdeiI2FWbO48tpWFNmczsa0DGpUKEFSQjl8MdpUQUQk2imAgyUjA+66CyZPdnN7334bzj4bH9CienlaVC/vdYUiIhJCNAE1GL74AurUgXffhSeegE8/hbPP9roqEREJYQrgQOTkwLPPwsUXQ1YWfPYZDB4MPp/XlYmISIjTEHRe/fQTdOsGc+e6ZSXHjYPTT/e6KhERCRPqAefFvHlQu7br8b7yitvHV+ErIiKnQAF8KjIz4YEH4KqroHRpd+73jjvA6KpmERE5NRqCzq3vv4cOHWDVKhe6I0ZAfLzXVYmISJhSAOfGO+/8r6f7/vtw441eVyQiImFOQ9B/5cAB6NHD9Xxr1nT79ip8RUQkCBTAf2bdOrdd4MSJbmrRZ59B5cpeVyUiIhFCAXwia+HFF6FRI9i71y2q8cQTUEij9SIiEjxKlePt3g3du8PMmdCqlev9li3rdVUiIhKB1AP+zWefubm9n3wCzz8Ps2YpfEVEJN8ogLOy4P/+D5o3d9OKVqyAfv00t1dERPJVdA9Bb9sGnTrBkiXQtas791u8uNdViYhIFIjeAP7wQ3e+NzMT3noLOnf2uiIREYki0TcEffgw3HMPtGsH1arBmjUKXxERKXDRFcCbNkHDhvDyy9C/PyxbBued53VVIiIShaJjCNpamDAB+vSBYsXg44/hmmu8rkpERKJY5PeA9+51S0nefjtcfLFb4UrhKyIiHovsAF6xAhITYepUePJJt4/vWWd5XZWIiEiEBnBODvz739CsmRt+XrIEHnoIYiKzuSIiEn4i7xzwzp1uTu+nn8JNN8HYsVCqlNdViYiI/E5kBfCcOS589+93wduzp1a0EhGRkBQZY7JHj8L997uLq8qXh5QUd9GVwldEREJUWPaAs3MsyZvT2ZCWQf3M3TR5tDcmJQXuvhuGD4e4OK9LFBER+UthF8DZOZYuE1aSum0PV679lB7zXuZgoULETv0A3z9v8Lo8ERGRXAm7IejkzemkbtvDoI9eYuTs59hQripter5Ico2mXpcmIiKSa2HXA96QlsGho9msPrs6e2OLM7ppe3JifGxMy6BF9fJelyciIpIrYRfAF1YoQVwRH7NqXHbsWHwRHzUqlPCwKhERkVMTdkPQSQnlSKxUivgiPgwufBMrlSIpoZzXpYmIiORa2PWAfTGGt3o0InlzOhvTMqhRoQRJCeXwxWjKkYiIhI+wC2BwIdyienmd8xURkbAVdkPQIiIikUABLCIi4gEFsIiIiAcUwCIiIh5QAIuIiHhAASwiIuKBgALYGHO1MWazMeZbY8yDwSpKREQk0uU5gI0xPuAl4BqgBtDBGFMjWIWJiIhEskB6wA2Bb62131trjwLvAG2DU5aIiEhkCySAzwa2HXd/u/+YiIiI/I18X4rSGNML6OW/u98YszmIL18G+DmIrxdK1LbwFKlti9R2gdoWrsKlbZX/7IFAAngHUOm4+xX9x37HWjsWGBvA+/wpY0yKtbZ+fry219S28BSpbYvUdoHaFq4ioW2BDEF/AZxnjKlqjCkCtAdmBqcsERGRyJbnHrC1NssYcy8wF/ABr1lrNwStMhERkQgW0Dlga+3HwMdBqiUv8mVoO0SobeEpUtsWqe0CtS1chX3bjLXW6xpERESijpaiFBER8UDYBnC4L4NpjHnNGJNujFl/3LEzjDHzjTFb/P+e7j9ujDGj/W390hhT17vK/5oxppIxZpExZqMxZoMxpq//eCS0LdYYs8oYs87ftqH+41WNMSv9bXjXf1EixpjT/Pe/9T9exdMG/A1jjM8Ys9YYM9t/P1LatdUY85UxJtUYk+I/FvbfjwDGmFLGmKnGmK+NMZuMMU0ioW3GmAT/1+u3jwxjTL9IaNvxwjKATWQsg/k6cPUJxx4EFlhrzwMW+O+Da+d5/o9ewJgCqjEvsoAB1toaQGPgHv/XJhLadgRobq2tDSQCVxtjGgNPA89ba88FfgV6+J/fA/jVf/x5//NCWV9g03H3I6VdAJdbaxOPm7YSCd+PAKOAOdbaC4DauK9f2LfNWrvZ//VKBOoBB4HpREDbfsdaG3YfQBNg7nH3HwIe8rquPLSjCrD+uPubgbP8t88CNvtvvwp0ONnzQv0DmAFcEWltA+KBNUAj3GIAhfzHj31v4mYINPHfLuR/nvG69j9pT0XcL7TmwGzAREK7/DVuBcqccCzsvx+BksB/Tvy/j4S2ndCeK4Glkdi2sOwBE7nLYJa31u703/4RKO+/HZbt9Q9N1gFWEiFt8w/TpgLpwHzgO2CPtTbL/5Tj6z/WNv/je4HSBVpw7o0EBgE5/vuliYx2AVhgnjFmtXEr80FkfD9WBXYBE/2nDsYbY4oSGW07Xntgiv92RLUtXAM44ln3Z1zYXqJujCkGfAD0s9ZmHP9YOLfNWptt3bBYRdyGJBd4W1HgjDHXAenW2tVe15JPLrHW1sUNU95jjLn0+AfD+PuxEFAXGGOtrQMc4H9DskBYtw0A/3UHbYD3T3ws3NsG4RvAuVoGMwz9ZIw5C8D/b7r/eFi11xhTGBe+k6y10/yHI6Jtv7HW7gEW4YZmSxljfptTf3z9x9rmf7wksLtgK82VpkAbY8xW3K5mzXHnFsO9XQBYa3f4/03HnUdsSGR8P24HtltrV/rvT8UFciS07TfXAGustT/570dS28I2gCN1GcyZQDf/7W6486e/He/qv9KvMbD3uGGYkGKMMcAEYJO1dsRxD0VC28oaY0r5b8fhzm1vwgXxjf6nndi239p8I7DQ/1d7SLHWPmStrWitrYL7WVpore1EmLcLwBhT1BhT/LfbuPOJ64mA70dr7Y/ANmNMgv9QC2AjEdC243Tgf8PPEFltC8+LsPw/69cC3+DOwQ32up481D8F2Alk4v6S7YE7j7YA2AJ8Cpzhf67BXfX9HfAVUN/r+v+iXZfghoW+BFL9H9dGSNsuAtb627YeGOI/Xg1YBXyLGyo7zX881n//W//j1bxuQy7amATMjpR2+duwzv+x4bffFZHw/eivNxFI8X9PfgicHkFtK4obWSl53LGIaNtvH1oJS0RExAPhOgQtIiIS1hTAIiIiHlAAi4iIeEABLCIi4gEFsIiIiAcUwCIhzhjzf8aY+/23HzfGtPyL517/VxuT+HeU6ZrHOu41xnTPy+eKyB8pgEXCiLV2iLX20794yvW4HcL+wL9qVXdgch7f/jWgdx4/V0ROoAAWCUHGmMHGmG+MMZ8DCccdf90Yc6P/9r+N23f5S2PMcGPMxbh1c5/176F6zgkv2xy3rF+W//OTjTGj/M9db4xp6D8+yhgzxH/7KmPMYmNMjLX2ILD1t+eJSGAK/f1TRKQgGWPq4ZaETMT9jK4BVp/wnNJAO+ACa601xpSy1u4xxszErWQ19SQv3fTE1wHirbWJ/g0KXgNq4rb3/MIYswQYDVxrrf1tl6QUoBluBSwRCYB6wCKhpxkw3Vp70LqdpE62zvle4DAwwRhzA27D8r9zFm77uuNNAbDWLgZK+IP8IHA7brvFF6213x33/HSgwim1RkROSgEsEob8w8gNcTvgXAfMycWnHcKt4/y7l/qT+7Vw6/CeGLax/tcRkQApgEVCz2LgemNMnH8nn9YnPsG/33JJa+3HwH1Abf9D+4Dif/K6m4BzTzh2i//1LsHtILPXGFMZGADUAa4xxjQ67vnn4zaiEJEAKYBFQoy1dg3wLm4Hn09w22+eqDgw2xjzJfA50N9//B1goDFm7UkuwvoEuPSEY4eNMWuBV4Aex20neb+1Ng23S9d4Y8xvPeemuKFpEQmQdkMSiSLGmOnAIGvtFmNMMi5oU3L5uXWA/tbaLvlZo0i0UA9YJLo8iLsYKy/KAI8GsRaRqKYesIiIiAfUAxYREfGAAlhERMQDCmAREREPKIBFREQ8oAAWERHxgAJYRETEA/8PK9KpleJY7eMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# aruco vs laser displacement sensor\n",
    "\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import cv2\n",
    "import glob\n",
    "\n",
    "#df = pd.DataFrame(\n",
    "#    {\n",
    "#        \"file name\": ['setup1_img1.jpg','setup1_img2.jpg','setup1_img3.jpg','setup1_img4.jpg','setup1_img5.jpg'\n",
    "#                     ,'setup1_img6.jpg','setup1_img7.jpg','setup1_img8.jpg','setup1_img9.jpg'],\n",
    "#        \"sensor (mm)\": [0.5064,1.5353,2.3412,3.4370,4.5723,5.8718,7.2475,8.3390,9.7316]\n",
    "#    }\n",
    "#)\n",
    "\n",
    "#df = pd.DataFrame(\n",
    "#    {\n",
    "#        \"file name\": ['setup2_img1.jpg','setup2_img2.jpg','setup2_img3.jpg','setup2_img4.jpg','setup2_img5.jpg'\n",
    "#                     ,'setup2_img6.jpg','setup2_img7.jpg','setup2_img8.jpg','setup2_img9.jpg'],\n",
    "#        \"sensor (mm)\": [0.5393,1.6483,2.8225,3.7052,4.8530,6.0415,7.2658,8.4812,9.6174]\n",
    "#    }\n",
    "#)\n",
    "\n",
    "#df = pd.DataFrame(\n",
    "#    {\n",
    "#        \"file name\": ['setup3_img1.jpg','setup3_img2.jpg','setup3_img3.jpg','setup3_img4.jpg','setup3_img5.jpg'\n",
    "#                     ,'setup3_img6.jpg','setup3_img7.jpg','setup3_img8.jpg','setup3_img9.jpg'],\n",
    "#        \"sensor (mm)\": [0.7023,2.1092,2.9332,4.2508,5.2229,6.3942,7.8878,8.8936,9.7511]\n",
    "#    }\n",
    "#)\n",
    "\n",
    "#df = pd.DataFrame(\n",
    "#    {\n",
    "#        \"file name\": ['setup4_img1.jpg','setup4_img2.jpg','setup4_img3.jpg','setup4_img4.jpg','setup4_img5.jpg'\n",
    "#                     ,'setup4_img6.jpg','setup4_img7.jpg','setup4_img8.jpg','setup4_img9.jpg'],\n",
    "#       \"sensor (mm)\": [0.2644,1.5617,3.0174,4.2728,5.6098,6.8563,7.5603,8.6505,9.9494]\n",
    "#    }\n",
    "#)\n",
    "\n",
    "df = pd.DataFrame(\n",
    "    {\n",
    "        \"file name\": ['setup5_img1.jpg','setup5_img2.jpg','setup5_img3.jpg','setup5_img4.jpg','setup5_img5.jpg'\n",
    "                     ,'setup5_img6.jpg','setup5_img7.jpg','setup5_img8.jpg','setup5_img9.jpg'],\n",
    "        \"sensor (mm)\": [0.4588,1.4462,2.6874,3.8038,5.7312,6.5623,7.2250,8.2916,9.7053]\n",
    "    }\n",
    ")\n",
    "\n",
    "df.insert(2,\"cX\",0)\n",
    "df.insert(3,\"cY\",0)\n",
    "df.insert(4,\"dist (px)\",0)\n",
    "\n",
    "aruco_dict = cv2.aruco.Dictionary_get(cv2.aruco.DICT_6X6_100)\n",
    "# marker created via cv2.aruco.drawMarker(aruco_dict, 1, 8)\n",
    "images = glob.glob('aruco/undistorted/setup5*.jpg')\n",
    "for fname in images:\n",
    "    img = cv2.imread(fname)\n",
    "    gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)\n",
    "    ret,gray_thresh = cv2.threshold(gray,50,255,cv2.THRESH_BINARY)\n",
    "    fname_short = fname.rsplit(\"/\")[2]\n",
    "    aruco_params = cv2.aruco.DetectorParameters_create()\n",
    "    corners = cv2.aruco.detectMarkers(gray_thresh, aruco_dict, parameters=aruco_params)\n",
    "    cX = (corners[0][0][0][0][0] + corners[0][0][0][2][0]) / 2.0\n",
    "    cY = (corners[0][0][0][0][1] + corners[0][0][0][2][1]) / 2.0\n",
    "    df.loc[df.index[df[\"file name\"].str.contains(fname_short)][0],\"cX\"] = cX\n",
    "    df.loc[df.index[df[\"file name\"].str.contains(fname_short)][0],\"cY\"] = cY\n",
    "    \n",
    "cX_zero = df.loc[0, \"cX\"]\n",
    "cY_zero = df.loc[0, \"cY\"]\n",
    "for fname in images:\n",
    "    fname_short = fname.rsplit(\"/\")[2]\n",
    "    cX = df.loc[df.index[df[\"file name\"].str.contains(fname_short)][0],\"cX\"] \n",
    "    cY = df.loc[df.index[df[\"file name\"].str.contains(fname_short)][0],\"cY\"] \n",
    "    df.loc[df.index[df[\"file name\"].str.contains(fname_short)][0],\"dist (px)\"] = ((cX - cX_zero)**2 + (cY - cY_zero)**2)**(1/2)\n",
    "    \n",
    "\n",
    "d = np.polyfit(df['dist (px)'],df['sensor (mm)'],1)\n",
    "f = np.poly1d(d)\n",
    "df.insert(5,'reg',f(df['dist (px)']))\n",
    "ax = df.plot(x='dist (px)',y='sensor (mm)',style='o',ms=5,figsize=(8,6))\n",
    "df.plot(x='dist (px)',y='reg',color='Red',ax=ax)\n",
    "\n",
    "df.insert(6,'error (mm)',df['sensor (mm)'] - df['reg'])\n",
    "\n",
    "print(df['error (mm)'])\n",
    "\n",
    "print((df['error (mm)']).std())\n",
    "\n",
    "print(d)\n",
    "plt.savefig('aruco_cal.png', dpi=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}