{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.fft as fp\n",
    "import matplotlib.pyplot as plt\n",
    "from PIL import Image\n",
    "import cv2\n",
    "%matplotlib inline"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Uploading image and searching for spectrum\n",
    "\n",
    "Uploading an image \"dog.jpg \" using the imread function and immediately bring it\n",
    "to the double type for correct calculations. We will find the image spectrum using\n",
    "the FFT. We will search for the spectrum of each color channel separately,\n",
    "and also output the logarithm of the spectrum module for clarity.\n",
    "We additionally convert the spectrum [0, 1] range.\n",
    "This is necessary for the correct display of the spectrum, otherwise there will\n",
    "be a white screen on the graph. Then, we concatenate the obtained spectra."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Utility functions\n",
    "\n",
    "remmax = lambda x: x/x.max()\n",
    "remmin = lambda x: x - np.amin(x, axis=(0, 1), keepdims=True)\n",
    "touint8 = lambda x: (remmax(remmin(x)) * (256-1e-4)).astype(int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\Maksim Vinogradov\\AppData\\Local\\Temp\\ipykernel_3084\\3252584836.py:8: ComplexWarning: Casting complex values to real discards the imaginary part\n",
      "  imfft[:, :, i] = fp.fft2(img[:, :, i])\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 493828 characters]",
      "text/plain": [
       "<Figure size 700x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Loading image\n",
    "im = Image.open('dog.jpg')\n",
    "c, r = im.size\n",
    "\n",
    "img = np.array(im)\n",
    "imfft = np.zeros((r, c, 3))\n",
    "for i in range(3):\n",
    "    imfft[:, :, i] = fp.fft2(img[:, :, i])\n",
    "imf = touint8(np.log(abs(fp.fftshift(imfft)) + 1))\n",
    "\n",
    "fig1, axes1 = plt.subplots(1, 2, figsize=(7, 5))\n",
    "axes1[0].imshow(im)\n",
    "axes1[0].set_title('Original image', fontsize=20)\n",
    "axes1[1].imshow(imf)\n",
    "axes1[1].set_title('Image spectrum', fontsize=20)\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Setting the low pass filter\n",
    "\n",
    "It is necessary to set a low-pass filter in the form of a Gaussian function.\n",
    "The cutoff frequency is 0.1 1/pixel. To set a Gaussian filter in matlab, you\n",
    "need to translate the cutoff frequency into image fractions.\n",
    "To do this, we find the sizes `r` and `c` of the image by the size function and\n",
    "multiply the cutoff frequency by the smaller side."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 34680 characters]",
      "text/plain": [
       "<Figure size 700x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# making gaussian low-pass filter mask\n",
    "s = 0.1  # sigma in 1/pix\n",
    "sigma = 0.1 * min(im.size)\n",
    "\n",
    "def gauss(row, column, std):  \n",
    "    x, y = np.meshgrid(np.arange(-column/2, column/2), np.arange(-row/2, row/2))\n",
    "    arg   = -(x ** 2 + y ** 2) / (2 * std ** 2)\n",
    "    h = np.exp(arg)\n",
    "    h[h < 2.2204e-16 * np.amax(h)] = 0\n",
    "    sumh = np.sum(h)\n",
    "    if sumh:\n",
    "        h = h / sumh\n",
    "    return h\n",
    "    \n",
    "h = gauss(r, c, sigma)  # gaussian mask\n",
    "f = np.linspace(-0.5, 0.5, c)  # frequency grid\n",
    "\n",
    "# graphics\n",
    "plt.figure(figsize=(7, 4))\n",
    "norm = np.amax(h[round(r / 2), :]) ** (-1)  # normalizing coefficient\n",
    "plt.plot(f, norm * h[round(r / 2), :], 'b')\n",
    "plt.axvline(x=s, color='r', lw=1)\n",
    "plt.axvline(x=-s, color='r', lw=1)\n",
    "plt.grid()\n",
    "plt.title('Gaussian filter with $\\sigma = %.1f$ pix$^{-1}$' % (s,), fontsize=20)\n",
    "plt.xlabel('Frequency, pix$^{-1}$', fontsize=18)\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Filtering in frequency domain\n",
    "\n",
    "Let's apply the resulting low-pass filter to the image spectrum.\n",
    "As mentioned above, the spectrum of the image is found using FFT.\n",
    "We multiply the obtained spectrum by the filter matrix, and then using the\n",
    "inverse Fourier transform we get the processed image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 44028 characters]",
      "text/plain": [
       "<Figure size 1000x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\Maksim Vinogradov\\AppData\\Local\\Temp\\ipykernel_3084\\457963194.py:10: ComplexWarning: Casting complex values to real discards the imaginary part\n",
      "  im1fft[:, :, i] = fp.fftshift(fp.fft2(img[:, :, i])) * h\n",
      "Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 231740 characters]",
      "text/plain": [
       "<Figure size 700x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(10, 6))\n",
    "plt.imshow(h, cmap='gray')\n",
    "plt.title('Gaussian mask in frequency domain', fontsize=20)\n",
    "plt.show()\n",
    "\n",
    "\n",
    "im1fft = np.zeros((r, c, 3))\n",
    "im1 = np.zeros((r, c, 3))\n",
    "for i in range(3):\n",
    "    im1fft[:, :, i] = fp.fftshift(fp.fft2(img[:, :, i])) * h\n",
    "    im1[:, :, i] = abs(fp.ifft2(fp.ifftshift(fp.fftshift(fp.fft2(img[:, :, i])) * h)))\n",
    "im1 = touint8(im1)\n",
    "\n",
    "fig2, axes2 = plt.subplots(1, 2, figsize=(7, 4))\n",
    "axes2[0].imshow(im1)\n",
    "axes2[0].set_title('Original Image', fontsize=20)\n",
    "axes2[1].imshow(np.log(abs(im1fft) + 1))\n",
    "axes2[1].set_title('Image spectrum', fontsize=20)\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It can be seen that mainly low frequencies remained in the spectrum.\n",
    "The image, in accordance with this, has become more blurry."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Setting the spatial filter\n",
    "\n",
    "We get a mask for spatial filtering. To do this, we will find the inverse\n",
    "Fourier transform from the low-pass filter obtained above.\n",
    "Let's cut the resulting mask to the size of 11x11 so that it is symmetrical with\n",
    "respect to the central element."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 16728 characters]",
      "text/plain": [
       "<Figure size 700x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "hs = fp.ifftshift(fp.ifft2(h))  # gaussian filter in frequency domain\n",
    "\n",
    "plt.figure(figsize=(7, 4))\n",
    "plt.imshow(np.log(abs(hs)+1), cmap='gray')\n",
    "plt.title('Gaussian mask in spatial domain')\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This mask has been normalized so that the sum of all its elements is equal to one."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Filtering in the spatial domain\n",
    "\n",
    "Let's apply the resulting mask to the original image.\n",
    "To do this, it is necessary to perform a convolution operation. By default,\n",
    "the image is padded with zeros at the edges."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\Maksim Vinogradov\\AppData\\Local\\Temp\\ipykernel_3084\\3083132412.py:12: ComplexWarning: Casting complex values to real discards the imaginary part\n",
      "  im2fft[:, :, i] = fp.fft2(im2[:, :, i])\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 398544 characters]",
      "text/plain": [
       "<Figure size 700x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "uncut = 11  # size of spatial domain mask\n",
    "cut = np.floor(uncut/2)  # half-size of spatial domain mask\n",
    "\n",
    "hs_cut = hs[int(r / 2 - cut) : int(r / 2 + cut + 1), int(c / 2 - cut) : int(c / 2 + cut + 1)]  # cut mask\n",
    "H = np.sum(abs(hs_cut)) ** (-1) * abs(hs_cut)  # normalized mask\n",
    "\n",
    "im2 = np.zeros((r,c,3))\n",
    "im2 = cv2.filter2D(img, ddepth=-1, kernel=H, borderType=cv2.BORDER_CONSTANT)\n",
    "\n",
    "im2fft = np.zeros((r,c,3))\n",
    "for i in range(3):\n",
    "    im2fft[:, :, i] = fp.fft2(im2[:, :, i])\n",
    "\n",
    "im2fft_shift = touint8(np.log(abs(fp.fftshift(im2fft))+1))\n",
    "\n",
    "fig3, axes3 = plt.subplots(1, 2, figsize=(7, 4))\n",
    "axes3[0].imshow(im2)\n",
    "axes3[0].set_title('Original image', fontsize=20)\n",
    "axes3[1].imshow(im2fft_shift)\n",
    "axes3[1].set_title('Image spectrum', fontsize=20)\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The image was also blurred, and high frequencies were suppressed in the spectrum."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Comparison of results\n",
    "\n",
    "As can be seen, the results obtained by filtering in the frequency domain with a\n",
    "Gaussian filter and by filtering in the spatial domain with a Fourier image of\n",
    "the same filter give identical results. The images visually look exactly the same.\n",
    "The spectra of these images also look the same, but the spectrum of the image\n",
    "processed in the spatial domain has pronounced vertical and horizontal bands\n",
    "passing through the center of the spectrum. This is due to the fact that we used\n",
    "a cropped mask, which does not exactly correspond to the Fourier image of the filter.\n",
    "This led to the appearance of artifacts in the spectrum.\n",
    "The visual differences between these images are only at their edges.\n",
    "The image processed in the spatial domain has darkened edges,\n",
    "because the addition of zeros was used. This phenomenon can be eliminated by\n",
    "using the addition of a periodic continuation of the image along the edges."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 47636 characters]",
      "text/plain": [
       "<Figure size 1000x600 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "im3 = cv2.filter2D(img, ddepth=-1, kernel=H, borderType=cv2.BORDER_REPLICATE)\n",
    "\n",
    "plt.figure(figsize=(10, 6))\n",
    "plt.subplot(131); plt.imshow(im1[0:50, 0:50, :]); plt.title('Frequency processed image')\n",
    "plt.subplot(132); plt.imshow(im2[0:50, 0:50, :])\n",
    "plt.title('Spatially processed image\\nZero boundary option')\n",
    "plt.subplot(133); plt.imshow(im3[0:50, 0:50, :])\n",
    "plt.title('Spatially processed image\\nPeriodic boundary option')\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It can be seen that the edges are darkened in the middle image corresponding to\n",
    "the addition of zeros, but when using periodic continuation, this effect is not\n",
    "observed, and it is visually quite difficult to find any differences from frequency processing."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Conclusion\n",
    "In this work, an algorithm for linear image filtering in the frequency domain\n",
    "and in the spatial domain was implemented, and these two filtering methods were compared.\n",
    "• application of a low-pass filter in the frequency domain\n",
    "(multiplication by the image spectrum) is identical to convolution of its Fourier\n",
    "image with an image in the spatial domain\n",
    "• however, the results are still different if, when processing in the spatial domain,\n",
    "not the full Fourier image of the filter is used, but only a fragment of it.\n",
    "This leads to the appearance of artifacts in the image spectrum, and hence in the image itself.\n",
    "• when using the addition of zeros, a \"frame\" appears at the edges of the image\n",
    "(darkened pixels). To avoid this, it is necessary to use addition by duplicating\n",
    "values at the edges, periodic repetition of the image and other methods of additions\n",
    "• calculation of the inverse Fourier transform of the pot filter and cropping,\n",
    "output of the resulting mask to the screen."
   ]
  }
 ],
 "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.10.5"
  },
  "vscode": {
   "interpreter": {
    "hash": "369f2c481f4da34e4445cda3fffd2e751bd1c4d706f27375911949ba6bb62e1c"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
