commit 8027439bc102b6e27189d00d19e07091ce3e7d36 Author: wuhanstudio Date: Mon Mar 11 13:00:35 2024 +0000 Initial Commit diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..e0f09ae --- /dev/null +++ b/.gitignore @@ -0,0 +1,4 @@ +.ipynb_checkpoints +data/EXCV10 +data/MaskedFace +__pycache__ diff --git a/ECMM426.pdf b/ECMM426.pdf new file mode 100644 index 0000000..9ec8bf5 Binary files /dev/null and b/ECMM426.pdf differ diff --git a/Question 1.ipynb b/Question 1.ipynb new file mode 100644 index 0000000..dc4574a --- /dev/null +++ b/Question 1.ipynb @@ -0,0 +1,676 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "f94476eb", + "metadata": {}, + "source": [ + "## Question 1 (14 marks)\n", + "\n", + "Write two functions compute_gradient_magnitude(gr_im, kx, ky) and\n", + "compute_gradient_direction(gr_im, kx, ky) to compute the magnitude and direction of\n", + "gradient of the grey image gr_im with the horizontal kernel kx and vertical kernel ky." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "edf40439", + "metadata": {}, + "outputs": [], + "source": [ + "import cv2\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "cbfcddf7", + "metadata": {}, + "outputs": [], + "source": [ + "image_path = \"data/shapes.png\"\n", + "gr_im = cv2.imread(image_path, cv2.IMREAD_GRAYSCALE)" + ] + }, + { + "cell_type": "markdown", + "id": "da7b5144", + "metadata": {}, + "source": [ + "This is the default kernel for Sobel Filter." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "7ba9e609", + "metadata": {}, + "outputs": [], + "source": [ + "# Sobel Filter\n", + "kx_conv = np.array([\n", + " [1, 0, -1], \n", + " [2, 0, -2], \n", + " [1, 0, -1]\n", + "])\n", + "\n", + "ky_conv = np.array([\n", + " [1, 2, 1], \n", + " [0, 0, 0],\n", + " [-1, -2, -1]\n", + "])" + ] + }, + { + "cell_type": "markdown", + "id": "8f3279b4", + "metadata": {}, + "source": [ + "**However, some students used the cross-correlation kernel, and thus produced different results.**" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "1316ddae", + "metadata": {}, + "outputs": [], + "source": [ + "kx_cross = np.array([\n", + " [-1, 0, 1], \n", + " [-2, 0, 2], \n", + " [-1, 0, 1]])\n", + "\n", + "ky_cross = np.array([\n", + " [-1, -2, -1], \n", + " [0, 0, 0],\n", + " [1, 2, 1]\n", + "])" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "6ccd97f3", + "metadata": {}, + "outputs": [], + "source": [ + "# kx_cross = np.flip(np.flip(kx_conv, 0), 1)\n", + "# ky_cross = np.flip(np.flip(ky_conv, 0), 1)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "cf6af87b", + "metadata": { + "scrolled": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAawAAAGiCAYAAAC7wvLcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABvo0lEQVR4nO3dd3gU5drH8e/MbEk2nUAKhIQAgST0TkQRJYKICooFDwgqNgwqoKh4PPYjvnosx4oVUEQUu6gggsJRQxeE0GuAkARCekjZ3Xn/wESjIElIdnd278917XXBzuzOPdlkfjvzPPM8iq7rOkIIIYSHU91dgBBCCFEXElhCCCEMQQJLCCGEIUhgCSGEMAQJLCGEEIYggSWEEMIQJLCEEEIYggSWEEIIQ5DAEkIIYQgSWEIIIQzBbYH18ssv06ZNG/z8/OjXrx+rV692VylCCCEMwC2B9cEHHzB16lQeeugh1q9fT7du3Rg6dCi5ubnuKEcIIYQBKO4Y/LZfv3706dOHl156CQCn00nr1q25/fbbue+++1xdjhBCCAMwuXqDlZWVrFu3junTp9c8p6oqqamppKenn/Q1FRUVVFRU1Pzf6XRy7NgxwsPDURSlyWsWQgjRuHRdp7i4mJYtW6KqdbvY5/LAOnr0KA6Hg8jIyFrPR0ZGsm3btpO+ZsaMGTzyyCOuKE8IIYQLHThwgJiYmDqt6/LAaojp06czderUmv8XFhYSGxvL2VyECbMbKxNCCNEQdqr4ka8JCgqq82tcHljNmzdH0zRycnJqPZ+Tk0NUVNRJX2O1WrFarX953oQZkyKBJYQQhvNb74n6NOu4vJegxWKhV69eLF26tOY5p9PJ0qVLSUlJcXU5QgghDMItlwSnTp3K+PHj6d27N3379uX555+ntLSU66+/3h3lCCGEMAC3BNbVV1/NkSNHePDBB8nOzqZ79+4sWrToLx0xhBBCiGpuuQ/rTBUVFRESEsIgRkgblhBCGJBdr+IHPqewsJDg4OA6vUbGEhRCCGEIElhCCCEMQQJLCCGEIRjixuHGogYFofj5ubsM71NViaOg0N1V+Bw1KAj9+HF0u93dpYh6UKxWVH8/HIVFYLwuBG7lU4FVPDSZ2Ck7CDBVursUr+HUFTa+2YXwN08+DqRoIqpG4UWdCNlagP7ryYc0Ex5EUdBCgqnq0pbDZ/lTmlBJwltVKOkb3V2ZofhUYFXZVF6N+5oQ1d/dpXiVrsFd3V2Cz9ES4mk2cT+7f4inzQ4/nOXl7i5JnIQaEADtWnOkbxhFg8u4vesyxgZvRVUUzml+A60eTsa5YYu7yzQMnwosIbyB6ufH3msiWNb2KWYEnseOLzvAugx3lyV+o1itaBEtKOrTioOpOsP7bOSViDeJMfljVjTABsB3vd5kwD1pJDzeAceWHe4t2iAksIQwGL1TOy66ZCXRpkCmR3zPuaN60y5DzrLcSTGZ0JqHc7xra7IGWOgwaA9Pt36FrhYHNtUCBP7lNRFaAF+c9QrD06aQ9HgU9sPZri/cYCSwhDAQ1WZj15XBzI34HxBAtCmQMcOX89NXfVF+2uDu8nyLoqAGBuJMakNOnyAYcoyHkudznv+R35odtN8ep5ZksbHgohe5puhO2j9TheNonktKNyoJLCEMxNGjA2OHLae5FlDz3OTwdbxzxUA6bAjAWVrqxup8g2K2oLVuSWHPSA4NdXJ5z3VMabGCSK36kl/92sh7WS28efWrTCy6jdiXNuEsLm6awr2ABJYQBqEGBbHjSj/mNVsD/B5YIao/04Z+yYKFF2Jaus59BXozVcMU0Zzy5BgOnm8h/qxM/tPmFTpbqghU/TjZJb/6GOgH/7rufZ7KH03kOxtxlpU1Tt1eRgJLCIOoSOnIvy78pNbZVbVxwXv5v9EqSWuCcRQVuaE6L6QoqDYbeoc2HB4Ygn1gIf/q/Ann+h8gQrOhKae/5FcfVwbmseu271lYMoiwjzZIm+RJSGAJYQBaaAg7roKrAg8Clr8st6kWZgz8iJnnXIH1qzWuL9CLqH5+KLGtyOsXQe6gKq7r/TPXhK4h3uT32yW/MzubOhVNUZkWvokdaREcONoFy+K1cmPxn0hgCeHpFIXSczry6qA5v/U4O7lLA3L45zV2Oq6NwJGT68ICvYCqYWoZRXmHKDIvtNBrwHb+Gb2A3lYHVsXMHy/BNiWrYuaV1t9yxT1BOEq7o/64QULrDySwhPBwpqhIsseWM9i/jL+7BGVTLTzTbwFPDxxD4EdH5EBXB1poCI4OsRwcHIQ1JY8Hkz7gbL8cwlR/NEXFHcOtBqp+fNDhI4ZOv5bQBzuhr8uQz/I3ElhCeDJF4dh5bXi7zyu/XY76e0P8j3HvFccJ/TFS7us5BdVmQ2ndktyBLTh2dgVTei9lRFAGrTTbbyHlmrOpvxOi+rOwyxwGTJtI28c64twsw2+BjNYuhEczxbTC/o9j9LfWbX2bauGNPu+Qe2E8KErTFmdAWvNwii7uyuGnTLw0/SW2pr7G7WH7iTUF/hZWnqO5FsCi/q+w7Y4gTPFx7i7HI3jWJySE+J2qkZvamjc6vVuvg+kAqxPr1TmY2sQ2YXHG5DiaR8jSHQS/Gcw1393KnKI4Spye2xsv3hzI+6kz2XlLS7QWLdxdjttJYAnhobS2sTQbe4AuFnP9XqeovJE4l6zhrUBtvG7X3sKRdwy/L1eTdO8u5t8+jK4fT+bpY+0odB53d2kn1d9P4/UrX2PvxAS05uHuLsetJLCE8ECKycShi6N5pf38Bl2qSrLYaHXFXrSObZugOu/gyM/H/N06Ev+1hcW3n0uvD6bwyJFk8h2ed9PuIH8n/xn3Nnvu6IAW3szd5biNBJYQHkhNiKfLVVtoZ274PT8vxi9g/8gWqDJp6d9yFBWhfb+ehIcz+DGtL/3m3sUjR5LJdXjWMFfDbeU8849ZHLwuES042N3luIUElhAeRjGZ2D+yBU/HLDyj94k12eg/4lf0Tu0aqTLv5iwuRv1xA+2e2MyPaX05e+7d3HaoP3urStxdWo0L/cu49vrFHLmi04m5tnyMBJYQHkbplMC5l68n2nRmIypoispDLRexb0Qwqs3WSNV5v5rgemwje25L4OKZ93DpzgvJqDyOQ3e6tTZNUbkjbBs9b91A8bDOPnf2LIElhAdRbTZ2Xx3Ko1FLG+X9Wmk2Lhy+BkePDo3yfr7EWVaGvmYTsf/dgP3WIK554S4u3DaCdRWVbguuKt3Bnqoqws2lHEvUUFtGuaUOd5Ebh4XwII7uCVwzfMVJB7htCE1ReShiBX0u702CTD/SIM6yMti6k5Y796IujOPGIZMJuTSLpxIW0MuiNfn9WyXOcrZXqczMPY/vd3UgcKU/EevLaLN9O/ZjBU26bU8jgSWEh1ADAth5mY0F4Wup75xKfydMszHpwkV8+dVgTMtk+pGG0u12HDt2E7lnP+q3cUw87w7US/J4Mfl9elmp00gk9TW/OIz7V11Gi++shG0ppsPufTgKi8DpwNHoW/N8cklQCA9R2S+RSRd989tstY3rxpBt7LlSQwsNafT39jXVwdXirTVE3F7BHY9MovvKcXxRaqNCr2rUbbUxH0UvNRG+7hj6+q048vPB6YtRdYIElhAeQAsOZs+VGjeGNM2YcYGqHw8M+oLi8xJlyKZGotvt2PfuJ+zd1bS5I58n/zWOLituZHZRRKONntHfT2PxRc9x+AmFiqE9fa6TxZ9JYAnhAcrO6ch/Bs//bfbapnFNUCZZoyoxRUY02TZ8ktOB/VAWQR+uIWFqDrOmjaTLV7fzf3kJHG2Ee7k6mAP4X6/ZdH/8Fw5O6okprnUjFG1MElhCuJkWFsaBa+xcGpDfpNuxqRae6/cB+YNkYNwm4XRgP5yN38I1JN2/m8V3nUvfT+7i8aOJZxxcgaofT0et4sWJM9nyzyjU7skoJt/rgiCBJYQ7qRqFF3Tk9ZR3mqTR/s+G2gopuKIEU6uWTb4tn6XrOPKOYV6yjsSHt7P0nnPos3AK7xWf2TiAZkVjkL+TZUOfo+jJcsou9r1LhBJYQriRqVU0RdcUMcivcRvrT8WqmHmxx/scSY2VgXGbmGq1QkQ4Be3NtIjNJ8HSOPOTxZsDWdx5Hmc/tJKsW3piivade7F875xSCE+hKCemD+n2EpoLzq6qDfKronxkAaal0dgPHHTZdn2FGhCAntyWA+cHETN0P6/F/5dulhNfFhpLoOrHIxG/kDgxi8eir6TDa1bse/c32vt7KgksIdzEFNOKihEF9LG6tj1JU1Se6/Ih08+/mbB3D/t0N+lGoyiYIiMo7hfH4bM0Bp67iScjl5JkURs1qP7IrGiMCcqlw1UvcU2LW0h4Owx11WZ0u71JtucJJLCEcAdVI/eC1rza9WW3zHQ7yK8KrjqKaXkr7PsyXb59b6FYrahtY8kd0JySISVM7bKQkYE7idACgDpOE30GNEWlvx8sS32em9tcQ/FrvQlZtAVHUVGTb9sdpA1LCDcwxcVgvSqH/k1/TDspTVF5Nfk9sobH+GRvs8aihoZQ0DUc6xU5fNj3DW4OyfotrFwr3hzI5x0/JvWf/yNzYmev7fougSWEiykmE1nDW/Fa4ntuObuq1stqIfqKfSjJ7d1Wg9E5cnIJ+TqDgEeDuHrWVP6x9zwyKt0zc7FNtfBg8008e9MbbJkeDf27opgtbqmlqUhgCeFiats42l29g06Wxh+Cqb6eb7uAAxc187nu0Y3JWVyM8vNG2jy7ibwpMYx+5S6Gbr2Y1RVVLh/VXVNUhtiq+G7Ys2gzjpJ/TS+vmuxRAksIF1L9/Nh/RSTPx33m7lKAE6ModL5kGyS2dXcphucsLobVm2j13Fq0223c8vSdDNh4FSvKT0wL4krtzIF83uFLrrn3G/ZO6YypbRuXbr+pSGAJ4UJ6UjvOGrGRmDOcnLExPRrzJfsvDpWzrEaiV1Xi2LKDyDfW0uwOB1P/PZHea8byRamNMmely+owKxqTw/Yxa/yLbHskDOc5PVCsbmo0bST1DqwVK1ZwySWX0LJlSxRF4bPPPqu1XNd1HnzwQaKjo/H39yc1NZWdO3fWWufYsWOMGTOG4OBgQkNDmTBhAiUlnjMNtRBNQfXzY+/lwfy75bfuLqWWDuYABo1Yj7O7TPLYmPSqShy79tJ8zhpi7izh34+Op+fPE3inqDmFTte1c/X300gf9BLaI7nkX9XT0JcI6x1YpaWldOvWjZdffvmky5966ileeOEFZs6cyapVqwgICGDo0KGUl/8+evGYMWPIyMhgyZIlLFy4kBUrVnDzzTc3fC+EMABHz46MvDjdLb3ITufhqKXsuSwALSzM3aV4Hd1ux77/AKHvrabd1GPMfPgKev5wG7OLIsh3lLmkhggtgE86fMrY+75m3x3GvURY78AaNmwYjz/+OJdddtlflum6zvPPP88DDzzAiBEj6Nq1K++88w5ZWVk1Z2Jbt25l0aJFvPnmm/Tr14+zzz6bF198kfnz55OVlXXGOySEJ1IDAthzuT93Nf/J3aWcVIQWQNrF31B6doIMjNtUnA7sBw8R9OEaEu/O4s1/XkbPRXfyfH6bRhnV/XSsiokRQRkkDtlJznnRhpwbrVHbsPbu3Ut2djapqak1z4WEhNCvXz/S09MBSE9PJzQ0lN69e9esk5qaiqqqrFq16qTvW1FRQVFRUa2HEEZS1bcjt134rUeeXVW7PmQrmZc70Jo3d3cp3s3pwJ6dQ8DHq0iavodP7x1Cn6+m8PSxdk1yqbBCr2JdRSVX7h7KBe9O49iTbYhYuPvEzMUG06h3DGZnnxjcMTIystbzkZGRNcuys7OJiKg9H4/JZKJZs2Y16/zZjBkzeOSRRxqzVCFcRrXZ2H2RletDNgM2d5dzSiGqP3f3+5YP+w/Db+FR0HV3l+T1HEfz8Pu6gORfIvnsrFRevfhcHuy3kMsC95/xzNOFzuN8U9qSZ3elUraiBa2+L6bdtgwcRUUYdTAuQ/QSnD59OoWFhTWPAwcOuLskIerM3qcjd170NWGa54ZVtbHBO8gc6USLaOHuUnzHbxNABn60msS7D/D2fSPp/s0dzCxoRW4DLxVW6FVMP3weTz87mhZpFbR6ehWs3mT4IZsaNbCiok4Mc5+Tk1Pr+ZycnJplUVFR5Obm1lput9s5duxYzTp/ZrVaCQ4OrvUQwgi04GB2jTYzIWTn6Vf2ACGqPw+ctZDiATLJo8vpOo4jR/D/fA1J0/fwwdSL6P/FVP4vL4GD9pJ63YRsVcyktfieYz0cJz5HLxnguFEDKz4+nqioKJYuXVrzXFFREatWrSIlJQWAlJQUCgoKWLduXc06y5Ytw+l00q9fv8YsRwi3O57SgUfO+wSbapwhcq4I2kvWZZWYIiNOv7JofLqO42gelsVrSXxgB99OHsi5H93Nv3K7s7eq7sHVyeLPB0NfZuvdUZji45q4aNeod2CVlJSwYcMGNmzYAJzoaLFhwwYyMzNRFIXJkyfz+OOP88UXX7Bp0ybGjRtHy5YtGTlyJABJSUlceOGF3HTTTaxevZqffvqJSZMmMXr0aFq2lFlQhffQwsLYP8bJ1UGH3V1KvYSo/jzV72PyBstZllvpOo78fExL19Hx4S2sntqbIe9P467svpQ4y0//eqCv1cz8i15i69QorxgQt96BtXbtWnr06EGPHj0AmDp1Kj169ODBBx8E4J577uH222/n5ptvpk+fPpSUlLBo0SL8/nAX/XvvvUdiYiKDBw/moosu4uyzz+b1119vpF0SwgMoCiUDE3g+ZX6TzYfUlC4OyKP0iiJMMa3cXYpvUxRUmw2leTOOt7Cgm3SOVgRS4Kz7nFd9rWbmDX+FbXe2wtTK2CcF9e4lOGjQIPS/6T2kKAqPPvoojz766CnXadasGfPmzavvpoUwDFPLaHLHHGeorRAwXmBZFTMzu83lzgvTCH/zkPQYdDHVZkONiqCoeySHU1Ta9TzAQ3Fv0dNaQHMtAKjf0F79/TRmjZzJjeUTSXjhRLd6I5KJcIRobKrG0fPjeL33q4Y8u6rW16pjHZWDtrwtjh273V2OT9BatKCiSywHhliI73OAf8W9TR9r4R96mDb8Pr6BfvD61a9xi+MW2r2k48jJPf2LPIwhurULYSSmVtE4r84jxWrsnllmReO5jh9w6KJIww+aahSKn5WKZiaIKyMtdhln+5U26u0Qg/ydvPaP19g9qR2mqMjTv8DDSGAJ0ZhUjaxLY3mn82zMiubuas5YD4tK65F7UWT6EZewHzhI4GfrSLjvGP9+7MRguR+XBNe5k0VdVIfWrrS2aC2Mdb+dBJYQjUhr34a4q3aTZPH8m4Trwqxo/F+bT8i8KEymH3GRmsFy554YLPe56dfQ5dtJvJgfR66jtFEmhawOrX23JKCFN2uEql1DAkuIRqKYTGReHskL8R+7u5RGlWi20uXibeid27u7FN/y22C5AR+vIumu3Xxy1xDOeedubjs0gL1VZz4d0yB/J0+Mf4fMmxINM0q/BJYQjURJbE/3S7cQ60GTMzYGTVGZEfMF+4cHy1mWmzjy87EsWkPbJzez984OXPjONG471J9M+5kF18iAEqZf/wGHrksyRGhJYAnRCFQ/P/aPaMbTMQvdXUqTiDcHMvTS1Th6dHR3KT7NWVyMkr6RtjM2sue2BC565R5G7z2fjMrjDb5UOCYoj3tu/YADE5I8fnJHCSwhGoGjR0eGXL6aaC87u/qjByKWs3eEDTXAc6dI8RXOsjL0NZuI+c9qiq4PY8yzdzFw0xWsKD8x8G19jQnK48lb3+bAzZ1Rg4KaoOLGIYElxBlSAwLYO8LG/RHL3V1Kk2quBTD+ou8pPydZhmzyELrdjmPHbiJfWkXoRAd3PTaRHunX80Wprd49C4fbyrnvpg/IGeu5oSWBJcQZquyXyMSLF3v05IyN5fZmv7BvhCLTj3gapwP73v00m72S+DuPMeOhcXT5dhIzC1rVa1LIMUF5XD/pa3LGeGZoSWAJcQa00BD2XKGRFrbd3aW4RIjqz/RBv00/ohr/PjOvo+vYD2UR/P4qku7Zy4LbL6THJ5N5/GhinacouS10L+Nu/4YjV3dGMXvWLAMSWEKcgbIBHXjk/E8MPQRTff0jaA+HLrVjipSzLI/12xQlpqXr6PjPLfzvxj4MnjONGw+cy9bKsr8NLk1RuT10D/+YvJj8a3p5VGhJYAnRQFp4M/Zf4eSaIGMOJNpQgaofj6V8Rt7gNigmGY7U0zmLi2H1JtrO2MjhibFc9crdXLbrIjZUVFCln3z4ME1RuTNsF5fc/T1Fl/f0mNCSwBKigYrOS+C5sz/wiiGY6uuqwFyOX1GImhDv7lJEHTnLytB/yaDVs6ux3xTAdc9OYciWy1lZ7jhpz0JNUZkWvokud22k+LKeHvHlRAJLiAYwxbTi6FVlXGIrcncpbmFWNN7uNocDl7SQm4kNpqZn4cur8J+ocufDk+i58joWlVn/0rPQqph5ruVyOt+9kZKRvdweWhJYQtSXqnFkcCwv95qHpvjun1B3i4mkS7fLkE1G5XTg2LWX0LmraXNHPjOmjKfLkrS/9Cy0qRaea7mcjndncHyYe8+0fPevTYgGMkU0p+ryfAb51f8GTW+iKSqPt/5tyCabdwz265OcDuyHsvBbuIakKbuYf/dFdF90O28VRtWccdlUCy/GLCPy3t1UndvNbT1EJbCEqA9VI2d4W17r+q5Pn11V62AO4OyLNuLonuDuUsSZ0nUcBYVYv15D0tQdvHf7cLp9NJlHjiST7yjDplqY3eYbzP/MxjHQPaElf3FC1IMpthVBo7Poa/Wdbuyn86/oxSeGbPLAG01FwziKijB/t44O/9zEj2l9SZl1F5MO9SPHUckHHT7C8lA2DjecaUlgCVFHisnEwRExvNnhPXeX4lFiTYFcmroKe3dpy/I2ztJS1B83EP/EL+wZH8elL97DhL2X8ESbT7E8mI1zQFeXDtMlgSVEHanxsbQbtZN2Zu8d4LahprRYwd6RVjnL8lLO8nIcW3bQ6vnVlN0QzIQnJ1NY4UfOXeXoZ3VzWWhJYAlRB4rJROaoKF5o86m7S/FIMaZArk1dwfGzE2VgXC+m2+04du6hxVtrCL2xgoAPQjjaxR+le7JLLg9KYAlRB0pye1Iu20iMF08fcqbuaLaW/SNBa97c3aWIJqbb7dgPHiL4/ZVELdiO08+EqVV0k39ZkcAS4jTUgAD2XBnGQ9GL3V2KRwvTbDx07ucUndtWzrJ8iCPvGEr6RhzZuU2+LQksIU6jsl8i40csk7OrOrgiMJOcyyswtYx2dynCxfSqStD1Jt2GBJYQf0MLDWHPlRppzTa4uxRDCFT9eLrPRxwdHCdnWaLRSWAJcSqKQsmgjjx+/seEqP7ursYwhtnyOX55AabYGHeXIryMBJYQp2CKiuTQFVVcEZjt7lIMxaqYearzJ+Skxsgkj6JRSWAJcTKKQsHZcTzTb4FPTc7YWM7zL8E54himODnLEo1HAkuIkzC1jKZgdAnDbYXuLsWQrIqZFzu/z6GLW8n0I6LRSGAJ8WeKwtHz43ih+3yfnJyxsfS3QuTITGjfxt2lCC8hgSXEn5jiWmO/Ko+BfpXuLsXQNEXl+XYfcuDiZnKWJRqF++c8diHFqbO3SiVELXF3KV7DCSgOd1fRiFSNI4NaMa3DfA7aj59+ffG3LAo0Pz8LFrWFDVvcXY4wOEXXm/hOryZQVFRESEgIgxiBqR4N4lqHdlS0Dm26wnyUNbsER8Z2d5fRKBSzBbp3pCrY4u5SvIbi0LEcKsCxc4+7SxEexK5X8QOfU1hYSHBwcJ1e41NnWI4duzHtcHcV3sebTrD0qkpYs8m3/jBcwJt+R4T7SBuWEEIIQ5DAEkIIYQgSWEIIIQxBAksIIYQhSGAJIYQwBAksIYQQhlCvwJoxYwZ9+vQhKCiIiIgIRo4cyfbtte+/KS8vJy0tjfDwcAIDAxk1ahQ5OTm11snMzGT48OHYbDYiIiKYNm0adrv9zPdGCCGE16pXYC1fvpy0tDRWrlzJkiVLqKqqYsiQIZSWltasM2XKFL788ksWLFjA8uXLycrK4vLLL69Z7nA4GD58OJWVlfz888/MmTOH2bNn8+CDDzbeXgkhhPA6ZzTSxZEjR4iIiGD58uUMHDiQwsJCWrRowbx587jiiisA2LZtG0lJSaSnp9O/f3+++eYbLr74YrKysoiMjARg5syZ3HvvvRw5cgSL5fQjDDR0pAshhBCeoSEjXZxRG1Zh4YmpF5o1awbAunXrqKqqIjU1tWadxMREYmNjSU9PByA9PZ0uXbrUhBXA0KFDKSoqIiMj46TbqaiooKioqNZDCCGEb2lwYDmdTiZPnsyAAQPo3LkzANnZ2VgsFkJDQ2utGxkZSXZ2ds06fwyr6uXVy05mxowZhISE1Dxat27d0LKFEEIYVIMDKy0tjc2bNzN//vzGrOekpk+fTmFhYc3jwIEDTb5NIYQQnqVBY3xOmjSJhQsXsmLFCmJifp8COyoqisrKSgoKCmqdZeXk5BAVFVWzzurVq2u9X3Uvwup1/sxqtWK1WhtSqhBCCC9RrzMsXdeZNGkSn376KcuWLSM+Pr7W8l69emE2m1m6dGnNc9u3byczM5OUlBQAUlJS2LRpE7m5uTXrLFmyhODgYJKTk89kX4QQQnixep1hpaWlMW/ePD7//HOCgoJq2pxCQkLw9/cnJCSECRMmMHXqVJo1a0ZwcDC33347KSkp9O/fH4AhQ4aQnJzMtddey1NPPUV2djYPPPAAaWlpchYlhBDilOrVrV1RlJM+P2vWLK677jrgxI3Dd911F++//z4VFRUMHTqUV155pdblvv379zNx4kR++OEHAgICGD9+PE8++SQmU93yU7q1CyGEsTWkW7tPzTgshBDCM7j8PiwhhBDCVSSwhBBCGIIElhBCCEOQwBJCCGEIElhCCCEMQQJLCCGEIUhgCSGEMAQJLCGEEIYggSWEEMIQJLCEEEIYggSWEEIIQ5DAEkIIYQgSWEIIIQxBAksIIYQhSGAJIYQwBAksIYQQhiCBJYQQwhAksIQQQhiCBJYQQghDkMASQghhCBJYQgghDEECSwghhCFIYAkhhDAECSwhhBCGIIElhBDCECSwhBBCGIIElhBCCEOQwBJCCGEIJncXIISnUKxW9IoKd5fhHoqCarOB6qLvsA4HzrIy12xLeA0JLCE4EVZlw7oR9Mth7PsPuLscl9OahbH/5kSOJ5WjqHqTbstZoRG11ETIh2vR7fYm3ZbwLhJYQgDOXok4bz3KwcWtafnqEZzl5e4uybUUBWfPYjanvIVNtTTppnZUlXLZ7mmENOlWhDeSNizh87QWLdg51sqC5HfoefUmHL0S3V2SWyhK055ZVdNwzXaE95HAEr5N1Tg2tB3/l/oB0aZAHmv5DbtGW9HCwtxdmRDiTySwhE8zxcVw/IpCRgQcBSDGFMiU8xdRcm4CKIqbqxNC/JEElvBZqp8fBy5vxdvd5mBVzDXPXx+yndx/HMfUJtaN1Qkh/kwCS/gsZ/cOpFz9C72stTsZBKp+vN77XQ5d3ArFanVTdUKIP5PAEj5Ja9GCnWP8eThqyUmXD7A66XzNFvQeHV1cmRDiVCSwhM9RrFaODm/PE0M+JNoUeNJ1NEXlsVYL2T0qAC1UOmAL4QkksITv6dqBoLGHGBV49G9XizcHcudFX1M0OBFUzUXFCSFORQJL+BSteTi7Rgcwq8M8zMrpQ+j6kO3kjS7D1Ka1C6oTQvydegXWq6++SteuXQkODiY4OJiUlBS++eabmuXl5eWkpaURHh5OYGAgo0aNIicnp9Z7ZGZmMnz4cGw2GxEREUybNg27DM8CgGJu2hEGfJ6qUTSoPVOGfUXsKS4F/lmg6sec3m+z7+qW0gFDeBxf+52sV2DFxMTw5JNPsm7dOtauXcv555/PiBEjyMjIAGDKlCl8+eWXLFiwgOXLl5OVlcXll19e83qHw8Hw4cOprKzk559/Zs6cOcyePZsHH3ywcffKgBSzhfILukH/rj73S+gqprgYCv5RwvXBu+v1ul4WjRFX/UjV2Z2bqDIh6ke12ahK7UXVgM4+9UW3XoF1ySWXcNFFF5GQkECHDh3497//TWBgICtXrqSwsJC33nqLZ599lvPPP59evXoxa9Ysfv75Z1auXAnAt99+y5YtW5g7dy7du3dn2LBhPPbYY7z88stUVlY2yQ4ahaKpZPcz0+LZTA5O7oXWsb3cuNqIVD8/DlzWije7v1PvsfI0RWVq83T2XKOghTdrogqFOD3FbEHtlsTuB7oR9/gO8jtYQPWd40SD27AcDgfz58+ntLSUlJQU1q1bR1VVFampqTXrJCYmEhsbS3p6OgDp6el06dKFyMjImnWGDh1KUVFRzVnayVRUVFBUVFTr4Y0c/jozYxcx99bnOPyUieKr+8kBspE4enak39Ub6WNt2B93cy2A/547j6MXd5QOGMLlFJMJU3wcuRN6YX4hn+/GPs2UqCXomu+EFTQgsDZt2kRgYCBWq5Vbb72VTz/9lOTkZLKzs7FYLISGhtZaPzIykuzsbACys7NrhVX18uplpzJjxgxCQkJqHq1be28DuIpKd6uVn3q9Q9ojC9j2YAf0s7rJZcIzoDUPZ9c//JjR8ls0peH9jIbZimk2PhMtIb4RqxPi72mhIRRf1ovMZwN4etrrfNr+6zq3wXqbev/1duzYkQ0bNrBq1SomTpzI+PHj2bJlS1PUVmP69OkUFhbWPA4c8P75imyqhTFBeay47D8EPplF7vU9MbVqKZcJ60vVKDyvPXed/zXNtYAzeiuzovFquw/Ye3WEfIEQTU4xmVA7J7L77mSufWQha/vOYbC/44y+dBldvefDslgstG/fHoBevXqxZs0a/vvf/3L11VdTWVlJQUFBrbOsnJwcoqKiAIiKimL16tW13q+6F2H1OidjtVqx+ugBIsYUyIJ2i1l69w/cNnAMsW9HY/0xw/fma2ogU1wMZWMLuTFkD2A+7fqnE28O5NLLfmbd/3qifb/+zAsU4iS0Fi3IHdGeltfuZWn8m8SYAmmM31+jO+OodjqdVFRU0KtXL8xmM0uXLq1Ztn37djIzM0lJSQEgJSWFTZs2kZubW7POkiVLCA4OJjk5+UxL8VqaojLEVsWvA19n2PM/sO+enqjdkqQt5TQUq5WDI1vxZtd3ag1ue6bub5HO7mtVtBYtGu09hYDfev8N6c22Z1rz7P2v8kXCot/CSkA9A2v69OmsWLGCffv2sWnTJqZPn84PP/zAmDFjCAkJYcKECUydOpXvv/+edevWcf3115OSkkL//v0BGDJkCMnJyVx77bVs3LiRxYsX88ADD5CWluazZ1D1YVMtTGu2m68nPIX1hWMcu64vpqhIuUx4MoqCo38yfUdvpLulcSfWDlH9eW3gHPKGtZcvDaJxKApaUgL7p3bn/P/8yNbBrzHQz91FeZ56/SXn5uYybtw4Dh8+TEhICF27dmXx4sVccMEFADz33HOoqsqoUaOoqKhg6NChvPLKKzWv1zSNhQsXMnHiRFJSUggICGD8+PE8+uijjbtXXi7eHMhnCYv59n4zE88bS5vZMVjSt+IsK3N3aR5Da96cHdeYeb/lt2jKmbVdncy5/mUEjjuEtq49joztjf7+wndozcM5dmECAddl8VXCU8Sb5fLfqdQrsN56662/Xe7n58fLL7/Myy+/fMp14uLi+Prrr+uzWXEKQ2xVbDn/dZ7p1pl3Px5M/LxsHLv3g9Ph7tLcS9U4NqQdDw365Iw7WpyKVTEzq8M8UsdOo90TQTiLi5tkO8J7KVYrzt5JbL3OzGvnvcVg/wo0RS7//R3f7W7iJayKmfubb2fhDU9R8pJO4T/6oEVGuLsst9I6tkUZe4QxQYebdDuxpkAmXPwdx89OlMuyou4UBa1jezLv7kXvl9aTMexlhtiqfLr3X13JT8hLtDMHsrTzR9z/0DtsfTQOfUB3VD/fuwiu2mzsv6wFryW9V6fBbc9UWtgmssZVnrjlQIjT0MKbUXx1P8petPPVzU/xROSv9R55xZdJYHkRs6JxaUAZm4a/SPILm8mc0hNTfJzvfPtXFCr7J3HBqNV0ctH4aoGqH7P6zeLQZXE+NaabqB/FbEFP6cbWf7fnwcdnsbTTJ7+1VYn6kMDyFGrjfRSBqh/PRK3mo1v+w8HnbJRc0ReteXijvb+nMkVGsGe0ygMRy116eaWvVaftlTtROrV32TaFQSgKpphWHLqzN4kvbmHDxf/lQluFXP5roMbt7+vhTPFxVMZ45th8FYEm7CGNN82KpqgkWWys7vMOnyZFMH3ZlSS8E4O6dit6lfcNNKz6+ZEzvC3/GfRek3W0OBWzovFym085e9zddNwfhiM/36XbF55JCw2hZFBHCq4r5sPu/yHJYgP8G+/90TkeAfazOqHY9UZ734ZQ7E60X7Y3+YAGPhVY+f2iueyBJQSpnjdKhAOVLn4HGv16tlUxMzoon2GX/Jfpvc/jp7m9ifloH/asw6C795e8MTm7JtB8TCaXBuQDrr83KtoUyL3DvmDu8ovx/2KNV/1sRf0oJhNqh7bsuKEZD178EdcE5WBWbI2+nXiTxlP/mE32VaE4dfde9n9lx0BiJthAAqvxOCwKN4duJkRtvG85RhGi+vNiy5/ZMWUZI8+6leh3WuH/fYZX3LulhYWxbayNRe1mYm6Ce67qamzQPl68roSgdS2xHzzktjqE+5iiIjk8si3tx+7gf7FvE20KpKm+QNlUC5cGlAHu/xte0SKbApo+NOVCqg+pvkz4y4A3ufLpRex+sBtq10QUk4G/t6gaBUM78sTQD+lgdl9YwYkDyMxuc8kaKR0wfI0aFETlhX3Y/mw0b9z7Xz5su/S3sBKNSQLLB9lUC2mhB/juH08T9MoR8sb1MWynDFOb1lSMOcZlgbmnX9kF+luh97Ub0Xt0dHcpwhVUDS0pgd33d+bqZ78m49w36GWVLytNRQLLh8WaAnk/fgkvPvASW5+Kpyq1l6Hu3VJtNvZd3ZJZXRp3cNszoSkqT7Vcws5xNrTQEHeXI5qQ1qIFeTf0xfZ6Pj+O/Q+3hh7ymN9DbyWB5eM0RaW/n8amIS+R+uyP7L+7J1qHdoYY1LWqbyIXX/Wzy+65qqswzcbzQ9+l4EIZUd8bKWYLznN7sPWpWP57/8t81O47IlzcM9VXGbjxQjSmQNWP+5tv5+qb1nHbeaMpmNOX8K934Dia5+7STkprHs6OMRrzm//UJIPbnqlhtmL+c30O2sa2OLbudNl2VT8/FL/6z3ygBAehqq7r2eg0gxoWBg24xcJZetw9t2YoClq7NuwbHc34q5fweViGjFLhYhJYopZ25kAWJy3kh4dUrh90AwlvtT5x71ZFhbtL+52qkTesAy+dN8vl91zVlVnReK3je4y6/G7injngsgk39aR27L4yGEds/banKDC67VqXDGcVpCrEnbufHR1i6tf7X4eADf60/mC/y3thasHBFAxLxv+mLL7t8NRvc1RJWLmaBJY4qUH+TjYNeYkne/bh0w/OIe6TXBw7dnvE/UVaYjuCrztIqn8xnjwNQ5LFxpVXLeenn/qh/eCa2YnV3QdA7cSPA19sYC+1pg+sCC2ARYlfQWL9XvdVmR+PLboOe1Z20xR2EorZgrNvMluvtfBq6iwutFUA0vvPXaQNS5xSoOrH4xGb+OLWpyh5wUHR6H5o4e4dKUQNCGDP6HBeb/++IRq47w3/hT3X47LZiR1FRbSbm8/t+0fi0J0u2aYrFDqPM+n7a2m2aIdrps9RNUxt23BoSm/6vryOTRe/8FtYCXeSwBKn1c4cyIoun3Lvo3PZ+mQ79JRuKO6YIVpRKD8nmWtGLDfMwKE21cK8c97gyCWum53YuWUn+2clsKGy8Yb6crf/O9KP9nPtOPKONfm2tNAQSi/rTfGrCl9NeorHIzYRqBqn96w3k8ASdTYyoISMYS+T/FIGWZN6oSW0delI8KbYGA5dV8Fd4a65vNZY+lgVWly7Hy3JRYPjOh20+Go3V/7vVsqcxh83MqPyOF++dzbayi1Nuh3FbEHp3ZmtT3bkniff5fvOHxMrN/96FAksUS821cLz0Wv54Pb/cOQ5E2WX9XXJ/UaKyUTOBTHM6jvbcN92NUXlrXYfsufqZqi2xh9T7mQcuUeIn63w7LGuLtleU3HoTibtHE3sJ1lN1/FHUTDFtSb71t50mLmdjItf4tKAMhlR3QPJJyIaJMli49Oub5M/toSqzvFNv8FuHWl57V76Wt3f6aMhok2B3DrqGypTklyzQV3H/NNm3vvofHZXlbhmm03gs9JQKt6Mxr7vQJNtQ7FYyDu7FWePW8eTUT9JV3UPJoEl6q1KdzC7KILUWfcQ+0AV6s+bmnR7WnAwu0YH81r8Ry7pdt1Ubg7Zwd5xOqaoSJdsT6+oIH7BEW7cMYYq3QUdFRpZofM4dy8bTejirU3a0UKvqCD0g7Vsu6cTnZfdyroK419G9VYSWKJeNlRU0PnH65lz56XEz1iPY0sT99pSFMrO6chtwxYbfjDR6g4YOZe4ru3PsX0P5W9H8+1xz7xf7e88cSSFDrPLcRQUNvm2dLsd7fv1JN51gFsev5Mrdqdy1FHa5NsV9SOBJeok31HG9ZnnMGHGZNrfmYNl8VqX3AxrimnFkRvKSAvb3uTbcoX+fhptxu1E7eSiwXGdDsK+3cGk5WMpdB53zTYbQUblcb5+7yyUddtcul3HkSOEz1pN2Y2hnPXe3bxXHE6FXuXSGsSpSWCJv1WlO/i4JJi+8+4i++ZWtHh7DfbsHJdsWzFbODQyllk9Zxvinqu6eqXN5+weG4YaFOSS7TnyjpHwVhX/yh5oiHuzypyVXLn2JmI/OuieIZicDhzbd9Hu0Y28MeVyuv44gQ2eNNKLD5ORLsQp7a0qYcQvNxH2VhDtl/zisuGFqqkd4insUckXhT35RqnfZcehQZvo79f07V0H7SXMPJaCuZ71BXXKo6pXAi4bAWPNVpa/14eMO/9HV4tn97KcWZBIxNv+2PdtdmsdzrIyrN+spf36CG4YOZlu123m2VaLCdNc09NT/JUElviLEmc59x0+lx/f6UXrzzKxH9yG0w1DMilFpST+V2WtWr+u2Y5ACwvv7czPPd5v8k4aHxV35sd/9scvq36zvkboOmrBUVx1a69eVUmrLw9zY+q1/NR9vsd2XjlsL2HmwqG0/zEDj+gmouvYs3OIeDOP7PQEet0wlScu/IArA/Ok27sbSGCJGg7dydLjVm5ZcjMdXy8latNq7Hb3jZZgP3AQGtCb2dw8nOOVUY1f0ElU6Rr+B4px/lq/thYdcPXFOceeTCyzevP+45GMCz7q4q2fnkN3krbvMtrNL8RRVOTucmrR7Xb0jVtJ/Fcwr3x/Jc/ffIz3O802zIgr3kK+IgjgxKWtARuv4tH7biBp+nb0XzLQ3RhWogk4HQR/t42Hvx3lkT3gfqpQ2T+3PfrmHe4u5ZQcRUX4f76G8LQqRvz3Hu7I6iOdMlxIAsvHVehV3J/TlQtfvofw2yoJ+Hi1S7oRC/dwFBTS8e1iph260KM6YJQ4y7kh/Xqivjng+V+UdB373v20fHEt225PJvGbifxwXPWon6e3ksDyYSvKIWnpLaxN60HM8+uw78v0iOlDRNPSN+9g/ftd2FTpOWcGLxzrRut3NJfPc3Um9KpKlPSNJE3bzX0P3szw7Zdw2G7cUUWMQALLBx11lHLF7lSmPTyRxKn7UX7e6FkTNIompdvtxHySydhfbvCIy1mZ9hLe+ex8/H7casgvTI78fEI/WIs6KYBBc6bxyJFkj/i5eiMJLB/i0J08e6wtA+bcTdmt4YTOXY3jaJ67yxJuYD9wkGZzAvigONqtdVTpDibsvIa2HxzDWep57Wp1pdvtOLbsoO2Tv/LjpH4kLp7ICtfeBeITfKqXoKXEyfDNYzBrHtFhthZV0ekXvo9HIn5pki7HK8sdXL/uOqLfsNL2x404yurXDVt4n4BlW3l04RVcePUzRGjuGbpp6XEbBe/FEL59jVu239icpaWo//uFpG0tuPuiiUTdsJeX4j9qkmlKjjpK+efhVHYVtcCpu26an5PZtzuSJHY1+XZ8KrCClm1DSfd3dxknZ7Xw4aSzeGD02kYNrKOOUqYeHMbWNzrR9us92LNzXN6dWngmZ3ExCe8UMDnlYt5ts9Tl9xUddZQycdlNJH25C4end7SoJ8eRI4TNOYJ9ZXuGjr6HW676momhOxt1xJYsh0b6/B60/vgAVLp3wN4k5z4c+U3fWcunAstRUAge2gNO9fNDcbZqtPer0h3MLGjLi59fRLt5+TTbshq7K6YWF4bi3LyT7e/0Zc2939HfxQNgPJ/Xn3bzHDiOHHHthl3IsX0Xbf7vIF+kp/LiuPN4f8Dr9LU2XmhpFeDIzvWZNmhpw/JCqyuq6PS/6/l0ygW0e3wjzs3bmnZEdWFcTgdRX+xl7MoJlDhd1+jya2U5n71/DubVrh3c1h2c5eWYv11Lx8mZTJxxB2P3DSLXA++DMwIJLC9y1FHKhMyzmTjjDhIm52D+di1OaasSp2E/nE3rOSZeLejkku1V6FXcmHEtcQuyfOr303E0j+Zvribv1igGvHc3bxVGSW/CepLA8gJVuoP5xWH0/+Ausm6NpfkbK102orrwDn7LN/P2R0PZUdX03/w/LYlAmxvepLMIeyynA+ev22j/6EbevfsSuv44gZXlcvWjriSwDO7XynJ6rBrHzMlXkPBIBs4NWwx5L4twL2d5OfHzchiXMb5Jv/Uftpdw/6KrCfu2iSf+9HDOsjKsX62h/Z053Pb0JP6x9zy5TFgHElgGVeg8zrTsHox9YSpxdxRg/WYNzuJid5clDMyxax/qnOZ8VNI0Awc7dCeP5aSS8H4ZjrxjTbINo7Fn5xDxxhry06JJWXAXs4siKHO6t8efJ5PAMpjqCRV7fD6ZX2/tTMsX1xpqOBvhwZwOQhdv5YFlTTM47poKnZ/e64n6a9Pfr2Mkut2Oc8MWOj66lbfvG0mv9AlkVBpndmhXOqPAevLJJ1EUhcmTJ9c8V15eTlpaGuHh4QQGBjJq1Chycmq3p2RmZjJ8+HBsNhsRERFMmzbNrdNYGMWOqlIGbBjNM//6B4n/3AarN7lnRlbhtRwFhbT7wM7juec26vuWOMu55dexxHx20Kc6WtSHo6AQ/8/X0PbuAq5+9S7uzelOvkN+Vn/U4MBas2YNr732Gl271p5cb8qUKXz55ZcsWLCA5cuXk5WVxeWXX16z3OFwMHz4cCorK/n555+ZM2cOs2fP5sEHH2z4Xni5Emc5D+R2YeSb02g+2UHQh6tkRHXRZMwrt/Ddh335tbLxurk/ebQPwXOCse/3wY4W9aHr2PcfoPUrm1g/uQc9F93JV2V+VOm+2973Rw0KrJKSEsaMGcMbb7xBWFhYzfOFhYW89dZbPPvss5x//vn06tWLWbNm8fPPP7Ny5UoAvv32W7Zs2cLcuXPp3r07w4YN47HHHuPll1+m0s13a3uaExMqanRdehsrJ/ch7un1OHbukU4Vokk5y8uJ/fgwt2wd0ygdMH6tLOeTBecQtNSYg9u6g7O4GHXFBpIe2M/D/76ewZuvYHeVjATfoMBKS0tj+PDhpKam1np+3bp1VFVV1Xo+MTGR2NhY0tPTAUhPT6dLly5ERkbWrDN06FCKiorIyMg46fYqKiooKiqq9fB2h+0lXLhtBNMfu5mkew+h/bAeZ7mMpilcw7EnE2V2C94vPrPRV6rvuWrz6VGPm0XY4+k6jpxcwt9dQ+AUM5e8cWLCyEKn77Zv1Tuw5s+fz/r165kxY8ZflmVnZ2OxWAgNDa31fGRkJNnZ2TXr/DGsqpdXLzuZGTNmEBISUvNo3bp1fcs2jBK9iv/LS2DQu9Mw3+5P2Dur5Z4q4Xq/dcB47JvLz2iOp/eLW6HNDcexbXcjFudbqkeCb/PCZrbdkUyPzyfzXnE45XrjD5Lt6eo1luCBAwe48847WbJkCX5+rht4bPr06UydOrXm/0VFRV4ZWpZ8hZQfJtFmtkrblb/iMPB0C8L4HAWFtPuonAfPGsrMmP/Ve3DcXEcpj353GYnf7sDhw/dcNRZHURHKzxtJ3N6MV5ZfSe7IcsJKfOsSa71+A9etW0dubi49e/bEZDJhMplYvnw5L7zwAiaTicjISCorKykoKKj1upycHKKiTtzbERUV9Zdeg9X/r17nz6xWK8HBwbUe3kZ3OGn1QymJ9+VgWrrO0HMDCe+hrt3K6ve7saai/gfG6YeG0v79CrnnqpE58o4R+NEaEh4qotnmInD4zpeBegXW4MGD2bRpExs2bKh59O7dmzFjxtT822w2s3Tp0prXbN++nczMTFJSUgBISUlh06ZN5Obm1qyzZMkSgoODSU5ObqTdMh69qhJl5a/YD2W5uxQhaugVFcR8dpAbN4yr1w2t6yoqWf1xV7R13j+4rVs4HTh27kFfvwXdh24JqtclwaCgIDp37lzruYCAAMLDw2uenzBhAlOnTqVZs2YEBwdz++23k5KSQv/+/QEYMmQIycnJXHvttTz11FNkZ2fzwAMPkJaWhtVqbaTdMijpQSU8kH3/AcLf7sPrHTswOWzfadev0Ku4fuN4Yj86hF06CjUtHztmNPpIF8899xwXX3wxo0aNYuDAgURFRfHJJ5/ULNc0jYULF6JpGikpKYwdO5Zx48bx6KOPNnYpQojGoOsErNjGy18NI7MOHTDeLGxL8LvB2PdluqA44UvOeALHH374odb//fz8ePnll3n55ZdP+Zq4uDi+/vrrM920EMJFHEVFtJ9XyB39RrGg/dennBX7sL2E574ZTodl23H42Ld/0fRkLEEhRJ3om3ew/4N2/FR+8hlzHbqTKQcupf38Uhz5+S6uTvgCCSwhRJ3odjstvznELWvHnnR24nWVDjI+S0TZtNMN1QlfIIElhKgz+75Moub6/WV24jJnJePX3kDsR4dkRBbRZCSwhBB1p+vYlmXw5udDao1t905RPM3n2aSjhWhSElhCiHpxlpbS9sMCrt92LYXO42TaS/jPV5cS9L9dPtfNWriWBJYQot70bXuomh3JzPyuTN53Ge3nF+E4mufusoSXO+Nu7UJ4osoKE3vt5Zhp2m/82RUh4GzSTXgkvaKCsCU7mdPuAvyP6LTYvN7dJQkfIIElvI5eUkrMLDPXfX5Xk29Lq9QJztvvi5mF42ge8e8dgvIK7BUV7i5H+AAJLOF1nOXlmL9di1lRXLI9uw+329j37nd3CcKHSGAJ7+XDQSKEN5JOF0IIIQxBAksIIYQhSGAJIYQwBAksIYQQhiCBJYQQwhAksIQQQhiCBJYQQghDkMASQghhCBJYQgghDEECSwghhCFIYAkhhDAECSwhhBCGIIElhBDCECSwhBBCGIIElhBCCEOQwBJCCGEIElhCCCEMQQJLCCGEIUhgCSGEMAQJLCGEEIYggSWEEMIQJLCEEEIYggSWEEIIQ5DAEkIIYQgSWEIIIQxBAksIIYQhSGAJIYQwBAksIYQQhiCBJYQQwhAksIQQQhiCBJYQQghDqFdgPfzwwyiKUuuRmJhYs7y8vJy0tDTCw8MJDAxk1KhR5OTk1HqPzMxMhg8fjs1mIyIigmnTpmG32xtnb4QQQngtU31f0KlTJ7777rvf38D0+1tMmTKFr776igULFhASEsKkSZO4/PLL+emnnwBwOBwMHz6cqKgofv75Zw4fPsy4ceMwm8088cQTjbA7QgghvFW9A8tkMhEVFfWX5wsLC3nrrbeYN28e559/PgCzZs0iKSmJlStX0r9/f7799lu2bNnCd999R2RkJN27d+exxx7j3nvv5eGHH8ZisZz5HgkhhPBK9W7D2rlzJy1btqRt27aMGTOGzMxMANatW0dVVRWpqak16yYmJhIbG0t6ejoA6enpdOnShcjIyJp1hg4dSlFRERkZGafcZkVFBUVFRbUeQgghfEu9Aqtfv37Mnj2bRYsW8eqrr7J3717OOecciouLyc7OxmKxEBoaWus1kZGRZGdnA5CdnV0rrKqXVy87lRkzZhASElLzaN26dX3KFkII4QXqdUlw2LBhNf/u2rUr/fr1Iy4ujg8//BB/f/9GL67a9OnTmTp1as3/i4qKJLSEEMLHnFG39tDQUDp06MCuXbuIioqisrKSgoKCWuvk5OTUtHlFRUX9pddg9f9P1i5WzWq1EhwcXOshhBDCt5xRYJWUlLB7926io6Pp1asXZrOZpUuX1izfvn07mZmZpKSkAJCSksKmTZvIzc2tWWfJkiUEBweTnJx8JqUIIYTwcvW6JHj33XdzySWXEBcXR1ZWFg899BCapnHNNdcQEhLChAkTmDp1Ks2aNSM4OJjbb7+dlJQU+vfvD8CQIUNITk7m2muv5amnniI7O5sHHniAtLQ0rFZrk+ygEEII71CvwDp48CDXXHMNeXl5tGjRgrPPPpuVK1fSokULAJ577jlUVWXUqFFUVFQwdOhQXnnllZrXa5rGwoULmThxIikpKQQEBDB+/HgeffTRxt0rIYQQXkfRdV13dxH1VVRUREhICIMYgUkxu7scIYQQ9WTXq/iBzyksLKxzvwQZS1AIIYQhSGAJIYQwBAksIYQQhiCBJYQQwhAksIQQQhiCBJYQQghDkMASQghhCPWeD0v8TvXzQ2nCQX+bil5ZibO01N1l+DTFakW12dxdRg1HYRE4He4uQ/yBGhSEs7RMPpc/kMBqKFWjeHg3/G/LItK/2N3V1FmlU2PX3K5EvLkG3W53dzk+SQ0IIGdcV2JH7yHQXOHuctiR3wLrWx2wfbZWDo6eQtWo6p2A5UA+jl173V2Nx5DAaiDVz8rhAQqrOsynuRbg7nLqzKE76ZIajfpBEI78fHeX43Oqw+q+KfO4KrDQ3eUAUBFXxaR7B7FV7U3Q15twlpW5uySfpwbY2HWphZDtUUTsP4ReVenukjyCtGE1kBoSTEj7fMJUY10S1BSVlJh9EB7q7lJ8jhoQQM61Xbln8nyPCSsAq2LmpZgf6HTPrxwb1Q3Vz8/dJfk8R9d2jDt/BfYLC9Aimru7HI8hgdVAVW2juLzNRjTFeD/C8RE/UtQtwt1l+BTVZiN3bFfumvIho4M878zWqph5odUK+t25lvxR3SW03EixWjl4vo1xoauZ1HE5JT1iQFHcXZZHMN7R1hMoCsVxflwctNHdlTRIsrmU4lgNVM3dpfiEE2HVjalTP2RMUJ67yzklq2Lm6eif6T35FwktN9JaRhF97kFiTTauDtrNwfNVtBCZtBYksBpE9ffnSE+IMRmz00KY6k/lWcXyR+AC1WF159QFHh1W1ayKmWda/iih5S6KQn6/aO5t8w2aohKi+pPSfxt6TLS7K/MIElgNoAYHEdLxmOHar6ppispZrfeihEpgNSXVZuPImBNhNS74qLvLqTMJLfdRAwM5fIGds/x+73k8Jfpbsgc2QzFJHzkJrAaoahfNxbEZhmy/qnZti58p7Bnl7jK8VvWZ1R13GSusqklouYejazuu6bmaQPX3n3dXi4ZjSD5aVKQbK/MMxj3iustv7VeXBa93dyVnpLOl+EQ7ljTmNro/XgY0YlhVk9ByLcVq5eB5Nm5q9nOt582Kxg0J6ZR1ivb5v1cJrHpS/f052l0xbPtVtTDVn/L+JWihoe4uxauoNhu51xo/rKpJaLmOFhVB83MOE2v66wgo14ds5cAFJp//e5XAqiclKJDQ5DzDtl9V0xSVs+P2oIQEubsUr1EdVrdP+dgrwqqaVTHzVPT/6HLHJgpHSmg1CUUhv38r7m33zUmbGkJUf/qftQ17YqxPn2VJYNWTvV00F7TaZuj2q2pXN19FUQ/pfdQY/hhW1wXnurucRmdTLfy31fckTd5MwWUSWo1NDQzk8GAH5/oVnHKdm6N+4NCgAFSr1XWFeRjjH3VdSVEoivfnqtA17q6kUXS3FlAYL+1YZ+rEcEvduHPKAq8Mq2o21cKLMcvoPHkTBZdLaDUmR+e2XNF7ba3OFn+WYnXQ+oL9KPGtXViZZ5HAqgfV35+8LsZvv6oWrvpzvHeZz18XPxPVYwNOnuwdbVanY1MtvBTzA10nb+TYVT3QguXWiDOlmC1kDQzgtvD//e16ZkXj9tilFHQN99kvmRJY9aAEBRLe9Yjh26+qaYrKoHY7pR2rgdSAALLHd+OuyR/6RFhVqx7GafCUnzg8trOE1hnSIlsQcl72STtb/NlZ1mMcHmJHa+6b4wtKYNWDvW0050fv8Ir2q2pXhq+muLvcj1VfakAA2dd149473zfECBaNzaqYeSTiF668dal3hZaiYIqKRA1y3Ze4/AGtub/913U6roRpNm7tu5yq5BgXVAYoClpwMKZozzhGyK3T9VDUzp/RYasB77l2391aQH57E/6KArru7nIMoSas7njfIweydRWzojEtfAvaRCfzlQuIfnczjqIid5d1ZnQdR+sIDg4Owj9HJ3xTMeq+wzgLCptk/jgtOJjDg5yc7ZcP1O3KzaVBG3nnrAtovdKKXtEE86kpClqzMJzxLTnSK4jKIIWY7wrgcHbjb6ueJLDqSLXZONZJoaXmXRPchav+HO9xHC0kGEeB50x54amqLwP6elhVMysaU5ttw3GLygJ9MNFzjR9aSsZuHMO789g/32ZFcUe+z06g5PuOtFhfgd+mAziPFTTa/FSO5DZc3nctIfVoZmhvthI1+CDqZ7E4tu488yIUBdXfHzUkmONdYshLtqCed4zRbVfR27aHO968BT1j15lvpxFIYNWREhRIs27e035VTVNUzmu/g6zQcJDA+ls1YXWnhNUfVZ9pcSsswPih5SwrI35eDjPPG8TH7b+BiA0UdS5ncVkrZu4/lyM/tCc8w07glqPoWTknJrxswNUJxWrl0MBAXmm+Agis8+vMisa0Not5rMv1BG7b1bArI4qCarOhtIqisFtzjvRUad9/PxNafsrFAQcJVKzYcXD2L2No80EWdg+ZQFICq47s7aI5N9qY81+dzqjwtczoMh6/fZnuLsVjVfcGvPvODySsTuKPofWRfj5R72UYOrQcO/dwaE5/1jyg099PI0yzMToonys7fUJRUjnLjkfxQW4f1mxNJup7jdBfC2DvgXqFlxbRguDzs4k11f9LcD9rPlmDHSQtDcORd6xOr1FMJlSbDWe71uR3CSb3LAfndNvGw1Hz6WzWsamW39Y8Uc9nJaFY5oRh37Oq3vU1FQmsOspPtHFV6GrActp1jaabJY9jiSZafunuSjyTGhREztjO3DXZs+ezcrfq0HJOVPiE84icuxlncfHpX+iJdJ0WX+zgH+fexObBM2sO5pqiEqbZGBVYxKjApVTEL2LPBVW8evRcFm7uSsR3ZoIOVGDedui0lw6PnteaJxJex6zUf166MM3GTSkrWNZlANoPpw4sxWxBbRZKZVIMR7v4UTGwmEvbbyIt/EciNStWxQyY/1qbo5Rpi24hcck2PKkRRAKrDtSAAPKTIM5UhTcGVoRmo6zrcbTQEGnH+hM1IICsCV2YPlEuA9aFWdG4NzwDJsIn+nlEvmfc0HIczaPd27H8X/cePNIi46TrWBUzSRYzL7Rcw9PRP7NnUBXLyjoyd39fCn9uT9SqSvwzsnDkHkW3V9WcfWlhYRw5v5Kz/co5WWDUxWXBv/B+38G0+snyezD+1h6ltG5JWdswcnubCU7J5Yb4RVwSsIPmmv9vAfn3lyAfzRlEh3fLPO54IIFVB2pwEB167/e69qtqmqIyKGEn2aHNpB3rTxzdExh742IJq3o4caa1ibJbLHxbNYAW7/2Cs7zc3WU1iGnVVj764FzG3LKaDuaAv123OrySLAe4NWQ/uZ3KWD62Nc/svIDi9XFErnYQsKcAfd9BHAkx3NDzp9/OcBqmvdlKwKBc1M9j0Q8ehnatKUwKJXuATudu+7m/1ed0seQTodl+a8qoWztZRuVxvv+gD61+Xdfg2pqKBFYdVLWN4qKIH7yy/araiPD1PJc4Bou0Y9XiNKv09N/n7jIMx6qYeShiHQsv6YTyRRAYNLCc5eXEfZTD+HPG8WPXBXU+BmiKSrQpkNFB+Yzq8T5Hux5n4z/CWXC0D8t+6QwKzA5dS306W/yZWdF4uMOX3H7DDdibh3J+p208EvkuCaYqglW/eoVUtQq9imt/vY7Y9/djb4ou82dIAqsO8jr5c1FgBvD337CMrKc1l6OdzbRc5O5KPIxvjoDTKKyKmVD/clCN/UN07NyD/l4/liZYGWKrqvfrzYpGtCmQaFMFF8b+SGHMEgqdDqJNDQ+raoP9y/h29NNEapbf2tlO3iZVV5+WROD/bij2g9vPuLam4L2nDI1Etdko7KATozX8l8AIIjV/SpIrvWfEgsYi91ILXafZ4t3csuw6Cp3Hz/jtQlR/YhshrOBEGMabA//Qw6/hch2l/PPrqwlZsq0RKmsaElinoYYE077nAayKd5+MmhWN85O3oYSHubsUITyO48gREmZX8vTRvu4upUk4dCf/OpxKwrxSj+to8UcSWKdhj4vgwsgMr26/qnZJ+AaOJ7RwdxlCeCR13TY+e/8ctlaWubuURve/chOr3u0Bv+5wdyl/y/uPwmfoaLcAhgZscXcZLtHbms2Rbt7XbV+IxqBXVBD7aQ7XbxmHQ3e6u5xGU+IsZ8LP42n10Z6mGZuwEUlg/Q3Vz4/C9hBn8u7LgdUiNX9KOko7lhCn4ti5B3VOcz4u9Z5L5y8c60bcuxp2Dxjc9nQksP6GGhZKOx9ov6pmVjSGdM2AiHB3lyKEZ9J1Qr/dzn1LRlPiNGZX/T/aXVXCux8PxvqjMa4iSWD9DXvrFlwctckn2q+qDQvbRHlbCawaxu6RLZqAIz+fhDll3J99jrtLOSNVuoNbd11D/IKjJ8ZANIB6H4kPHTrE2LFjCQ8Px9/fny5durB27dqa5bqu8+CDDxIdHY2/vz+pqans3Fl7CPxjx44xZswYgoODCQ0NZcKECZSUlJz53jSyoz0COT/Ac7t4NoWz/HLI7mvx2Sm4hagLZeMOlr/bh9UV9b8vy1N8VRZC8axWOLbtdncpdVavwMrPz2fAgAGYzWa++eYbtmzZwjPPPENY2O/Xc5966ileeOEFZs6cyapVqwgICGDo0KGU/+FO9zFjxpCRkcGSJUtYuHAhK1as4Oabb268vWoMqkZxPMSb6j8wpZGFq/4cb1eBGtg494kI4Y30igpafXWYW3691pAdMEqc5UxNv4pmi3eD05OGt/179Wqc+b//+z9at27NrFmzap6Lj4+v+beu6zz//PM88MADjBgxAoB33nmHyMhIPvvsM0aPHs3WrVtZtGgRa9asoXfv3gC8+OKLXHTRRfznP/+hZcuWjbFfZ8wU2YKoHtk+035VTVNULunyK7uiWoNBBy1tVHLjsDgFx+59BLzbl/c6RjAu+Ki7y6mXJ470pc07Co6jxqq7XmdYX3zxBb179+bKK68kIiKCHj168MYbb9Qs37t3L9nZ2aSmptY8FxISQr9+/UhPTwcgPT2d0NDQmrACSE1NRVVVVq06+bwrFRUVFBUV1Xo0NWd4KFfGrPep9qtqF4X+SkWs9/SCOiNyZVSciq4T/N02Hv3qCvIdxmgDAthaWcZnn5yN9ceMhk3+6Eb1Ohrv2bOHV199lYSEBBYvXszEiRO54447mDNnDgDZ2Se6RUZGRtZ6XWRkZM2y7OxsIiIiai03mUw0a9asZp0/mzFjBiEhITWP1q1b16fsBjnSL4yzbI0w/bQB9bQeI7ufVdqxQM6wxN9yFBTSfn4p/8w+392l1EmFXsUt28fQ5pM8Q46gX6/Acjqd9OzZkyeeeIIePXpw8803c9NNNzFz5symqg+A6dOnU1hYWPM4cOBAk24PVaMwAZK8e/jAUwpT/ShLkHYsQM6wxGmp2/fz/Zc92VvleR3H/uyXCpWihdHoe5v4GNpE6hVY0dHRJCcn13ouKSmJzMwTU1JERUUBkJOTU2udnJycmmVRUVHk5ubWWm632zl27FjNOn9mtVoJDg6u9WhKpsgWRHbP8bn2q2pmReOSLr+iREecfmUhfJWiYGrVkuwxnRh48S+0NFndXdFpdTRXEDkyk7yrumGKijz9CzxMvQJrwIABbN9ee9j5HTt2EBcXB5zogBEVFcXSpUtrlhcVFbFq1SpSUlIASElJoaCggHXrfp8cbNmyZTidTvr169fgHWlMzvBQrm69zifbr6pdFPorFTGh7i5DCI+k+vnhHNidrTOiePru13ml1ZlNxugqYZqNLxI/Zer0+Wy9vw3074pi9fygrVavI/KUKVNYuXIlTzzxBLt27WLevHm8/vrrpKWlAaAoCpMnT+bxxx/niy++YNOmTYwbN46WLVsycuRI4MQZ2YUXXshNN93E6tWr+emnn5g0aRKjR4/2mB6CR/v6bvtVtZ7WY2T3t4LqW936/0LasMQfKQqm1jFkTexJwn+2su68lxjs7zDUl1urYmZ0UD4rRj5D0NNZ5F7fE1Mrzzj2nk69rnn16dOHTz/9lOnTp/Poo48SHx/P888/z5gxY2rWueeeeygtLeXmm2+moKCAs88+m0WLFuHn51ezznvvvcekSZMYPHgwqqoyatQoXnjhhcbbqzOhahR0hI5m491b0ZjCVD+Od6hA9ffDWVrq7nKEcDvVZuP4oE7svbaSd/q9SB+rgqbY3F1Wg8WYAnmv7Td8d/ePTOo9lvZzIjGt2erRnTEUXTdYv0ZOXGYMCQlhECMwNfJpuCk6iry3Avmp24eG+tbUFCYd6seesa1xbN/l7lLcxjGoJ9PfmsNgf+PcXOlJBm66jMBxpThyck+/sqdSFLT28ey7Oop/XLWM25v9Qojq7+6qGlWuo5Qncs9l6ft9if3oIPZ9mU2+TbtexQ98TmFhYZ37Jfj2EfkknOGhXNtmlc+HFcClYeupaBXi7jKEcBs1KIjyi/uw/0l/Ftz4DA803+Z1YQUQoQXwTNRqZt72EjufDMU+uBdqQIC7y/oLOSr/kaKQ1yuMPv573F2JR+huLThxP5avt2MJ36MoaEkJ7J3WmREzvmNVv7foZPG+oPojTVEZ4Key5uzX6PefNeyf2g0toa1H3Y/pm/22T0VRKUiCjmY74Pk9fppamOpHWftKFE1DN9B4Y43Kc/5WhYuoAQGUn5NM1vUVzOv7X3pZLYDfaV/nLUJUfx6L2MCg8Vu5JX48beeGYVm51SNGdJczrD/QWoRj7ViITZFZd+HE/VgpibtR41q5uxQhmt5vbVUHbu/GgCdXsnrAa7+Fle/RFJUhtipWXfBf2jyxnYNp3THFx7m7LAmsP9JbNGNcwirMilwCqzYmMp2q6FB3lyFEk1KDgjh+aR8OPO3PnFue5/GITV7ZVlVfEVoAM1sv5+WJr7Dj32FUDO/j1rYtCaxqikJhp1Cfv//qz3pajpLb099327EM14dW1IuioHZOZPf9nRn57+9Y1We2z55VnYpZ0RjoB6vPeYV+j68h845uaEkJKGbX/5ykDauaonKkh0KyuRww7r0Vja2ZZqWoewXRZhN6hQ+2Y0kbltfSQkMoGJqE6YYcvkt6mlhTICBhdSphmo0nItYz6uY13JoyFv93exCybCeOvGMuq0EC6zdai3D8EwsIVn2ncbUurIqZlA57KGwVjX3PPneX43pyhuV1FLMFunVgx+hApg//jHHBhzArMtBzXWiKSl+ryo893+X1dh14qf+FtJ9XDL/uRK+qbPLtyyXB3+jRzbmuvdx/dTJjItOpiGvm7jLcQ86wvIoaFETF+V3ZMcnKE5e+T3//veQ4jlPmbPqDrbcocZZToldxQcBWpl/0Gdtv86dqYBdUW9NfmZIzLDjRftUx+Lf2Kx9tq/kbPS1HOdrFj8jlmqGm0xbiL5xO/A8UETc/hOcXXUNFiEJJDFTFVRDerITooCIGNNtND/99tDIVEaQ6CFJUAlWrz3TGqtIdlDgrKNadFDhNHLCHsrq0HRsLYsgsDCU/LxBzjoXgnWAtctKmyIF1/zGcjqY/NkhgAYqmkdtboaO5Amm/+qvmmj9FvcqJtphxlktgCeNylpbClh1YtpxorQo2mYiwWlEsFtBU7M1CWdLqHL5oMZjKAIXjkQqlcXaCo4uJCiqmY0gOZwfvINGSQ5BiJ0hVCFIthhip/Y/KnJWU6VUUO3WOOS1sKI9lY2lrthdGklUUTElOILb9JgIP6ZjLnPgfqcScVUjUsSNEObLRK6tq3ZflqqOCBBaghoVhS5D2q1MxKxr92++lMDoS59797i7HtaQNy6vpdju63Q7VAzwfzUPbAdUtWqqf34kws1pRAvzZFp3M+pY9qQxSqQhRKI11YokroXVYATEBBfQM3k+iNYs2pkKiNQs21b2dOCr0KvbbK8l2BLC2rC0bimM4VBrKwbxQHAdtBO5XsRboWAsd2LKOYzp8jNYVR6DiIHplZa2BcD3hq6oEFuCMjeC69unSfvU3xkX+zFPtxmL2tcASPs1ZXg7VB+0joOzLJAAIAFA11AAbip8fisVMVvPW7G6VxPFwEyUxCtPGf8R1we4d9HdWYRteeWsE/rk6ATlVWLNK8DuaTzvHYaiowFFSWnOZXwfsbq329CSwgLLWAdJ+dRrdLHkUxVsIVxQw3gD/QjQ+pwNncTEUF5/4/6EsrBvBCoSc1Q3b9RVuLQ/gqD2IVt/l49y4FQDnbw+j8vlTCsVkIruvRluz584B4wmaaVbyUqpQDTQ7aWNQq5xsKY9xdxnCYHRVIUIrdncZNDcVe9TgtWfK5wNLDQpCSyghXIZh+VtWxUyntodQw32re7tpWyavzhvOxyV1m69H/C7XUcqBrGZg9/QLTY2vKsSMn1Ll7jII0sqxh3jPl0yfDyw9LppRCRuk/aoObmz1P8oTo91dhks58o4RP3sfD78xlvnFYe4uxzByHaUMXT+Bjq9UunQkBI+gKJRGmghS3X9vV6hWSnm4xWvOsnz+KF0aF8ilIevdXYYh9LZmU9jGe37568p+KIvWb2TwxOvX8GGJTGh5OrmOUi785QYiHzejr93s7nJcTtE0yqIUglT396uL0ooojVTBS76Qe8deNJBitpDTV6ONyf3fhIwgUvMnb0AVqr/vXT51FBQS81YGj70xRi4P/o1cRynDNlxPi8ctJ8LKFzvoKCoVYTo2D/hiF6pWUt5MQVHdX0tj8OnAUgP80ToWS/tVHZkV7UQ7VohvHrAdBYW0fnMrD78xVkLrJI46Shm+8Xqa/9sP1vhoWAGoCo4QO7YzvJm40HmcCv3M2sFCVAV7kPd8Dj4dWM74GEa2+1Xar+rhxlb/43gn353Q0ZGffyK03hzLZ6UyYGq1o45Shm28nrB/+8OqTb4bVr9R/BxnNJTT3qoSUlbexDW7LyLf0fCZfm2KGYefLpcEvUFZXAAjQ9e5uwxD6W3NpjjW99qx/siRn0/rN7byrzfGSWhxIqwu+vU6wv7tj7LyV58PKzUoEFtgBWoDR07+qdxJ6sd30+ZfFRQ+HMuwTeM4bC9p0HtpioIz0I4a6L5JFxuTzwaWYrWSLe1X9Rap+XP07CqXjMzsyf4YWr58ebA6rEIft0lY/UYJsBEdUtSg164ohxvmpZH41F4cW3di+mEDYdPNDPxxEjuqSuv9fioqAc2Oo9i8o9nDZwNLtdnwSy6Q9qt6MisaXdodRA0OcncpblcdWr7apiVhdXJ6gD9tg47Wq6nBoTv5uCSYm9+9jfYv7cGenXNigdOBc+NWEh4vY+i3k/mpvH7jVKgotA3PQ7d5xzipPhtYzvYxjIrfKO1XDXBDyx853kVGf4A/tGm9PtanurwftJcwcNUtElYn4Qj2o6Mtp+7r607+Ly+JJ54ZQ9tnNv8eVn9cZ8sOkv7vKDfMS6v371m4tRTd3ztmUvbZo3VpaxvXhK5xdxmG1NuaTXGM2afbsf7IkZ9P7NzdPDJ7DN+WGWuaiYY46ihlyOpbiXvMIWF1Ek6rRke/rDqt69CdvFIQzycvnk/EvM04ik59KdGxay/tX9rDI7PG8GFJCA799GdbmqKSYMtFN3vHOKk+GViK1UpuL5VIzSd3/4xFav4cTbH7fDvWH9mzc2gzZx+TZ93EojLvGQrnz3IdpQxcdQtxjzlw/rpNwuokdFWhRR3GESxzVnJnVgrvzRhG83fWnxhI9zTs2TnEvbGdJ14aw8zCuDqFVqS5EN1LjnXesRf1pPr7EdL9KIGK9x5YmpJZ0ejS4QBqkPSQ+yP7oSziZ+/jrrcm8EWp94V5rqOUc1eeOLNybtwqYXUyikJFmBnLacZEL3NW8o/dl7Dln10ImbcGvaLuI7s7juYR/dZG5vzfxdyX0+u092qFamVUhlm94oqITwaWs10Mw2K2SPvVGRgbnc7xztKO9Wf2Q1nEvbGdf752nVeF1mF7CeeuvJU2j1XVTFUhTkJROR6uYlNPPeBvofM4V+wcSfEjMZi/+6VmPqr6cJaW0mzeOn76v37ceuB8ypyn7u0cqpVS3sw7ZpLyySN2cdtAxoSudncZhnaW3yGK2nhHQ25jcxzNo/WbGfzzteu84j6tw/YSzls58URY/brN3eV4NEXTKG2lEHqKI+tBewm9VkzEeXcYph82NCisqulVlQR/vJbdTyYxds9wCp3HT7pelFZCWZR3jCdo/D2oJ8Vq5Wg3ab86U9Gajby+dtQA77ghsbGdGMYpgwdnGvvm4uqwin9UwqpOVIWKZg6sJwmHg/YSzll6Jx0eKUFfl3FGYVVNt9uxLVxP4SOxXLhpLEcdf71XK0h1UBGme8V4gj531FZtNsJ6HpH2qzOkKSrdOmZKO9bfcBQUEvPGZsOG1mF7CYPSbzsRVpslrOpCURT0IDt+Su1LcOsqKhn45V0kPZGPY/uuRt2mbrdj+n4DQY8GMnDVLRz806gYQYpKVZB3DM9k/D2oJ2fblgxttVXarxrBP6JXcrxra3eX4dEcRUWGDK3qsGr7WKWEVT2ZrA5M/N6NfEU5jJ0zmcSH9+DYuadpNup0oKzaTJuHqjhn6Z1kVP5+edCmmnH6O0HOsIynqH2QtF81krP8DpHfXtqxTqc6tP712jhDdMSovgwoYVV/SmAANtuJHn8O3cnS4xq3vHMb8S/vxHHkSNNu3OnAkbGd5MfyuGRh7VExFJsdNdj4o7H4VGApVit5nRVaat5xE527RWs2CnpUoQbJME2n4ygqOtER4/Xr+KrMc4fJyXWUnmizelwuAzZIWAjxYcdwovNsfgJ3/fcW4l/c1vRh9Qf2PftI+s9hJsz9fVQM/4AKFJPxj3s+FViqzUbzPjn4K3JW0Bg0RSUleZeMK1hHjoJCWr+RwX2v3+CRoVXddT3+celg0VDOQH9iA47xxNEuLHh6CNGvr8eRd8zlddj3ZdL2hR3MeGEMrxe0Jz78GHqQ8TtI+VRg6XHRnBe1U9qvGtHIFus5nhTt7jIM40RobeXeN2/gh+Oe83uY6yjl3J/SaPO4XcLqDDgCzHy3ryOLnhxIs/nrcZaXu6+Wo3lEz81g7v9dxJbMaK8YT9A77iaro/JIGwu292Dp4Q7uLsVrlFZYCAoz4f0j6DUeR34+sa9tZZJ+K0GD6j5IalPK3h9O0gsFOLbscHcphqZVOIh5VkNdsx69yv1TFzkKCgn7YD1BBzqh2Os30rsnUnTdeOOrFBUVERISwiBGYKrHNNSK1YpiMf63DI9TVeXWb5JG5VG/jw4HzrKGz2wrTlADAnAeL2+Ue6wak2IyoVgsHvUZ2/UqfuBzCgsLCa5jhxCfOsPSKyrqNWaXEE1Jfh+9j7O0/pMsuoJut6PbTz1clFHU6yJ6mzZtUBTlL4+0tDQAysvLSUtLIzw8nMDAQEaNGkVOTu1LHpmZmQwfPhybzUZERATTpk3D7gU/SCGEEE2rXoG1Zs0aDh8+XPNYsmQJAFdeeSUAU6ZM4csvv2TBggUsX76crKwsLr/88prXOxwOhg8fTmVlJT///DNz5sxh9uzZPPjgg424S0IIIbzRGbVhTZ48mYULF7Jz506Kiopo0aIF8+bN44orrgBg27ZtJCUlkZ6eTv/+/fnmm2+4+OKLycrKIjIyEoCZM2dy7733cuTIESx1vJ7f0DYsIYQQnqEhbVgN7ldbWVnJ3LlzueGGG1AUhXXr1lFVVUVqamrNOomJicTGxpKeng5Aeno6Xbp0qQkrgKFDh1JUVERGRsYpt1VRUUFRUVGthxBCCN/S4MD67LPPKCgo4LrrrgMgOzsbi8VCaGhorfUiIyPJzs6uWeePYVW9vHrZqcyYMYOQkJCaR+vWMn6dEEL4mgYH1ltvvcWwYcNo2bJlY9ZzUtOnT6ewsLDmceDAgSbfphBCCM/SoG7t+/fv57vvvuOTTz6peS4qKorKykoKCgpqnWXl5OQQFRVVs87q1bUHnq3uRVi9zslYrVasVpkORAghfFmDzrBmzZpFREQEw4cPr3muV69emM1mli5dWvPc9u3byczMJCUlBYCUlBQ2bdpEbm5uzTpLliwhODiY5OTkhu6DEEIIH1DvMyyn08msWbMYP348JtPvLw8JCWHChAlMnTqVZs2aERwczO23305KSgr9+/cHYMiQISQnJ3Pttdfy1FNPkZ2dzQMPPEBaWpqcQQkhhPhb9Q6s7777jszMTG644Ya/LHvuuedQVZVRo0ZRUVHB0KFDeeWVV2qWa5rGwoULmThxIikpKQQEBDB+/HgeffTRM9sLIYQQXs+nxhIUQgjhGVx6H5YQQgjhShJYQgghDMGQo7VXX8W0UwWGu6AphBDCThXw+/G8LgwZWHl5eQD8yNdurkQIIcSZKC4uJiQkpE7rGjKwmjVrBpyYqqSuO+rJioqKaN26NQcOHKhz46Onkn3xXN60P7Ivnquu+6PrOsXFxfUaLcmQgaWqJ5reQkJCvOIDrhYcHOw1+yP74rm8aX9kXzxXXfanvicc0ulCCCGEIUhgCSGEMARDBpbVauWhhx7ymuGcvGl/ZF88lzftj+yL52rK/THkSBdCCCF8jyHPsIQQQvgeCSwhhBCGIIElhBDCECSwhBBCGIIElhBCCEMwZGC9/PLLtGnTBj8/P/r168fq1avdXdJfrFixgksuuYSWLVuiKAqfffZZreW6rvPggw8SHR2Nv78/qamp7Ny5s9Y6x44dY8yYMQQHBxMaGsqECRMoKSlx4V6cMGPGDPr06UNQUBARERGMHDmS7du311qnvLyctLQ0wsPDCQwMZNSoUeTk5NRaJzMzk+HDh2Oz2YiIiGDatGnY7XZX7gqvvvoqXbt2rbkLPyUlhW+++cZw+3EyTz75JIqiMHny5JrnjLQ/Dz/8MIqi1HokJibWLDfSvgAcOnSIsWPHEh4ejr+/P126dGHt2rU1y410DGjTps1fPhtFUUhLSwNc+NnoBjN//nzdYrHob7/9tp6RkaHfdNNNemhoqJ6Tk+Pu0mr5+uuv9X/+85/6J598ogP6p59+Wmv5k08+qYeEhOifffaZvnHjRv3SSy/V4+Pj9ePHj9esc+GFF+rdunXTV65cqf/vf//T27dvr19zzTUu3hNdHzp0qD5r1ix98+bN+oYNG/SLLrpIj42N1UtKSmrWufXWW/XWrVvrS5cu1deuXav3799fP+uss2qW2+12vXPnznpqaqr+yy+/6F9//bXevHlzffr06S7dly+++EL/6quv9B07dujbt2/X77//ft1sNuubN2821H782erVq/U2bdroXbt21e+8886a5420Pw899JDeqVMn/fDhwzWPI0eOGHJfjh07psfFxenXXXedvmrVKn3Pnj364sWL9V27dtWsY6RjQG5ubq3PZcmSJTqgf//997quu+6zMVxg9e3bV09LS6v5v8Ph0Fu2bKnPmDHDjVX9vT8HltPp1KOiovSnn3665rmCggLdarXq77//vq7rur5lyxYd0NesWVOzzjfffKMriqIfOnTIZbWfTG5urg7oy5cv13X9RO1ms1lfsGBBzTpbt27VAT09PV3X9RMBrqqqnp2dXbPOq6++qgcHB+sVFRWu3YE/CQsL0998803D7kdxcbGekJCgL1myRD/33HNrAsto+/PQQw/p3bp1O+kyo+3Lvffeq5999tmnXG70Y8Cdd96pt2vXTnc6nS79bAx1SbCyspJ169aRmppa85yqqqSmppKenu7Gyupn7969ZGdn19qPkJAQ+vXrV7Mf6enphIaG0rt375p1UlNTUVWVVatWubzmPyosLAR+HzV/3bp1VFVV1dqfxMREYmNja+1Ply5diIyMrFln6NChFBUVkZGR4cLqf+dwOJg/fz6lpaWkpKQYdj/S0tIYPnx4rbrBmJ/Lzp07admyJW3btmXMmDFkZmYCxtuXL774gt69e3PllVcSERFBjx49eOONN2qWG/kYUFlZydy5c7nhhhtQFMWln42hAuvo0aM4HI5aOw0QGRlJdna2m6qqv+pa/24/srOziYiIqLXcZDLRrFkzt+6r0+lk8uTJDBgwgM6dOwMnarVYLISGhtZa98/7c7L9rV7mSps2bSIwMBCr1cqtt97Kp59+SnJysuH2A2D+/PmsX7+eGTNm/GWZ0fanX79+zJ49m0WLFvHqq6+yd+9ezjnnHIqLiw23L3v27OHVV18lISGBxYsXM3HiRO644w7mzJlTqx4jHgM+++wzCgoKuO666wDX/p4ZcnoR4T5paWls3ryZH3/80d2lNFjHjh3ZsGEDhYWFfPTRR4wfP57ly5e7u6x6O3DgAHfeeSdLlizBz8/P3eWcsWHDhtX8u2vXrvTr14+4uDg+/PBD/P393VhZ/TmdTnr37s0TTzwBQI8ePdi8eTMzZ85k/Pjxbq7uzLz11lsMGzasXvNYNRZDnWE1b94cTdP+0vskJyeHqKgoN1VVf9W1/t1+REVFkZubW2u53W7n2LFjbtvXSZMmsXDhQr7//ntiYmJqno+KiqKyspKCgoJa6/95f062v9XLXMlisdC+fXt69erFjBkz6NatG//9738Ntx/r1q0jNzeXnj17YjKZMJlMLF++nBdeeAGTyURkZKSh9ufPQkND6dChA7t27TLcZxMdHU1ycnKt55KSkmoucRr1GLB//36+++47brzxxprnXPnZGCqwLBYLvXr1YunSpTXPOZ1Oli5dSkpKihsrq5/4+HiioqJq7UdRURGrVq2q2Y+UlBQKCgpYt25dzTrLli3D6XTSr18/l9ar6zqTJk3i008/ZdmyZcTHx9da3qtXL8xmc6392b59O5mZmbX2Z9OmTbX+AJcsWUJwcPBf/rBdzel0UlFRYbj9GDx4MJs2bWLDhg01j969ezNmzJiafxtpf/6spKSE3bt3Ex0dbbjPZsCAAX+59WPHjh3ExcUBxjsGVJs1axYREREMHz685jmXfjaN1m3ERebPn69brVZ99uzZ+pYtW/Sbb75ZDw0NrdX7xBMUFxfrv/zyi/7LL7/ogP7ss8/qv/zyi75//35d1090aQ0NDdU///xz/ddff9VHjBhx0i6tPXr00FetWqX/+OOPekJCglu6tE6cOFEPCQnRf/jhh1pdW8vKymrWufXWW/XY2Fh92bJl+tq1a/WUlBQ9JSWlZnl1t9YhQ4boGzZs0BctWqS3aNHC5V2O77vvPn358uX63r179V9//VW/7777dEVR9G+//dZQ+3Eqf+wlqOvG2p+77rpL/+GHH/S9e/fqP/30k56amqo3b95cz83NNdy+rF69WjeZTPq///1vfefOnfp7772n22w2fe7cuTXrGOkYoOsnemTHxsbq995771+WueqzMVxg6bquv/jii3psbKxusVj0vn376itXrnR3SX/x/fff68BfHuPHj9d1/US31n/96196ZGSkbrVa9cGDB+vbt2+v9R55eXn6NddcowcGBurBwcH69ddfrxcXF7t8X062H4A+a9asmnWOHz+u33bbbXpYWJhus9n0yy67TD98+HCt99m3b58+bNgw3d/fX2/evLl+11136VVVVS7dlxtuuEGPi4vTLRaL3qJFC33w4ME1YWWk/TiVPweWkfbn6quv1qOjo3WLxaK3atVKv/rqq2vdt2SkfdF1Xf/yyy/1zp0761arVU9MTNRff/31WsuNdAzQdV1fvHixDvylRl133Wcj82EJIYQwBEO1YQkhhPBdElhCCCEMQQJLCCGEIUhgCSGEMAQJLCGEEIYggSWEEMIQJLCEEEIYggSWEEIIQ5DAEkIIYQgSWEIIIQxBAksIIYQh/D8vrQsWTaAfZAAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plt.imshow(gr_im)" + ] + }, + { + "cell_type": "markdown", + "id": "535502fe", + "metadata": {}, + "source": [ + "## Method 1: cv2.Sobel (Convolution)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "36c99d87", + "metadata": {}, + "outputs": [], + "source": [ + "def m1_compute_gradient_magnitude(gr_im, kx, ky):\n", + "\n", + " Gx = cv2.Sobel(gr_im, cv2.CV_64F, 1, 0, ksize=kx.shape[0])\n", + " Gy = cv2.Sobel(gr_im, cv2.CV_64F, 0, 1, ksize=ky.shape[0])\n", + " \n", + " # Compute the magnitude of gradients\n", + " magnitude = np.sqrt(Gx**2 + Gy**2)\n", + " \n", + " return magnitude\n", + "\n", + "def m1_compute_gradient_direction(gr_im, kx, ky):\n", + "\n", + " Gx = cv2.Sobel(gr_im, cv2.CV_64F, 1, 0, ksize=kx.shape[0])\n", + " Gy = cv2.Sobel(gr_im, cv2.CV_64F, 0, 1, ksize=ky.shape[0])\n", + " \n", + " # Compute the direction of gradients\n", + " direction = np.arctan2(Gy, Gx)\n", + "\n", + " return direction" + ] + }, + { + "cell_type": "markdown", + "id": "a351fc50", + "metadata": {}, + "source": [ + "**Here, we use k_conv**." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "9523e11f", + "metadata": {}, + "outputs": [], + "source": [ + "m1_magnitude = m1_compute_gradient_magnitude(gr_im, kx_conv, ky_conv)\n", + "m1_direction = m1_compute_gradient_direction(gr_im, kx_conv, ky_conv)" + ] + }, + { + "cell_type": "markdown", + "id": "26e30167", + "metadata": {}, + "source": [ + "## Method 2: cv2.filter2D (Cross-Correlation)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "3224665e", + "metadata": {}, + "outputs": [], + "source": [ + "def m2_compute_gradient_magnitude(gr_im, kx, ky):\n", + "\n", + " Gx = cv2.filter2D(gr_im, cv2.CV_64F, kx)\n", + " Gy = cv2.filter2D(gr_im, cv2.CV_64F, ky)\n", + "\n", + " # Compute the magnitude of gradients\n", + " magnitude = np.sqrt(Gx**2 + Gy**2)\n", + " \n", + " return magnitude\n", + "\n", + "\n", + "def m2_compute_gradient_direction(gr_im, kx, ky):\n", + " Gx = cv2.filter2D(gr_im, cv2.CV_64F, kx)\n", + " Gy = cv2.filter2D(gr_im, cv2.CV_64F, ky)\n", + "\n", + " # Compute the direction of gradients\n", + " direction = np.arctan2(Gy, Gx)\n", + "\n", + " return direction" + ] + }, + { + "cell_type": "markdown", + "id": "ab719729", + "metadata": {}, + "source": [ + "**Here, we use k_cross**." + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "561a90aa", + "metadata": {}, + "outputs": [], + "source": [ + "m2_magnitude = m2_compute_gradient_magnitude(gr_im, kx_cross, ky_cross)\n", + "m2_direction = m2_compute_gradient_direction(gr_im, kx_cross, ky_cross)" + ] + }, + { + "cell_type": "markdown", + "id": "a84d8da8", + "metadata": {}, + "source": [ + "## Method 3: scipy.signal.convolve2d (Convolution)" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "e5ef9996", + "metadata": {}, + "outputs": [], + "source": [ + "from scipy.signal import convolve2d" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "0351c477", + "metadata": {}, + "outputs": [], + "source": [ + "def m3_compute_gradient_magnitude(gr_im, kx, ky):\n", + "\n", + " Gx = convolve2d(gr_im, kx, mode='same')\n", + " Gy = convolve2d(gr_im, ky, mode='same')\n", + " \n", + " # Compute the magnitude of gradients\n", + " magnitude = np.sqrt(Gx**2 + Gy**2).astype(np.float64)\n", + " \n", + " return magnitude\n", + "\n", + "def m3_compute_gradient_direction(gr_im, kx, ky):\n", + "\n", + " Gx = convolve2d(gr_im, kx, mode='same')\n", + " Gy = convolve2d(gr_im, ky, mode='same')\n", + " \n", + " # Compute the direction of gradients\n", + " direction = np.arctan2(Gy, Gx).astype(np.float64)\n", + " \n", + " return direction" + ] + }, + { + "cell_type": "markdown", + "id": "2393ad09", + "metadata": {}, + "source": [ + "**Here, we use k_conv**." + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "93365cc4", + "metadata": {}, + "outputs": [], + "source": [ + "# Compute the gradient magnitude and direction\n", + "m3_magnitude = m3_compute_gradient_magnitude(gr_im, kx_conv, ky_conv)\n", + "m3_direction = m3_compute_gradient_direction(gr_im, kx_conv, ky_conv)" + ] + }, + { + "cell_type": "markdown", + "id": "17c87bab", + "metadata": {}, + "source": [ + "## Method 4: scipy.ndimage.convolve (Convolution)" + ] + }, + { + "cell_type": "markdown", + "id": "ec74ef97", + "metadata": {}, + "source": [ + "This is what ChatGPT returns. But some students forget to convert the data type to float, causing errors." + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "d8f7c62a", + "metadata": {}, + "outputs": [], + "source": [ + "from scipy import ndimage" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "id": "88e3a971", + "metadata": {}, + "outputs": [], + "source": [ + "def m4_compute_gradient_magnitude(gr_im, kx, ky):\n", + "\n", + " Gx = ndimage.convolve(gr_im.astype(float), kx)\n", + " Gy = ndimage.convolve(gr_im.astype(float), ky)\n", + "\n", + " # Compute the magnitude of gradients\n", + " magnitude = np.sqrt(Gx**2 + Gy**2).astype(np.float64)\n", + "\n", + " return magnitude\n", + "\n", + "def m4_compute_gradient_direction(gr_im, kx, ky):\n", + "\n", + " Gx = ndimage.convolve(gr_im.astype(float), kx)\n", + " Gy = ndimage.convolve(gr_im.astype(float), ky)\n", + " \n", + " # Compute the direction of gradients\n", + " direction = np.arctan2(Gy, Gx).astype(np.float64)\n", + " \n", + " return direction" + ] + }, + { + "cell_type": "markdown", + "id": "ed9fad3f", + "metadata": {}, + "source": [ + "**Here, we use k_conv**." + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "b15ee435", + "metadata": {}, + "outputs": [], + "source": [ + "m4_magnitude = m4_compute_gradient_magnitude(gr_im, kx_conv, ky_conv)\n", + "m4_direction = m4_compute_gradient_direction(gr_im, kx_conv, ky_conv)" + ] + }, + { + "cell_type": "markdown", + "id": "3be8af33", + "metadata": {}, + "source": [ + "## Method 5: cv2.addWeighted / cv2.phase (Wrong Sum)" + ] + }, + { + "cell_type": "markdown", + "id": "6a9b84f8", + "metadata": {}, + "source": [ + "Some students followed the OpenCV documentation: https://docs.opencv.org/3.4/d2/d2c/tutorial_sobel_derivatives.html" + ] + }, + { + "cell_type": "markdown", + "id": "48ffd1e3", + "metadata": {}, + "source": [ + "They used $G = |G_x| + |G_y|$, rather than $G = \\sqrt{G_x^2 + G_y^2}$" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "a2b9510e", + "metadata": {}, + "outputs": [], + "source": [ + "def m5_compute_gradient_magnitude(gr_im, kx, ky):\n", + " x,y = gr_im.shape\n", + " # calculate the derivatives in x and y directions\n", + " grad_x = cv2.filter2D(gr_im, cv2.CV_64F, kx)\n", + " grad_y = cv2.filter2D(gr_im, cv2.CV_64F, ky)\n", + "\n", + " gradSum = cv2.addWeighted(grad_x, 0.5, grad_y, 0.5, 0)\n", + "\n", + " return gradSum\n", + "\n", + "def m5_compute_gradient_direction(gr_im, kx, ky):\n", + " # calculate the derivatives in x and y directions\n", + " grad_x = cv2.filter2D(gr_im, cv2.CV_64F, kx)\n", + " grad_y = cv2.filter2D(gr_im, cv2.CV_64F, ky)\n", + "\n", + " direction = cv2.phase(grad_x, grad_y, angleInDegrees=False)\n", + " return direction" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "id": "9ce074c8", + "metadata": {}, + "outputs": [], + "source": [ + "m5_magnitude = m5_compute_gradient_magnitude(gr_im, kx_cross, ky_cross)\n", + "m5_direction = m5_compute_gradient_magnitude(gr_im, kx_cross, ky_cross)" + ] + }, + { + "cell_type": "markdown", + "id": "1cd12219", + "metadata": {}, + "source": [ + "## Method 6: cv2.filter2D (Wrong, parameter -1)" + ] + }, + { + "cell_type": "markdown", + "id": "a01fc4bc", + "metadata": {}, + "source": [ + "Some students added an extra parameter -1 in the wrong place." + ] + }, + { + "cell_type": "markdown", + "id": "67cf59d3", + "metadata": {}, + "source": [ + "It's either -1 or cv2.CV_64F. Can't use both.\n", + "```\n", + "grad_x = cv2.filter2D(gr_im, -1, kx)\n", + "grad_x = cv2.filter2D(gr_im, cv2.CV_64F, kx)\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "id": "2003be8a", + "metadata": {}, + "outputs": [], + "source": [ + "def m6_compute_gradient_magnitude(gr_im, kx, ky):\n", + " grad_x = cv2.filter2D(gr_im, -1, cv2.CV_64F, kx)\n", + " grad_y = cv2.filter2D(gr_im, -1, cv2.CV_64F, ky)\n", + " \n", + " return np.sqrt(grad_x**2 + grad_y**2).astype(np.float64)\n", + "\n", + "def m6_compute_gradient_direction(gr_im, kx, ky):\n", + " grad_x = cv2.filter2D(gr_im, -1, cv2.CV_64F, kx)\n", + " grad_y = cv2.filter2D(gr_im, -1, cv2.CV_64F, ky)\n", + " \n", + " return np.arctan2(grad_y, grad_x).astype(np.float64)" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "id": "d647ef2b", + "metadata": {}, + "outputs": [], + "source": [ + "m6_magnitude = m6_compute_gradient_magnitude(gr_im, kx_cross, ky_cross)\n", + "m6_direction = m6_compute_gradient_direction(gr_im, kx_cross, ky_cross)" + ] + }, + { + "cell_type": "markdown", + "id": "710f2e14", + "metadata": {}, + "source": [ + "## Save outputs" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "id": "4e98ca9d", + "metadata": {}, + "outputs": [], + "source": [ + "np.save('data/question1_magnitude.npy', m1_magnitude)\n", + "np.save('data/question1_direction.npy', m1_direction)" + ] + }, + { + "cell_type": "markdown", + "id": "600779b5", + "metadata": {}, + "source": [ + "## Put students' implementations here" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "id": "a953ff6a", + "metadata": {}, + "outputs": [], + "source": [ + "def compute_gradient_magnitude(gr_im, kx, ky):\n", + " # Ensure the image is a float64 for computation\n", + " gr_im_float64 = gr_im.astype(np.float64)\n", + "\n", + " # Compute gradients in x and y direction\n", + " grad_x = cv2.filter2D(gr_im_float64, -1, kx.astype(np.float64))\n", + " grad_y = cv2.filter2D(gr_im_float64, -1, ky.astype(np.float64))\n", + "\n", + " # Compute gradient magnitude\n", + " magnitude = np.sqrt(grad_x**2 + grad_y**2)\n", + " \n", + " return magnitude\n", + "\n", + "def compute_gradient_direction(gr_im, kx, ky):\n", + " # Ensure the image is a float64 for computation\n", + " gr_im_float64 = gr_im.astype(np.float64)\n", + " \n", + " # Compute gradients in x and y direction\n", + " grad_x = cv2.filter2D(gr_im_float64, -1, kx.astype(np.float64))\n", + " grad_y = cv2.filter2D(gr_im_float64, -1, ky.astype(np.float64))\n", + "\n", + " # Compute gradient direction\n", + " direction = np.arctan2(grad_y, grad_x)\n", + " \n", + " return direction" + ] + }, + { + "cell_type": "markdown", + "id": "a02a91e2", + "metadata": {}, + "source": [ + "Different APIs use different kernels:\n", + "- cv2.Sobel(gr_im, cv2.CV_64F, 1, 0, ksize=kx.shape[0]) ==> k_conv\n", + "- cv2.filter2D(gr_im, cv2.CV_64F, kx) ==> k_cross\n", + "- scipy.signal.convolve2d(gr_im, kx, mode='same') ==> k_conv\n", + "- ndimage.convolve(gr_im.astype(float), kx) ==> k_conv" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "id": "63663950", + "metadata": {}, + "outputs": [], + "source": [ + "# For convolution\n", + "# magnitude = compute_gradient_magnitude(gr_im, kx_conv, ky_conv)\n", + "# direction = compute_gradient_direction(gr_im, kx_conv, ky_conv)\n", + "\n", + "# For Cross-Correlation\n", + "magnitude = compute_gradient_magnitude(gr_im, kx_cross, ky_cross)\n", + "direction = compute_gradient_direction(gr_im, kx_cross, ky_cross)" + ] + }, + { + "cell_type": "markdown", + "id": "88424988", + "metadata": {}, + "source": [ + "## Test (Should output ALL PASS)" + ] + }, + { + "cell_type": "markdown", + "id": "3a0e2419", + "metadata": {}, + "source": [ + "Restart and Run ALL for each submission" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "id": "166652ce", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "PASS: Method 1\n", + "PASS: Method 2\n", + "PASS: Method 3\n", + "PASS: Method 4\n", + "ALL PASS\n" + ] + } + ], + "source": [ + "assert np.allclose(m1_magnitude, magnitude), np.allclose(m1_direction, direction)\n", + "print (\"PASS: Method 1\")\n", + "\n", + "assert np.allclose(m2_magnitude, magnitude), np.allclose(m2_direction, direction)\n", + "print (\"PASS: Method 2\")\n", + "\n", + "assert np.allclose(m3_magnitude, magnitude), np.allclose(m3_direction, direction)\n", + "print (\"PASS: Method 3\")\n", + "\n", + "assert np.allclose(m4_magnitude, magnitude), np.allclose(m4_direction, direction)\n", + "print (\"PASS: Method 4\")\n", + "\n", + "print (\"ALL PASS\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fb616c86", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "what", + "language": "python", + "name": "what" + }, + "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.16" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/Question 2.ipynb b/Question 2.ipynb new file mode 100644 index 0000000..6f01333 --- /dev/null +++ b/Question 2.ipynb @@ -0,0 +1,419 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "f96d2a83", + "metadata": {}, + "source": [ + "## Question 2 (20 marks)" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "6f53891a", + "metadata": {}, + "outputs": [], + "source": [ + "import cv2\n", + "import numpy as np" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "d3bbb31a", + "metadata": {}, + "outputs": [], + "source": [ + "n_clusters = 100" + ] + }, + { + "cell_type": "markdown", + "id": "28068e50", + "metadata": {}, + "source": [ + "Read images" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "87dd5c72", + "metadata": {}, + "outputs": [], + "source": [ + "im_book = cv2.imread('data/books.jpg', cv2.IMREAD_GRAYSCALE)\n", + "im_mount = cv2.imread('data/mount_rushmore_1.jpg', cv2.IMREAD_GRAYSCALE)\n", + "im_notre = cv2.imread('data/notre_dame_1.jpg', cv2.IMREAD_GRAYSCALE)" + ] + }, + { + "cell_type": "markdown", + "id": "79c57454", + "metadata": {}, + "source": [ + "## Generate Clusters" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "2e571277", + "metadata": {}, + "outputs": [], + "source": [ + "def get_features(image):\n", + " image = image[:, :, np.newaxis]\n", + "\n", + " # Initialize a SIFT detector\n", + " sift = cv2.SIFT_create()\n", + "\n", + " # Detect keypoints and compute descriptors\n", + " keypoints, descriptors = sift.detectAndCompute(image, None)\n", + "\n", + " return keypoints, descriptors" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "12d042c8", + "metadata": {}, + "outputs": [], + "source": [ + "from sklearn.cluster import KMeans\n", + "\n", + "def get_clusters(keypoints, descriptors, n_clusters=100):\n", + "\n", + " # Perform k-means clustering\n", + " kmeans = KMeans(n_clusters=n_clusters, random_state=42, n_init=10)\n", + " kmeans.fit(descriptors)\n", + "\n", + " # Assign descriptors to clusters\n", + " clusters = kmeans.predict(descriptors)\n", + "\n", + " # Convert keypoints to locations (x, y coordinates)\n", + " locations = np.array([kp.pt for kp in keypoints], dtype=np.int64)\n", + " \n", + " return clusters, locations" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "b04681a0", + "metadata": {}, + "outputs": [], + "source": [ + "kpts_book, des_book = get_features(im_book)\n", + "kpts_mount, des_mount = get_features(im_mount)\n", + "kpts_notre, des_notre = get_features(im_notre)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "70d23c90", + "metadata": {}, + "outputs": [], + "source": [ + "clusters_book, locations_book = get_clusters(kpts_book, des_book, n_clusters=n_clusters)\n", + "clusters_mount, locations_mount = get_clusters(kpts_mount, des_mount, n_clusters=n_clusters)\n", + "clusters_notre, locations_notre = get_clusters(kpts_notre, des_notre, n_clusters=n_clusters)" + ] + }, + { + "cell_type": "markdown", + "id": "1518e19f", + "metadata": {}, + "source": [ + "## Method 1 (Two FOR loops)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "04f3901b", + "metadata": {}, + "outputs": [], + "source": [ + "def m1_generate_bovw_spatial_histogram(im, locations, clusters, division):\n", + " \"\"\"\n", + " Create bag of visual words representation of an image based on the division information.\n", + " \n", + " Parameters:\n", + " im (numpy.ndarray): Image array of data type uint8.\n", + " locations (numpy.ndarray): Array of shape (N, 2) with Cartesian coordinates (x, y).\n", + " clusters (numpy.ndarray): Array of shape (N,) with quantised cluster id.\n", + " division (list): List of integers of length 2 indicating division along Y and X axes.\n", + " \n", + " Returns:\n", + " numpy.ndarray: 1-dimensional array representing the BoVW spatial histogram.\n", + " \"\"\"\n", + "\n", + " # Determine the size of each division\n", + " div_height = im.shape[0] // division[0]\n", + " div_width = im.shape[1] // division[1]\n", + " \n", + " # Initialize the histogram\n", + " num_clusters = np.unique(clusters).size\n", + " histogram = np.zeros((division[0] * division[1] * num_clusters,), dtype=np.int64)\n", + "\n", + " # Two FOR loops\n", + " for div_y in range(division[0]):\n", + " for div_x in range(division[1]):\n", + " # Define the bounds of the current division\n", + " y_start = div_y * div_height\n", + " y_end = (div_y + 1) * div_height\n", + " x_start = div_x * div_width\n", + " x_end = (div_x + 1) * div_width\n", + "\n", + " # Find features within the current division\n", + " div_mask = (locations[:, 1] >= y_start) & (locations[:, 1] < y_end) & \\\n", + " (locations[:, 0] >= x_start) & (locations[:, 0] < x_end)\n", + " div_locations = locations[div_mask]\n", + " div_clusters = clusters[div_mask]\n", + "\n", + " # Calculate the histogram for the current division\n", + " for i in range(num_clusters):\n", + " cluster_mask = (div_clusters == i)\n", + " histogram[div_y * division[1] * num_clusters + div_x * num_clusters + i] = np.sum(cluster_mask)\n", + " \n", + " return histogram" + ] + }, + { + "cell_type": "markdown", + "id": "0f64e74a", + "metadata": {}, + "source": [ + "## Method 2 (One FOR loop)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "a597d45d", + "metadata": {}, + "outputs": [], + "source": [ + "def m2_generate_bovw_spatial_histogram(im, locations, clusters, division):\n", + "\n", + " img_shape = np.shape(im)\n", + "\n", + " height = img_shape[0]\n", + " width = img_shape[1]\n", + "\n", + " ## Possible Mistakes: Some students swapped x and y\n", + " div_x = division[1]\n", + " div_y = division[0]\n", + "\n", + " x_size = width / div_x\n", + " y_size = height / div_y\n", + "\n", + " num_divisions = division[0] * division[1]\n", + "\n", + " num_clusters = np.max(clusters) + 1\n", + "\n", + " histogram = np.zeros(num_clusters * num_divisions)\n", + "\n", + " # One FOR loop\n", + " for i in range(len(locations)):\n", + " point = locations[i]\n", + " cluster = clusters[i]\n", + "\n", + " x_div = np.ceil((point[0] + 1) / x_size).astype(np.int64) - 1\n", + " y_div = np.ceil((point[1] + 1) / y_size).astype(np.int64) - 1\n", + "\n", + " # Possible Mistakes: Some students miscalculated the boundary condition\n", + " # x_div = np.ceil(point[0] / x_size).astype(np.int64) - 1\n", + " # y_div = np.ceil(point[1] / y_size).astype(np.int64) - 1\n", + "\n", + " # Calculate the array position\n", + " div = x_div + (y_div * div_x)\n", + " array_pos = (div * num_clusters) + cluster\n", + " \n", + " # Update the histogram\n", + " histogram[array_pos] = histogram[array_pos] + 1\n", + "\n", + " return histogram.astype(int)" + ] + }, + { + "cell_type": "markdown", + "id": "36f2fe5e", + "metadata": {}, + "source": [ + "## Put students' implementations here" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "3804d695", + "metadata": {}, + "outputs": [], + "source": [ + "# Be careful, some students used a different function name (e.g. bowv rather than bovw)\n", + "def generate_bovw_spatial_histogram(im, locations, clusters, division):\n", + " # Determine the number of clusters\n", + " num_clusters = np.unique(clusters).size\n", + "\n", + " # Initialize histogram\n", + " spatial_histogram = np.zeros(num_clusters * np.prod(division), dtype=np.int64)\n", + "\n", + " div_size_y = im.shape[0] // division[0]\n", + " div_size_x = im.shape[1] // division[1]\n", + "\n", + " for div_y in range(division[0]):\n", + " for div_x in range(division[1]):\n", + " start_y = div_y * div_size_y\n", + " end_y = (div_y + 1) * div_size_y if div_y < division[0] - 1 else im.shape[0]\n", + " start_x = div_x * div_size_x\n", + " end_x = (div_x + 1) * div_size_x if div_x < division[1] - 1 else im.shape[1]\n", + " for loc, cluster_id in zip(locations, clusters):\n", + " x, y = loc\n", + " if start_y <= y < end_y and start_x <= x < end_x:\n", + " index = (div_y * division[1] + div_x) * num_clusters + cluster_id\n", + " spatial_histogram[index] += 1\n", + " return spatial_histogram" + ] + }, + { + "cell_type": "markdown", + "id": "0d516e6a", + "metadata": {}, + "source": [ + "## Test (Should output ALL PASS)" + ] + }, + { + "cell_type": "markdown", + "id": "a269b35f", + "metadata": {}, + "source": [ + "Restart and Run ALL for each submission" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "95d72a5c", + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Testing division: [1, 1]\n", + "PASS: Book\n", + "PASS: Mount\n", + "PASS: Notre\n", + "Testing division: [2, 2]\n", + "PASS: Book\n", + "PASS: Mount\n", + "PASS: Notre\n", + "Testing division: [2, 3]\n", + "PASS: Book\n", + "PASS: Mount\n", + "PASS: Notre\n", + "ALL PASS\n" + ] + } + ], + "source": [ + "histograms = []\n", + "for division in [ [1, 1], [2, 2], [2, 3] ]:\n", + " print('Testing division:', division)\n", + "\n", + " m1_histogram_book = m1_generate_bovw_spatial_histogram(im_book, locations_book, clusters_book, division)\n", + " m1_histogram_mount = m1_generate_bovw_spatial_histogram(im_mount, locations_mount, clusters_mount, division)\n", + " m1_histogram_notre = m1_generate_bovw_spatial_histogram(im_notre, locations_notre, clusters_notre, division)\n", + "\n", + " m2_histogram_book = m2_generate_bovw_spatial_histogram(im_book, locations_book, clusters_book, division)\n", + " m2_histogram_mount = m2_generate_bovw_spatial_histogram(im_mount, locations_mount, clusters_mount, division)\n", + " m2_histogram_notre = m2_generate_bovw_spatial_histogram(im_notre, locations_notre, clusters_notre, division)\n", + "\n", + " # Students' implementations\n", + " histogram_book = generate_bovw_spatial_histogram(im_book, locations_book, clusters_book, division)\n", + " histogram_mount = generate_bovw_spatial_histogram(im_mount, locations_mount, clusters_mount, division)\n", + " histogram_notre = generate_bovw_spatial_histogram(im_notre, locations_notre, clusters_notre, division)\n", + " \n", + " assert np.allclose(m1_histogram_book, m2_histogram_book)\n", + " assert np.allclose(m1_histogram_book, histogram_book)\n", + " print(\"PASS: Book\")\n", + "\n", + " assert np.allclose(m1_histogram_mount, m2_histogram_mount)\n", + " assert np.allclose(m1_histogram_mount, histogram_mount)\n", + " print(\"PASS: Mount\")\n", + "\n", + " assert np.allclose(m1_histogram_notre, m2_histogram_notre)\n", + " assert np.allclose(m1_histogram_notre, histogram_notre)\n", + " print(\"PASS: Notre\")\n", + "\n", + " histograms.append( [m1_histogram_book, m1_histogram_mount, m1_histogram_notre] )\n", + "print(\"ALL PASS\")" + ] + }, + { + "cell_type": "markdown", + "id": "00a11df8", + "metadata": {}, + "source": [ + "## Save Output" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "ea1fadb1", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "D:\\Anaconda3\\envs\\what\\lib\\site-packages\\numpy\\lib\\npyio.py:521: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.\n", + " arr = np.asanyarray(arr)\n" + ] + } + ], + "source": [ + "np.save('data/question_3_histogram.npy', histograms)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ddc3036c", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "what", + "language": "python", + "name": "what" + }, + "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.16" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/Question 3.ipynb b/Question 3.ipynb new file mode 100644 index 0000000..e07d119 --- /dev/null +++ b/Question 3.ipynb @@ -0,0 +1,224 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "924a2a50", + "metadata": {}, + "source": [ + "## Question 3 (10 marks)" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "2b3d1ba2", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "478bad7a", + "metadata": {}, + "outputs": [], + "source": [ + "points = np.load('data/points.npy').astype(np.uint8)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "e07b0285", + "metadata": {}, + "outputs": [], + "source": [ + "def compute_rotation_matrix(points, theta):\n", + " \"\"\"\n", + " Write a function compute_rotation_matrix(points, theta) to compute the rotation matrix in\n", + " homogeneous coordinate system to rotate a shape depicted with 2-dimensional (x,y) coordinates\n", + " points with an angle 𝜃 (theta in the definition) in the anticlockwise direction about the centre of the shape.\n", + "\n", + " Parameters:\n", + " points: a 2-dimensional numpy array of data type uint8 with shape 𝑘 × 2. Each row\n", + " of points is a Cartesian coordinate (x, y).\n", + " \n", + " theta: a floating-point number denoting the angle of rotation in degree.\n", + " \n", + " Returns:\n", + " The expected output is a 2-dimensional numpy array of data type float64 with shape 3 × 3.\n", + " \"\"\"\n", + "\n", + " # Convert theta from degrees to radians\n", + " theta_rad = np.radians(theta)\n", + "\n", + " # Calculate the centre of the shape\n", + " centre = np.mean(points, axis=0)\n", + "\n", + " # Define the translation matrices to move the centre of the shape to the origin and back\n", + " translation_to_origin = np.array([[1, 0, -centre[0]],\n", + " [0, 1, -centre[1]],\n", + " [0, 0, 1]], dtype=np.float64)\n", + "\n", + " translation_back = np.array([[1, 0, centre[0]],\n", + " [0, 1, centre[1]],\n", + " [0, 0, 1]], dtype=np.float64)\n", + "\n", + " # Define the rotation matrix about the origin\n", + " rotation = np.array([[np.cos(theta_rad), -np.sin(theta_rad), 0],\n", + " [np.sin(theta_rad), np.cos(theta_rad), 0],\n", + " [0, 0, 1]], dtype=np.float64)\n", + "\n", + " # Combine the translation and rotation into a single transformation matrix\n", + " rotation_matrix = np.dot(np.dot(translation_back, rotation), translation_to_origin)\n", + " \n", + " return rotation_matrix" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "16ef4247", + "metadata": {}, + "outputs": [], + "source": [ + "rotation_matrices = []\n", + "\n", + "for t in range(0, 365, 5):\n", + " rotation_matrices.append( compute_rotation_matrix(points, t) )" + ] + }, + { + "cell_type": "markdown", + "id": "a130c201", + "metadata": {}, + "source": [ + "## Save Output" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "23967b5d", + "metadata": {}, + "outputs": [], + "source": [ + "np.save('data/question_3_rotation_matrices.npy', rotation_matrices)" + ] + }, + { + "cell_type": "markdown", + "id": "a808d8a4", + "metadata": {}, + "source": [ + "## Put students' implementations here" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "b0836f8d", + "metadata": {}, + "outputs": [], + "source": [ + "def compute_rotation_matrix(points, theta):\n", + " # Convert points to float64\n", + " points = points.astype(np.float64)\n", + " # Calculate centre\n", + " centre = np.mean(points, axis=0)\n", + " # Compute rotation matrix\n", + " rotation_matrix = np.array([[np.cos(np.radians(theta)), -np.sin(np.radians(theta)), 0],\n", + " [np.sin(np.radians(theta)), np.cos(np.radians(theta)), 0],\n", + " [0, 0, 1]])\n", + " # Translation matrix to origin\n", + " translation_to_origin = np.array([[1, 0, -centre[0]],\n", + " [0, 1, -centre[1]],\n", + " [0, 0, 1]])\n", + " # Translation matrix to original position\n", + " translation_to_centre = np.array([[1, 0, centre[0]],\n", + " [0, 1, centre[1]],\n", + " [0, 0, 1]])\n", + " # Combine transformations with data type float64 \n", + " combined_matrix = np.dot(np.dot(translation_to_centre, rotation_matrix), translation_to_origin).astype(np.float64)\n", + " return combined_matrix\n", + "\n", + " return 0" + ] + }, + { + "cell_type": "markdown", + "id": "73b68192", + "metadata": {}, + "source": [ + "## Test (Should output ALL PASS)" + ] + }, + { + "cell_type": "markdown", + "id": "1c0a88a6", + "metadata": {}, + "source": [ + "Restart and Run ALL for each submission" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "132d734b", + "metadata": { + "scrolled": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "73 PASS\n", + "ALL PASS\n" + ] + } + ], + "source": [ + "n_pass = 0\n", + "for t in range(0, 365, 5):\n", + " if np.allclose(compute_rotation_matrix(points, t), rotation_matrices[int(t / 5)]):\n", + " n_pass = n_pass + 1\n", + "\n", + "print(n_pass, \"PASS\")\n", + "assert n_pass == len(rotation_matrices)\n", + "print(\"ALL PASS\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0fac308a", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "what", + "language": "python", + "name": "what" + }, + "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.16" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/README.md b/README.md new file mode 100644 index 0000000..9479457 --- /dev/null +++ b/README.md @@ -0,0 +1,2 @@ +## ECMM426 Template + diff --git a/ca_utils.py b/ca_utils.py new file mode 100644 index 0000000..aeb2f5a --- /dev/null +++ b/ca_utils.py @@ -0,0 +1,187 @@ +import cv2 +import math +import torch +import pickle +import numpy as np +import torch.nn as nn +import torch.nn.functional as F + +def im2single(im): + im = im.astype(np.float32) / 255 + return im + +def single2im(im): + im *= 255 + im = im.astype(np.uint8) + return im + +def load_interest_points(eval_file): + """ + This function is provided for development and debugging but cannot be used in + the final handin. It 'cheats' by generating interest points from known + correspondences. It will only work for the 3 image pairs with known + correspondences. + + Args: + - eval_file: string representing the file path to the list of known correspondences + - scale_factor: Python float representing the scale needed to map from the original + image coordinates to the resolution being used for the current experiment. + + Returns: + - x1: A numpy array of shape (k,) containing ground truth x-coordinates of imgA correspondence pts + - y1: A numpy array of shape (k,) containing ground truth y-coordinates of imgA correspondence pts + - x2: A numpy array of shape (k,) containing ground truth x-coordinates of imgB correspondence pts + - y2: A numpy array of shape (k,) containing ground truth y-coordinates of imgB correspondence pts + """ + with open(eval_file, 'rb') as f: + d = pickle.load(f, encoding='latin1') + scale_factor = 1.0 + return d['x1'] * scale_factor, d['y1'] * scale_factor, d['x2'] * scale_factor, d['y2'] * scale_factor + +def show_interest_points(img, X, Y): + """ + Visualized interest points on an image with random colors + + Args: + - img: A numpy array of shape (M,N,C) + - X: A numpy array of shape (k,) containing x-locations of interest points + - Y: A numpy array of shape (k,) containing y-locations of interest points + + Returns: + - newImg: A numpy array of shape (M,N,C) showing the original image with + colored circles at keypoints plotted on top of it + """ + newImg = img.copy() + for x, y in zip(X.astype(int), Y.astype(int)): + cur_color = np.random.rand(3) + newImg = cv2.circle(newImg, (int(x), int(y)), 10, cur_color, -1) + return newImg + +def conv3x3(in_planes, out_planes, stride=1): + """ + 3x3 convolution with padding + """ + return nn.Conv2d(in_planes, out_planes, kernel_size=3, stride=stride, padding=1, bias=False) + +class BasicBlock(nn.Module): + expansion = 1 + + def __init__(self, inplanes, planes, stride=1, downsample=None): + super(BasicBlock, self).__init__() + self.conv1 = conv3x3(inplanes, planes, stride) + self.bn1 = nn.BatchNorm2d(planes) + self.relu = nn.ReLU(inplace=True) + self.conv2 = conv3x3(planes, planes) + self.bn2 = nn.BatchNorm2d(planes) + self.downsample = downsample + self.stride = stride + + def forward(self, x): + out = self.conv1(x) + out = self.bn1(out) + out = self.relu(out) + out = self.conv2(out) + out = self.bn2(out) + + if self.downsample is not None: + residual = self.downsample(x) + else: + residual = x + + out += residual + out = self.relu(out) + return out + +class Bottleneck(nn.Module): + expansion = 4 + + def __init__(self, inplanes, planes, stride=1, downsample=None): + super(Bottleneck, self).__init__() + self.conv1 = nn.Conv2d(inplanes, planes, kernel_size=1, bias=False) + self.bn1 = nn.BatchNorm2d(planes) + self.conv2 = nn.Conv2d(planes, planes, kernel_size=3, stride=stride, padding=1, bias=False) + self.bn2 = nn.BatchNorm2d(planes) + self.conv3 = nn.Conv2d(planes, planes*4, kernel_size=1, bias=False) + self.bn3 = nn.BatchNorm2d(planes*4) + self.relu = nn.ReLU(inplace=True) + self.downsample = downsample + self.stride = stride + + def forward(self, x): + out = self.conv1(x) + out = self.bn1(out) + out = self.relu(out) + + out = self.conv2(out) + out = self.bn2(out) + out = self.relu(out) + + out = self.conv3(out) + out = self.bn3(out) + + if self.downsample is not None: + residual = self.downsample(x) + else: + residual = x + + out += residual + out = self.relu(out) + return out + +class ResNet(nn.Module): + + def __init__(self, block, layers, in_channels=3, channels=[16, 32, 64], num_classes=10, flatten=True): + super(ResNet, self).__init__() + self.name = "resnet" + self.flatten = flatten + self.channels = channels + self.inplanes = channels[0] + self.conv1 = nn.Conv2d(in_channels, channels[0], kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) + self.bn1 = nn.BatchNorm2d(channels[0]) + self.relu = nn.ReLU(inplace=True) + self.layer1 = self._make_layer(block, channels[0], layers[0]) + self.layer2 = self._make_layer(block, channels[1], layers[1], stride=2) + self.layer3 = self._make_layer(block, channels[2], layers[2], stride=2) + self.avgpool = nn.AdaptiveAvgPool2d(1) # global pooling + self.fc = nn.Linear(channels[2], num_classes) # global pooling + if flatten: + self.feature_size = channels[2]*block.expansion + + for m in self.modules(): + if isinstance(m, nn.Conv2d): + n = m.kernel_size[0] * m.kernel_size[1] * m.out_channels + m.weight.data.normal_(0, math.sqrt(2. / n)) + elif isinstance(m, nn.BatchNorm2d): + m.weight.data.fill_(1) + m.bias.data.zero_() + + def _make_layer(self, block, planes, blocks, stride=1): + downsample = None + if stride != 1 or self.inplanes != planes * block.expansion: + downsample = nn.Sequential( + nn.Conv2d(self.inplanes, planes * block.expansion, kernel_size=1, stride=stride, bias=False), + nn.BatchNorm2d(planes * block.expansion) + ) + + layers = [] + layers.append(block(self.inplanes, planes, stride, downsample)) + self.inplanes = planes * block.expansion + for _ in range(1, blocks): + layers.append(block(self.inplanes, planes)) + + return nn.Sequential(*layers) + + def forward(self, x): + x = self.conv1(x) + x = self.bn1(x) + x = self.relu(x) + + x = self.layer1(x) + x = self.layer2(x) + x = self.layer3(x) + + if self.flatten: + x = self.avgpool(x) + x = torch.flatten(x, 1) + x = self.fc(x) + return x diff --git a/data/books.jpg b/data/books.jpg new file mode 100644 index 0000000..ba9ecc5 Binary files /dev/null and b/data/books.jpg differ diff --git a/data/mask.png b/data/mask.png new file mode 100644 index 0000000..3cc0044 Binary files /dev/null and b/data/mask.png differ diff --git a/data/mount_rushmore_1.jpg b/data/mount_rushmore_1.jpg new file mode 100644 index 0000000..9e6e2b9 Binary files /dev/null and b/data/mount_rushmore_1.jpg differ diff --git a/data/notre_dame_1.jpg b/data/notre_dame_1.jpg new file mode 100644 index 0000000..da3924b Binary files /dev/null and b/data/notre_dame_1.jpg differ diff --git a/data/points.npy b/data/points.npy new file mode 100644 index 0000000..a243090 Binary files /dev/null and b/data/points.npy differ diff --git a/data/question1_direction.npy b/data/question1_direction.npy new file mode 100644 index 0000000..27666ca Binary files /dev/null and b/data/question1_direction.npy differ diff --git a/data/question1_magnitude.npy b/data/question1_magnitude.npy new file mode 100644 index 0000000..a731eab Binary files /dev/null and b/data/question1_magnitude.npy differ diff --git a/data/question_3.npy b/data/question_3.npy new file mode 100644 index 0000000..89ef3f2 Binary files /dev/null and b/data/question_3.npy differ diff --git a/data/question_3_histogram.npy b/data/question_3_histogram.npy new file mode 100644 index 0000000..7f6f492 Binary files /dev/null and b/data/question_3_histogram.npy differ diff --git a/data/question_3_rotation_matrices.npy b/data/question_3_rotation_matrices.npy new file mode 100644 index 0000000..df8ab58 Binary files /dev/null and b/data/question_3_rotation_matrices.npy differ diff --git a/data/shapes.png b/data/shapes.png new file mode 100644 index 0000000..2bad356 Binary files /dev/null and b/data/shapes.png differ