zadanieanalizalab5/Lab_5.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c36cf7d5-c332-4482-af0e-3ba30ea9bb8e",
   "metadata": {},
   "source": [
    "<b>Uwaga:</b> Należy rozwiązać oba poniższe zadania. Na zaliczenie potrzeba co najmniej 3 pkt. \n",
    "<p></p>\n",
    "\n",
    "1. <b>Zadanie 1.</b>(2x1.5 pkt): Zreplikować wyniki  związane z  materiałem na stronie <a href=\"https://thepythoncodingbook.com/2021/08/30/2d-fourier-transform-in-python-and-fourier-synthesis-of-images/\">2d fourier transform in python and fourier synthesis of images</a>. Istotą tego zadania jest poprawne zrozumienie działania transformaty Fouriera na obrazach  - jego zaliczenie będzie polegało na odpowiedzi na dwa pytania teoretyczne.\n",
    "\n",
    "2. <b>Zadanie 2.</b>(2x1 pkt): Napisać kod dwóch filtrów opartych o transformatę Fouriera: dolno i górnoprzepustowego. W pierwszym przypadku przykładowy obraz dobrej jakości zaburzyć kolejno trzema różnymi szumami dostępnyhmi w programie Gimp. W każdym przypadku zastosować filtr dolnoprzepustowy możliwie dokładnie przywracający obraz przed zaszumieniem. Następnie wygenerować filtr górnoprzepustowy wizualizujący w trzech powyższych przypadkach postać dodanego szumu.  Materiał referencyjny znajduje się m.in. na stronie\n",
    "<a href=\"https://fairyonice.github.io/Low-and-High-pass-filtering-experiments.html\">Low and High pass filtering experiments.html</a>. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03b53c2c",
   "metadata": {},
   "source": [
    "## Zadanie 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "9cd51714",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgMAAAClCAYAAADBAf6NAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA5MklEQVR4nO1dSXNkV9W8rcEaS7OHIBzYH4NxgA3G9pIdv4C/w54FWzZs+BEs2PIP7PYUGCKw8cI2BqyhVF0aWt2WvgWcclbq3PvOfVV66lZlRnSEhCTfeu/dSk7lOZn33tXV1VUSBEEQBGFmMXfbL0AQBEEQhNuFigFBEARBmHGoGBAEQRCEGYeKAUEQBEGYcagYEARBEIQZh4oBQRAEQZhxqBgQBEEQhBmHigFBEARBmHEsRH/x+9//fvrqq6/S+fl56vV66eWXX04/+clP0ttvv53eeOON9Morr6Tnn38+zc3NpaOjo/TJJ5+kDz/8ML377rvpgw8+SJ9++mk6PDxM8/Pz6bnnnkuvvPJKeuONN9Lbb7+dXn/99fTSSy+lXq+Xzs/P0z//+c/08ccfp/v376f79++njz/+OH355Zfp/Pw8ra+vp5deeim99tpr6a233ko///nP049+9KP03HPPpYWFhXR0dJQ+/fTT0drvv//+aO25ubn07LPPjtZ+66230uuvv55efvnltLGxkR4+fDha+/3330/vvvtu+vjjj9MXX3yRzs7O0traWvrud797be0XXnghLSwspH6/n/7xj3+kDz/8ML3zzjvpww8/TJ988kna399P9+7dS88++2z64Q9/OLb2//3f/6WNjY306NGj9NVXX6W//vWv6b333kv3799Pf/nLX9Lnn3+eTk9P09raWnrxxRfTa6+9lt5888305ptvpldffTW98MILaXFxMR0fH6fPPvts7Lpt7ZRS2t3dvbb29773vbS5uZkePXqU/vWvf6W//e1vo7U/+uij9MUXX6STk5O0urqaXnzxxfTjH/84vfnmm+mtt95Kr776avrOd76TFhcX02AwSJ999ln66KOP0v3799N7772X/v73v6evv/46XV1dpd3d3fSDH/wg/fSnP02///3vb2wz53Dv3r3O18ytPT8/P/p6bm68Fn/8+PHY95gHdpeywUr3JKXxa728vMz+rGvcxtriXfFuV7wbVgZOTk7So0eP0vz8fFpdXU0bGxtpd3c37ezspK2trbS2tpbm5+fT48eP08nJSTo+Pk6Hh4fp8PAwHR8fp/Pz83R1dZWeeeaZ1Ov10tbWVtrb20vb29tpfX09LS0tpaurq/Tw4cM0GAzS0dFROjw8TP1+f7T23NxcWllZSZubm2l7ezvt7Oykzc3NtLq6mhYWFty1Hzx4kM7Pz9Pl5WVaXFxM6+vraWtra/S6e71eWlpaSiml9PDhw/TgwYPU7/fTwcFB6vf7aTgcjq29sbGRdnZ2Rmuvra2lhYWF9M0336TT09PR2kdHR+n4+DidnZ2N1l5bWxutvb29nTY2NtLS0lK6d+/eaG277sPDwzQcDtPFxUWam5tLy8vLaXNzM7s2XvfR0VEaDAbp7OwsffPNN9eum9e+uLgYW/vo6CidnJyki4uLdO/evbS8vJx6vV7a3d1Nu7u7Y2tfXl6OXffh4eG1tVdXV0fPTBCEOMS74t2ueDdcDJyeno4WWVtbS9vb26NNubGxkZaXl8du7uHhYTo4OEhHR0ejm2sXyBt6fX3dfbC2oXFtu7l7e3vu2rahbe0HDx6khw8fjq1tr93WXlxcdNfu9/ujtRcWFkbXbf9ym6q09s7OzmjttbU1d+2Dg4N0fHycTk5Oxta2TWWbcnl5Oc3NzaWHDx+mk5OTsQ09GAzS+fl5unfvXlpaWkq9Xm+0thHBM888ky4vL9PZ2Vnq9/vp8PAw7e/vj4jg8ePHY9dta29sbKSVlZU0NzeXLi4u0nA4HP09rp1SSktLS2lzczPt7u6mvb29FnT4dOHevXtj/6b1u3cZug95iHfFu13xbrgYsE2FF4jV6eLi4tgFHhwcjG7u6enp2AVubW2NNuXm5uboAh89epSGw6H7YO0Ce73e6AbZpsKba28IrLRya29tbYXWvrq6SktLS2ljY8OtrC8vL9P5+XkaDAZjG9rWtqreq6zn5+fTo0ePxjaVVZhc1dubCTeVVfW4ofG6bW3e0Lg2bqp+v39t7fX19WskZJ8o8LpxQ+OnGaysBUGIQ7wr3u2Kd8PFgC3C1Yp3c/v9fjo6Orp2c3lT8c3FC8TqFi+QN9Xq6mqam5sbSVVHR0dpf38/HRwcjKQqu7lYWW9vb6der+eubRuD126S6GxDo0R3eXk52lQ7OzujDY0P1qusUaLjDW1ro0THFSLKZLipmiQ6vO4miY7XNonu9PR0TKLD/SIIQhziXfFuV7wbLgaurq6u9UCwWkkpjR4s3lzsv1jfiSWb0gXazbVNhWv3er2RVHV+fj5WnVrfCde2v/XkIuz9WN+Kez9NEh1WmMPhcCRV8YPF60apCjc0y4M1Ep0nD+J9Y4nONmRJosNPFNjz8qp6XNveTEYEQjs87VL60/zabxPiXfFuV7wbLgbs5mIPhB/s2dlZse+E0kVpmMJuED5YHuTgvpOtbRsD1+bez/b29qj3Yw+Wby6vvb6+PnbdtnZJorO1PbkIZTLc0Pb3KNHhhvb6ThGJDjd0qe+EVb239tra2phEx70+lKrW1tbGqvrNzc0WdCgIswvxrni3K94NFwN2c5sGGvDBegMNOEyBF3h2dtbY+8EKs6nnxTIZVpg8SRvpO6HEZj2vnETX7/ddic6TyfC6cUN7vR+7bpPosOdlbwiW6FgezEl0OMWL143DNyZVmUTHG9o+UWDPCzd0r9drQYeCMLsQ74p3u+LdcDHg3VweaMAKEfsv8/PzI6mKp1mx94OyifVfvGEKu0DeVGwPwbWxQrSeV9PaJpOxRIe9H0+iM3sIS3S4tkl03hQvS3S8qbwJ4pxE1zRJi5vKk+hwU+UkOk8m47WNCARBiEO8K97tinfDxUCb/gsPU3gP1ru5tqnwAvnBel5P3lQmF9naTZYcezNx/6Uk0eU2FUt0uPbKyoor0Xm2GB5gYYnOKuMmic7kIpboeIAlZweyTYUSHd5z7nmxRGfXPevg3r966UIJ4l3xble8Gy4GTLKxG9QU3NDk9cz5TA8ODtJgMBh5PTk0gns/TV5Pb1Nx78ezaLDXkzeVJ9Hh2imNez1tU/V6vWsPFitEswPxFG9OouPKukai8/pObMnBXh9KdJ5UFZHoBEGIQ7wr3u2Kd8PFQNTraT2QSbyeXnJWjdfTe7CeVGU3Fx9MyevpSXRRr6e9ISIyGW8qHvyJWHJMouP0Ku63YYXZJNEtLy+P1kYSsuvOSXSYViYIQhziXfFuV7wbLgZQ7jHJJndzu/Z6epGMN+H19Na2Kq3J68lTvHjd+/v7o6rei8FkO5BJdLb2/v5+GgwGbmKY13eqicG0Z8YxmDlbjK2NlblVt5LCBaEO4l3xble826oYaOr95Hog0Z5XrdcTq1NvbXs4JZnM6794UhV7PXFDc8+Lez+5nhdu6GgMpq1tlfFgMBhddy4Gky05TRIdVqdNthhOK7PKGNe+S4ftRKG5AGESiHfFu13xbrgY8AY5Il7PaAymbSy+QK4QJ/F6bkEMptd3ysVgWpXHEl2k91OS6HJ9p4hEZ29Gs9W0icE0aTGXVhaV6DgxzJOqbG1BEOIQ74p3u+LdcDFgQyyTxmCy15MPichdIPZ+puH1xAdrm8qLwbQNzV5PT6LD3g9LdGjJ8dYuSXS5KErPklMj0eEnCuw7eRnokbQyPkyFJbqLi4sWdCgIswvxrni3K94NFwMoF5W8nth/8WIwnxSvZzQGEy0aXhQlP9ic1zN35CZX9fhgm2IwvbVxgtj+RWIwLy8vxwJOShJdUwymrYtSla0965iWtfAutR/u0rVMG+Jd8W5XvFvVJuAHG/V64oOdptczGoOJvZ9oDKYX3OB5PW1DR7yeXgymvXaT6HI+U8+SwydledddE4PJG9reTJhWxm9kLwYT7zknhgm
      "text/plain": [
       "<Figure size 640x480 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# gratings.py\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "x = np.arange(-500, 501, 1)\n",
    "\n",
    "X, Y = np.meshgrid(x, x)\n",
    "\n",
    "wavelength = 100\n",
    "angle = np.pi/9\n",
    "grating = np.sin(\n",
    "    2*np.pi*(X*np.cos(angle) + Y*np.sin(angle)) / wavelength\n",
    ")\n",
    "\n",
    "plt.set_cmap(\"gray\")\n",
    "\n",
    "plt.subplot(131)\n",
    "plt.imshow(grating)\n",
    "plt.axis(\"off\")\n",
    "\n",
    "# Calculate the Fourier transform of the grating\n",
    "ft = np.fft.ifftshift(grating)\n",
    "ft = np.fft.fft2(ft)\n",
    "ft = np.fft.fftshift(ft)\n",
    "\n",
    "plt.subplot(132)\n",
    "plt.imshow(abs(ft))\n",
    "plt.axis(\"off\")\n",
    "plt.xlim([480, 520])\n",
    "plt.ylim([520, 480])\n",
    "\n",
    "# Calculate the inverse Fourier transform of \n",
    "# the Fourier transform\n",
    "ift = np.fft.ifftshift(ft)\n",
    "ift = np.fft.ifft2(ift)\n",
    "ift = np.fft.fftshift(ift)\n",
    "ift = ift.real  # Take only the real part\n",
    "\n",
    "plt.subplot(133)\n",
    "plt.imshow(ift)\n",
    "plt.axis(\"off\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "778ad991",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgMAAAD1CAYAAADNj/Z6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOz9eZDsZ3Xfj797757pfZt97r1zV6EdIWRFYASkDGY3gYDjGJw4hSlju+LETpwEx3HVt4KdsuO4wAabFJjYjtkMtmww2BhsExBgISEkXenuc2ef3pfp6e7p7ffH/F5nnr4YSU4MSGieKpWke2e6P59nOed93ud9zuMZjUYjHY7DcTgOx+E4HIfjGTu83+kHOByH43AcjsNxOA7Hd3YcgoHDcTgOx+E4HIfjGT4OwcDhOByH43AcjsPxDB+HYOBwHI7DcTgOx+F4ho9DMHA4DsfhOByH43A8w8chGDgch+NwHI7DcTie4eMQDByOw3E4DsfhOBzP8HEIBg7H4Tgch+NwHI5n+DgEA4fjcByOw3E4DsczfPif7A96PJ5v5XMcjsNxOJ7EeDo2DH3rW9+qVqulVqulcDiswWCgZrOpZrOp6elpBQIBjUYj+Xw+FQoFBYNBhcNh+f1+xeNx9Xo91et15XI57e3tye/3azQaaXd3VxMTE2o2m5qYmFA4HJbH49He3p729vbU6XQUj8fl8XhUq9U0GAw0MTGheDyujY0NhcNh1Wo1TU1NqdvtamJiQoFAQJ1OR9vb2wqFQkokEvL7/QqHwxoOh2o0GpqcnNTe3p58Pp8kKRwOa2NjQ9lsVr1eTxMTE9rZ2VEsFlOj0ZAk+Xw+BYNBxWIxra2taXZ2Vh6PR48++qiOHj2q3d1dDYdDJZNJ1Wo17e7uKpfLyev1qt/vKxwO27PVajUFg0GFQiEVi0UlEgl5PB75fD4NBgPVajUlk0nF43HVajWNRiObX55jNBrZvITDYYVCIXW7Xfl8Pvn9+24hFAppOBxqd3dXgUBAkjQcDu29+YzRaKRKpaLBYCCfz6fRaKRUKiVJajQa8vl8thbRaFQej0ftdls7Ozs2l81m09a80+nYnLnP2+/31e12FQwGtbu7q9FoZGuaSCR0+fJlLS0tyefzaWNjQ0tLS9rd3VUoFFKv17O5rNVqikQi6na7ikaj2t3dld/vtzVlDiKRiAaDgXZ2drS7u6vJyUlJUiAQ0O7urrxer+0LSdrb21O/3x+bl+FwqMFgoHQ6Lb/fr36/bz/r8XgUDofVbDY1Go0UiUQUCoXsWVqtliKRiAqFgvL5vHw+n0qlkrxer62L3+9XKBTS3t6ehsOhJicn7e/5vMFgoMFgoPe9732Pe04PmYHDcTgOx7d01Go1tdtthcNheb1eBQIBJRIJ5XI5tVotSfvBhsfjUSaTUb1e18TEhCKRiIbDoRnWXq+ndrstr9crj8ejwWBgQKHf76vf76vX62k4HJqDw3lFIhEDIp1OR8Ph0MCG1+tVKpUyJ7e3t6d4PK54PC6fz2eOsN/vKxqNyuv1mnOZmJiQtO8Yu92uRqORGWhARDgc1uTkpPr9vgaDgfx+v4bDoUajkaLRqL2bz+fTcDhUMBiU1+tVp9Ox92y1WvaOGHeAB0CCd+S5O52O+v2+Op2Odnd31el01Ov11O/3NRwOVa1WbY12d3ft2SqVinq9ngGuSqWiZrOpbrdra9Dv9+Xz+WwtfD6frVM0GpUktdtt9Xo9SbJ3dn9/cnJStVpNPp/P1sbj8cjr9ZpDDQaD8vl8KhaLBqwqlYrtD4/HM7ZOfE4gEDCneO07er1ec/qADeZ7MBgoGAzK4/HYXAGS6vW6pH1AFAqFFA6Hbc12dnbU6/UUCAQ0OTmpWCymyclJRSIRRaNRdTod7e3tqVgsGtgNh8OSpGAwaOvPGvP+Xq/XQMne3p6CwaABAdad547FYkokEvL5fIrFYopEIqrX6yoWiwawHm88aWbgcByOw3E4/m9Gu92WJHNsktTr9ZRMJhUKhdRut9Xtdi3SOnbsmEVR3W5XnU5Ho9FIoVBI0n7ECMsQDAY1Pz+ver1uTsXv95shxQjiNIPBoAaDgaLRqLrdriYnJ83Y7+3taTQaKRAIWATb7/fHolYcvvtswWBQnU5HqVTKDPNwODTAEQwGJcn+XJJFjel02t6J58UxBINB7e3tKRwOmzPZ29vTysqKTp8+bRE84KPf79tch0IhtVotc9R8HlGyJDWbTXOYExMT5mhisZgCgYA5nZ2dHXuXdDptAI13l6RYLKZyuWzAie/r9XpqNpuKxWKKRqMaDAbq9/sW9U5NTanX66lWqymTyYwBvUAgYIAum83K6/WaI+/3+4rH42N74ujRo/L5fAqFQga+WKe9vT0NBgNzwsFgUN1u1xiRyclJNRoNTUxMGLAgig8GgwoEAorH4+p2uwoEArp69aqOHDmiXq+nQqGgUCikSCRiwI51AfywD/v9vvb29uT1eg38SFImk7E5Zm5brZZCoZD9bLfbNbADkAiHw8ZwhEIhdTodY5HY74DqJxqHzMDhOByH41s6IpGIfD6fut2uCoWC2u22AoGAOcd6vW5RXSwWkySL4IPBoNLptOLxuEajkTlEInPpIFIjRSDJnByG0ev1GgUNIPF4PBZZNRoN9Xo9dTodi1x7vZ7K5bK2traMIvb5fOasW62WRdCJRMIiVb/fb8/GcxChE9UDKPr9/thz8HeTk5PmCDqdjiYmJjQcDrW3t6eZmRn73K2tLXNWnU5H9XrdQBURL8+Ic8UhTk9Pm2PH0UqyqBuaOZFIKJ1OK5vNKhQKmeOBnSAlw+eTLoBpkPYdKhExz8W8dLtdSTIgxPq7z8qzsW9Ie+zu7tpeYO1J5zD/zE0oFLIUgN/vN4fP9wKcAKKdTkc7Ozu2TtD5Pp9P+XzegE65XB5bT5gN3i0QCCiZTCqTyWhxcdH2Is82MTFhIEfaB8qhUMj29HA4VLfbNZAEYGFeWC+YKADHzs6O2u22Wq3WITNwOA7H4fjOD4z+5OSk2u22OYOdnR2LNolwut2udnZ2lEgkNDExoU6no1gsJr/fPwYQJiYmzNFAkZMvT6fT5gCh3aFeyccSCUsHEZwkxeNxNZtN+Xw+7ezsKJPJaGJiwtIAOItrc+tEeUTee3t7ikQi2tnZkdfrVSwWM1ahUqkok8mYQ/J6vYrH42q328YQMFwwEA6HLU9MmoHoOxQKaTQamfNkTsiv8/foD4gUc7mcRcBExoALom8cTSQSGZvXeDyuVqtlkffk5KSxDICPYrGoyclJi+YBS5FIRI1Gw+j4cDhsqRBJRtXz96SZer2esTvNZtP+nznHITNPu7u7NicwJAARNAI45mQyqcnJSUs14UBJf0j7YIQ5cEEtrEE0GrWUCcwD68geYp6YD/QD7FWPx2OOfDQa2R6AsQDkkmLjz/lemAXA2GAwMDD7eOOQGTgch+NwfEsHefvd3V1NTU2p1WrpwoUL6vf7RqP3ej3Lr2azWaNPI5GIJGl9fV2NRkPxeNycFYYPgyzJ8uySTA9AJClJOzs7ZrBhAdycq3RA4Xs8HkWjUSWTSYtsMa6DwUCpVErRaNQiYgy3x+NRJBKxqM3r9arVapnzQZ+Ag5AOUigAl93dXQMPRN68UyqVMrp4YmJC7XbbIunJyUmNRiMDTXzH2tqaOUtYAFgZGBq0FYFAwIAZjA3fJ+3n3plj1mB3d9cAEPNEPr/ZbBqLwHu7oIXInd9vNpsG1FzGpNvt2s8DCCcnJy2dA0gLBAKanZ1Vr9dTLBYzfUOn07HUjssawSwRcZMuYq6g2tPptAKBgKWLAIOzs7MGvojiq9WqzRPgAuEn+0SSaV5gQfg+5sfj8ajValkqCoYNoSYg201bwZC47AMg6/HGITNwOA7H4fiWDhwGiufFxUWL7DF2RFI4WxxWrVbT5OSkEomEOWmcar/ftxw11OvExIRarZZFX6QGcH7JZNKeC8NKnnd3d1flctlEgeSAcSTQ0RhYBhQ1ETJODkdMxJZKpTQcDhWNRk1sRpSL6hxwQwSNw69UKpbz39nZ0WAwMMcJY4J
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# fourier_synthesis.py\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "image_filename = \"Earth.png\"\n",
    "\n",
    "def calculate_2dft(input):\n",
    "    ft = np.fft.ifftshift(input)\n",
    "    ft = np.fft.fft2(ft)\n",
    "    return np.fft.fftshift(ft)\n",
    "\n",
    "# Read and process image\n",
    "image = plt.imread(image_filename)\n",
    "image = image[:, :, :3].mean(axis=2)  # Convert to grayscale\n",
    "\n",
    "plt.set_cmap(\"gray\")\n",
    "\n",
    "ft = calculate_2dft(image)\n",
    "\n",
    "plt.subplot(121)\n",
    "plt.imshow(image)\n",
    "plt.axis(\"off\")\n",
    "plt.subplot(122)\n",
    "plt.imshow(np.log(abs(ft)))\n",
    "plt.axis(\"off\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "672ca19f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# fourier_synthesis.py\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "image_filename = \"Elizabeth_Tower_London.jpg\"\n",
    "\n",
    "def calculate_2dft(input):\n",
    "    ft = np.fft.ifftshift(input)\n",
    "    ft = np.fft.fft2(ft)\n",
    "    return np.fft.fftshift(ft)\n",
    "\n",
    "def calculate_2dift(input):\n",
    "    ift = np.fft.ifftshift(input)\n",
    "    ift = np.fft.ifft2(ift)\n",
    "    ift = np.fft.fftshift(ift)\n",
    "    return ift.real\n",
    "\n",
    "def calculate_distance_from_centre(coords, centre):\n",
    "    # Distance from centre is √(x^2 + y^2)\n",
    "    return np.sqrt(\n",
    "        (coords[0] - centre) ** 2 + (coords[1] - centre) ** 2\n",
    "    )\n",
    "\n",
    "def find_symmetric_coordinates(coords, centre):\n",
    "    return (centre + (centre - coords[0]),\n",
    "            centre + (centre - coords[1]))\n",
    "\n",
    "def display_plots(individual_grating, reconstruction, idx):\n",
    "    plt.subplot(121)\n",
    "    plt.imshow(individual_grating)\n",
    "    plt.axis(\"off\")\n",
    "    plt.subplot(122)\n",
    "    plt.imshow(reconstruction)\n",
    "    plt.axis(\"off\")\n",
    "    plt.suptitle(f\"Terms: {idx}\")\n",
    "    plt.pause(0.01)\n",
    "\n",
    "# Read and process image\n",
    "image = plt.imread(image_filename)\n",
    "image = image[:, :, :3].mean(axis=2)  # Convert to grayscale\n",
    "\n",
    "# Array dimensions (array is square) and centre pixel\n",
    "# Use smallest of the dimensions and ensure it's odd\n",
    "array_size = min(image.shape) - 1 + min(image.shape) % 2\n",
    "\n",
    "# Crop image so it's a square image\n",
    "image = image[:array_size, :array_size]\n",
    "centre = int((array_size - 1) / 2)\n",
    "\n",
    "# Get all coordinate pairs in the left half of the array,\n",
    "# including the column at the centre of the array (which\n",
    "# includes the centre pixel)\n",
    "coords_left_half = (\n",
    "    (x, y) for x in range(array_size) for y in range(centre+1)\n",
    ")\n",
    "\n",
    "# Sort points based on distance from centre\n",
    "coords_left_half = sorted(\n",
    "    coords_left_half,\n",
    "    key=lambda x: calculate_distance_from_centre(x, centre)\n",
    ")\n",
    "\n",
    "plt.set_cmap(\"gray\")\n",
    "\n",
    "ft = calculate_2dft(image)\n",
    "\n",
    "# Show grayscale image and its Fourier transform\n",
    "plt.subplot(121)\n",
    "plt.imshow(image)\n",
    "plt.axis(\"off\")\n",
    "plt.subplot(122)\n",
    "plt.imshow(np.log(abs(ft)))\n",
    "plt.axis(\"off\")\n",
    "plt.pause(2)\n",
    "\n",
    "# Reconstruct image\n",
    "fig = plt.figure()\n",
    "# Step 1\n",
    "# Set up empty arrays for final image and\n",
    "# individual gratings\n",
    "rec_image = np.zeros(image.shape)\n",
    "individual_grating = np.zeros(\n",
    "    image.shape, dtype=\"complex\"\n",
    ")\n",
    "idx = 0\n",
    "\n",
    "# All steps are displayed until display_all_until value\n",
    "display_all_until = 200\n",
    "# After this, skip which steps to display using the\n",
    "# display_step value\n",
    "display_step = 10\n",
    "# Work out index of next step to display\n",
    "next_display = display_all_until + display_step\n",
    "\n",
    "# Step 2\n",
    "for coords in coords_left_half:\n",
    "    # Central column: only include if points in top half of\n",
    "    # the central column\n",
    "    if not (coords[1] == centre and coords[0] > centre):\n",
    "        idx += 1\n",
    "        symm_coords = find_symmetric_coordinates(\n",
    "            coords, centre\n",
    "        )\n",
    "        # Step 3\n",
    "        # Copy values from Fourier transform into\n",
    "        # individual_grating for the pair of points in\n",
    "        # current iteration\n",
    "        individual_grating[coords] = ft[coords]\n",
    "        individual_grating[symm_coords] = ft[symm_coords]\n",
    "\n",
    "        # Step 4\n",
    "        # Calculate inverse Fourier transform to give the\n",
    "        # reconstructed grating. Add this reconstructed\n",
    "        # grating to the reconstructed image\n",
    "        rec_grating = calculate_2dift(individual_grating)\n",
    "        rec_image += rec_grating\n",
    "\n",
    "        # Clear individual_grating array, ready for\n",
    "        # next iteration\n",
    "        individual_grating[coords] = 0\n",
    "        individual_grating[symm_coords] = 0\n",
    "\n",
    "        # Don't display every step\n",
    "        if idx < display_all_until or idx == next_display:\n",
    "            if idx > display_all_until:\n",
    "                next_display += display_step\n",
    "                # Accelerate animation the further the\n",
    "                # iteration runs by increasing\n",
    "                # display_step\n",
    "                display_step += 10\n",
    "            display_plots(rec_grating, rec_image, idx)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d60f08d5",
   "metadata": {},
   "source": [
    "## Zadanie 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "7b39ab4f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAF7CAYAAADMhrFdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOz9O69t3ZOnBT4RMcaYc661z/u+WUmrEALENwCJSwm/JIw228BpgfgAOGlRDggLFwM+BB+gJdRSqtUWEi2w8XGqIMn8n3P2WnOOS0S0Md8uCVE0SQOdF+3HO/ucvbXOXnONEZdf/EIyM/niiy+++OKLL774G4T+Vb+AL7744osvvvjii/+1fAUwX3zxxRdffPHF3zi+Apgvvvjiiy+++OJvHF8BzBdffPHFF1988TeOrwDmiy+++OKLL774G8dXAPPFF1988cUXX/yN4yuA+eKLL7744osv/sbxFcB88cUXX3zxxRd/4/gKYL744osvvvjii79xfAUwX3zxxRdffPHF3zj+Wgcw/+l/+p/yL/wL/wL7vvP3/t7f47/8L//Lv+qX9MUXX3zxxRdf/DXgr20A85/9Z/8Zf/Inf8J/8B/8B/zX//V/zb/4L/6L/Bv/xr/Bf/ff/Xd/1S/tiy+++OKLL774K0b+ui5z/Ht/7+/xr/6r/yr/yX/ynwAQEfxz/9w/x7/77/67/Hv/3r/3V/zqvvjiiy+++OKLv0rKX/UL+CcxxuC/+q/+K/7BP/gH//hrqsrf//t/n//iv/gv/onf03un9/6P/xwR/Pmf/zl//Md/jIj8H/6av/jiiy+++OKL/+1kJj9//uSf+Wf+GVT/5xtFfy0DmD/7sz/D3fm7f/fv/o++/nf/7t/lv/lv/pt/4vf8R//Rf8R/+B/+h///eHlffPHFF1988cX/wfy3/+1/yz/7z/6z/7N//9cygPn/hX/wD/4Bf/Inf/KP//z9+3f++X/+n+f//n/7f/DLLx+EJDYKshtndx7FiG2Qb8M2kAVXDUwUvcAx6kMQT5Y4uRQkSXVUjQghxsTMWAS1GLsLpwhzOYdVJpNljvYNtU7BmOYUNTQrdb0Y1hAFqaBLWAH4Rs/OlkZkktWJEKYFx6q4OmSSYnACJTA3RpyEGp6GNMVPAZ2kLKYK22qgwUpHzCgGtgrTXywDsUYOwSvsLnRbqDtmBddCzkGVDY+gRoA4eIUiTN5E7mQm+s3xH8Jmg9AH4Y5sibrRL2WWpMWFWMWXEgUyFrssVlRWh8MGA7Dd0FehtOS1kqFOVaF7QXDqcPISekm0CHsmUxWLoOyBvwRvjWu72K+CakPLZM2F0pCHUhxKERIhlqIoGQFrMM2Zy0hVfILmYjZwKTy5OEPZ88TkV/IVpCSUCppsm5AySVdmdZo15J14KGETAZpWgo47qDWKDkA4I1AFT2A505zqD2iFcr5RMU4Xxp606ZQlfLakbcL16cgmbKvQRWlr0HJh0lhSMRl0q1hblFXITOpeaXNQbDG0wRLeNdmzoCiuga5kFYFMRATpjlajTOdqhSZKySBNeXenikMaHhNSiahMDX5J4fLAPy5SHmgGEvfntoShJVhRGHaxrYJn4gp1Je5GPSaS0AX2vpF6YnUnpGNTOVMIklaDTKWsBfGNkReiiluiqqgnxKBEIx+L+EPAcbBiUTe5X9NSdC+UNZgW6ABqUDkoVlnRkYQ3wS+aDAEzRd2RJcih2GX0EqQGvpTawTehBkgBb8nDDF9B+gJdrFbJK5GslLl4iYEGoQUGiPxkrEowEBJRxXIDK2g/wSpRIJajKXQPyl5AYJnQCPxTUU36AvkA64KH4AQftphMMoXwjerKaJNcSasgJDaVn8AjldOCyEnrDSnGoZMrkpDGUqGoQF6ogC6I94JjI8vFksYlQcsCTLayyP6BmKO58MOwa/L2jSGdDYGYqBdqANV4i0MoGwtV0FC8KhcXh1VaPnAcEnIZdeukFMocDN2o5c0aT0ZC8YVL3OdTHywTHjPpx8S00VJwXyxz2qoEhhN4AUvwnrhN5F2wXxprvdGu9I+Czca1TejOcxWIQR6NGoG+Lq5y3xu1NCIHWj5QH/SZOMExldkqWya5OdkDLIhURCaxFK+JDbh048FEYvDSX7BycvRC042lE6/JSkMvQ7dO+kBqhQb2Q1mbc6SxRvA2YWuFfiWPXehXZxOFJqQvph2Qg3Rokogl8VJWKRQzrFyECCcDyZ0iQemLUTdWLpS76lJECSmM73/Ov/Vv/Z/59u3b/9d7/69lAPNP/VP/FGbGP/pH/+h/9PV/9I/+Ef/0P/1P/xO/Z9s2tm37n379tz/i+OVJdrA9QBX9gM2NIYP1MHZXQgVdE1wp3wTZkhRFxDnEmG9lUydr4qKUWMh8kFlwm6xvynEm3uHXhxHhsA4OcXJXHg1cN4oLUhUvjYhfkJ5oDWY4Wzk41k88Nmpxlhm6NkQVMecZjVlABmxyEWdj/easLdkXbH7gBLmEF4WjTtwO3AdFjXolsUOuYBU4fGNuTuaDhwzECxdgm7MPZVenb4nmwmQjy05kp8aT1EUVQ0k+zdn9gcwHqk4UQ75NzlV4RMX3gkghLVleaM+J9Q+iOWYnawiWv2HjxJvxy6Z8XxtPMUY02DrnlhxdMYxa3uxy4N9g/pzYYWycoIU5Fu1hZBSKCFWN2Av0yXEkEWCPHXMHlK1tTILiiUjl3d7IMgBsfWPzN59HsgM/A2R/sn0o5e2IGjU3JL4xqtC2CaocW8FXghk2fmU+Fm7Bc044Kq7B8AdFFvWXDf9e2f/I+D6EbS+MUbEY/BaJ+cH3j5NfwpBVeJ+NX74VznC+ISyBVQeHPNmjkHvwdzJYxVEM1zfiH7R1seoHTR1tB6sVypqkGuGLoUJZG20povfvrdTkyEJpi5FBRqGJ3AF9cUwMLZ2pD3YLqJXaJu9T+WgLkQPZF2Mmczl7hWawpvCIQHJHoxF1Mmpj0851Obt949KFWUPnoqZR2iKmYKlkKVQxHiOprfDePihdiHZiCRWjLqAYowefx+CbVYoYewg/1NmmgCzMn2Re6PZB/GboTycT2ODh93vYBUQLVZRfxmRtSvpGaTBoPN6Lb89AX0p5Otu7ogcIF2IH8wn7pyHbYpZB/ijwW2LT0JLsdrDWhDaRKEQuWlV+mhFjgFT+yB4MW0yfrI87cTlmoTRjrUHNCsdEliIfjVxKrAF1Q3A2Km05+QFjOlsxcgvmz8T2SsqF1A2VRbmC3HbUA3tAY+dawiEDyaBrMiXIAY80qgqmSs4gH5UiSWsXTMPS6QtqBFKeZFSuZ+fbbxs6J5dvzKPyR5kMDbAncl3onsTa0P6gKvQP4dclSDR+irPLgf1o8CtoFR4J3WEPgQ2GF5oO9vUrUoN9KO9t0FRJEcSFfCl8fONpsPoHz187JQvLoJ3CtgW975wFci2+IXAkNStrBeJOrYUlylpgPlilwKMSZ8f/SFAK274xfnX+DoYtYVoDg1Sl5kDmRcgTfju4dOPXvNCpqHbWNMrHwS9RcSDVWaOicqGykwSrBNuZXO1X4mPgo1Bl8UsxVjUe/gtSleKF9otjSzniYKXSVbEPYUeYgK4nGoXxa1DzJC2p9cHjGBA7v2xAgW3f4ExkO4n8hq5gbx/kDFw7j6bE40kXY6zBw77hc7AZ+O7oWZASPFP5XiqPOsFPJjsfHf5Q74zmf0n+8ddyCqm1xr/8L//L/Omf/uk//lpE8Kd/+qf86//6v/6/6mfJKOhyTCfbKqg0jmxYDPZwjjOYY5GeZBZs2whZhBpyTXgnHgMpA9cdnUpcwYzKFCMIjrKxfSYJhA2iBBBog6hCeQQvKmMEg4rHYo2T5IWp4S40U2Y6ncrSRWAUV6QkWQerB8wL5kRscE3jks6Vib6FHJ2/sGBFo6+JrpNXEa4VdJRcf6Dron+CX7BNZfXFijtz6hO8dySDbb655M1nCNc70HMhK5GerLHRM0ktnLPzXkH7XmCC28XIJK+LNTqiMFRIT/QK/CysovhlZOn4Cf7ziYbQ5+uuMvWJnsavEWRNTAclAyXhgC0GxRtrW4zPi6VBbxc/pfF
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Got range [5.757753222571058e-06..1.1031931990186616].\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAB9IAAAG9CAYAAABJdz36AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOy9e7gkVXX3/1l776rq7nNmhotcBBRQkIgSjShGuWrUQUVDVBCNEUxiUBREf9GAb6LiBSJqREWRiNFXQeMrQXNRETGaYETjk2A0XggawQSjKDAwc0531d57rd8fu+e8HAd1QBR83Z/nmQe6TnXXrurqqlVrfdd3i5kZlUqlUqlUKpVKpVKpVCqVSqVSqVQqlUqlUqlUAHB39gAqlUqlUqlUKpVKpVKpVCqVSqVSqVQqlUqlUrkrUQvplUqlUqlUKpVKpVKpVCqVSqVSqVQqlUqlUqncglpIr1QqlUqlUqlUKpVKpVKpVCqVSqVSqVQqlUrlFtRCeqVSqVQqlUqlUqlUKpVKpVKpVCqVSqVSqVQqt6AW0iuVSqVSqVQqlUqlUqlUKpVKpVKpVCqVSqVSuQW1kF6pVCqVSqVSqVQqlUqlUqlUKpVKpVKpVCqVyi2ohfRKpVKpVCqVSqVSqVQqlUqlUqlUKpVKpVKpVG5BLaRXKpVKpVKpVCqVSqVSqVQqlUqlUqlUKpVKpXILaiG9UqlUKpVKpVKpVCqVSqVSqVQqlUqlUqlUKpVbUAvplUrlLs+73/1uRISrr776Nr/305/+NCLCpz/96Tt8XLdERHjFK17xY9e5+uqrERHe/e53/0zHUqlUKpVKpfKLzk8T//0y84UvfIGHP/zhLCwsICJ88Ytf5BWveAUismq9PfbYg+OOO+7OGeSc+h1XKpVKpVK5K7A5VvrBD37wE9f9aWKoPfbYgyOOOOJ2vfcXkfe+9738yq/8Ck3TsM022wBw2GGHcdhhh62s84uWK93a8f688tGVSuXnQ7izB1CpVCqVSqVSqVQqlUrlpyPGyFFHHcVoNOKNb3wjk8mE3Xfffave+9WvfpX/83/+D8cddxx77LHHz3aglUqlUqlUKpX/p/n617/Occcdx+GHH84pp5zCZDLZ6vd+9KMf5Z//+Z9/YsNSpVKp/LyohfRKpXKX53d+53c45phj6LruNr/3kEMOYTqd0rbtz2Bkt43dd9+d6XRK0zR39lAqlUqlUqlUKv+P8c1vfpNrrrmGd7zjHfz+7//+yvI//uM/5pRTTvmx7/3qV7/KaaedxmGHHfZzK6T/NDF+pVKpVCqVyp3BlVdeiXPV5Pcn8elPfxpV5U1vehN77bXXyvJLLrnkJ773ox/9KG9961t/oQvpd6V8dKVS+emphfRKpXKXZWlpiYWFBbz3eO9v12c45xiNRnfwyG4fInKXGUulUqlUKpVK5f8trrvuOoAV68zNhBAI4c559N8cz98aP02MX6lUKpVKpXJnUAWAW8ePikvvrMKymTGbzRiPxz+X7d2V8tGVSuWnp8qnKpXKz5wrrriCxz72saxdu5bFxUV+4zd+g8997nOr1tk8R+I//MM/cMIJJ7Djjjuy2267rfrbLedPVFVe8YpXsMsuuzCZTHjEIx7BV7/61S3mKrq1OWkOO+ww7n//+/PVr36VRzziEUwmE3bddVfOPPPMVWMahoGXvexl7L///qxbt46FhQUOPvhgPvWpT92u43Br8+gcd9xxLC4u8u1vf5sjjjiCxcVFdt11V9761rcC8OUvf5lHPvKRLCwssPvuu/O+971v1WfecMMN/OEf/iH77bcfi4uLrF27lsc+9rH827/92xbbv+aaa3jiE5/IwsICO+64Iy984Qv5+Mc/fqtz9nz+85/n8MMPZ926dUwmEw499FD+6Z/+6Xbtd6VSqVQqlcodxdve9jbud7/70XUdu+yyC8973vPYsGHDyt/f/OY3471ftewNb3gDIsKLXvSilWU5Z9asWcMf/dEf/djtbZ7L8pJLLuGBD3wgo9GIfffdl4suumjVerclJnvLW97C/e53PyaTCdtuuy0PfvCDV8V4Gzdu5OSTT2aPPfag6zp23HFHHv3oR/Ov//qvP3Kcxx13HIceeigARx11FCKyMv/krc2Rfkve/e53c9RRRwHwiEc8AhHZIj782Mc+xsEHH8zCwgJr1qzh8Y9/PF/5yle2GMPi4iLf/OY3edzjHseaNWv47d/+7R+73R+O8Tcf709/+tM8+MEPZjwes99++62M5aKLLmK//fZjNBqx//77c8UVV6z6zC996Uscd9xx3Ote92I0GrHzzjvzu7/7u1x//fVbbH/zNkajEfe+970599xzf+SxOv/889l///0Zj8dst912HHPMMfzXf/3Xj9y3SqVSqVQqv3hs2LCB4447jm222YZ169bxrGc9i+Xl5VXr3Noc6V/60pc49NBDGY/H7Lbbbrz61a/mXe961xZxzmY+85nPcMABBzAajbjXve7Fe97znp84ts05xde//vW88Y1vZPfdd2c8HnPooYfy7//+71uMZ2vioa2JOa+66iqe/OQns/POOzMajdhtt9045phjuOmmm37kWPfYYw9e/vKXA7DDDjsgIivd5T88R/oPc9xxx63kRDfHpLeMzVSVs846i/vd736MRiN22mknjj/+eG688cYtxnDEEUfw8Y9/fCWmPPfcc4HyPZ988snc4x73oOs69tprL1772teiqqs+Y/P5sG7dOrbZZhuOPfbYVc8YP44fl4/efL5MJhP22msvLrzwQgD+4R/+gYc+9KGMx2P22WcfLr300lWfec0113DCCSewzz77MB6P2X777TnqqKNu9Ry7Lefk1sT5lcovO7UjvVKp/Ez5yle+wsEHH8zatWt5yUteQtM0nHvuuRx22GErAcItOeGEE9hhhx142ctextLS0o/83FNPPZUzzzyTJzzhCaxfv55/+7d/Y/369cxms60a14033sjhhx/Ok570JI4++mguvPBC/uiP/oj99tuPxz72sQDcfPPNnHfeeTztaU/j2c9+Nhs3buSd73wn69ev55//+Z954AMfeLuPyy3JOfPYxz6WQw45hDPPPJMLLriA5z//+SwsLPC//tf/4rd/+7d50pOexNvf/nae+cxn8rCHPYw999wTgP/8z//kwx/+MEcddRR77rkn3/ve9zj33HM59NBD+epXv8ouu+wClG6gRz7ykfzP//wPL3jBC9h555153/ved6uigL//+7/nsY99LPvvvz8vf/nLcc7xrne9i0c+8pFcdtllHHDAAXfIflcqlUqlUqncFl7xildw2mmn8ahHPYrnPve5XHnllZxzzjl84Qtf4J/+6Z9omoaDDz4YVeUzn/kMRxxxBACXXXYZzjkuu+yylc+64oor2LRpE4cccshP3O5VV13FU5/6VJ7znOdw7LHH8q53vYujjjqKiy++mEc/+tHA1sdk73jHOzjppJN4ylOewgte8AJmsxlf+tKX+PznP8/Tn/50AJ7znOdw4YUX8vznP599992X66+/ns985jN87Wtf40EPetCtjvH4449n11135fTTT+ekk07iIQ95CDvttNNWHddDDjmEk046iTe/+c289KUv5b73vS/Ayn/f+973cuyxx7J+/Xpe+9rXsry8zDnnnMNBBx3EFVdcscoKPqXE+vXrOeigg3j9619/m+bD3Mw3vvENnv70p3P88cfzjGc8g9e//vU84QlP4O1vfzsvfelLOeGEEwA444wzOProo1dZrH7iE5/gP//zP3nWs57FzjvvzFe+8hX+/M//nK985St87nOfW0nEXnHFFRx++OHc/e5357TTTiPnzCtf+Up22GGHLcbzmte8hj/5kz/h6KOP5vd///f5/ve/z1ve8hYOOeQQrrjiii06rSqVSqVSqfxicvTRR7Pnnntyxhln8K//+q+cd9557Ljjjrz2ta/9ke+59tprV4SIp556KgsLC5x33nk/snP9G9/4Bk95ylP4vd/7PY499lj+4i/+guOOO47999+f+93vfj9xjO95z3vYuHEjz3ve85jNZrzpTW/ikY98JF/+8pdXYr+
      "text/plain": [
       "<Figure size 2500x1800 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "img = plt.imread(\"wmi1.jpg\")/float(2**8)\n",
    "plt.imshow(img)\n",
    "plt.show()\n",
    "\n",
    "shape = img.shape[:2]\n",
    "\n",
    "def draw_cicle(shape,diamiter):\n",
    "    assert len(shape) == 2\n",
    "    TF = np.zeros(shape,dtype=np.bool)\n",
    "    center = np.array(TF.shape)/2.0\n",
    "\n",
    "    for iy in range(shape[0]):\n",
    "        for ix in range(shape[1]):\n",
    "            TF[iy,ix] = (iy- center[0])**2 + (ix - center[1])**2 < diamiter **2\n",
    "    return(TF)\n",
    "\n",
    "\n",
    "TFcircleIN   = draw_cicle(shape=img.shape[:2],diamiter=150)\n",
    "TFcircleOUT  = ~TFcircleIN\n",
    "\n",
    "\n",
    "\n",
    "fft_img = np.zeros_like(img,dtype=complex)\n",
    "for ichannel in range(fft_img.shape[2]):\n",
    "    fft_img[:,:,ichannel] = np.fft.fftshift(np.fft.fft2(img[:,:,ichannel]))\n",
    "\n",
    "\n",
    "\n",
    "def filter_circle(TFcircleIN,fft_img_channel):\n",
    "    temp = np.zeros(fft_img_channel.shape[:2],dtype=complex)\n",
    "    temp[TFcircleIN] = fft_img_channel[TFcircleIN]\n",
    "    return(temp)\n",
    "\n",
    "fft_img_filtered_IN = []\n",
    "fft_img_filtered_OUT = []\n",
    "## for each channel, pass filter\n",
    "for ichannel in range(fft_img.shape[2]):\n",
    "    fft_img_channel  = fft_img[:,:,ichannel]\n",
    "    ## circle IN\n",
    "    temp = filter_circle(TFcircleIN,fft_img_channel)\n",
    "    fft_img_filtered_IN.append(temp)\n",
    "    ## circle OUT\n",
    "    temp = filter_circle(TFcircleOUT,fft_img_channel)\n",
    "    fft_img_filtered_OUT.append(temp) \n",
    "    \n",
    "fft_img_filtered_IN = np.array(fft_img_filtered_IN)\n",
    "fft_img_filtered_IN = np.transpose(fft_img_filtered_IN,(1,2,0))\n",
    "fft_img_filtered_OUT = np.array(fft_img_filtered_OUT)\n",
    "fft_img_filtered_OUT = np.transpose(fft_img_filtered_OUT,(1,2,0))\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "def inv_FFT_all_channel(fft_img):\n",
    "    img_reco = []\n",
    "    for ichannel in range(fft_img.shape[2]):\n",
    "        img_reco.append(np.fft.ifft2(np.fft.ifftshift(fft_img[:,:,ichannel])))\n",
    "    img_reco = np.array(img_reco)\n",
    "    img_reco = np.transpose(img_reco,(1,2,0))\n",
    "    return(img_reco)\n",
    "\n",
    "\n",
    "img_reco              = inv_FFT_all_channel(fft_img)\n",
    "img_reco_filtered_IN  = inv_FFT_all_channel(fft_img_filtered_IN)\n",
    "img_reco_filtered_OUT = inv_FFT_all_channel(fft_img_filtered_OUT)\n",
    "\n",
    "fig = plt.figure(figsize=(25,18))\n",
    "ax  = fig.add_subplot(1,3,1)\n",
    "ax.imshow(np.abs(img_reco))\n",
    "ax.set_title(\"original image\")\n",
    "\n",
    "ax  = fig.add_subplot(1,3,2)\n",
    "ax.imshow(np.abs(img_reco_filtered_IN))\n",
    "ax.set_title(\"low pass filter image\")\n",
    "\n",
    "\n",
    "ax  = fig.add_subplot(1,3,3)\n",
    "ax.imshow(np.abs(img_reco_filtered_OUT))\n",
    "ax.set_title(\"high pass filtered image\")\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}