classes/3_Lecture_Filters.ipynb

{"cells":[{"cell_type":"markdown","metadata":{"id":"4cpfMcRjPZM1"},"source":["<Body>   \n","<img src=\"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAVEAAAB+CAYAAACd+yIVAAABOGlDQ1BrQ0dDb2xvclNwYWNlQWRvYmVSR0IxOTk4AAAokWNgYFJILCjIYRJgYMjNKykKcndSiIiMUmB/xsDEwMLAySDMwJyYXFzgGBDgwwAEMBoVfLvGwAiiL+uCzJpivO1cjMWGdoHpqRe3pXw2xVSPArhSUouTgfQfIM5OLigqYWBgzACylctLCkDsHiBbJCkbzF4AYhcBHQhkbwGx0yHsE2A1EPYdsJqQIGcg+wOQzZcEZjOB7OJLh7AFQGyovSAg6JiSn5SqAPK9hqGlpYUmATeTDEpSK0pAtHN+QWVRZnpGiYIjMKRSFTzzkvV0FIwMjIwYGEDhDlH9ORAcnoxiZxBiCIAQmyPBwOC/lIGB5Q9CzKSXgWGBDgMD/1SEmJohA4OAPgPDvjnJpUVlUGMYmYwZGAjxATPHUlQo3ou4AAAAOGVYSWZNTQAqAAAACAABh2kABAAAAAEAAAAaAAAAAAACoAIABAAAAAEAAAFRoAMABAAAAAEAAAB+AAAAAAQtgrIAAD7VSURBVHgB7Z0JvFVTF8B3lCJKZoVGUpKKFFIRKUolMjUoZQgRiagkklSmMqRS0WAeylR9IiRjZMwQGSqhEJWI963/7u1r3/POOXe+77739vr93jvn7rPHdfZZe+211l6r1OjRo/M++vBD5SA+DGxburQadPXVqtb++8dXwOVyGHAYKNYYKL1gwQL16SefFOtBpnNw5cqVU2vXrVO10lmpq8thwGGgyGJgmyLb80LqeJkyZQqpZdesw4DDQC5iwBHRXHwrrk8OAw4DRQYDjogWmVflOuow4DCQixhwRDQX34rrk8OAw0CRwYAjokXmVbmOOgw4DOQiBhwRzcW34vrkMOAwUGQw4IhokXlVrqMOAw4DuYgBR0Rz8a24PjkMOAwUGQw4IlpkXlVud/Tff/9Vf/zxh9q8eXNud9T1LusY+Pvvv9Xvv/+u/vnnn6y3nY0GS2ejkZLWxvr169Vff/2VlmGXlmOmGPjvsMMOqlSpUmmpM9VKNm3apF5//XX1ysKF6uOPP1bLly9Xv/76q+5fXl6e2m677dS+++2nDpCjsUcccYQ65thj1T777JNqs5Hya9euVbSTLgDH4Jd+pxvoJ/3NVWBO7brrrmnr3po1a9TLL7+sFi1apL74/HP19ddf64WVdsDFTjvtpGrWrKnq1KmjmjdvrpodfbTacccd09Z+YVSUEBHt1q2bGn7D8LT2s80JbdRnn30WV51NmjRRDz38UFx5TaYXX3xR9T63t/mZletVAwfqieRt7F+ZRH8FcGoQym233TZSBM7OJsTbbLONqlSpkiZG1atXV/Xq1VONDz9cHXzwwVkjrt9/952aPHmyevzxxzVnwUfQ+dRTVaNGjdR+QjS33357zW2sE6Lx+Rdf6A/p3nvvVUOHDlXNmjVTPXv1UscKQU0V2rRpo/4QzsYLcDpwPX4AgQSHNpDX5o4gpnzkVatVUzVr1FAH1K6t+123bl27WEL3GzZsUI0POyyhMmSmr/SHazzjYr4wnkQXl+3Klo37+wsbxGuvvaam3H+/eumll1T58uVVu3btVP/+/VXdgw5Se+21lx4Lu5TVq1erjz76SL0kx815DnHt1KmTOrd3b1VDcF4UoVTbtm3z4j07X7FiRbXHHntEjZOPY+h1Q6PS/H789ONP6uyzzy7w6JtvvokiFgUyWAmcW99333014itUqKDan9xede/e3cqh1JYtW9TEiRPV/Pnz5UP7QzGJV61aFZUnlR98ZBMnTVIQ9EThnXfeUacJ0fGDIUOGqF7nnhv1iEn322+/qZXff6/elrKLFy9WC2WVtz+U3XbbTbVv317jtmatzJzo37hxo7r9ttvUtGnT9LtqKETzxhtvVPEQF97Hww89pG6++Wa93T9cCP91w4bFVTYKGXH8mDF9uho8eLBvzukzZqijjjqqwLPNf/6p1vz4o/ryyy/V22+9pZ577jn17bffRuXbfffdZa6drHrJIlClSpWoZ7F+IOI4WBY8P4A7byqc+kFCaCAgEBu4QjizskLcDDz6yCNqoCzMfnC/EC44fQPMGbbOP//8s1ot855dwodCtBbLzuGnn34y2SLXVInol7JYDh8+XL366qu6Tr7Hy6+4QkErYgF9HDN6tHr44Yc1AwF9uGLAAMW3XZQgIU6UD5o/G6rXqG7/DLz/e8vf6gtBeCrwp0x4uw6I0pa/twjx6RWpdszoMWrChAmR30X5hg+JRYs/CNd5552nt0d3jR+vuUHGxkScMmWK/mP1v3rQoIQ/9DAcvf/+++qySy9VLHbAGWecoW4cMSKKaw4rDzd1dteu6kghYD3kA3tLCFXHDh3UgCuvVH369MkaFx3Ux7KyMMNF8weXPPCqqzQhHSaE/ud8ogPxuV848GlTp6pOp5yirrnmGr0rCKozVvoJwkn37dtX1a9fP1bWhJ8zZ/hjcT3wwAMjBJaFl0X4zjvvVG++8UbC9foVYFEdedNNersOp3/7HXcoYcr8svqm0cebR43SC8mVQjwfeOABxc7xtttvV40bN/Ytk4uJ2+RipxLpE5MCbgeAyPJiizOwlR8zdqz+g4uw4ZlnnlEntG6tnn32WTs56fs5s2er07t0iRDQdsLx3jRyZNwE1G6Yfs+YOVPtsssuetvJx3fJxRfnnCKK7eVJJ52kwGXVqlXtIeht9WOPPqray2K1YsWKqGfx/oBrQ8SRCQIa1gfGdeSRR6pZs2apCy64ICxrzGeIDQYI0Rt23XWR98ecTISA2o107NhRjZSdCrBy5Up1pizUcKdFBYo8EYUz/vTTTzW+kbVASEsCdO7cWY2TBcQLiC8uvugi9UiKk5DFqF+/fhFRC2KUkUJA+RiTBV1H/sdCHRD77iJnR1GVa7DnnnuqccLxe+Wo9JMPvec552jxRCL97iAceDeP+CmR8unIy/uD205GHEX7yOl7i/zy8ccei3TnrLPO0iKlSEISN6eKmOsU4fIBZMBXSx/vvuuuJGrKfpEiT0RB2Xei8ABWrUyf7FNXmOP/Wp9wguoinKIfXHvtteq9997zexQzDQIMl2HD2FtvTYsWtbVwyiijDLC97yMfZS6aRqG088rcTb/hRMePG2d+xnW99LLL4sqX6UwQ0stEqZMosOO7WHYPWGUYYGG8NkAObfLEex12/fVaLmzyi8N4LUYxv3P1WiyI6MYNGzV+N27aes1VZGeiX3wMyB29wIQfIpPbVkJ58/j9xnRpkMhVbUAbnk4Z1ZUiD7UVJ5jDEC0gZUiBSw5q+5yePYMeqemiyELpFg+wfUekkSsAJ4oiKxEYIcrE+fPmRRW5QpRImIelA1Daeon7DTfcoOa+8EI6qs9YHcWCiGLeAZhrxrCVgxXvvffeqkWLFr49w4bTaE19M3gSsWJAFODFI/KvdAJbZS+H9+STTyo0zSmBKE/SDchFGzRs6FstohObK/PNlJ+IvWwugZGRxtsn3s9UUazZUFtMwE4WEUU6gW19tWrVoqqEUGMFkKtQLIhoriI3W/2yTVy8bc4W5VC8gH3rL7/8EpUdZUQmTKfQ2HthlMhLbesL7/PC+n3ooYcGNv3Ou+8GPrMfHBRg5mTnyfZ9vH3CtnOomOB5oau8w1Rk5N76+I2t9JkiY7WBxQqbUqNAtp/lwr0jornwFlLsQ8MGDQJreEPMWuIBbBExmPZCkMzVmy/R33B42IzagNICUxcvJ2znKYz7A4XjCoKvv/oq6FFUulfTH/WwkH54Ob6gbgwW+Tr2rjZgGXKy2M5mAjqJtt6r0ENpPPG++zLRXMp1OiKaMgoLv4IacowuCNAkx9J+Y9EwZsyYAlXAFbQ67rgC6elKOM6n7qVLl6qnn346uSYyIBOlIxV33jmwP+s8nHtQRgz2cw3i6RNG+gSz9ALiiQpxGNR7y8Xze3exiz7EhzEYL9YS2EXnGjgimmtvJIn+cJILoXwQ/CgncsIAI2e/PAeLMiST55qbNG3q261bxeYwqa1bBmSidHBHOcYYBPFaFcRzgieojUyl7xyyOJg2bxENuR9kWsbb1GduoMTLRbMnR0T9ZkgRTAvTkIZpkLHJ4ySOHzTxbLf98qSSxnFHW0tv6vpejrnOnTvX/Cz0K8cog2AH8RcQD/iNM55ymcwTy+EKJnLvB5jJeUUx6e5nkBwaI/yw95HufsRTnyOi8WCpCOTxM3My3Q77WHAYEeRb4EDxtJNJQFwQpLR68MEHM9l0QnX/Ih6qggA7yTAA99js2s5lwvJn8xnc8VVi1B4EM8TfgB+gTOJIaSah9gEH+FYPQ/DEE0/4PiusREdECwvzaW43TO4Z5ursWTneGAQHBEzkoPzJpO8f4DSF890/xRBDJNNeMmWWi3OSIDg8hiMaiGhv8RGQi4AYCKcqfsDRznkBuwEWDjx2ZRKqiHMW+ucHHEfOJXBENJfeRpJ9QX7oNU0yVUFAg2RfbOXhRIMgG8bhOP4Igvn/+1/Qo6ymLxLlih8gQkn2vLhffbmUxkmyoG1z9Sy4rIPbDfJBu2TJkpxSMDkimkszN8m+4Pg26GRSmAJgmfgc8HrlMl3AUUimuQ3a2rtyZdNkgSua4cIGuNBP5NCCH1wkBxOCFii//EUpLcw0rkrIO0vnGCsHtMNcf/PNN9PZVEp1OSKaEvpyo/Dbb78d2BH8YAbB+2JOFARBEzgof7Lpe8nppSB4T9zwFTYEaadxZ3d+it6QCntsYe1jahYE2Zobe4YcS8VFY66AI6K58iZS6EeQXBOZpp8tpmlq2bJl5rbAFU40GxDGyeGMOmhLmY2+YfrllQuyzUTGeZd4GMpFZVG68PJp2NxIYziRsP7uHGKHGq8j+bD60/WsoOeKdNXs6skKBuAYcBriBT5wfH96T37Y+SBSQRBmYB5UJpn0WAbbHBbItCbY22+UKhPE5+dYsVe1Aa9Z/S65RMV7XNIuW5TusX01Dqn9+h228PnlTzYtzLYWM7hcAUdEc+VNJNEPThoFeT8aLt5vgmztTFOciQ6CMOP9oDLJpMcy5sf8KltElDPa88RLEQbdhAsB8L6Okw08+u8vgfdKAvzwww+hw8zW3CgfEsAubO6Gdj4DDx0RzQBSs1ElH3zfCy+MOKQ2bWIWAgdK8K9YQB1BUFZMc7IBBOgLgw2eM9thecOeETfJmHpxNh+TMILdcWzzc4lKiWNvzKrKyLjriA3kpRIS5QhxvnKYBJkrztt2P5zFwnmsd+ZXZzJpYe3gZ4EdQ1ieZNpMpowjoslgrZDLvCEf+0DxyWmcUZvutGjZUvsQDTJgN/nMdVNIFIBsTc6wQwL0M8z+1Ywjnutt4lTae2qID/FXiYxgIrASmZLInPuLLPmQQw7RYX1LGgHVOA+ZFzzfLsbCF8/7iCdPGR8/uXY55ka25qndrvfeEVEvRgrp9ySJIFpaJichJPxkQXzwBBrD56bXhyVn0AkmR+TIRKBUSOYgk6mQIsk9inHePV2u1gh+5hftk07jN+A9sT0krhInpYwXKbatLEydJWwF13T1JTlEZa9U2LzIXi+USr932Mz03hHRzOA14VpryxbyCYnnfqPIMglFTAhdZELr169Xa0RGhUmHcXaBkTfnzvEjCtFN1uQkzA40KQcgCY9aYvbIliwMymX4ZAxtE00VkyX++nzwgTpXQlejWMEy4Jk5c/QfclnkzOn08B827sJ8Fgvnf+cHhsx0H7fEmBth8zfTfbPrd0TUxkYh3h/drJmOO4/xO8bdaB8hoHCgDcWzepfTT9fxuDlFxF+Y1j3eYewY4vkpFnGLt41Y+RhfGIR5UAorl+wzwnjcc889Osqp4UipC3OwM+QdjJIQv6eedlqy1ReJcrEUR8giswGbQ+YG/kxzYSsPHhwRzcZsSKANtvIoNLIBnDwJOo0TS7mQrv7FsgOtXKVKupqKux6USQTTw1G1DRDVqyUWVC3R0jfw8Xdp5y3K98ReQnQRJNJBIZcNgIkIgsoSFidXwBnb58qbKIR+BJ1NpitBx0HT3c1fQzwk0VaVQiCitBvkmAN/A9dLVMriDHB4xMEKgjDiFlQmmfTfQuZG2NxNpq1Uyjgimgr2inhZZK9BEOb+LahMMunr1q0LLMaHEsuONLBwig+QgeJJyA/wsbkkzthKfuWLQhoy9yDI1txYGzI36mTYTWPQ2P3SHRH1w0oJSQvbkqLMygaskhNJQRDWv6Ay6UxHFh0EaPKLM/iF5zDj/SHkkIbJk45r2NwI61862k6kjqwR0ZJob5fIiyiMvNiTopn2gzVr1ii2rpmGsON72ZINB42xZkjsqkRCUQfVn8vpzY46KrB7q7JERDny6wfIa5vE8OPqVy5TaSkTUWN2E6uDmdSkGfu9f7Zk/qOPNc6i9By8Hd+6tW+XUaIETWLfAkkmfiVu/IIgzHlKUJl0pod5redYqDcCZjrbLuy6GggXvltAcD1Of2Ua8GAfdPyUHQLHcXMFUiaiGzdsjGssscwm4qokIJM5ibLpz00BOVxyEAbCwt5mIwb855995ts1YvgEccm+BUiURSGdEMtv5vLly9PZXE7VxQLbvl073z4xbtv8yzdTiokcxQ2Cdu3bBz0qlPSUieiq1avi6jicaKYIacWdK+o+rFsbrKSIq5MlMBPEqnZAXPWwiZwOVKGZD+I2unbrlngTMU4/JVphmMNo6vr2m28SrbJI5Q96BxyT/SbDY8eXgR/gG6Jz585+jwotLWUiunrVaoU3oXggU+YqVSpvtSUMk6/F07+SmicoBtA7Ic6e04ErHH74AfOkjZweKmzYO4Yt4krxMFWcgVNzrVq18h1iYc2N0+SgQ4UKFXz7VFiJKRNR2Pp4vUzX2r9W2se5XdntVPUa1XW9HweEcUh7o8WswlPkbLifmzc85mdy2+bnBxXUXnbZZcmdRknzdp4gc2GEdHUxJ6K8iysHDvT1GZDp8Bz4ifACXOgl4s811yBlIsqAFr68MK5x1T+4flz5Esl0eOPDtasynEisWLEikaIubz4GOEI6eMiQAvjgNJHfZC6QMYkEiDO+O72AfWInIepJQZq38/ShWrVqgV3JJZ+WgZ1M8QGinjPPPLNALQsWLNCu6Ao8SEMCNrh8z164QFw/7h5gTeLNm83faSGic8RJQzzmMEFedFIZ8IknnaiLz58/P+5q8Mzdo0cPNXz4cDXomkGqeYvmcZctrhmbN2+uunTpUmB43qOPBTIkmfDqK68UkIfiFo+YRrlkDufHoZshF/ftvBnnoGuuKeDkhuiyC1580WRJ6/XRRx8tUB/EvG/fvgXScyEhLUQUU5jnnnsu5njqHlRXhYXIjVmBJ0OlSpWU0S7PmjnL89T/ZzNx9LFw4UI17Pphqlv3buq8885T06ZNU9MemKbKly/vX6iEpMKNIgez4YUXXshI/PepgnMvXDFggPZg5U2P+3eat/O0i4/RICgJ23nGzqmxW2+7rcDiRgyqdAMn2GZ74sqzjceVYSbNJFMZR1qIKB0YO2as2vzn5ph96dq1a8w88Wbo37+/JnwQxXjkoVWrVlUT7pugKlQsKJiGE7tl9C3xNl0s82E9cd9990UdtcQO+M5x49I63nfeeUe97Il3f+KJJ6oLUo2emYHtfNipKbixeJWqaUVgIVSGcfuQoUOjWkamHSTXjsqYwI+7775bYSNqA56zcumYp9037tNGRDF5GD1mtLf+Ar+7duuauP1fgVqU1hpSFx/5jTfc6JOjYFLv3r0VvjiDgA+5Zq2aQY9LRDqnmCZOnBjlBf6hWbN0CI10IACxz4gbo98XZlZjxfN8LgIRU8MCsxEDqqQAIrDzzz8/argjRoxQ6fI9+9VXX2mn2HYDA0WxRYyrXIa0EVEGOXnSZPXUU0+FjhdHqhDbVPxhNmnaRLijO7XW8KYRN0WCioU2LA8bNGwQK4tqcEjsPDErKeIZ8JA/YcKECCHlI0FjHsv3ZzzDhtOwrTkaNmqkJou3frZsuQgYnbds2TKwa99l4fROYOOF8ODqQYNUz549Iy3jSvEO2WqnCsyx/swxYYoM9OvXT12Yo3JQ00euaSWiVDjgigFqzuw53AYCW+dbhfPAhCRRQPkxbeo0zVFOmjhJJSKXKaVKxW4ujiyxKyn6OQiHMeuhh5SJP0+cbxRxqQBbP/uDIwTxzJkzo8QHqdSfqbIdO3YMrPrDDz8MfFZcHwy97joFMTXAwohILRW46aab1AcSVQCAwbpRdiv9L788lSqzVjbtRJTtGivIrWNvDdXYd+jYQT319FOBcW+8GID7nD5juhp1yyiJyFhGjRYtLluJRGDpB0tjZv/wg5L3UQQhhTPKc8RbETGcgBnTp2tCmoztKA47zuvTR88JFARXiXNjPMjnKgdq46R5ixaBYZvjUajadRWXe7b1KAc5X898uFDk2Zg9JQo4fh51881qiuxGAELdzBTx0dlp1J0k2qdE82fMs/04UUa8JMoDuJeGjfxdiiEshjB+LU4oWMmWfbpMrV27Vm/T0QhWrlJZ1T6gtnzETSJyVMI0DL52sHo3CX+OkydPlqBjnVXZcmV98YSZVDqOOhI+wY/QhMWM2SKLT5AzF+MbwLfTGU7Uk1q4RQjo7bJtY7IvlXhP18t7rVevXszW10u4k/HjxysC8fHBICoYKgqKVBQFiBX8vK6DwyAA9374Zbsea0dEHmS2Z4m9pNdZNccTx4wZoy4XrsmIqHBMskmUI4naNPr1j/GEjYt4R0HlMj1vWsjiMnfuXFEqj1EPya7l3F69tOlgPwmaaHYwQe+D9M/Eb8J1wtVycg3ztm5y1PcyURbn2omksDHwrFTbtm3z2KplErDDRCt/zDHHaGQl2hYECaPvmTNmqueff973A4q3To6x3XHnHQXMmd566y2F4un39b+HVoUGe6IQhDBXXL1EZsQCki5YLJOMkA2FDRjfs/2eLhExOWJL0LbW4gUK3444UMZEDNnWTxLkDccirwj3+byYvqFtPVaC6vU45xx19NFHpzyMemJ2tGHDhpTroYKd5Aih2UbGqhDjesQRcJ/esCZ4fCI+E2NHyXqJ7MbOPvvsWFVGnqOgOirNYWE+EycesRaISAdSvPnyiy/UlClT1JNPPqm/T0Q1iO3qiHPrPcRLPv3AkgFzSPA9Xw5aLFq0SOP/VDkLf458M+k0f0xxOAkVzwoRNT2CAGFw30iUCXXq1tGsO556UDZhYA0Ht2njJj0RV65aKR/i52rp0qUa2ZiTpAt2ly0ITgzQxNPeotcXqXlz58VFnOMhoog0/DilZPsfKzZ7svUmW46xIQt8+eWX1UdyxS0cOwgIG1t1AuBBVPYXTT+iAD6mdLouS5c22Iw/UfzSPh6uUCrhRAVOEa5vF7Fbrlqtmg4kCPeaKBT2uBLtr19+YsEj+1702mvqU9k1rpBdJuFEIKBYxlSUgy41JNAikQOOlnmBZUamOWa/fqYzLatENJ0dL6y64iGihdU3167DgMNA9jGQdsVS9ofgWnQYcBhwGCg8DDgiWni4dy07DDgMFAMMOCJaDF6iG4LDgMNA4WHAEdHCw71r2WHAYaAYYMAR0WLwEt0QHAYcBgoPA46IFh7uXcsOAw4DxQADjogWg5fohuAw4DBQeBhwRLTwcO9adhhwGCgGGHBENMGXyKkqBw4DDgMOAwYDpVsff7zKVChj00hxupaW46m77bprcRqSG4vDgMNAChgoJeeg81Io74o6DDgMOAyUaAy47XyJfv1u8A4DDgOpYsAR0VQx6Mo7DDgMlGgMOCJaol+/G7zDgMNAqhhwRDRVDLryDgMOAyUaA46IlujX7wbvMOAwkCoGHBFNFYOuvMOAw0CJxoAjoiX69bvBOww4DKSKAUdEU8WgK+8w4DBQojHgiGiJfv1u8A4DDgOpYsAR0VQx6Mo7DDgMlGgMOCJaol+/G7zDgMNAqhhwRDRVDLryDgMOAyUaA46IlujX7wbvMOAwkCoGHBFNFYOuvMOAw0CJxoAjoiX69bvBOww4DKSKAUdEU8WgK+8w4DBQojHgiGiJfv1u8A4DDgOpYsAR0VQx6Mo7DDgMlGgMOCJaol+/G7zDgMNAqhgonWoFrnxiGPjtt98CC5QrW1aVLVcu6vnGjRsVEUbLSXpZeR4GhMv68ssv1U8//aS2lYB6NWrUULvvvnukyKZNm9Rff/2lSpcurcqXLx9JNzcbNmxQW7ZsUdttt53afvvtdTL5KVdKflWoWNFkjbqaPHaiXYed7ndvygf1iz7RNxuC8tp5wnDN+OijHxg82c/C8pPPjMEuY99XDMAdeX7//Xf177//2tkj97xz3j3wxx9/qH/++UeVKlVKVahQIZLH3Pz5559q8+bNqkyZMmqHHXYwyVFXyi9fvlytW7dOz4M99thD7bvvvrpOO2NYn+x3m2yf1q9fr5ivQfOaOc/cB4JwZ+rYcccd9Xy3+8+9mc8mHbyBF+aOF9bLd0mwuaD3bOaE/T4idRCozkF2MCCTPK9a1aqBf2NGjy7QkWOPPVbnv/jiiws8Mwny4eTdcccdeY0aNixQ96mdO+e9//77OuvgwYP18+7dupmiUdc+vXvr5/0vuyyS/vjjj+u0A2vXjqR5b5568skC7TLO2gcckNfltNPyZs2alSdExlss8nv48OG6fMMGDXzzLV682Lf+/WvVyuvYoUPepIkT82SSR+rjJhauH3rooaj89o+rr7rKt71D6tfPA0evvPKKnV3fP/Lww75lzPv+09M/u4JjjzkmsOyw666LZG130kmRfG+88UYk3dwYPF526aUmKXL9cc2avEFXX513UN26kTpM3xgXz8CZgeNatSqQz+QfIvPIQIeTT47kW/z66yY5ch1x4436+SWXXBJJ46bJ4Yfr9Pvvvz8q3fyY+8ILkXqF8JvkyHXVqlV51atV03lmzpwZSbdvzHw2/eZKmbZt2uS9tGCBnTWvzoEH6roee+yxqHTzg3lG+e7du5ukyFWT5C3CbSyaNilCWBO5qbjX3qpB+066yPI3X1fff/B+IsUjeQ9u007tsu9+6q9NG9Xi6VMj6Yee0kXtuOtukd/x3ix99mn166qVOvu+DRqpGo2bRhV97+nH1fof10Sl2T9KbbONKivc2g47V1J7H1hX7Vq1eoHV2s6f6P1OwkmU9XBCO3i4ww8++EB9JVwDMH/ePM2JsOrawAopRFG98847keTddttNr/Jr167V6d9+84065JBDIs8zeQMHDMDxwBW99dZb+u/xxx5TkydPLsDNkm/27Nm6zC+//KJefvlldbyE8Q4Cu364FVkg9N+jjz6qpk6bpvbaa68CRf1wDdcfD9jtwdnOnz9f/8nHpK4fPty3CvBfAIQLigVwSV4OsrznfZs6rhk0SD33/PMxdyfkX7ZsmerWtav6+eefdXE41f2qVtW7CyFGinGB95t8cOLXpx132sl0I+o6SPr0wgsvFNhNRWVK04+nn35az3Gqe/KJJ9SZZ54ZWPN2Mq4K0mfmIxz2p59+qvr06aOelnlXt27dwHLmwddff63nGL9fe/VVvdOzd3iaiP4jRHThpLtNmYSu+9ZvGCGiX721WAjglITKm8xV6tXXRPRvIQp2X2q3ODYpIvr+M0+pb5a8ras/qvu5BYjoktmPq5UffWCaj3llsWh6ZnfV+LSz1Db5hCJmoZAMt99+uxIuMySHUsLh6edMeibAC/LRnHraaVFlhgwZEiGg7dq3V1dddZXaZ599dB4+EOEk1YknnRRVJpM/3nr7bbXLLrtoIvqNEG/av+fuu3Ufr7jiCjVxUvRi/dprr6mfRfzAGCGKjDmMiM6eM0dPfLa/q1auVM88+6y67dZb1WeffaYuOP989YSU30YWQBviwbWd39zXqlVLzf/f//TPzbJV/lSI0b333quES1IPPPCA2m+//dS5vXub7PrKIve2taBFPYzxo2/fvuqiiy+OkWvr46+++kqNHzdOXTFgQGh+tts9zzlHE1C2qgOuvFKdccYZEWLNIrZUFqIf1vgzFOdfcIHq169faBvm4YoVK9Qdd96pBg4caJIydoVwAsybt2XOff/995F57230pBNPVLfedptOZkGR3Zne6j8hdcRDRO3vkDk6R4hvr3PPjTQTPdsiye7Gi4Hfflit5t42Sk3v10f9/ecm7+O0/0YOaDi0vhddpOvnpduAbAsODzhRJso4+agMASWtcuXKSrZRvvIinmcSjEz2Svlo4VCA/wlBgjO1wXwMfc47T8uq/vfiiwpZVyyAUO4jsrwL5CMfM2aMzr506VL17DPPxCqa1HNk1Q0aNFD33HOPOu6443QdIkLRu4OkKkyikOwfdalq1arpKwSdxSMMpk6Zon744QedhZ1Ar169IgSURN5To0MP1fMnrJ5Yz0yf7pswQXN6sfKn8vyTjz9Wn3/+uZb3n3rqqboqQ+hi1XvggQeq2rVr62ws3vHAk/nMTNB3qDlRb0VtB1yrdqtWw5vs+7tcAGu/V+066vhLwldJu8K9DjjQ/pnV+/onnqwOObFDVJv/bBHB9q+/qNXLPlEfzn1WbfxlnX6+4p231DM3D1edho2Myp/uH6++8opiO169enV1nhAYPl6Rg6nVq1ervffeWze3QAiOgcsvv9zc5ty1h3BCcIMoCuCmDz/8cN1Hfs8TMQXQpUsX9dGHHyqRN6rnnntOc0v6QRz/2p98shp1yy1qpXAjz0v9/LYBomy2sqSjVIMrSwZQTlxw4YV6QWBr+PqiRar1CSdEqoLQ2W3xAM7cyx1HClg3GwQfdlmUGDtZ3xdtA2eddZZ68MEH1Xfffaeuvvpqze0H1f+8cM0AhP+II4/U94n82xSjT6auM2Q7PWvmTMXuQ2TK6smnnoprzKtlt/TRRx+ZaiJX6gkC6gbA+wnyJzJ33d7FwjDEgpWye/niiy90tspVqsTKrt59912NZ7bv7BSm3H+/+liI+JdSR63999flfYkoMsB9Dk5Nhrb9ThVU9cZNYnYyFzJUqrxPYF/rt22vmp97gZpx6flq1SdbX/aHz89RTc/opmWlyfb/QuGg7InfVeSa1157baQ6M1FOaNNGcw4tmjfXsjhkQXBfwLfyEQHIfGrK1jNXAW0uiwGT72vZ8hlAfoZMF+6gqsjoGCtEFK6CLWciIIovTUTZUnpBFGVRSQPlI79QCGGyQH8N2OMhDY1w48MOM4/19fXFiyMLX9QDzw/EHvwZOLlDBwW36wU02rfIooEc8P333tOihXNkofKDb/LxARedDMDt8megXbt2atz48eZn5Kr7NHq0Ol0WQ2T5cMD2ljeS0XMzceJExV+8gPjh6XwiCgE94ogjtKUCugParV+/foGqXhWR0dmy8GyUufbpJ59o0RjfzGke0ViBgpJgdn/Ht26trTkQNSGiQmxkxBZuO++HOU/a9hUqqs43jlalLFkoiqtUAJMYTFLM39/y2wByrHlz5+qfTBSgTdu2+mq2v/wwqootIqcJMpHRhXLgn+GiTJ/pktmCmTG2lolKvjfffFPLuBLptqlfKihQrFKlSlrhhNKJv50ClDUFCgYkRNqS5/a9yW7aMVejnDLPg65wnaYM151DzKKaCvFg4QVGC0FF/h0GzDc/gEt/VuTK/HnNyMiPjDeqT4LLIGCH0aNHD/0YEQtyyliA0s+u39zvvPPOvkUXCeePCR8mXhBQzJWMDN3+NuzCbNtff/11veCgW6CcWGfohd3O570HZ0Y81EYWeICFHoCQG/GKLyeqc7l/URioVGVfrZxa/sYinf7t0iVRzxP9MVm2BUGKJba8vGzgYuShQhgMkUUWhEyo7kEHafs+8kBA3xOO5FCRbYWBIS8Ix/3gb5HDAtv62NH55Y83DaUM8lugmnCkwI8//qj4IAC2po/J6g5AlJicTNJ4lSzk/0Q4DKB6tWr6av8bM3ZsIK7tfPHew1EbqOZpD6KzWMQuycD5ohiLd8zUz1b+pZde0hz4YNnFVBe7YC+g/EIbDTfW30fkg3hIzzEp+OKCBdq22K4DWXW8iiXKweUvkHoQNbCzEjM3u7oC9/3791c9e/YskA4TAT68YAglzEfLli3149/zZehzROl47eDBBexAm8surrdo41E8IkJqIeUaNmzorbrA75cFt8bW+CpRlqFQ/if/G2HRYrFv2rSpcpxoAdQFJ1SuWy/y8Peffozcp/vGbCGoFxkOsj6IjgGz1T/mmGNMklausNUJA2Mug2bXj3M15lR+htxh9cZ6dp9s19i2A23zOWqIpOkDht+MkT+TZsYYq26ePybmTYYTMxx7POWSyUP/xokGGoCLOuqoo5KpJqkyhvMxhZHtjho1Sv+EmL7tUdrxwCjB3luyRIkNpCmasSsmUYgagFcWLtRy/HQ1BgGcm79Dg0s0c8YoItEhoEvwwq677qqOPvpoLdPkGe8PbjYWsGU3wGJDe0ZJR7oh6I6IGizFcS0ncl4DW/I5RfM7XVdeEgokYIJoOt+UD8P8DbrmGp0OAYJgIthGRgW8IXI3NK8IvA2Qh7rQ9ANGK7lGzFlG3XyzPmVDOoThLpFzffvtt/yMy+xDZ4zxD2ItxuLqVuEEAbZdjRs31veGSLL9M+PjOidfu75cTl59KIomDfla6a0//vvPdpG6r8nHC3K/kzJkzgUOUTKcI/19VWwFATFq9z359V8P03vnJzqAiKNoAiL4sprFBAsiAqDwQcFnb9u9hNkqmvStLWrw61OyFUNAWYxRCsJZ2/PGbOltwudtp7fgooookxi/IfTePOY3HKhR3LKTsdsacdNNOhtKTHZZbjtvsBbH9c/f/zO9wQg/E/BUvqyFiX+cEB1b+dSxY0c1Ul4gqyhbYbYpvFC2ymzZWPmPlz+MveFSIMiIBSBM9erV01zgeCGWEKj77rtPKyQwg6I+NM0A21M/QsT2qbms5jawdcXg2wYUYHzsTHZDvHmOqMGYImGWg4Af6Nipk+LooQHuDxJRBVtmVvqDDz5YizPM89PEpAUZI+Oy5Xz7y4JyjyhAbHyZMtjObp9/dNKkIeO6xlLkmXT7Cl7rS/sQGj48m+B0E3mkn+IEebYXT9Q5VmxZzQJit2Hf804efvhhO0kxLkQ/YcDiirG84cbtvByZnCSmTdiK/vrrr+oOIaJYesgJHVVOiJERs9hl7PvJYtcLp29DjZo11dSpU+2kAve2qKHAwyQTDOeHGMw25aM65qw+CBFwKIU8mKldLeZ2l4gtLmPiHfoposiLjBixF5y1nBSLOkDQXuyxYQ74ZjDJc5woGIsTjHae7JyuygSYiYIA20sQIDCNGjXSzZp8bL2RJ/YULtScscZMBhMRCA3EmG0QgMkMAnVWbQgdhBFu0RDQliIeePiRRyL16ELWP+Rc9p+f4gAiQn0QUPrPSakbbrhBzRQzFHP23vQdrsBPa2yE+Mi4bEJMVxgL9RsCWqdOHXWVyAafEqsFlBJ+gGLB7jf3iBBiAUSTthgT93BAzVu0UFNE8zxcxhQE3rb4bcQZQWVIZ1vqLbs638YzrByL2c3523q/fOCYxa6zGJljKfGXzAtsat+UXQq44T2xyFKPF/z69INsbWOBLWqIlTee57YM3e/wyLGtWmnDe+Y8OoUgYOd2WL71xPXDhgVli2zVW4lpmNefBQrAI/PFOMzlUjI58jbLJBnVqmmkwl6TZiRl4jR/3NjIiaXqhzVR3e6aHKkz3psN69aqsW1bRLKf9+BjKhkb0mkX9ow6sdTqov6ROrmZfO5ZkRNLLXr3VS369I167v3xy8rv1fjTTlJ5+XLH1pcOVE3P6u7NFvO3kalUEu2j9+VQ2DyHOHqPAPIcRwmYasCN2UfPeMbW4iPh4H4S+SlCcFZrzHG8xJi8cJ9wrxCJHYVrrSPH32yOkDwAH78Rrm9N+e8/9VIGYgyXYwMfK1yQn2YaAgYRhKijOfcC4/glvz7GiFjCS/TQyqLB9XMmYeozuDS/7SsEMdCxRT6O7fw7SH5koH5bavKF4YnnjJPx+gHvIkiezfjMMVIWRxaVoLmBwnGz4BW7VEzG/IB3hV3mWqkLBeKee+6pDhBu1zsXTVt+ddh9Qg4JxxbYp3yTIsbPLscAc/QfESMFlbPnAIsjxJFjwQD48HvvBo98N9TLnGS83nfNroJ5DzB/mccQacRazCmYEUReLJwQTBYEL1CeeraVsr7b+Xm336KCjOhNZbvXqKWO7zfA/CxwXSFHLke23Cr/KvDQSugy6g5Vs8mRVkru3f4jk2T2jUMiBLSMfFAHt90qi0y0t0Hckqkn1nO4OcPRmTLmyocQS0Nv8kKcvETYPLOvTED+woBJF6vfdnk+8jBgHHZ9EGL7d1hZ+1kyZSgfhmO7fvs+HjzZ+e37eN4D+Q0xtcva91hsxALeleHEwvLGasuUNfJW89t7DTpWubsQrzDwzgEWoFjv04vHIDMpiKKXMHoZCBaXMIC48gf4EtHvP1oaVl4/2yxUOAzyhKpzDj4W/JvP2cXKV1jPObH0wtiR6rsP3ot0oVmPPqp8pXBCEMnsbhwGHAaKNQY0Ec3TnvTSO87thKXeZZ/YcsNyPnKY9PYkdm1vPTJDYf+Jsoh+g4/Nwq7//M3XEU9Qppa6rVqrZuf0MT/d1WHAYaCEY0AT0VKqVBQazrztHlW5TvjWIJZBdpW6ByclE43qSJZ+bFr/m4rFfXNaqVmP3qpln4sC5WJZ6q5rxmHAYSCHMOC7nefce0narpYpt71C7vnvP1vtKe33g+y31hHN1OFdzla4w3PgMOAw4DBgY8CXiNoZCvse2WoykEi5I7v2VM17X6jgSPHchCE9DpnLV9pVb++Tad+VcRhwGCgZGMg5Ilq2fLSt2t+b/0zqTdgnirx1+lWI6coOFXfWf37PXZrDgMOAw4AfBnLO2L60mDJgQmRgQ74fT/M73uuGX9ZGsu7gY4sYeehuHAYcBhwGUsBAznGijGX3ajXVqk+3+u4kZlOdlsclNERiJ+GJ3sBuVWuY20K/clwTI/ggY2i/55wuwbiX45B+xukYPGM4z8kU+9QJ8Ye2k/AJ2A9ipM05eo55EgXUDzDk/0DOq3N6BRu77+WkzYoA57h45zG2fpx6so8clsauUxxH28bVdnsYMhMWhH5jFE2/bR+d5MUXwB5iq2f3FeN5opkaYPeAwTxhPMxpLfPM7zpRjlV+Lr4F8C7kZ0OIgTnhIwxQP/2jD167QpOHEz34sfQDbA05thkEvDNwYAC8YZ+Jb9ggo36Tlyt9XSKORTgMge0tx0ptfNl5/e6JMmBOfvEcW9eacqTTDzemPIb6hBNZL/Nxd+lrE/FiRLRQP+BIKY47DGAgj0Nxe+57cWDy2lcM2r0Opek3YUHwaoaxP6ff8CMQ1ne7znTe5yQRrXXU0REi+t6cJxQyy/K7bHWiEM/gFz0wKZKt7I47qX3qFXTUGsmQ5ZuuZ5+t3Z0NCIiNQ9A54toYh690b9iwYdoXIh55vDGKeP6uxPTBbZg5I08aMFTiL1WSj2uaBHBjAkvkRX2CiLhEfjBDPJPjoMJ4A+Ic/9h85yHe/MSs6STn3oEZM2ao++V8thfw24hDX/twAPlGjhypiTongP6Qj5GTIfgBmCZxiwzgzf+UU05Rw66/3iRpt28EaPMCDnY534wbtiAjaRYhHElw5LGBHEU9WwK3eYGPsq+Ps2Zwh69Tzl17CQYOXwgC5wfEwxotjoqDAIfLxGvyAgTxHHEPh/NtYgh5YcWKFQrXbBBBCB+ElxM3nOrBsxeencwC5y1r/+4nnuBZ0Gzg9A5e6iVyaNRpM/wd0CYLOgsKfaQsxIzQNPhw8BKwaXK+HjeHXsCj0l3ifBpj9TsFBzjnDgPG+IksOAZwk3ednF1nUaUfPGcRLyOn5M4Tl3eXiXs9v5Nypny6rzm3nWeAh5x4sto2f/L8KSv9o4P6q7/yz3/HQgDOkt9+dFYk2yEndVDbyEdQ1IGTIcQo4jx5stDl9NO1Yw+b27Lr4hwwxMImejwnQNtXEvHQ/jME1JSHEzTPyY/7Oz5yPkYDcGyco8dZxPvyMfL7Q+FsINzHxAjaZ+rgigcf2oKrXCiuz/BJifMJCCmciR8QcgQCCndsYub45SPtGXHYQv1fCNf7muwccFQiYZt1/XD3fiChfyPjN3gII6CmDrhokx+8EYgP3wZ4p2LB5diiDXD9ncQRDRwci8JSweEr4lUKfLLgLBLnw7wbjlXGA53FoQvtL5d63xWuliiYM2VRJBCfgY/lHXWWBY1z/ePvuivS5hLxYcvChXclzuWzk/ECnLw9PpyfsNu6acQInfVu+W2ec2Xh4fSQnfZxvrMaCkhoas0w7CLfA57A6DPvCB+uBKTDwc6l/fpFOYvx9indv3OSiOIAmciaBr59f4maen53He/IpHmvhFpecM8d6unh10Ye4ZG+xbkXRn4X5Ru2u63EyQLeY8wZ4kTHgxcoOBvCC3uBbRqxZ0zgL/s5W0vvn/3c3Js8EFQIwZES0wfiZgCHFwCcgjm3DlfTQfp1jngZihfglmiLseB0GI4NH6yILOBE7C2qqRPCidgA57wmbo555r2accCBsk3EaTDEjTPZ1I9DEi+YMvbVmyfotykD3vBahTMRCKRehCxOHI4d7hEnLA8IhwfhMpwqZXEreK8QJXzQXplAxE3aB6dwdYabNxwyOL1YvB7R9kNCwFgADZfHuyPMCrsVfNEOHTrUd4j2+PD32lg84MMQAOaZffVLJw2nOoPF6TIuHfF2xfs0wA6E3dFZsvDggUli0ZtHGb9mjIji+f329sfF/cfRShuOuaCfqnPs8ZGkHz5fpib26KIeuKiXWvTAZPXRvOfUspdfVEueekzNHjFU3dnxBPXa1ImR/Jw84tDA9iIzKy6A5yC2bGHeZ8LGijwV4kZYDq93e4gQZ5PZaqULIEJb41NurdFshY2fxnS1Qz34Vh0mYg8+NPyt2gBRgYDjdg8iQL9icaN2ee6RY98sPlhxckFAtkwDBJJFD2Jh5IpwcLgIxG+sTUDsvsDRw4kuFNd4sSKB2uXMPcSU92QcHbPVRnxAELgg+S6LIO7piE5r+mrq87uWYV4IUU4UEAUxb28S0YEt+7frIYQ44o0JVlwo+3km7jNGRDFeX//jD3H/bfwt2gsQXohOueEWdZScEiolL9YA0TZfvOs29cSQgeqRqy5Vz4wUeeHsJ7R9p8lDpNIe90zNKVmo6VsqV5RC+MYkWB0hGJIBomriEeklqzzcxmyp8xT5cPmI0gEQLjgp2/s+0RmbNWumt/SDxH2d1zNTqu3iIg1uke2tDRBVuCfcoCG3a9myZSS+k50v1j1hJdhqGqfMsfKn+px3hWchiCewSMQYQCfZWoeBeU5coUQBIoXIYL98L1BEMwWQT4cBCxSEEbFHGOD4m50AbhcTBUQV7DwI8RwEcORET0D04OeqMahcKulaWAiRqrh35Ug9mBklA2yf7XoSqcPvhNS2pcuoVn0vUwcd10Zzn58smBfxpORXNyeKOFnU+NQzVawx7LTbHpG+xvJY5ddWYaURnEzHkhFHvPNkq2w8ycTbn6NFgQPHyZbehPo14Zn9tvLUq6NXypbPwECJJW+CpJk0Pj62WsAvQqQhNHyIcAYGINCTxMnvLaJsIRok2y5iCrFdxnVeqsDWlu33Go8PTrydM27jcYgommyLkW8GcXRBfakmMlU/F3uE0y2dL8enLIT6zvwwIkF1xUqPxKPKV/6YdulDGJg4Uz96lEZhZXiGa7fhIj5gcTszP9oqbaK4CVLYmTprVK+ub72KKmS6Zl6sE0sEorliAWBHtjV1xLryXsMIqClv8EZfvM6bTZ50XjURZet76VPzUq6Xs+X8pRvwJ0q0zRPXD1aICX74bJlwnuvU37K11YR7z73Uvoc0UnvtXzuKaw3rR5dRt4c9ztlnyI5GyrYSjSjbyxH5Avp4Owwhg1gi4Mekh60PW3lc6FXP/xC8dSFGYAtsAFMrL8CFoCHmw0WDi7d15Jxe4oiLMwgrIZGJUnmzaOrZBuLo2OuOzNtGPL/Rwtsu05D14smfsB7G5ymu4HCthlgjUSLKNhduxws9ZKw2Xuw+ePPG+5uxAKY9g0tkvubery58igJheUw5lIkLJRoCW2xEFcwP3g9cN4DlAwsk75e5FwSImYCynsUQTpp5gSNnQoUMEbkpHuWNLDeoPr90xuMn7/bmNXniGb+3bDK///sykimd5TIQzNpHH6P/stx02ppjIgY54KURJl3YZCUPKzkaSLS/hCpIFE6TbSJaTD4gzFlQ/lxvadG99bE9CnIobPIyYQlvAedxgmj4kXui6Q0C5Gto8FFgwMWhlLDjmweVC0tHHgrn1EE4TQMm6B+EweaKeQ5HP1jS7QXClPO7wqGh/feLV95ElCXJbFH92jFpJtYW3vsB7EcBOGjEIkFgLAjisRlt0qSJlqFuEGUVnDrKQMOxUz91zBVxD7JY7IeDAOsAwNsm4hXmBcqw1iKPZ14g000GqJt+QNTDiDDjZzGw7VGTaS/eMtvEm9HlSw8GMPkwQntvjSadPLHgPLELxRs9wce8ZjCxyiJXwkiaLT2mPxBtFC7pALgmOGTCyfrZCHrbIAxKS+F6jG2q93kivxEVACg6ABarOcLloiybIB+y/XeliCQgimwv44UpYsbEIof8L9MAoZg2daoWT6DNBkxYl8k+NrmmP4wZUQlbcEKZxILKIv5gUWXXwIJsE1DKmuisk2TBCwK40AfFJAp5M/PKDyCm2JIi350loWKSAfpCeORHJIRNEKAEe1E0/xjex/MdBdWTSLojoolgKw15OVUyX4Jpma2lXeWj+ZNDyyDtBz6aTLgnIhYiQJ8iH1uicLrYjGLShEkMooEgbWei9ZKf0yUoRYgoyqknA5hm+YUagZgZkyeTN9Er8tXpYvbDuAznhqwXkUV32cpj/2r/sQhhrRCvlh5iC6cM0SfWeCYBQn2dcOYoeJAnGk4ZccuZEtmTeOh3i72mH9wgOwpEGBcKd+93us2vTFgaJlco5FBmzpg+vUBW+soBCIgXce0h3kHQQog6lgPYiMajxffWw9jhLin/ntioeoFvil0NTMEAWSSzBUVqO58tpGSyHSbaqaIFRy6JgqaeTFJkX5wOwizjePnYCwjPA2RRbK/YMhNaOVFgVedDZfJjMhIGaLdtpQl5OfZ5UMj2DsXBAvnYsSbgJBTAiSmsCvgYkEUiZ+P0Cdsvgs3FCwQi21kIIFER1goBRoMMkUOOZ59wwhgbboTTUF6AMLWVxeNxicWO3ae9iLwofVwmJ3QgECjJ4KrpN3LjIAN6NMfr8mMAmbZ2kT7G2uKz+zAiB7hPYpsTipfjrZxaw67SBsQenM6hHxz5ZLHiiC1abxYR+oom/aKLLrKLpXSPDB65JgQd7Tvaf2TpK8Q4fqos4Lw/bFS7d//PtjuoQWSiyGAHCeGlbCIAgWa3QTvI1LkeK1p+dj8cOuD7YdEkqmxQFM9E2os3ryOi8WIqTfmIOcM2etiwYfrkhamWCcJRR4isF9gmESjNDzitw4kRiKEX4ESCODwmHts3iA9xwv1gBzGm5mPxIxxs/wwR3UmiE5DPBk49wRUhh+QEENwMRvYQJrabJgIpfYRYcGzVBraVO+bHsDHp4Ih2UIppkMWFtuHQCEmMXSWyMIBAZ0vElMZEuNSJnn/kh4hjSWDkvtT/gByT1SD1E8QPKwO2oshCDVdoqmJxoQxKKv5sAD9hRBQCDz4Igw1ghsX7IkIqRzcLLKaSB9k0hARbVY7z2njDR4J9HFdXGvKPUz8VPDj2y84C86BwoYgXMGKH6zdAX3kfXmLPc96fd17wvpn7nGQD98ZCxNRHW16RgnnGFXtgjjePHzdOPSFRbo2YAZk8YhvmWJjs1q4rXfelhBtI3Oo1Xa2X8HowwWA7DkHbX5QGaK6TAT5E/rwfeKy62EbDCVWLYTITq55En2OXiiaYbRemM7EUaYnWX5Ly8/44bglnbnPTmcQBIhl2T5xwQtZZWADpgvOEg4dYZ0sb7x3v/wFsy/ORf1iKRgAAAABJRU5ErkJggg==\" width=\"200\" align=\"right\" >   \n","<h1> <b>Computer Vision 2022/2023 </b> </h1>\n","<p><b>Andry Maykol Pinto </b> (amgp@fe.up.pt)</p>\n","<p>Department of Electrical and Computer Engineering</p>\n","</Body>  \n","FEUP\n","\n","---"]},{"cell_type":"markdown","metadata":{"id":"6VnUBWsOPf2E"},"source":["# **Image Filters**"]},{"cell_type":"markdown","metadata":{"id":"36GvUEzZsklu"},"source":["<h2>Table of contents</h2>\n","<div class=\"alert alert-block alert-info\" style=\"margin-top: 20px\">\n","<ol>\n","    <li><a href=\"#ref1\">An introduction to convolution kernels in OpenCV</a></li>\n","    <li><a href=\"#ref2\">Spatial linear filter: blurring</a></li>\n","    <li><a href=\"#ref3\">Spatial linear filter: Gaussian blur </a></li>\n","    <li><a href=\"#ref4\">Spatial non-linear filter: Median blur</a></li>\n","    <li><a href=\"#ref5\">Spatial non-linear filter: Bilateral filter</a></li>\n","    <li><a href=\"#ref6\">Similarity metrics</a></li>\n","    <li><a href=\"#ref7\">Sharpening filter</a></li>\n","    <li><a href=\"#ref8\">Exploring the Fourier transform</a></li>\n","    <li><a href=\"#ref9\">Low Pass Filter (LPF)</a></li>\n","    <li><a href=\"#ref10\">High Pass Filter (HPF)</a></li>\n","    <li><a href=\"#ref11\">Edge detectors</a></li>\n","</ol>\n","<p></p>\n","</div>\n","<hr>\n","\n"]},{"cell_type":"markdown","metadata":{"id":"pnCWYcbDgi0t"},"source":["<a id=\"ref1\"></a>\n","# An introduction to convolution kernels in OpenCV\n","\n","OpenCV provides a function to convolve a kernel with an image, named as `cv2.filter2D()`. It is possible to use 2D-convolution kernels and the OpenCV library to apply different blurring and sharpening techniques to an image.\n","\n","**Example**:\n","\n","\n","$$\n","kernel= \\frac{1}{25}\n","\\begin{bmatrix}\n","     1          & 1    & 1 & 1 & 1      \\\\\n","     1          & 1    & 1 & 1 & 1      \\\\\n","     1          & 1    & 1 & 1 & 1      \\\\\n","     1          & 1    & 1 & 1 & 1      \\\\\n","     1          & 1    & 1 & 1 & 1      \n","\\end{bmatrix}\n","\\\\\n","$$ "]},{"cell_type":"markdown","metadata":{"id":"XGGsLDOTySB3"},"source":["Filtering with the above kernel results in the following procedure: for a sliding 5x5 window being centered on each pixel, all pixels falling within this window are summed up and the result is then divided by 25. This division aims to normalize the result, keeping the same range of values. "]},{"cell_type":"code","execution_count":1,"metadata":{"id":"oVrvVX3nYSnE"},"outputs":[{"ename":"ModuleNotFoundError","evalue":"No module named 'google'","output_type":"error","traceback":["\u001b[0;31m---------------------------------------------------------------------------\u001b[0m","\u001b[0;31mModuleNotFoundError\u001b[0m                       Traceback (most recent call last)","Cell \u001b[0;32mIn[1], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m \u001b[39mfrom\u001b[39;00m \u001b[39mgoogle\u001b[39;00m\u001b[39m.\u001b[39;00m\u001b[39mcolab\u001b[39;00m \u001b[39mimport\u001b[39;00m files\n\u001b[1;32m      2\u001b[0m uploaded \u001b[39m=\u001b[39m files\u001b[39m.\u001b[39mupload()\n","\u001b[0;31mModuleNotFoundError\u001b[0m: No module named 'google'"]}],"source":["from google.colab import files\n","uploaded = files.upload()"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"W4_hGVbyuBZP"},"outputs":[],"source":["import cv2\n","import numpy as np\n","from matplotlib import pyplot as plt\n","\n","img = cv2.imread('img_example.jpg')\n","\n","if img is None:\n","  print('Could not read image')\n"," \n","# define the kernel\n","b_kernel = np.ones((5,5), np.float32)/25\n","\n","# converting from BGR to RGB for displaying in plt.\n","cv2.cvtColor(img, cv2.COLOR_BGR2RGB, img); \n","\n","# apply the convolution\n","dst = cv2.filter2D(src=img, ddepth=-1, kernel=b_kernel)\n","\n","# displaying original image and final\n","plt.figure(figsize=(18,19))\n","plt.subplot(121)\n","plt.imshow(img)\n","plt.title('Original')\n","\n","plt.subplot(122)\n","plt.imshow(dst)\n","plt.title('Averaging 5x5')"]},{"cell_type":"markdown","metadata":{"id":"a3baMGcWc64c"},"source":["Images also can be filtered with various low-pass filters (LPF), high-pass filters (HPF).\n","\n","*   A LPF removes noise, or blurring the image. \n","*   A HPF helps in finding edges in an image.\n","\n","In this context, **frequency refers to the rate of change of pixel values**. Sharp edges would be high frequency content because the pixel values change rapidly in that region.\n","\n"]},{"cell_type":"markdown","metadata":{"id":"CyZyoUAOyyIm"},"source":["<a id=\"ref2\"></a>\n","# Spatial linear filter: blurring\n","\n","aims to averaging the pixel values within a neighborhood. This is also called a low pass filter. A low pass filter is a filter that allows low frequencies and blocks higher frequencies.  The smoothening effect is controlled by the **size of the kernel**, depending on how much we want to smoothen the image. A bigger size leads to averaging over a larger area. "]},{"cell_type":"code","execution_count":null,"metadata":{"id":"FXl361yfc4xd"},"outputs":[],"source":["import cv2\n","import matplotlib.pyplot as plt\n","import numpy as np\n","\n","#load image\n","img = cv2.imread('img_example.jpg')\n","\n","# converting from BGR to RGB for displaying in plt.\n","cv2.cvtColor(img, cv2.COLOR_BGR2RGB, img); \n","\n","#define kernel\n","kernel_3x3 = np.ones((3,3), np.float32) / 9.0\n","#b_kernel = np.array([[0, 0, 0],\n","#                   [0, 1, 0],\n","#                   [0, 0, 0]])\n","\n","#kernel_5x5 = np.ones((5,5), np.float32) / 25.0\n","print(kernel_3x3)\n","\n","# apply kernel\n","dst = cv2.filter2D(src=img, ddepth=-1, kernel=kernel_3x3)\n","\n","# displaying original image and final\n","plt.figure(figsize=(18,19))\n","plt.subplot(121)\n","plt.imshow(img)\n","plt.title('Original')\n","\n","plt.subplot(122)\n","plt.imshow(dst)\n","plt.title('Blurring 3x3')\n"]},{"cell_type":"markdown","metadata":{"id":"J9kDtKWthxO_"},"source":["OpenCV has provided ` blur()` method to blur an image.\n","\n","The syntax of this function is `dest = blur(source, kernel_size)`"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"REAixFKy0ld-"},"outputs":[],"source":["import cv2\n","import numpy as np\n","from matplotlib import pyplot as plt\n","\n","#Load image\n","img = cv2.imread('img_example.jpg')\n","\n","# converting from BGR to RGB for displaying in plt.\n","cv2.cvtColor(img, cv2.COLOR_BGR2RGB, img); \n","\n","#Blurring\n","blur = cv2.blur(img,(7,7))\n","\n","# displaying original image and final\n","plt.figure(figsize=(18,19))\n","plt.subplot(121)\n","plt.imshow(img)\n","plt.title('Original')\n","\n","plt.subplot(122)\n","plt.imshow(blur)\n","plt.title('Blurring 7x7')"]},{"cell_type":"markdown","metadata":{"id":"tg1YDeSv2Ahd"},"source":["<a id=\"ref3\"></a>\n","# Spatial linear filter: Gaussian blur\n","\n","In this type, the image is blurred and the contrast it reduced.\n","\n","The syntax of this function is\n","`dst\t=\tcv2.GaussianBlur(src, kernel_size, sigmaX, sigmaY)`\n","\n","The final two arguments are sigmaX and sigmaY, meaning the Gaussian kernel standard deviations, in the X (horizontal) and Y (vertical) direction. The default setting of sigmaY is zero. If you simply  set sigmaX to zero, then the standard deviations are computed from the kernel size (width and height respectively). \n","Gaussian filtering is highly effective in removing Gaussian noise from the image."]},{"cell_type":"code","execution_count":null,"metadata":{"id":"n0V0RnbD2HjZ"},"outputs":[],"source":["import cv2\n","import numpy as np\n","from matplotlib import pyplot as plt\n","\n","#Load image\n","img = cv2.imread('img_example.jpg')\n","\n","# converting from BGR to RGB for displaying in plt.\n","cv2.cvtColor(img, cv2.COLOR_BGR2RGB, img); \n","\n","#Gaussian blur\n","blur = cv2.GaussianBlur(src=img,ksize=(5,5),sigmaX=0)\n","\n","plt.figure(figsize=(18,19))\n","plt.subplot(121)\n","plt.imshow(img)\n","plt.title('Original')\n","\n","plt.subplot(122)\n","plt.imshow(blur)\n","plt.title('Gaussian blur 5x5')"]},{"cell_type":"markdown","metadata":{"id":"A4qmvhXx24U9"},"source":["<a id=\"ref4\"></a>\n","# Spatial non-linear filter: Median blur\n","\n","It's a non-linear filter since the central element of an image is replaced by median of all pixels.\n","\n","The syntax of this function is `dest = cv2.medianBlur(src, kernel_size)`.\n","\n","This is highly effective in removing salt-and-pepper noise. One interesting thing to note is that, in the Gaussian and box filters, the filtered value for the central element can be a value which may not exist in the original image. However this is not the case in median filtering, since the central element is always replaced by some pixel value in the image.\n","\n","The kernel size must be **a positive odd integer**."]},{"cell_type":"code","execution_count":null,"metadata":{"id":"CM5raL2F38nn"},"outputs":[],"source":["import cv2\n","import numpy as np\n","import os\n","from matplotlib import pyplot as plt\n","\n","def noisy_sp(image):\n","      row,col,ch = image.shape\n","      s_vs_p = 0.5\n","      amount = 0.15\n","      out = np.copy(image)\n","      # Salt mode\n","      num_salt = np.ceil(amount * image.size * s_vs_p)\n","      coords = [np.random.randint(0, i - 1, int(num_salt))\n","              for i in image.shape]\n","      out[coords] = 1\n","\n","      # Pepper mode\n","      num_pepper = np.ceil(amount* image.size * (1. - s_vs_p))\n","      coords = [np.random.randint(0, i - 1, int(num_pepper))\n","              for i in image.shape]\n","      out[coords] = 0\n","      return out\n","\n","#Load image\n","img = cv2.imread('img_example.jpg')\n","\n","# converting from BGR to RGB for displaying in plt.\n","cv2.cvtColor(img, cv2.COLOR_BGR2RGB, img); \n","\n","# Add S%P noise\n","img_noisy = noisy_sp(img)\n","\n","#Median blur\n","median = cv2.medianBlur(src=img_noisy, ksize=5)\n","\n","# Plotting\n","plt.figure(figsize=(18,19))\n","plt.subplot(121)\n","plt.imshow(img_noisy)\n","plt.title('Original')\n","\n","plt.subplot(122)\n","plt.imshow(median)\n","plt.title('Median blur')"]},{"cell_type":"markdown","metadata":{"id":"6G4aapEFrZmc"},"source":["> *Discuss this outcome.*"]},{"cell_type":"markdown","metadata":{"id":"CmT6HGUU7sS2"},"source":["<a id=\"ref5\"></a>\n","# Spatial non-linear filter: Bilateral filter\n","\n","\n","The bilateral filter `cv2.bilateralFilter()` is highly effective at noise removal while preserving edges. The operation is slower compared to other filters. The last two arguments here specify the color and space neighborhood. \n","\n","The bilateral filter also uses a Gaussian filter in the space domain, but it also uses one more (multiplicative) Gaussian filter component which is a function of pixel intensity differences. \n","\n","The bilateral filter performs the averaging only at pixels that are similar to each other in intensity. It also takes a color neighborhood metric to see if it can replace the current pixel that is similar in intensity as well.\n"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"1ASry3Zm8k40"},"outputs":[],"source":["import cv2\n","import numpy as np\n","import os\n","from matplotlib import pyplot as plt\n","\n","def noisy_gauss(image, mean=0, var=0.001):\n","    image = np.array(image/255, dtype=float)\n","    noise = np.random.normal(mean, var ** 0.5, image.shape)\n","    out = image + noise\n","    if out.min() < 0:\n","        low_clip = -1.\n","    else:\n","        low_clip = 0.\n","    out = np.clip(out, low_clip, 1.0)\n","    out = np.uint8(out*255)\n","    return out\n","\n","\n","#Load image\n","img = cv2.imread('img_example.jpg')\n","\n","# converting from BGR to RGB for displaying in plt.\n","cv2.cvtColor(img, cv2.COLOR_BGR2RGB, img); \n","\n","# Add S&P noise\n","img_noisy = noisy_sp(img)\n","\n","#Bilateral filtering\n","bf = cv2.bilateralFilter(src=img_noisy,d=11,sigmaColor=175,sigmaSpace=175)\n","\n","# Plotting\n","plt.figure(figsize=(18,19))\n","plt.subplot(121)\n","plt.imshow(img_noisy)\n","plt.title('Original')\n","\n","plt.subplot(122)\n","plt.imshow(bf)\n","plt.title('Bilateral filtering')"]},{"cell_type":"markdown","metadata":{"id":"iiNQb_weqU27"},"source":["Change the parameters of the filter.\n","\n","> *Discuss the outcomes.*"]},{"cell_type":"markdown","metadata":{"id":"-Y3YbkJxJL1m"},"source":["<a id=\"ref6\"></a>\n","# Similarity metrics\n","\n","The similarity metrics will quantify the degree of visual and semantic similarity of a pair of images. In other words, they compare two images and return a value that tells you how visually similar they are. There multiple similarity measures (proposed by researchers), namely:\n","1. Mean Squared Error (MSE)\n","2. Root Mean Squared Error (RMSE)\n","3. Peak Signal-to-Noise Ratio (PSNR)\n","4. Structural Similarity Index (SSIM)\n","5. Universal Quality Image Index (UQI)\n","6. Multi-scale Structural Similarity Index (MS-SSIM)\n","7. Erreur Relative Globale Adimensionnelle de Synthèse (ERGAS)\n","8. Spatial Correlation Coefficient (SCC)\n","9. Relative Average Spectral Error (RASE)\n","10. Spectral Angle Mapper (SAM)\n","11. Visual Information Fidelity (VIF)\n","\n","\n","The MSE, RMSE, <a href=\"https://ieeexplore.ieee.org/abstract/document/1284395/\" class=\"cw kn\" target=\"_blank\" rel=\"noopener ugc nofollow\">PSNR</a> and the <a href=\"https://ieeexplore.ieee.org/abstract/document/1284395/\" class=\"cw kn\" target=\"_blank\" rel=\"noopener ugc nofollow\">SSIM</a> are the most frequently used. \n","\n","There is a library in python called **sewar** that implements some similarity measures. This library is not installed in colaboratory by default and so, we need to do it.\n","\n","```\n","!pip install sewar\n","```"]},{"cell_type":"markdown","metadata":{"id":"4PgHF2w4p_yQ"},"source":["----\n","**Exercise 1:** \n","\n","Measure the impact of the Average, Gaussian, Median and Bilateral filters for image corrupted by Gaussian and salt-and-pepper noise.\n","\n","a) Load the image \"img_example.jpg\" and show two images corrupted with  salt-and-pepper and gaussian noise (mean = 0, variance = 0.01), respectively.\n","\n","b) Determine the PSNR and SSIM for each noisy image, considering the original one. \n","\n","c) Apply the Average (3x3), Gaussian (3x3), Median(3x3) and Bilateral (d=3) filters to each noisy image.\n","\n","d) Determine the PSNR and SSIM for the result of each filter, in respect to the original image.\n","\n","e) Repeate points 2 to 4 for kernels with different sizes: 3x3, 5x5, 7x7 and 11x11. Show a graphical representation for the PSNR and SSIM for different kernels. \n","> application notes: see https://matplotlib.org/stable/tutorials/introductory/sample_plots.html\n"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"pZ5wrLzGgSax"},"outputs":[],"source":["!pip install sewar"]},{"cell_type":"markdown","metadata":{"id":"-s2DkWuXN4yX"},"source":["Solution a)."]},{"cell_type":"code","execution_count":null,"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"elapsed":7167,"status":"ok","timestamp":1634020806926,"user":{"displayName":"Andry Pinto","photoUrl":"https://lh3.googleusercontent.com/a/default-user=s64","userId":"12885272950202163621"},"user_tz":-60},"id":"lZB70RqmHCTi","outputId":"7e55ab4a-d3c3-4f5a-a74c-ef77b3fe7510"},"outputs":[{"name":"stdout","output_type":"stream","text":["SSIM:  (1.0, 1.0)\n"]}],"source":["import cv2\n","import numpy as np\n","import os\n","from matplotlib import pyplot as plt\n","from sewar.full_ref import mse, rmse, psnr, uqi, ssim, ergas, scc, rase, sam, msssim, vifp\n","\n","#def noisy_gauss... TODO\n","#def noisy_sp...  TODO\n","\n","#Load image\n","img = cv2.imread('img_example.jpg')\n","\n","# converting from BGR to RGB for displaying in plt.\n","cv2.cvtColor(img, cv2.COLOR_BGR2RGB, img); \n","\n","## TODO\n","#print(\"PSNR: \", psnr(blur, org))\n","#print(\"SSIM: \", ssim(blur, org))"]},{"cell_type":"markdown","metadata":{"id":"Mu075tddOywA"},"source":["Solution b)"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"wYvmdTISO0cF"},"outputs":[],"source":["#TODO"]},{"cell_type":"markdown","metadata":{"id":"k3qpzKFgO4uz"},"source":["Solution c)"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"7a4lCKlNO2pA"},"outputs":[],"source":["#TODO"]},{"cell_type":"markdown","metadata":{"id":"eJLgdyUxO1dX"},"source":["Solution d)"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"_CoFM2k9O9oS"},"outputs":[],"source":["#TODO"]},{"cell_type":"markdown","metadata":{"id":"lN2SlYMIn7Bp"},"source":["\n","_________\n"]},{"cell_type":"markdown","metadata":{"id":"3Sim7bmGsfBw"},"source":["<a id=\"ref7\"></a>\n","# Sharpening filter\n","\n","A sharpening filter will sharpen the edges (high frequency components) in the image.\n","\n","**Example:**\n","\n","$$\n","kernel= \n","\\begin{bmatrix}\n","     -1          & -1    & -1       \\\\\n","     -1          & 8     & -1       \\\\\n","     -1          & -1    & -1         \n","\\end{bmatrix}\n","\\\\\n","$$ \n","\n"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"QflWyDlLRPh8"},"outputs":[],"source":["import cv2\n","import numpy as np\n","from matplotlib import pyplot as plt\n","\n","#Load image\n","img = cv2.imread('img_example.jpg')\n","\n","# converting from BGR to RGB for displaying in plt.\n","cv2.cvtColor(img, cv2.COLOR_BGR2RGB, img); \n","\n","#Define HF kernels\n","kernel_sharpen_1 = np.array([[0,-1,0],\n","                             [-1, 4,-1], \n","                             [0,-1,0]])\n","\n","kernel_sharpen_2 = np.array([[-1,-1,-1],\n","                             [-1, 8,-1], \n","                             [-1,-1,-1]])\n","\n","\n","kernel_sharpen_3 = np.array([[-1,-1,-1,-1,-1], \n","                             [-1,-2,-2,-2,-1],\n","                             [-1,-2, 32,-2,-1],\n","                             [-1,-2,-2,-2,-1],\n","                             [-1,-1,-1,-1,-1]])\n","\n","# Applying different kernels to the input\n","output_1 = cv2.filter2D(img, -1, kernel_sharpen_1)\n","output_2 = cv2.filter2D(img, -1, kernel_sharpen_2)\n","output_3 = cv2.filter2D(img, -1, kernel_sharpen_3)\n","\n","#Plots\n","plt.figure(figsize=(18,19))\n","plt.subplot(121)\n","plt.imshow(img)\n","plt.title('Original')\n","\n","plt.subplot(122)\n","plt.imshow(output_1)\n","plt.title('HFC (-1, 4, -1)')\n","\n","plt.figure(figsize=(18,19))\n","plt.subplot(121)\n","plt.imshow(img)\n","plt.title('Original')\n","\n","plt.subplot(122)\n","plt.imshow(output_2)\n","plt.title('HFC (-1, 8, -1)')\n","\n","\n","\n","plt.figure(figsize=(18,19))\n","plt.subplot(121)\n","plt.imshow(img)\n","plt.title('Original')\n","\n","plt.subplot(122)\n","plt.imshow(output_3)\n","plt.title('HFC (-1,-2, 32,-2,-1)')\n"]},{"cell_type":"markdown","metadata":{"id":"fuutiI9IezwD"},"source":["Let's add the high frequency components:"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"n2bcJAJ-eyKT"},"outputs":[],"source":["hfc_1 = img + output_1\n","hfc_2 = img + output_2\n","hfc_3 = img + output_3\n","\n","#Plots\n","plt.figure(figsize=(18,19))\n","plt.subplot(121)\n","plt.imshow(img)\n","plt.title('Original')\n","\n","plt.subplot(122)\n","plt.imshow(hfc_1)\n","plt.title('Sharpen (-1, 4, -1)')\n","\n","\n","plt.figure(figsize=(18,19))\n","plt.subplot(121)\n","plt.imshow(img)\n","plt.title('Original')\n","\n","plt.subplot(122)\n","plt.imshow(hfc_2)\n","plt.title('Sharpen (-1, 8, -1)')\n","\n","\n","plt.figure(figsize=(18,19))\n","plt.subplot(121)\n","plt.imshow(img)\n","plt.title('Original')\n","\n","plt.subplot(122)\n","plt.imshow(hfc_3)\n","plt.title('Sharpen (-1,-2, 32,-2,-1)')\n"]},{"cell_type":"markdown","metadata":{"id":"1MJKCcZcT326"},"source":["> *Discuss this outcome.*"]},{"cell_type":"markdown","metadata":{"id":"JLuoNmkO2CMR"},"source":["____________"]},{"cell_type":"markdown","metadata":{"id":"kyFN6e7iiSY0"},"source":["<a id=\"ref8\"></a>\n","# Exploring the Fourier transform\n","\n","This concept is useful when manipulating images because it allows us to identify regions in images where a signal (such as the values of image pixels) changes a lot, and also regions where the change is less dramatic.\n","\n","The waveforms you observe all around you are the sum of other waveforms. \n","\n","The **magnitude spectrum** of an image is another image that provides a representation of the original image in terms of its changes.\n","\n","This forms the foundation of the aforementioned processing operations: HPFs and LPFs.\n","\n","\n","*   `cv2.dft()` performs a forward or inverse Discrete Fourier transform of a 1D  or 2D floating-point array\n","*    `cv2.idft()` calculates the inverse Discrete Fourier Transform of a 1D or 2D array.\n","\n","\n","\n","\n","\n"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"BWfinKmSjTy2"},"outputs":[],"source":["import cv2\n","import numpy as np\n","from matplotlib import pyplot as plt\n","\n","#Load image\n","img = cv2.imread('img_example.jpg', cv2.IMREAD_GRAYSCALE)\n","\n","# converting from BGR to RGB for displaying in plt.\n","cv2.cvtColor(img, cv2.COLOR_BGR2RGB, img); \n","\n","# Convert to  floating\n","img_float32 = np.float32(img)\n","\n","# Foward FT (the result, though being a complex array, has complex-conjugate symmetry)\n","dft = cv2.dft(img_float32, flags = cv2.DFT_COMPLEX_OUTPUT) \n","\n","#Shift the zero-frequency component to the center of the spectrum.\n","dft_shift = np.fft.fftshift(dft)\n","\n","# Determine the magnitude spectrum.\n","magnitude_spectrum = 20*np.log(cv2.magnitude(dft_shift[:,:,0],dft_shift[:,:,1]))\n","\n","# Plot\n","plt.figure(figsize=(18,19))\n","plt.subplot(121)\n","plt.imshow(img, cmap = 'gray')\n","plt.title('Original')\n","\n","plt.subplot(122)\n","plt.imshow(magnitude_spectrum, cmap = 'gray')\n","plt.title('Magnitude Spectrum')\n","\t"]},{"cell_type":"markdown","metadata":{"id":"nY3IX7c_mQpA"},"source":["<a id=\"ref9\"></a>\n","## Low Pass Filter (LPF)\n","\n","Create a mask with high value (1) at low frequencies and 0 at HF region."]},{"cell_type":"code","execution_count":null,"metadata":{"id":"c2GKVJZPmP-5"},"outputs":[],"source":["rows, cols = img.shape\n","crow, ccol = int(rows/2) , int(cols/2)     # center\n","\n","# create a mask first, center square is 1, remaining all zeros\n","mask = np.zeros((rows, cols, 2), np.uint8) \n","# mask.shape is (rows, cols, 2) which matches the returned np.ndarray from cv2.dft().\n","\n","mask[crow-30 : crow+30, ccol-30 : ccol+30] = 1  \n","\n","\n","# apply mask and inverse DFT\n","fshift = dft_shift*mask\n","f_ishift = np.fft.ifftshift(fshift)\n","img_back = cv2.idft(f_ishift)\n","img_back = cv2.magnitude(img_back[:,:,0], img_back[:,:,1])\n","\n","\n","#Plot\n","plt.figure(figsize=(18,19))\n","plt.subplot(121)\n","plt.imshow(img, cmap = 'gray')\n","plt.title('Original')\n","\n","plt.subplot(122)\n","plt.imshow(img_back, cmap = 'gray')\n","plt.title('Magnitude Spectrum')\n","plt.show()     \n","\n"]},{"cell_type":"markdown","metadata":{"id":"mIJ9vTMeovY2"},"source":["<a id=\"ref10\"></a>\n","## High Pass Filter (HPF)"]},{"cell_type":"markdown","metadata":{"id":"5XXORMDvom5C"},"source":["----\n","**Exercise 2:** \n","\n","Implement a HPF using the Fourier Transform based on the example above."]},{"cell_type":"code","execution_count":null,"metadata":{"id":"s0NE7cRDlCKM"},"outputs":[],"source":["#TODO"]},{"cell_type":"markdown","metadata":{"id":"cGwHFCCSaNRW"},"source":["> *Discuss this outcome: the difference between LPF and HPF*"]},{"cell_type":"markdown","metadata":{"id":"PrSYocp0aqdu"},"source":["------------------"]},{"cell_type":"markdown","metadata":{"id":"Egig1Xxnq944"},"source":["<a id=\"ref11\"></a>\n","# Edge detectors\n","\n","It is used to identify the boundaries (edges) of objects or regions within an image. Edges are characterized by sudden changes in pixel intensity and so, the there is a need to analyse the neighboring pixels.\n","\n","In computer vision, the most well known edge detectors are:\n","\n","*   Sobel Edge Detection\n","*   Canny Edge Detection\n","\n"]},{"cell_type":"markdown","metadata":{"id":"QyPxCJozuT7q"},"source":["The Sobel Operators are the following:\n","\n","\n","$$\n","Sobel_X = \n","\\begin{bmatrix}\n","     -1          & 0    & 1       \\\\\n","     -2          & 0    & 2       \\\\\n","     -1          & 0    & 1         \n","\\end{bmatrix}\n","\\\\\n","$$ \n","\n","$$\n","Sobel_Y = \n","\\begin{bmatrix}\n","     1         & 2    & 1       \\\\\n","     0         & 0     & 0       \\\\\n","     -1        & -2    & -1         \n","\\end{bmatrix}\n","\\\\\n","$$ \n"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"OhCmYjkJriob"},"outputs":[],"source":["import cv2\n","import numpy as np\n","from matplotlib import pyplot as plt\n","\n","#Load image\n","img = cv2.imread('img_example.jpg')\n","\n","# converting from BGR to RGB for displaying in plt.\n","cv2.cvtColor(img, cv2.COLOR_BGR2RGB, img); \n","\n","# Convert to graycsale\n","img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)\n","\n","# Blur the image for better edge detection\n","img_blur = cv2.GaussianBlur(img_gray, (3,3), 0) \n","\n","# Sobel Edge Detection\n","sobelx = cv2.Sobel(src=img_blur, ddepth=cv2.CV_64F, dx=1, dy=0, ksize=5) # Sobel Edge Detection on the X axis\n","sobely = cv2.Sobel(src=img_blur, ddepth=cv2.CV_64F, dx=0, dy=1, ksize=5) # Sobel Edge Detection on the Y axis\n","sobelxy = cv2.Sobel(src=img_blur, ddepth=cv2.CV_64F, dx=1, dy=1, ksize=5) # Combined X and Y Sobel Edge Detection\n","\n","\n","# Display Sobel Edge Detection Images\n","plt.figure(figsize=(18,19))\n","plt.subplot(221)\n","plt.imshow(img)\n","plt.title('Original')\n","\n","plt.subplot(222)\n","plt.imshow(sobelx, cmap = 'gray')\n","plt.title('Sobel X')\n","\n","plt.subplot(223)\n","plt.imshow(sobely, cmap = 'gray')\n","plt.title('Sobel Y')\n","\n","plt.subplot(224)\n","plt.imshow(sobelxy, cmap = 'gray')\n","plt.title('Sobel XY')\n","plt.show()    "]},{"cell_type":"markdown","metadata":{"id":"2VXlbxfHuh4n"},"source":["----\n","**Exercise 3:** \n","\n","Use the Canny Edge detector to obtain a binary image capturing the contours.\n","\n","This detector is formed by 4 steps:\n","1. Noise Reduction\n","2. Calculating Intensity Gradient of the Image\n","3. Suppression of False Edges\n","4. Hysteresis Thresholding\n","\n","> application notes: search for Canny Edge detector. See the https://docs.opencv.org/4.5.3/da/d5c/tutorial_canny_detector.html to be inspired. \n","> application notes: search for Otsu’s Binarization (`cv2.threshold`)\n"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"uJhS06CasQc4"},"outputs":[],"source":["#TODO"]},{"cell_type":"markdown","metadata":{"id":"kMywrnS9wcmN"},"source":["______________________________"]},{"cell_type":"markdown","metadata":{"id":"LrF2bdLJRSlK"},"source":["# References\n","Learning OpenCV 4 Computer Vision with Python 3, 3rd edition, Joseph and Joe, Packt"]}],"metadata":{"colab":{"collapsed_sections":[],"provenance":[]},"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.16"},"vscode":{"interpreter":{"hash":"9ac03a0a6051494cc606d484d27d20fce22fb7b4d169f583271e11d5ba46a56e"}}},"nbformat":4,"nbformat_minor":0}