codes/07_fourier-transform.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Dependencies"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\n",
        "import numpy as np\n",
        "import scipy as sp\n",
        "from matplotlib.gridspec import GridSpec\n",
        "\n",
        "# local dependencies\n",
        "from utils.dft import (\n",
        "    fourier_basis_1d,\n",
        "    fourier_basis_2d,\n",
        "    dft,\n",
        "    idft,\n",
        "    dft2,\n",
        "    idft2\n",
        ")\n",
        "\n",
        "from utils.filters import (\n",
        "    ideal as ideal_filter,\n",
        "    sinc as sinc_filter,\n",
        "    gaussian as gaussian_filter\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Discrete Fourier Transform (DFT)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 1-dimensional signal"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# create basis vectors with period(N) = 5\n",
        "signal_1d_length = 5\n",
        "basis_vectors = fourier_basis_1d(signal_1d_length)\n",
        "\n",
        "# plot\n",
        "fig, axs = plt.subplots(nrows=3, ncols=signal_1d_length, figsize=(16, 4), layout='compressed')\n",
        "fig.suptitle(\"1D DFT basis vectors\")\n",
        "\n",
        "for i in range(signal_1d_length):\n",
        "    axs[0, i].imshow(basis_vectors[i, :].real.reshape(1, -1), cmap='gray', vmin=basis_vectors.real.min(), vmax=basis_vectors.real.max())\n",
        "    axs[0, i].set_title(f\"basis: {i} [real]\")\n",
        "    axs[0, i].set_yticks([])\n",
        "    axs[0, i].set_xticks(range(signal_1d_length))\n",
        "    axs[1, i].imshow(basis_vectors[i, :].imag.reshape(1, -1), cmap='gray', vmin=basis_vectors.imag.min(), vmax=basis_vectors.imag.max())\n",
        "    axs[1, i].set_title(f\"basis: {i} [imaginary]\")\n",
        "    axs[1, i].set_yticks([])\n",
        "    axs[1, i].set_xticks(range(signal_1d_length))\n",
        "    axs[2, i].imshow(np.abs(basis_vectors[i, :]).reshape(1, -1), cmap='gray', vmin=np.abs(basis_vectors).min(), vmax=np.abs(basis_vectors).max())\n",
        "    axs[2, i].set_title(f\"basis: {i} [magnitude]\")\n",
        "    axs[2, i].set_yticks([])\n",
        "    axs[2, i].set_xticks(range(signal_1d_length))\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "real part of basis vectors:\n",
            "[[ 1.          1.          1.          1.          1.        ]\n",
            " [ 1.          0.30901699 -0.80901699 -0.80901699  0.30901699]\n",
            " [ 1.         -0.80901699  0.30901699  0.30901699 -0.80901699]\n",
            " [ 1.         -0.80901699  0.30901699  0.30901699 -0.80901699]\n",
            " [ 1.          0.30901699 -0.80901699 -0.80901699  0.30901699]]\n"
          ]
        }
      ],
      "source": [
        "print(f\"real part of basis vectors:\\n{basis_vectors.real}\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "imaginary part of basis vectors:\n",
            "[[ 0.          0.          0.          0.          0.        ]\n",
            " [ 0.         -0.95105652 -0.58778525  0.58778525  0.95105652]\n",
            " [ 0.         -0.58778525  0.95105652 -0.95105652  0.58778525]\n",
            " [ 0.          0.58778525 -0.95105652  0.95105652 -0.58778525]\n",
            " [ 0.          0.95105652  0.58778525 -0.58778525 -0.95105652]]\n"
          ]
        }
      ],
      "source": [
        "print(f\"imaginary part of basis vectors:\\n{basis_vectors.imag}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 2-dimensional signal"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "metadata": {},
      "outputs": [],
      "source": [
        "# create basis images with period(N, N) = (8, 8)\n",
        "signal_2d_length = 8\n",
        "basis_images = fourier_basis_2d(signal_2d_length)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# plot\n",
        "fig, axs = plt.subplots(nrows=signal_2d_length, ncols=signal_2d_length, figsize=(16, 18), layout='compressed')\n",
        "fig.suptitle(\"2D DFT basis images [real part]\")\n",
        "\n",
        "for i in range(signal_2d_length):\n",
        "    for j in range(signal_2d_length):\n",
        "        axs[i, j].imshow(basis_images[i, j].real.reshape(signal_2d_length, signal_2d_length), cmap='gray', vmin=basis_images.real.min(), vmax=basis_images.real.max())\n",
        "        axs[i, j].set_title(f\"basis: {i},{j}\")\n",
        "        axs[i, j].set_yticks([])\n",
        "        axs[i, j].set_xticks([])\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# plot\n",
        "fig, axs = plt.subplots(nrows=signal_2d_length, ncols=signal_2d_length, figsize=(16, 18), layout='compressed')\n",
        "fig.suptitle(\"1d basis vectors [imaginary part]\")\n",
        "\n",
        "for i in range(signal_2d_length):\n",
        "    for j in range(signal_2d_length):\n",
        "        axs[i, j].imshow(basis_images[i, j].imag.reshape(signal_2d_length, signal_2d_length), cmap='gray', vmin=basis_images.imag.min(), vmax=basis_images.imag.max())\n",
        "        axs[i, j].set_title(f\"basis: {i},{j}\")\n",
        "        axs[i, j].set_yticks([])\n",
        "        axs[i, j].set_xticks([])\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Reconstructing a 1d signal"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# 1d array of length 5 (N=5 which is considered as the period of the signal)\n",
        "arr1 = np.array([3, 0, 2, 1, 5])\n",
        "\n",
        "# DFT : transform signal from spatial domain to frequency domain\n",
        "dft_arr1 = dft(arr1)\n",
        "\n",
        "# IDFT: transform signal from frequency domain to spatial domain\n",
        "idft_arr1 = idft(dft_arr1)\n",
        "\n",
        "# take only the real part because our original signal is a real part\n",
        "idft_arr1 = idft_arr1.real\n",
        "\n",
        "# clip the signal values in range(0, 5) [in this case we know that for example our original signal must have values in this range]\n",
        "idft_arr1 = idft_arr1.clip(arr1.min(), arr1.max())\n",
        "\n",
        "# change dtype from float to int (original signal had integer values)\n",
        "idft_arr1 = np.round(idft_arr1).astype(np.int32)\n",
        "\n",
        "# plot\n",
        "fig, axs = plt.subplots(nrows=1, ncols=3, figsize=(12, 2), layout='compressed')\n",
        "fig.suptitle(\"Reconstructing 1D signal\")\n",
        "\n",
        "axs[0].imshow(arr1.reshape(1, -1), cmap='gray')\n",
        "axs[0].set_title(\"Original signal\")\n",
        "axs[1].imshow(np.abs(dft_arr1).reshape(1, -1), cmap='gray')\n",
        "axs[1].set_title(\"Magnitude in frequency domain\")\n",
        "axs[2].imshow(idft_arr1.reshape(1, -1), cmap='gray')\n",
        "axs[2].set_title(\"Reconstructed signal\")\n",
        "\n",
        "for ax in fig.axes:\n",
        "    ax.set_yticks([])\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Reconstructing a 2d signal"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# 2d array of length 3x3\n",
        "arr2 = np.array([[1, 1, 3], [2, 1, 2], [2, 3, 3]])\n",
        "\n",
        "# DFT : transform signal from spatial domain to frequency domain\n",
        "dft_arr2 = dft2(arr2)\n",
        "\n",
        "# IDFT: transform signal from frequency domain to spatial domain\n",
        "idft_arr2 = idft2(dft_arr2)\n",
        "\n",
        "# take only the real part because our original signal is a real part\n",
        "idft_arr2 = idft_arr2.real\n",
        "\n",
        "# clip the signal values in range(1, 3) [in this case we know that for example our original signal must have values in this range]\n",
        "idft_arr2 = idft_arr2.clip(arr2.min(), arr2.max())\n",
        "\n",
        "# change dtype from float to int (original signal had integer values)\n",
        "idft_arr2 = np.round(idft_arr2).astype(np.int32)\n",
        "\n",
        "# plot\n",
        "fig, axs = plt.subplots(nrows=1, ncols=3, figsize=(12, 5), layout='compressed')\n",
        "fig.suptitle(\"Reconstructing 2D signal\")\n",
        "\n",
        "axs[0].imshow(arr2, cmap='gray')\n",
        "axs[0].set_title(\"Original signal\")\n",
        "axs[1].imshow(np.abs(dft_arr2), cmap='gray')\n",
        "axs[1].set_title(\"Magnitude in frequency domain\")\n",
        "axs[2].imshow(idft_arr2, cmap='gray')\n",
        "axs[2].set_title(\"Reconstructed signal\")\n",
        "\n",
        "for ax in fig.axes:\n",
        "    ax.set_xticks(range(3))\n",
        "    ax.set_yticks(range(3))\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Discrete Fourier Transform Using numpy\n",
        "   - Discrete Fourier Transform (DFT) is calculated using the efficient Fast Fourier Transform (FFT) algorithm."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## f(x) = 1\n",
        "   - DFT(f(x)) = impulse signal"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "metadata": {},
      "outputs": [],
      "source": [
        "# ideal signal\n",
        "signal_length = 200\n",
        "signal = np.ones(shape=(signal_length,))\n",
        "\n",
        "# fft\n",
        "fft_signal = np.fft.fft(signal)\n",
        "\n",
        "# shift the low frequencies to the center\n",
        "fftshift_signal = np.fft.fftshift(fft_signal)\n",
        "\n",
        "# magnitude of the signal\n",
        "mag_fft_signal = np.abs(fft_signal)\n",
        "mag_fftshift_signal = np.abs(fftshift_signal)\n",
        "\n",
        "# phase of the signal\n",
        "phase_fft_signal = np.angle(fft_signal)\n",
        "phase_fftshift_signal = np.angle(fftshift_signal)\n",
        "\n",
        "# reconstruct the signal\n",
        "reconstructed_signal = np.fft.ifft(np.fft.ifftshift(fftshift_signal)).real"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# plot\n",
        "fig, axs = plt.subplots(nrows=1, ncols=3, figsize=(16, 4), layout='compressed')\n",
        "\n",
        "axs[0].plot(signal)\n",
        "axs[0].set_title('Original signal')\n",
        "axs[1].plot(fftshift_signal.real)\n",
        "axs[1].set_title('fft [real]')\n",
        "axs[2].plot(fftshift_signal.imag)\n",
        "axs[2].set_title('fft [imaginary]')\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Square 1D Signal"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 14,
      "metadata": {},
      "outputs": [],
      "source": [
        "# ideal signal\n",
        "signal_length = 200\n",
        "ideal_signal = np.zeros(shape=(signal_length,))\n",
        "ideal_signal[:signal_length // 10] = 1\n",
        "\n",
        "# fft\n",
        "fft_ideal_signal = np.fft.fft(ideal_signal)\n",
        "\n",
        "# shift the low frequencies to the center\n",
        "fftshift_ideal_signal = np.fft.fftshift(fft_ideal_signal)\n",
        "\n",
        "# magnitude of the signal\n",
        "mag_fft_ideal_signal = np.abs(fft_ideal_signal)\n",
        "mag_fftshift_ideal_signal = np.abs(fftshift_ideal_signal)\n",
        "\n",
        "# phase of the signal\n",
        "phase_fft_ideal_signal = np.angle(fft_ideal_signal)\n",
        "phase_fftshift_ideal_signal = np.angle(fftshift_ideal_signal)\n",
        "\n",
        "# reconstruct the signal\n",
        "reconstructed_ideal_signal = np.fft.ifft(np.fft.ifftshift(fftshift_ideal_signal)).real"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# plot\n",
        "fig, axs = plt.subplots(nrows=1, ncols=3, figsize=(16, 4), layout='compressed')\n",
        "\n",
        "axs[0].plot(ideal_signal)\n",
        "axs[0].set_title('Original signal')\n",
        "axs[1].plot(fftshift_ideal_signal.real)\n",
        "axs[1].set_title('fft [real]')\n",
        "axs[2].plot(fftshift_ideal_signal.imag)\n",
        "axs[2].set_title('fft [imaginary]')\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# plot\n",
        "fig, axs = plt.subplots(nrows=1, ncols=3, figsize=(16, 4), layout='compressed')\n",
        "\n",
        "axs[0].plot(ideal_signal)\n",
        "axs[0].set_title('Original signal')\n",
        "axs[1].plot(mag_fft_ideal_signal)\n",
        "axs[1].set_title('Magnitude')\n",
        "axs[2].plot(mag_fftshift_ideal_signal)\n",
        "axs[2].set_title('Magnitude [fftshift]')\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# plot\n",
        "fig, axs = plt.subplots(nrows=1, ncols=3, figsize=(16, 4), layout='compressed')\n",
        "\n",
        "axs[0].plot(ideal_signal)\n",
        "axs[0].set_title('Original signal')\n",
        "axs[1].plot(mag_fftshift_ideal_signal)\n",
        "axs[1].set_title('Magnitude [shifted version]')\n",
        "axs[2].plot(phase_fftshift_ideal_signal)\n",
        "axs[2].set_title('Phase [shifted version]')\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# plot\n",
        "fig, axs = plt.subplots(nrows=1, ncols=3, figsize=(16, 4), layout='compressed')\n",
        "\n",
        "axs[0].plot(ideal_signal)\n",
        "axs[0].set_title('Original signal')\n",
        "axs[1].plot(mag_fftshift_ideal_signal)\n",
        "axs[1].set_title('Magnitude [shifted version]')\n",
        "axs[2].plot(reconstructed_ideal_signal)\n",
        "axs[2].set_title('Reconstructed signal')\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 19,
      "metadata": {},
      "outputs": [],
      "source": [
        "# sinc signal\n",
        "signal_length = 200\n",
        "t = np.linspace(-20, 20, signal_length)\n",
        "sinc_signal = np.sinc(t)\n",
        "\n",
        "# fft\n",
        "fft_sinc_signal = np.fft.fft(sinc_signal)\n",
        "\n",
        "# shift the low frequencies to the center\n",
        "fftshift_sinc_signal = np.fft.fftshift(fft_sinc_signal)\n",
        "\n",
        "# magnitude of the signal\n",
        "mag_fftshift_sinc_signal = np.abs(fftshift_sinc_signal)\n",
        "\n",
        "# phase of the signal\n",
        "phase_fftshift_sinc_signal = np.angle(fftshift_sinc_signal)\n",
        "\n",
        "# reconstruct the signal\n",
        "reconstructed_sinc_signal = np.fft.ifft(fft_sinc_signal).real"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# plot\n",
        "fig, axs = plt.subplots(nrows=1, ncols=3, figsize=(16, 4), layout='compressed')\n",
        "\n",
        "axs[0].plot(sinc_signal)\n",
        "axs[0].set_title('Original signal')\n",
        "axs[1].plot(mag_fftshift_sinc_signal)\n",
        "axs[1].set_title('Magnitude [shifted version]')\n",
        "axs[2].plot(reconstructed_sinc_signal)\n",
        "axs[2].set_title('Reconstructed signal')\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 2D Block signal"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# create a 2d square ideal low-pass filter\n",
        "image_size = (256, 256)\n",
        "block_size = 8\n",
        "block_2d = ideal_filter(image_size, block_size, shape='square')\n",
        "\n",
        "# apply 2d fast-fourier transform\n",
        "fft_block_2d = np.fft.fft2(block_2d)\n",
        "fftshift_block_2d = np.fft.fftshift(fft_block_2d)\n",
        "\n",
        "# magnitude & phase\n",
        "abs_fft_block_2d = np.abs(fft_block_2d)\n",
        "phase_fft_block_2d = np.angle(fft_block_2d)\n",
        "\n",
        "# magnitude & phase [shifted version]\n",
        "abs_fftshift_block_2d = np.abs(fftshift_block_2d)\n",
        "phase_fftshift_block_2d = np.angle(fftshift_block_2d)\n",
        "\n",
        "# reconstruct the image\n",
        "ifft_block_2d = np.fft.ifft2(fft_block_2d).real.clip(0, 1)\n",
        "\n",
        "# plot\n",
        "fig, axs = plt.subplots(nrows=2, ncols=3, figsize=(12, 6), layout='compressed')\n",
        "\n",
        "axs[0, 0].imshow(block_2d, cmap='gray')\n",
        "axs[0, 0].set_title(\"Original image\")\n",
        "axs[0, 0].axis('off')\n",
        "axs[0, 1].imshow(np.log2(abs_fft_block_2d + 1), cmap='gray')\n",
        "axs[0, 1].set_title(\"log2(magnitude)\")\n",
        "axs[0, 1].axis('off')\n",
        "axs[0, 2].imshow(phase_fft_block_2d, cmap='gray')\n",
        "axs[0, 2].set_title(\"phase\")\n",
        "axs[0, 2].axis('off')\n",
        "axs[1, 0].imshow(ifft_block_2d, cmap='gray')\n",
        "axs[1, 0].set_title(\"Reconstructed image\")\n",
        "axs[1, 0].axis('off')\n",
        "axs[1, 1].imshow(np.log2(abs_fftshift_block_2d + 1), cmap='gray')\n",
        "axs[1, 1].set_title(\"log2(magnitude) [shift]\")\n",
        "axs[1, 1].axis('off')\n",
        "axs[1, 2].imshow(phase_fftshift_block_2d, cmap='gray')\n",
        "axs[1, 2].set_title(\"phase [shift]\")\n",
        "axs[1, 2].axis('off')\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Ideal low-pass & cardinal sine (sinc) filters"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 22,
      "metadata": {},
      "outputs": [],
      "source": [
        "# create a 2d circular ideal low-pass filter\n",
        "image_size = (256, 256)\n",
        "cutoff_frequency = 16\n",
        "ilpf = ideal_filter(image_size, cutoff_frequency)\n",
        "\n",
        "# apply 2d fast-fourier transform\n",
        "fft_ilpf = np.fft.fft2(ilpf)\n",
        "fftshift_ilpf = np.fft.fftshift(fft_ilpf)\n",
        "\n",
        "# magnitude & phase\n",
        "abs_fft_ilpf = np.abs(fft_ilpf)\n",
        "phase_fft_ilpf = np.angle(fft_ilpf)\n",
        "\n",
        "# magnitude & phase [shifted version]\n",
        "abs_fftshift_ilpf = np.abs(fftshift_ilpf)\n",
        "phase_fftshift_ilpf = np.angle(fftshift_ilpf)\n",
        "\n",
        "# reconstruct the image\n",
        "ifft_ilpf = np.fft.ifft2(fft_ilpf).real.clip(0, 1)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# plot\n",
        "fig, axs = plt.subplots(nrows=1, ncols=3, figsize=(12, 3), layout='compressed')\n",
        "\n",
        "axs[0].imshow(ilpf, cmap='gray')\n",
        "axs[0].set_title(\"Original image\")\n",
        "axs[0].axis('off')\n",
        "axs[1].imshow(np.abs(fftshift_ilpf.real), cmap='gray')\n",
        "axs[1].set_title(\"fft [real]\")\n",
        "axs[1].axis('off')\n",
        "axs[2].imshow(np.abs(fftshift_ilpf.imag), cmap='gray')\n",
        "axs[2].set_title(\"fft [imaginary]\")\n",
        "axs[2].axis('off')\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# plot\n",
        "fig, axs = plt.subplots(nrows=2, ncols=3, figsize=(12, 6), layout='compressed')\n",
        "\n",
        "axs[0, 0].imshow(ilpf, cmap='gray')\n",
        "axs[0, 0].set_title(\"Original image\")\n",
        "axs[0, 0].axis('off')\n",
        "axs[0, 1].imshow(np.log2(abs_fft_ilpf + 1), cmap='gray')\n",
        "axs[0, 1].set_title(\"log2(magnitude)\")\n",
        "axs[0, 1].axis('off')\n",
        "axs[0, 2].imshow(phase_fft_ilpf, cmap='gray')\n",
        "axs[0, 2].set_title(\"phase\")\n",
        "axs[0, 2].axis('off')\n",
        "axs[1, 0].imshow(ifft_ilpf, cmap='gray')\n",
        "axs[1, 0].set_title(\"Reconstructed image\")\n",
        "axs[1, 0].axis('off')\n",
        "axs[1, 1].imshow(np.log2(abs_fftshift_ilpf + 1), cmap='gray')\n",
        "axs[1, 1].set_title(\"log2(magnitude) [shifted version]\")\n",
        "axs[1, 1].axis('off')\n",
        "axs[1, 2].imshow(phase_fftshift_ilpf, cmap='gray')\n",
        "axs[1, 2].set_title(\"phase [shifted version]\")\n",
        "axs[1, 2].axis('off')\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# create a 2d sinc filter\n",
        "image_size = (256, 256)\n",
        "sigma = 2\n",
        "sinc = sinc_filter(image_size, sigma)\n",
        "\n",
        "# apply 2d fast-fourier transform\n",
        "fft_sinc = np.fft.fft2(sinc)\n",
        "fftshift_sinc = np.fft.fftshift(fft_sinc)\n",
        "\n",
        "# magnitude & phase [shifted version]\n",
        "abs_fftshift_sinc = np.abs(fftshift_sinc)\n",
        "phase_fftshift_sinc = np.angle(fftshift_sinc)\n",
        "\n",
        "# reconstruct the image\n",
        "ifft_sinc = np.fft.ifft2(fft_sinc).real.clip(0, 1)\n",
        "\n",
        "# plot\n",
        "fig, axs = plt.subplots(nrows=1, ncols=3, figsize=(12, 3), layout='compressed')\n",
        "\n",
        "axs[0].imshow(sinc, cmap='gray')\n",
        "axs[0].set_title(\"Original image\")\n",
        "axs[0].axis('off')\n",
        "axs[1].imshow(np.log2(abs_fftshift_sinc + 1), cmap='gray')\n",
        "axs[1].set_title(\"magnitude [shifted version]\")\n",
        "axs[1].axis('off')\n",
        "axs[2].imshow(ifft_sinc, cmap='gray')\n",
        "axs[2].set_title(\"Reconstructed image\")\n",
        "axs[2].axis('off')\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Signal shifting effect on magnitude & phase"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# create two separaet vertical lines\n",
        "image_size = 256\n",
        "vertical_1 = np.zeros(shape=(image_size, image_size))\n",
        "vertical_1[:, int(1.5 * image_size // 4): int(2.5 * image_size // 4)] = 1\n",
        "\n",
        "vertical_2 = np.zeros(shape=(image_size, image_size))\n",
        "vertical_2[:, int(.5 * image_size // 4): int(1.5 * image_size // 4)] = 1\n",
        "\n",
        "# fft2\n",
        "fftshift_vertical_1 = np.fft.fftshift(np.fft.fft2(vertical_1))\n",
        "fftshift_vertical_2 = np.fft.fftshift(np.fft.fft2(vertical_2))\n",
        "\n",
        "# magnitude\n",
        "abs_fftshift_vertical_1 = np.abs(fftshift_vertical_1)\n",
        "abs_fftshift_vertical_2 = np.abs(fftshift_vertical_2)\n",
        "\n",
        "# phase\n",
        "phase_fftshift_vertical_1 = np.angle(fftshift_vertical_1)\n",
        "phase_fftshift_vertical_2 = np.angle(fftshift_vertical_2)\n",
        "\n",
        "# plot\n",
        "fig, axs = plt.subplots(nrows=2, ncols=3, figsize=(12, 6), layout='compressed')\n",
        "\n",
        "axs[0, 0].imshow(vertical_1, cmap='gray')\n",
        "axs[0, 0].axis('off')\n",
        "axs[0, 0].set_title(\"vertical_1\")\n",
        "axs[0, 1].imshow(np.log2(abs_fftshift_vertical_1 + 1), cmap='gray')\n",
        "axs[0, 1].axis('off')\n",
        "axs[0, 1].set_title(\"Magnitude_1\")\n",
        "axs[0, 2].imshow(phase_fftshift_vertical_1, cmap='gray')\n",
        "axs[0, 2].axis('off')\n",
        "axs[0, 2].set_title(\"Phase_1\")\n",
        "axs[1, 0].imshow(vertical_2, cmap='gray')\n",
        "axs[1, 0].axis('off')\n",
        "axs[1, 0].set_title(\"vertical_2\")\n",
        "axs[1, 1].imshow(np.log2(abs_fftshift_vertical_2 + 1), cmap='gray')\n",
        "axs[1, 1].axis('off')\n",
        "axs[1, 1].set_title(\"Magnitude_2\")\n",
        "axs[1, 2].imshow(phase_fftshift_vertical_2, cmap='gray')\n",
        "axs[1, 2].axis('off')\n",
        "axs[1, 2].set_title(\"Phase_2\")\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# create two separate horizontal lines\n",
        "image_size = 256\n",
        "horizontal_1 = np.zeros(shape=(image_size, image_size))\n",
        "horizontal_1[int(1.5 * image_size // 4): int(2.5 * image_size // 4), :] = 1\n",
        "\n",
        "horizontal_2 = np.zeros(shape=(image_size, image_size))\n",
        "horizontal_2[int(.5 * image_size // 4): int(1.5 * image_size // 4), :] = 1\n",
        "\n",
        "# fft2\n",
        "fftshift_horizontal_1 = np.fft.fftshift(np.fft.fft2(horizontal_1))\n",
        "fftshift_horizontal_2 = np.fft.fftshift(np.fft.fft2(horizontal_2))\n",
        "\n",
        "# magnitude\n",
        "abs_fftshift_horizontal_1 = np.abs(fftshift_horizontal_1)\n",
        "abs_fftshift_horizontal_2 = np.abs(fftshift_horizontal_2)\n",
        "\n",
        "# phase\n",
        "phase_fftshift_horizontal_1 = np.angle(fftshift_horizontal_1)\n",
        "phase_fftshift_horizontal_2 = np.angle(fftshift_horizontal_2)\n",
        "\n",
        "# plot\n",
        "fig, axs = plt.subplots(nrows=2, ncols=3, figsize=(12, 6), layout='compressed')\n",
        "\n",
        "axs[0, 0].imshow(horizontal_1, cmap='gray')\n",
        "axs[0, 0].axis('off')\n",
        "axs[0, 0].set_title(\"horizontal_1\")\n",
        "axs[0, 1].imshow(np.log2(abs_fftshift_horizontal_1 + 1), cmap='gray')\n",
        "axs[0, 1].axis('off')\n",
        "axs[0, 1].set_title(\"Magnitude_1\")\n",
        "axs[0, 2].imshow(phase_fftshift_horizontal_1, cmap='gray')\n",
        "axs[0, 2].axis('off')\n",
        "axs[0, 2].set_title(\"Phase_1\")\n",
        "axs[1, 0].imshow(horizontal_2, cmap='gray')\n",
        "axs[1, 0].axis('off')\n",
        "axs[1, 0].set_title(\"horizontal_2\")\n",
        "axs[1, 1].imshow(np.log2(abs_fftshift_horizontal_2 + 1), cmap='gray')\n",
        "axs[1, 1].axis('off')\n",
        "axs[1, 1].set_title(\"Magnitude_2\")\n",
        "axs[1, 2].imshow(phase_fftshift_horizontal_2, cmap='gray')\n",
        "axs[1, 2].axis('off')\n",
        "axs[1, 2].set_title(\"Phase_2\")\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Signal rotation effect on magnitude & phase"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# load lenna.jpg\n",
        "lenna_1 = plt.imread('../assets/images/misc/lenna.jpg')\n",
        "lenna_2 = np.rot90(lenna_1)\n",
        "\n",
        "# fft2\n",
        "fftshift_lenna_1 = np.fft.fftshift(np.fft.fft2(lenna_1))\n",
        "fftshift_lenna_2 = np.fft.fftshift(np.fft.fft2(lenna_2))\n",
        "\n",
        "# magnitude\n",
        "abs_fftshift_lenna_1 = np.abs(fftshift_lenna_1)\n",
        "abs_fftshift_lenna_2 = np.abs(fftshift_lenna_2)\n",
        "\n",
        "# phase\n",
        "phase_fftshift_lenna_1 = np.angle(fftshift_lenna_1)\n",
        "phase_fftshift_lenna_2 = np.angle(fftshift_lenna_2)\n",
        "\n",
        "# plot\n",
        "fig, axs = plt.subplots(nrows=2, ncols=3, figsize=(12, 6), layout='compressed')\n",
        "\n",
        "axs[0, 0].imshow(lenna_1, cmap='gray')\n",
        "axs[0, 0].axis('off')\n",
        "axs[0, 0].set_title(\"lenna_1\")\n",
        "axs[0, 1].imshow(np.log2(abs_fftshift_lenna_1 + 1), cmap='gray')\n",
        "axs[0, 1].axis('off')\n",
        "axs[0, 1].set_title(\"magnitude_1\")\n",
        "axs[0, 2].imshow(phase_fftshift_lenna_1, cmap='gray')\n",
        "axs[0, 2].axis('off')\n",
        "axs[0, 2].set_title(\"phase_1\")\n",
        "axs[1, 0].imshow(lenna_2, cmap='gray')\n",
        "axs[1, 0].axis('off')\n",
        "axs[1, 0].set_title(\"lenna_2\")\n",
        "axs[1, 1].imshow(np.log2(abs_fftshift_lenna_2 + 1), cmap='gray')\n",
        "axs[1, 1].axis('off')\n",
        "axs[1, 1].set_title(\"magnitude_2\")\n",
        "axs[1, 2].imshow(phase_fftshift_lenna_2, cmap='gray')\n",
        "axs[1, 2].axis('off')\n",
        "axs[1, 2].set_title(\"phase_2\")\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Signal flip effect on magnitude & phase"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# load lenna.jpg\n",
        "cm_1 = plt.imread(\"../assets/images/dip_3rd/CH02_Fig0222(b)(cameraman).tif\")\n",
        "cm_2 = np.fliplr(cm_1)\n",
        "\n",
        "# fft2\n",
        "fftshift_cm_1 = np.fft.fftshift(np.fft.fft2(cm_1))\n",
        "fftshift_cm_2 = np.fft.fftshift(np.fft.fft2(cm_2))\n",
        "\n",
        "# magnitude\n",
        "abs_fftshift_cm_1 = np.abs(fftshift_cm_1)\n",
        "abs_fftshift_cm_2 = np.abs(fftshift_cm_2)\n",
        "\n",
        "# phase\n",
        "phase_fftshift_cm_1 = np.angle(fftshift_cm_1)\n",
        "phase_fftshift_cm_2 = np.angle(fftshift_cm_2)\n",
        "\n",
        "# plot\n",
        "fig, axs = plt.subplots(nrows=2, ncols=3, figsize=(12, 6), layout='compressed')\n",
        "\n",
        "axs[0, 0].imshow(cm_1, cmap='gray')\n",
        "axs[0, 0].axis('off')\n",
        "axs[0, 0].set_title(\"cameraman_1\")\n",
        "axs[0, 1].imshow(np.log2(abs_fftshift_cm_1 + 1), cmap='gray')\n",
        "axs[0, 1].axis('off')\n",
        "axs[0, 1].set_title(\"magnitude_1\")\n",
        "axs[0, 2].imshow(phase_fftshift_cm_1, cmap='gray')\n",
        "axs[0, 2].axis('off')\n",
        "axs[0, 2].set_title(\"phase_1\")\n",
        "axs[1, 0].imshow(cm_2, cmap='gray')\n",
        "axs[1, 0].axis('off')\n",
        "axs[1, 0].set_title(\"cameraman_2\")\n",
        "axs[1, 1].imshow(np.log2(abs_fftshift_cm_2 + 1), cmap='gray')\n",
        "axs[1, 1].axis('off')\n",
        "axs[1, 1].set_title(\"magnitude_2\")\n",
        "axs[1, 2].imshow(phase_fftshift_cm_2, cmap='gray')\n",
        "axs[1, 2].axis('off')\n",
        "axs[1, 2].set_title(\"phase_2\")\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Energy equality of a signal in spatial & frequency domain"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# load lenna.jpg\n",
        "lenna = plt.imread('../assets/images/misc/lenna.jpg')\n",
        "\n",
        "# fft2\n",
        "fft_lenna = np.fft.fft2(lenna, norm='ortho')\n",
        "fftshift_lenna = np.fft.fftshift(fft_lenna)\n",
        "\n",
        "# reconstructed lenna\n",
        "reconstructed_lenna = np.fft.ifft2(fft_lenna, norm='ortho').real.clip(0, 255).astype(np.uint8)\n",
        "\n",
        "# enery calculation\n",
        "energy_1 = np.sum(np.abs(lenna.astype(np.float64))**2)\n",
        "energy_2 = np.sum(np.abs(fft_lenna) ** 2)\n",
        "energy_3 = np.sum(np.abs(reconstructed_lenna.astype(np.float64)) ** 2)\n",
        "\n",
        "# plot\n",
        "fig, axs = plt.subplots(nrows=1, ncols=3, figsize=(12, 3), layout='compressed')\n",
        "\n",
        "axs[0].imshow(lenna, cmap='gray')\n",
        "axs[0].axis('off')\n",
        "axs[0].set_title(f\"original [E= {round(energy_1)}]\")\n",
        "axs[1].imshow(np.log2(np.abs(fftshift_lenna) + 1), cmap='gray')\n",
        "axs[1].axis('off')\n",
        "axs[1].set_title(f\"frequency [E= {round(energy_2)}]\")\n",
        "axs[2].imshow(reconstructed_lenna, cmap='gray')\n",
        "axs[2].axis('off')\n",
        "axs[2].set_title(f\"reconstructed [E= {round(energy_3)}]\")\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Reconstruct from only magnitude / phase"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# load lenna.jpg\n",
        "lenna = plt.imread('../assets/images/misc/lenna.jpg')\n",
        "\n",
        "# fft2\n",
        "fftshift_lenna = np.fft.fftshift(np.fft.fft2(lenna))\n",
        "abs_fftshift_lenna = np.abs(fftshift_lenna)\n",
        "phase_fftshift_lenna = np.angle(fftshift_lenna)\n",
        "\n",
        "# reconstructed lenna [only magnitude]\n",
        "reconstructed_abs_lenna = np.fft.ifft2(np.fft.ifftshift(abs_fftshift_lenna)).real.clip(0, 255).astype(np.uint8)\n",
        "\n",
        "# reconstructed lenna [only phase]\n",
        "reconstructed_phase_lenna = np.fft.ifft2(np.fft.ifftshift(np.exp(1j * phase_fftshift_lenna))).real\n",
        "\n",
        "# plot\n",
        "fig, axs = plt.subplots(nrows=1, ncols=3, figsize=(16, 6), layout='compressed')\n",
        "\n",
        "axs[0].imshow(lenna, cmap='gray')\n",
        "axs[0].axis('off')\n",
        "axs[0].set_title(f\"original\")\n",
        "axs[1].imshow(reconstructed_abs_lenna, cmap='gray')\n",
        "axs[1].axis('off')\n",
        "axs[1].set_title(f\"reconstructed [only magnitude]\")\n",
        "axs[2].imshow(reconstructed_phase_lenna, cmap='gray')\n",
        "axs[2].axis('off')\n",
        "axs[2].set_title(f\"reconstructed [only phase]\")\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Ringing effect when using ideal filters"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# load cameraman.tif\n",
        "cm = plt.imread('../assets/images/dip_3rd/CH02_Fig0222(b)(cameraman).tif')\n",
        "\n",
        "# create several ideal low-pass filters\n",
        "ilpf_1 = ideal_filter(cm.shape, 32)\n",
        "ilpf_2 = ideal_filter(cm.shape, 16)\n",
        "ilpf_3 = ideal_filter(cm.shape, 8)\n",
        "\n",
        "# fft2\n",
        "fft_cm = np.fft.fftshift(np.fft.fft2(cm))\n",
        "\n",
        "# applying filters\n",
        "cm_ilpf_1 = fft_cm * ilpf_1\n",
        "cm_ilpf_2 = fft_cm * ilpf_2\n",
        "cm_ilpf_3 = fft_cm * ilpf_3\n",
        "\n",
        "# reconstruct cm\n",
        "ifft_cm_ilpf_1 = np.fft.ifft2(np.fft.ifftshift(cm_ilpf_1)).real.clip(0, 255).astype(np.uint8)\n",
        "ifft_cm_ilpf_2 = np.fft.ifft2(np.fft.ifftshift(cm_ilpf_2)).real.clip(0, 255).astype(np.uint8)\n",
        "ifft_cm_ilpf_3 = np.fft.ifft2(np.fft.ifftshift(cm_ilpf_3)).real.clip(0, 255).astype(np.uint8)\n",
        "\n",
        "# plot\n",
        "fig, axs = plt.subplots(nrows=4, ncols=3, figsize=(12, 16), layout='compressed')\n",
        "\n",
        "axs[0, 0].imshow(ilpf_1, cmap='gray')\n",
        "axs[0, 0].axis('off')\n",
        "axs[0, 0].set_title(\"ilpf_1\")\n",
        "axs[0, 1].imshow(ilpf_2, cmap='gray')\n",
        "axs[0, 1].axis('off')\n",
        "axs[0, 1].set_title(\"ilpf_2\")\n",
        "axs[0, 2].imshow(ilpf_3, cmap='gray')\n",
        "axs[0, 2].axis('off')\n",
        "axs[0, 2].set_title(\"ilpf_3\")\n",
        "axs[1, 0].imshow(np.log2(np.abs(cm_ilpf_1) + 1), cmap='gray')\n",
        "axs[1, 0].axis('off')\n",
        "axs[1, 0].set_title(\"Magnitude 1\")\n",
        "axs[1, 1].imshow(np.log2(np.abs(cm_ilpf_2) + 1), cmap='gray')\n",
        "axs[1, 1].axis('off')\n",
        "axs[1, 1].set_title(\"Magnitude 2\")\n",
        "axs[1, 2].imshow(np.log2(np.abs(cm_ilpf_3) + 1), cmap='gray')\n",
        "axs[1, 2].axis('off')\n",
        "axs[1, 2].set_title(\"Magnitude 3\")\n",
        "axs[2, 0].imshow(np.angle(cm_ilpf_1), cmap='gray')\n",
        "axs[2, 0].axis('off')\n",
        "axs[2, 0].set_title(\"Phase 1\")\n",
        "axs[2, 1].imshow(np.angle(cm_ilpf_2), cmap='gray')\n",
        "axs[2, 1].axis('off')\n",
        "axs[2, 1].set_title(\"Phase 2\")\n",
        "axs[2, 2].imshow(np.angle(cm_ilpf_3), cmap='gray')\n",
        "axs[2, 2].axis('off')\n",
        "axs[2, 2].set_title(\"Phase 3\")\n",
        "axs[3, 0].imshow(ifft_cm_ilpf_1, cmap='gray')\n",
        "axs[3, 0].axis('off')\n",
        "axs[3, 0].set_title(\"Reconstruction 1\")\n",
        "axs[3, 1].imshow(ifft_cm_ilpf_2, cmap='gray')\n",
        "axs[3, 1].axis('off')\n",
        "axs[3, 1].set_title(\"Reconstruction 2\")\n",
        "axs[3, 2].imshow(ifft_cm_ilpf_3, cmap='gray')\n",
        "axs[3, 2].axis('off')\n",
        "axs[3, 2].set_title(\"Reconstruction 3\")\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# load cameraman.tif\n",
        "cm = plt.imread('../assets/images/dip_3rd/CH02_Fig0222(b)(cameraman).tif')\n",
        "\n",
        "# create several Gaussian low-pass filters\n",
        "glpf_1 = gaussian_filter(cm.shape, 32)\n",
        "glpf_2 = gaussian_filter(cm.shape, 16)\n",
        "glpf_3 = gaussian_filter(cm.shape, 8)\n",
        "\n",
        "# fft2\n",
        "fft_cm = np.fft.fftshift(np.fft.fft2(cm))\n",
        "\n",
        "# applying filters\n",
        "cm_glpf_1 = fft_cm * glpf_1\n",
        "cm_glpf_2 = fft_cm * glpf_2\n",
        "cm_glpf_3 = fft_cm * glpf_3\n",
        "\n",
        "# reconstruct cm\n",
        "ifft_cm_glpf_1 = np.fft.ifft2(np.fft.ifftshift(cm_glpf_1)).real.clip(0, 255).astype(np.uint8)\n",
        "ifft_cm_glpf_2 = np.fft.ifft2(np.fft.ifftshift(cm_glpf_2)).real.clip(0, 255).astype(np.uint8)\n",
        "ifft_cm_glpf_3 = np.fft.ifft2(np.fft.ifftshift(cm_glpf_3)).real.clip(0, 255).astype(np.uint8)\n",
        "\n",
        "# plot\n",
        "fig, axs = plt.subplots(nrows=4, ncols=3, figsize=(12, 16), layout='compressed')\n",
        "\n",
        "axs[0, 0].imshow(glpf_1, cmap='gray')\n",
        "axs[0, 0].axis('off')\n",
        "axs[0, 0].set_title(\"glpf_1\")\n",
        "axs[0, 1].imshow(glpf_2, cmap='gray')\n",
        "axs[0, 1].axis('off')\n",
        "axs[0, 1].set_title(\"glpf_2\")\n",
        "axs[0, 2].imshow(glpf_3, cmap='gray')\n",
        "axs[0, 2].axis('off')\n",
        "axs[0, 2].set_title(\"glpf_3\")\n",
        "axs[1, 0].imshow(np.log2(np.abs(cm_glpf_1) + 1), cmap='gray')\n",
        "axs[1, 0].axis('off')\n",
        "axs[1, 0].set_title(\"Magnitude 1\")\n",
        "axs[1, 1].imshow(np.log2(np.abs(cm_glpf_2) + 1), cmap='gray')\n",
        "axs[1, 1].axis('off')\n",
        "axs[1, 1].set_title(\"Magnitude 2\")\n",
        "axs[1, 2].imshow(np.log2(np.abs(cm_glpf_3) + 1), cmap='gray')\n",
        "axs[1, 2].axis('off')\n",
        "axs[1, 2].set_title(\"Magnitude 3\")\n",
        "axs[2, 0].imshow(np.angle(cm_glpf_1), cmap='gray')\n",
        "axs[2, 0].axis('off')\n",
        "axs[2, 0].set_title(\"Phase 1\")\n",
        "axs[2, 1].imshow(np.angle(cm_glpf_2), cmap='gray')\n",
        "axs[2, 1].axis('off')\n",
        "axs[2, 1].set_title(\"Phase 2\")\n",
        "axs[2, 2].imshow(np.angle(cm_glpf_3), cmap='gray')\n",
        "axs[2, 2].axis('off')\n",
        "axs[2, 2].set_title(\"Phase 3\")\n",
        "axs[3, 0].imshow(ifft_cm_glpf_1, cmap='gray')\n",
        "axs[3, 0].axis('off')\n",
        "axs[3, 0].set_title(\"Reconstruction 1\")\n",
        "axs[3, 1].imshow(ifft_cm_glpf_2, cmap='gray')\n",
        "axs[3, 1].axis('off')\n",
        "axs[3, 1].set_title(\"Reconstruction 2\")\n",
        "axs[3, 2].imshow(ifft_cm_glpf_3, cmap='gray')\n",
        "axs[3, 2].axis('off')\n",
        "axs[3, 2].set_title(\"Reconstruction 3\")\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Image sharpening using gaussian high-pass filter"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# load cameraman.tif\n",
        "cm = plt.imread('../assets/images/dip_3rd/CH02_Fig0222(b)(cameraman).tif')\n",
        "\n",
        "# create several Gaussian low-pass filters\n",
        "ghpf_1 = 1 - gaussian_filter(cm.shape, 64)\n",
        "ghpf_2 = 1 - gaussian_filter(cm.shape, 32)\n",
        "ghpf_3 = 1 - gaussian_filter(cm.shape, 16)\n",
        "\n",
        "# fft2\n",
        "fft_cm = np.fft.fftshift(np.fft.fft2(cm))\n",
        "\n",
        "# applying filters\n",
        "cm_ghpf_1 = fft_cm * ghpf_1\n",
        "cm_ghpf_2 = fft_cm * ghpf_2\n",
        "cm_ghpf_3 = fft_cm * ghpf_3\n",
        "\n",
        "# reconstruct cm\n",
        "ifft_cm_ghpf_1 = np.fft.ifft2(np.fft.ifftshift(cm_ghpf_1)).real.clip(0, 255).astype(np.uint8)\n",
        "ifft_cm_ghpf_2 = np.fft.ifft2(np.fft.ifftshift(cm_ghpf_2)).real.clip(0, 255).astype(np.uint8)\n",
        "ifft_cm_ghpf_3 = np.fft.ifft2(np.fft.ifftshift(cm_ghpf_3)).real.clip(0, 255).astype(np.uint8)\n",
        "\n",
        "# plot\n",
        "fig, axs = plt.subplots(nrows=2, ncols=3, figsize=(12, 6), layout='compressed')\n",
        "\n",
        "axs[0, 0].imshow(ghpf_1, cmap='gray')\n",
        "axs[0, 0].axis('off')\n",
        "axs[0, 0].set_title(\"ghpf_1\")\n",
        "axs[0, 1].imshow(ghpf_2, cmap='gray')\n",
        "axs[0, 1].axis('off')\n",
        "axs[0, 1].set_title(\"ghpf_2\")\n",
        "axs[0, 2].imshow(ghpf_3, cmap='gray')\n",
        "axs[0, 2].axis('off')\n",
        "axs[0, 2].set_title(\"ghpf_3\")\n",
        "axs[1, 0].imshow(ifft_cm_ghpf_1, cmap='gray')\n",
        "axs[1, 0].axis('off')\n",
        "axs[1, 0].set_title(\"Reconstruction 1\")\n",
        "axs[1, 1].imshow(ifft_cm_ghpf_2, cmap='gray')\n",
        "axs[1, 1].axis('off')\n",
        "axs[1, 1].set_title(\"Reconstruction 2\")\n",
        "axs[1, 2].imshow(ifft_cm_ghpf_3, cmap='gray')\n",
        "axs[1, 2].axis('off')\n",
        "axs[1, 2].set_title(\"Reconstruction 3\")\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# sharpening cameraman\n",
        "sharp_cm_1 = (cm.astype(np.float32) + ifft_cm_ghpf_1.astype(np.float32)).clip(0, 255).astype(np.uint8)\n",
        "sharp_cm_2 = (cm.astype(np.float32) + ifft_cm_ghpf_2.astype(np.float32)).clip(0, 255).astype(np.uint8)\n",
        "sharp_cm_3 = (cm.astype(np.float32) + ifft_cm_ghpf_3.astype(np.float32)).clip(0, 255).astype(np.uint8)\n",
        "\n",
        "# plot\n",
        "fig, axs = plt.subplots(nrows=1, ncols=4, figsize=(16, 4), layout='compressed')\n",
        "\n",
        "axs[0].imshow(cm, cmap='gray')\n",
        "axs[0].axis('off')\n",
        "axs[0].set_title(\"Original\")\n",
        "axs[1].imshow(sharp_cm_1, cmap='gray')\n",
        "axs[1].axis('off')\n",
        "axs[1].set_title(\"Sharpening 1\")\n",
        "axs[2].imshow(sharp_cm_2, cmap='gray')\n",
        "axs[2].axis('off')\n",
        "axs[2].set_title(\"Sharpening 2\")\n",
        "axs[3].imshow(sharp_cm_3, cmap='gray')\n",
        "axs[3].axis('off')\n",
        "axs[3].set_title(\"Sharpening 3\")\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Periodic noise removal"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# load cameraman.tif\n",
        "lenna = plt.imread('../assets/images/misc/lenna.jpg').astype(np.float64)\n",
        "\n",
        "# add priodic noise to lenna\n",
        "v, h = np.meshgrid(np.arange(lenna.shape[0]), np.arange(lenna.shape[1]))\n",
        "v_period = 16\n",
        "h_period = 32\n",
        "\n",
        "noisy_lenna = lenna + 20 * np.sin(2 * np.pi / v_period * v) + 50 * np.cos(2 * np.pi / h_period * h)\n",
        "\n",
        "# fft2\n",
        "fft_lenna = np.fft.fftshift(np.fft.fft2(noisy_lenna))\n",
        "\n",
        "# plot\n",
        "fig = plt.figure(figsize=(12, 8), layout='constrained')\n",
        "\n",
        "gs = GridSpec(2, 3, figure=fig)\n",
        "ax1 = fig.add_subplot(gs[0, 0])\n",
        "ax2 = fig.add_subplot(gs[1, 0])\n",
        "ax3 = fig.add_subplot(gs[:, 1:])\n",
        "\n",
        "ax1.imshow(lenna, cmap='gray')\n",
        "ax1.set_title(\"Original\")\n",
        "ax2.imshow(noisy_lenna, cmap='gray')\n",
        "ax2.set_title(\"Periodic noise\")\n",
        "ax3.imshow(np.log2(np.abs(fft_lenna) + 1)[206:306, 206:306], cmap='gray', vmin=0)\n",
        "ax3.add_artist(plt.Circle((18, 50), 3, color='red', fill=False))\n",
        "ax3.add_artist(plt.Circle((82, 50), 3, color='red', fill=False))\n",
        "ax3.add_artist(plt.Circle((50, 34), 3, color='red', fill=False))\n",
        "ax3.add_artist(plt.Circle((50, 66), 3, color='red', fill=False))\n",
        "ax3.set_title(\"Magnitude of noisy image [zoomed in center [206:306, 206:306]]\")\n",
        "\n",
        "\n",
        "for ax in fig.axes:\n",
        "    ax.axis('off')\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# destroy some frequencies\n",
        "fft_lenna[256 - 16, 256] = 0\n",
        "fft_lenna[256 + 16, 256] = 0\n",
        "fft_lenna[256, 256 - 32] = 0\n",
        "fft_lenna[256, 256 + 32] = 0\n",
        "\n",
        "# reconstruct lenna\n",
        "ifft_lenna = np.fft.ifft2(np.fft.ifftshift(fft_lenna)).real.clip(0, 255).astype(np.uint8)\n",
        "\n",
        "# plot\n",
        "fig = plt.figure(figsize=(12, 8), layout='constrained')\n",
        "\n",
        "gs = GridSpec(2, 3, figure=fig)\n",
        "ax1 = fig.add_subplot(gs[0, 0])\n",
        "ax2 = fig.add_subplot(gs[1, 0])\n",
        "ax3 = fig.add_subplot(gs[:, 1:])\n",
        "\n",
        "ax1.imshow(noisy_lenna, cmap='gray')\n",
        "ax1.set_title(\"Noisy lenna\")\n",
        "ax2.imshow(ifft_lenna, cmap='gray')\n",
        "ax2.set_title(\"Denoised lenna\")\n",
        "ax3.imshow(np.log2(np.abs(fft_lenna) + 1)[206:306, 206:306], cmap='gray')\n",
        "ax3.add_artist(plt.Circle((18, 50), 3, color='red', fill=False))\n",
        "ax3.add_artist(plt.Circle((82, 50), 3, color='red', fill=False))\n",
        "ax3.add_artist(plt.Circle((50, 34), 3, color='red', fill=False))\n",
        "ax3.add_artist(plt.Circle((50, 66), 3, color='red', fill=False))\n",
        "ax3.set_title(\"Magnitude of noisy image [zoomed in center [206:306, 206:306]]\")\n",
        "\n",
        "for ax in fig.axes:\n",
        "    ax.axis('off')\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## FFT instead of Circular Convolution"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# load cameraman.tif\n",
        "cm = plt.imread('../assets/images/dip_3rd/CH02_Fig0222(b)(cameraman).tif')\n",
        "\n",
        "# create a gaussian noise\n",
        "gaussian_noise = np.random.normal(loc=0, scale=10, size=cm.shape)\n",
        "\n",
        "# apply gaussian noise to the image\n",
        "noisy_cm = (cm + gaussian_noise).clip(0, 255).astype(np.uint8)\n",
        "\n",
        "# plot\n",
        "fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(8, 4), layout='compressed')\n",
        "\n",
        "axs[0].imshow(cm, cmap='gray')\n",
        "axs[0].set_title('Original image')\n",
        "axs[1].imshow(noisy_cm, cmap='gray')\n",
        "axs[1].set_title('Noisy image')\n",
        "\n",
        "for ax in fig.axes:\n",
        "    ax.axis('off')\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# smooth in spatial domain [using convolution]\n",
        "glpf = gaussian_filter(size=(5, 5), sigma=1, norm=True)\n",
        "\n",
        "# convolution\n",
        "denoised_cm_1 = sp.signal.convolve2d(noisy_cm, glpf, mode='same', boundary='fill')\n",
        "\n",
        "# clip and dtype\n",
        "denoised_cm_1 = denoised_cm_1.clip(0, 255).astype(np.uint8)\n",
        "\n",
        "# plot\n",
        "fig, axs = plt.subplots(nrows=1, ncols=3, figsize=(16, 8), layout='compressed')\n",
        "\n",
        "axs[0].imshow(cm, cmap='gray')\n",
        "axs[0].set_title('Original image')\n",
        "axs[1].imshow(noisy_cm, cmap='gray')\n",
        "axs[1].set_title('Noisy image')\n",
        "axs[2].imshow(denoised_cm_1, cmap='gray')\n",
        "axs[2].set_title('Denoised image')\n",
        "\n",
        "for ax in fig.axes:\n",
        "    ax.axis('off')\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# smooth in frequency domain [using fft instead of circular convolution]\n",
        "\n",
        "# image size : MxM, filter size : NxN -> upscale to (M+N-1, M+N-1)\n",
        "size = (cm.shape[0] + glpf.shape[1] - 1, cm.shape[0] + glpf.shape[1] - 1)\n",
        "\n",
        "# pad noisy image to size M+N-1\n",
        "padded_noisy_cm = np.zeros(shape=size)\n",
        "padded_noisy_cm[:cm.shape[0], :cm.shape[1]] = noisy_cm\n",
        "\n",
        "# pad filter to size M+N-1\n",
        "padded_glpf = np.zeros(shape=size)\n",
        "padded_glpf[:glpf.shape[0], :glpf.shape[1]] = glpf\n",
        "\n",
        "# plot\n",
        "fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(8, 4), layout='compressed')\n",
        "\n",
        "axs[0].imshow(padded_noisy_cm, cmap='gray')\n",
        "axs[0].set_title('padded noisy image')\n",
        "axs[1].imshow(padded_glpf, cmap='gray')\n",
        "axs[1].set_title('padded filter')\n",
        "axs[1].add_artist(plt.arrow(40, 40, -25, -25, color='red', head_width=10, head_length=10))\n",
        "\n",
        "for ax in fig.axes:\n",
        "    ax.axis('off')\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# fft\n",
        "fft_padded_noisy_cm = np.fft.fftshift(np.fft.fft2(padded_noisy_cm))\n",
        "abs_fft_padded_noisy_cm = np.abs(fft_padded_noisy_cm)\n",
        "\n",
        "fft_padded_glpf = np.fft.fftshift(np.fft.fft2(padded_glpf))\n",
        "abs_fft_padded_glpf = np.abs(fft_padded_glpf)\n",
        "\n",
        "fft_multiplication = np.multiply(fft_padded_noisy_cm, fft_padded_glpf)\n",
        "abs_fft_multiplication = np.abs(fft_multiplication)\n",
        "\n",
        "# plot\n",
        "fig, axs = plt.subplots(nrows=1, ncols=3, figsize=(16, 8), layout='compressed')\n",
        "\n",
        "axs[0].imshow(np.log2(abs_fft_padded_noisy_cm + 1), cmap='gray')\n",
        "axs[0].set_title('Magnitude(noisy_cm)')\n",
        "axs[1].imshow(np.log2(abs_fft_padded_glpf + 1), cmap='gray')\n",
        "axs[1].set_title('Magnitude(filter)')\n",
        "axs[2].imshow(np.log2(abs_fft_multiplication + 1), cmap='gray')\n",
        "axs[2].set_title('Magnitude(noisy_cm * filter)')\n",
        "\n",
        "for ax in fig.axes:\n",
        "    ax.axis('off')\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# reconstruction\n",
        "ifft_noisy_cm = np.fft.ifft2(np.fft.ifftshift(fft_multiplication)).real\n",
        "\n",
        "# downscale to size: cm.shape\n",
        "ifft_noisy_cm = ifft_noisy_cm[glpf.shape[0] // 2:cm.shape[0] + glpf.shape[0] // 2, glpf.shape[1] // 2:cm.shape[1] + glpf.shape[1] // 2]\n",
        "\n",
        "# clip & dtype\n",
        "denoised_cm_2 = ifft_noisy_cm.clip(0, 255).astype(np.uint8)\n",
        "\n",
        "# plot\n",
        "fig, axs = plt.subplots(nrows=1, ncols=3, figsize=(16, 8), layout='compressed')\n",
        "\n",
        "axs[0].imshow(noisy_cm, cmap='gray')\n",
        "axs[0].set_title('Noisy image')\n",
        "axs[1].imshow(denoised_cm_1, cmap='gray')\n",
        "axs[1].set_title('Denoised [convolution]')\n",
        "axs[2].imshow(denoised_cm_2, cmap='gray')\n",
        "axs[2].set_title('Denoised [fft]')\n",
        "\n",
        "for ax in fig.axes:\n",
        "    ax.axis('off')\n",
        "\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.12.3"
    },
    "author_name": "Amirhossein Heydari",
    "author_email": "AmirhosseinHeydari78@gmail.com",
    "author_github": "https://github.com/mr-pylin",
    "origin_repo": "https://github.com/mr-pylin/media-processing-workshop"
  },
  "nbformat": 4,
  "nbformat_minor": 2
}