{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "\n", "import matplotlib\n", "matplotlib.rcParams['figure.figsize'] = (6, 6)\n", "\n", "import math\n", "import cmath # math functions for complex numbers\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "import ipywidgets\n", "from ipywidgets import interact" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**TODO**:\n", "- ajouter notation simplifiée\n", " - fourier-transform-introduction\n", " - fourier-transform-slides\n", " - slides de Bologne\n", "- ajouter un plot interactif de f(t) (fonction 1D) pour les vecteurs (numpy) de coefs a et b donnés\n", "- ajouter un plot interactif pour calculer les vecteurs a et b d'une fonction 1D donnée\n", "- finir notation JL\n", " - sine transform\n", " - cosine transform\n", "- implementations python avec scipy" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The Fourier Transform" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following notations come from the following book: *Toutes les mathématiques et les bases de l'informatique*, Horst Stöcker (Dunod, 2002)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "TODO..." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def plot(a, b):\n", " plt.axis('equal')\n", " \n", " plt.axis([-2.*np.pi, 2.*np.pi, -2.*np.pi, 2.*np.pi])\n", " plt.xticks(np.arange(-2.*np.pi, 2.*np.pi, 1))\n", " plt.yticks(np.arange(-2.*np.pi, 2.*np.pi, 1))\n", " \n", " plt.axhline(y=0, color='k')\n", " plt.axvline(x=0, color='k')\n", " \n", " t = np.arange(-math.pi, math.pi, 0.1)\n", " y = a * np.cos(t) + b * np.sin(t)\n", " plt.plot(t, y, \"-r\",\n", " t, a*np.cos(t), \"--k\",\n", " t, b*np.sin(t), \"--k\",\n", " a*np.cos(t), b*np.sin(t), \"--b\")\n", " \n", " plt.grid()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "@interact(a=(-5., 5.), b=(-5., 5.))\n", "def cos_sin(a, b):\n", " plot(a, b)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The Fast Fourier Transform (J.-L. Stark's notation)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following notations come from the following book: *Astronomical Image and Data Analysis*, J.-L. Stark and F. Murtagh (Springer, 2002)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### One dimension case" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Fourier transform of a continuous function f(t) is defined by:\n", "\n", "$$\n", "\\hat{f}(\\nu) = \\int_{-\\infty}^{+\\infty} f(t) e^{-i 2\\pi\\nu t} dt\n", "$$\n", "\n", "the inverse Fourier transform is:\n", "\n", "$$\n", "f(t) = \\int_{-\\infty}^{+\\infty} \\hat{f}(\\nu) e^{i 2\\pi\\nu t} du\n", "$$\n", "\n", "The discrete Fourier transform is givent by:\n", "\n", "$$\n", "\\hat{f}(u) = \\frac{1}{N} \\sum_{k=-\\infty}^{+\\infty} f(k) e^{-i 2\\pi \\frac{uk}{N}}\n", "$$\n", "\n", "and the inverse discrete Fourier transform is:\n", "\n", "$$\n", "f(k) = \\sum_{u=-\\infty}^{+\\infty} \\hat{f}(u) e^{i 2\\pi \\frac{uk}{N}}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Two dimention case" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In case of images (2 variables), the Fourier Transform is:\n", "\n", "$$\n", "\\hat{f}(u, v) = \\frac{1}{MN} \\sum_{l=-\\infty}^{+\\infty} \\sum_{k=-\\infty}^{+\\infty} f(k,l) e^{-i 2\\pi \\left( \\frac{uk}{M} + \\frac{vl}{N} \\right) }\n", "$$\n", "\n", "and the inverse discrete Fourier transform is:\n", "\n", "$$\n", "f(k,l) = \\sum_{u=-\\infty}^{+\\infty} \\sum_{v=-\\infty}^{+\\infty} \\hat{f}(u,v) e^{i 2\\pi \\left( \\frac{uk}{M} + \\frac{vl}{N} \\right) }\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Complex notation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since $\\hat{f}(u, v)$ is generally complex, this can be written using the complex its real and imaginary parts:\n", "\n", "$$\n", "\\hat{f}(u, v) = Re[\\hat{f}(u, v)] + i Im[\\hat{f}(u, v)]\n", "$$\n", "\n", "with:\n", "\n", "$$\n", "Re[\\hat{f}(u, v)] = \\frac{1}{MN} \\sum_{l=-\\infty}^{+\\infty} \\sum_{k=-\\infty}^{+\\infty} f(k,l) \\cos\\left( 2\\pi \\left( \\frac{uk}{M} + \\frac{vl}{N} \\right) \\right)\n", "$$\n", "\n", "$$\n", "Im[\\hat{f}(u, v)] = \\frac{-1}{MN} \\sum_{l=-\\infty}^{+\\infty} \\sum_{k=-\\infty}^{+\\infty} f(k,l) \\sin\\left( 2\\pi \\left( \\frac{uk}{M} + \\frac{vl}{N} \\right) \\right)\n", "$$\n", "\n", "It can also be written using modulus and argument:\n", "\n", "$$\n", "\\hat{f}(u, v) = |\\hat{f}(u, v)| + e^{i ~ \\text{arg} ~ \\hat{f}(u, v)}\n", "$$\n", "\n", "$|\\hat{f}(u, v)|^2$ is called the power spectrum, and $\\Theta(u,v) = \\text{arg} ~ \\hat{f}(u, v)$ the phase." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Fourier Transform with Python" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import math\n", "\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from matplotlib import cm\n", "\n", "import PIL.Image as pil_img # PIL.Image is a module not a class..." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# %load %load https://raw.githubusercontent.com/jeremiedecock/snippets/master/python/numpy/fft_transform/fft2.py\n", "\n", "# Fast Fourier Transform snippet\n", "#\n", "# Additional documentation:\n", "# - Numpy implementation: http://docs.scipy.org/doc/numpy/reference/routines.fft.html\n", "# - Scipy implementation: http://docs.scipy.org/doc/scipy/reference/fftpack.html\n", "\n", "# SET OPTIONS ########################################################\n", "\n", "shift = True # Shift the zero to the center\n", "threshold = 0.0001 # The threshold value (between 0 and 1)\n", "file_path = \"./sample-images/julie_lebrun.jpeg\" # The file image to filter\n", "\n", "# GET DATA ###########################################################\n", "\n", "# Open the image and convert it to grayscale\n", "signal = np.array(pil_img.open(file_path).convert('L'))\n", "\n", "# Init plot\n", "fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(14, 8))\n", "\n", "# Plot\n", "ax1.imshow(signal, interpolation='nearest', cmap=cm.gray)\n", "ax1.set_title(\"Original image\")\n", "\n", "\n", "# FOURIER TRANSFORM WITH NUMPY #######################################\n", "\n", "# Do the fourier transform #############\n", "\n", "transformed_signal = np.fft.fft2(signal)\n", "\n", "if shift:\n", " transformed_signal = np.fft.fftshift(transformed_signal)\n", "\n", "ax2.imshow(np.log10(abs(transformed_signal)),\n", " interpolation='nearest',\n", " cmap=cm.gray)\n", "ax2.set_title(\"Fourier coefficients before filtering\")\n", "\n", "\n", "# Filter ###############################\n", "\n", "max_value = np.max(abs(transformed_signal))\n", "filtered_transformed_signal = transformed_signal * (abs(transformed_signal) > max_value*threshold)\n", "\n", "ax3.imshow(np.log10(abs(filtered_transformed_signal)),\n", " interpolation='nearest',\n", " cmap=cm.gray)\n", "ax3.set_title(\"Fourier coefficients after filtering\")\n", "\n", "\n", "# Do the reverse transform #############\n", "\n", "if shift:\n", " filtered_transformed_signal = np.fft.ifftshift(filtered_transformed_signal)\n", "\n", "filtered_signal = np.fft.ifft2(filtered_transformed_signal)\n", "\n", "ax4.imshow(abs(filtered_signal), interpolation='nearest', cmap=cm.gray)\n", "ax4.set_title(\"Filtered image\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# %load %load https://raw.githubusercontent.com/jeremiedecock/snippets/master/python/numpy/fft_transform/fft2_basic_test.py\n", "\n", "# Fast Fourier Transform snippets\n", "# \n", "# Documentation:\n", "# - Numpy implementation: http://docs.scipy.org/doc/numpy/reference/routines.fft.html\n", "# - Scipy implementation: http://docs.scipy.org/doc/scipy/reference/fftpack.html\n", "\n", "fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(14, 8))\n", "\n", "# MAKE DATA ##########################################################\n", "\n", "pattern = np.zeros((4, 4))\n", "pattern[1:3,1:3] = 1\n", "\n", "signal = np.tile(pattern, (2, 2))\n", "\n", "ax1.imshow(signal, interpolation='nearest', cmap=cm.gray)\n", "ax1.set_title(\"Original image\")\n", "\n", "\n", "# FOURIER TRANSFORM WITH NUMPY #######################################\n", " \n", "\n", "# Do the fourier transform #############\n", "\n", "transformed_signal = np.fft.fft2(signal)\n", "\n", "ax2.imshow(abs(transformed_signal), interpolation='nearest', cmap=cm.gray)\n", "ax2.set_title(\"Fourier coefficients before filtering\")\n", "\n", "#shifted_transformed_signal = np.fft.fftshift(transformed_signal)\n", "#shifted_filtered_signal = np.fft.ifftshift(transformed_signal)\n", "\n", "\n", "# Filter ###############################\n", "\n", "max_value = np.max(abs(transformed_signal))\n", "filtered_transformed_signal = transformed_signal * (abs(transformed_signal) > max_value*0.5)\n", "\n", "ax3.imshow(abs(filtered_transformed_signal), interpolation='nearest', cmap=cm.gray)\n", "ax3.set_title(\"Fourier coefficients after filtering\")\n", "\n", "\n", "# Do the reverse transform #############\n", "\n", "filtered_signal = np.fft.ifft2(filtered_transformed_signal)\n", "\n", "ax4.imshow(abs(filtered_signal), interpolation='nearest', cmap=cm.gray)\n", "ax4.set_title(\"Filtered image\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# %load %load https://raw.githubusercontent.com/jeremiedecock/snippets/master/python/numpy/fft_transform/fft2_with_noise.py\n", "\n", "# Fast Fourier Transform snippet\n", "# \n", "# Additional documentation:\n", "# - Numpy implementation: http://docs.scipy.org/doc/numpy/reference/routines.fft.html\n", "# - Scipy implementation: http://docs.scipy.org/doc/scipy/reference/fftpack.html\n", "\n", "# SET OPTIONS ########################################################\n", "\n", "shift = True # Shift the zero to the center\n", "threshold = 0.0001 # The threshold value (between 0 and 1)\n", "noise_coefficient = 0.5 # The noise coefficient (between 0 and 1)\n", "file_path = \"./sample-images/julie_lebrun.jpeg\" # The file image to filter\n", "\n", "# GET DATA ###########################################################\n", "\n", "# Open the image and convert it to grayscale\n", "img = np.array(pil_img.open(file_path).convert('L'))\n", "\n", "# Add noise\n", "# TODO: lets the user choose the noise method : uniform, normal, poisson, ...\n", "#noise = np.random.rand(*img.shape) * 255. * noise_coefficient\n", "noise = np.random.poisson(255. * noise_coefficient, size=img.shape)\n", "\n", "noised_img = img + noise\n", "\n", "# Init plots\n", "fig, ((ax1, ax4, ax7, ax10), (ax2, ax5, ax8, ax11), (ax3, ax6, ax9, ax12)) = plt.subplots(3, 4, figsize=(16, 10))\n", "\n", "# Display images\n", "ax1.imshow(img, interpolation='nearest', cmap=cm.gray)\n", "ax1.set_title(\"Original image\")\n", "\n", "ax2.imshow(noise, interpolation='nearest', cmap=cm.gray)\n", "ax2.set_title(\"Noise\")\n", "\n", "ax3.imshow(noised_img, interpolation='nearest', cmap=cm.gray)\n", "ax3.set_title(\"Noised image\")\n", "\n", "\n", "# FOURIER TRANSFORM WITH NUMPY #######################################\n", "\n", "# Do the fourier transform #############\n", "\n", "fourier_img = np.fft.fft2(img)\n", "fourier_noise = np.fft.fft2(noise)\n", "fourier_noised_img = np.fft.fft2(noised_img)\n", "\n", "if shift:\n", " fourier_img = np.fft.fftshift(fourier_img)\n", " fourier_noise = np.fft.fftshift(fourier_noise)\n", " fourier_noised_img = np.fft.fftshift(fourier_noised_img)\n", "\n", "ax4.imshow(np.log10(abs(fourier_img)),\n", " interpolation='nearest',\n", " cmap=cm.gray)\n", "ax4.set_title(\"Fourier coefficients\\nbefore filtering (image)\")\n", "\n", "ax5.imshow(np.log10(abs(fourier_noise)),\n", " interpolation='nearest',\n", " cmap=cm.gray)\n", "ax5.set_title(\"Fourier coefficients\\nbefore filtering (noise)\")\n", "\n", "ax6.imshow(np.log10(abs(fourier_noised_img)),\n", " interpolation='nearest',\n", " cmap=cm.gray)\n", "ax6.set_title(\"Fourier coefficients\\nbefore filtering (noised image)\")\n", "\n", "\n", "# Filter ###############################\n", "\n", "max_value = np.max(abs(fourier_img))\n", "\n", "# TODO: lets the user choose between a threshold based on max, mean, median, ... or a geometric mask (square or circle to the center)\n", "filtered_fourier_img = fourier_img * (abs(fourier_img) > max_value*threshold)\n", "filtered_fourier_noise = fourier_noise * (abs(fourier_noise) > max_value*threshold)\n", "filtered_fourier_noised_img = fourier_noised_img * (abs(fourier_noised_img) > max_value*threshold)\n", "\n", "ax7.imshow(np.log10(abs(filtered_fourier_img)),\n", " interpolation='nearest',\n", " cmap=cm.gray)\n", "ax7.set_title(\"Fourier coefficients\\nafter filtering (image)\")\n", "\n", "ax8.imshow(np.log10(abs(filtered_fourier_noise)),\n", " interpolation='nearest',\n", " cmap=cm.gray)\n", "ax8.set_title(\"Fourier coefficients\\nafter filtering (noise)\")\n", "\n", "ax9.imshow(np.log10(abs(filtered_fourier_noised_img)),\n", " interpolation='nearest',\n", " cmap=cm.gray)\n", "ax9.set_title(\"Fourier coefficients\\nafter filtering (noised image)\")\n", "\n", "\n", "# Do the reverse transform #############\n", "\n", "if shift:\n", " filtered_fourier_img = np.fft.ifftshift(filtered_fourier_img)\n", " filtered_fourier_noise = np.fft.ifftshift(filtered_fourier_noise)\n", " filtered_fourier_noised_img = np.fft.ifftshift(filtered_fourier_noised_img)\n", "\n", "filtered_img = np.fft.ifft2(filtered_fourier_img)\n", "filtered_noise = np.fft.ifft2(filtered_fourier_noise)\n", "filtered_noised_img = np.fft.ifft2(filtered_fourier_noised_img)\n", "\n", "ax10.imshow(abs(filtered_img), interpolation='nearest', cmap=cm.gray)\n", "ax10.set_title(\"Filtered image\")\n", "\n", "ax11.imshow(abs(filtered_noise), interpolation='nearest', cmap=cm.gray)\n", "ax11.set_title(\"Filtered noise\")\n", "\n", "ax12.imshow(abs(filtered_noised_img), interpolation='nearest', cmap=cm.gray)\n", "ax12.set_title(\"Filtered noised image\")\n", "\n", "\n", "# SAVE FILES ######################\n", "\n", "ax1.set_axis_off()\n", "ax2.set_axis_off()\n", "ax3.set_axis_off()\n", "ax4.set_axis_off()\n", "ax5.set_axis_off()\n", "ax6.set_axis_off()\n", "\n", "ax7.set_axis_off()\n", "ax8.set_axis_off()\n", "ax9.set_axis_off()\n", "\n", "ax10.set_axis_off()\n", "ax11.set_axis_off()\n", "ax12.set_axis_off()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# %load %load https://raw.githubusercontent.com/jeremiedecock/snippets/master/python/numpy/fft_transform/fft_with_noise.py\n", "\n", "# Fast Fourier Transform snippet\n", "# \n", "# Additional documentation:\n", "# - Numpy implementation: http://docs.scipy.org/doc/numpy/reference/routines.fft.html\n", "# - Scipy implementation: http://docs.scipy.org/doc/scipy/reference/fftpack.html\n", "\n", "# SET OPTIONS ########################################################\n", "\n", "shift = True # Shift the zero to the center\n", "threshold = 0.0001 # The threshold value (between 0 and 1)\n", "noise_coefficient = 0.5 # The noise coefficient (between 0 and 1)\n", "file_path = \"./sample-images/julie_lebrun.jpeg\" # The file image to filter\n", "\n", "# GET DATA ###########################################################\n", "\n", "# Make the signal\n", "t = np.arange(5. * -math.pi, 5. * math.pi, 0.01)\n", "#t = np.arange(0., 10. * math.pi, 0.01)\n", "#sig = np.sin(t) + 2. * np.sin(2. * t) + 4. * np.sin(3. * t)\n", "sig = np.cos(t) + 2. * np.cos(2. * t) + 4. * np.cos(3. * t)\n", "\n", "# Add noise\n", "# TODO: lets the user choose the noise method : uniform, normal, poisson, ...\n", "#noise = np.random.poisson(1. * noise_coefficient, size=sig.shape)\n", "noise = np.random.rand(*sig.shape) * noise_coefficient\n", "\n", "noised_sig = sig + noise\n", "\n", "# Init plots\n", "fig, ((ax1, ax4, ax7, ax10), (ax2, ax5, ax8, ax11), (ax3, ax6, ax9, ax12)) = plt.subplots(3, 4, figsize=(16, 10))\n", "\n", "# Display images\n", "ax1.plot(t, sig)\n", "ax1.set_title(\"Original signal\")\n", "\n", "ax2.plot(t, noise)\n", "ax2.set_title(\"Noise\")\n", "\n", "ax3.plot(t, noised_sig)\n", "ax3.set_title(\"Noised signal\")\n", "\n", "\n", "# FOURIER TRANSFORM WITH NUMPY #######################################\n", "\n", "# Do the fourier transform #############\n", "\n", "print(np.fft.rfft(sig))\n", "print(type(np.fft.rfft(sig)[0]))\n", "\n", "fourier_sig = np.fft.fft(sig)\n", "fourier_noise = np.fft.fft(noise)\n", "fourier_noised_sig = np.fft.fft(noised_sig)\n", "\n", "if shift:\n", " fourier_sig = np.fft.fftshift(fourier_sig)\n", " fourier_noise = np.fft.fftshift(fourier_noise)\n", " fourier_noised_sig = np.fft.fftshift(fourier_noised_sig)\n", "\n", "ax4.plot(t, abs(fourier_sig))\n", "ax4.set_title(\"Fourier coefficients\\nbefore filtering (signal)\")\n", "\n", "ax5.plot(t, abs(fourier_noise))\n", "ax5.set_title(\"Fourier coefficients\\nbefore filtering (noise)\")\n", "\n", "ax6.plot(t, abs(fourier_noised_sig))\n", "ax6.set_title(\"Fourier coefficients\\nbefore filtering (noised signal)\")\n", "\n", "\n", "# Filter ###############################\n", "\n", "max_value = np.max(abs(fourier_sig))\n", "\n", "# TODO: lets the user choose between a threshold based on max, mean, median, ... or a geometric mask (square or circle to the center)\n", "filtered_fourier_sig = fourier_sig * (abs(fourier_sig) > max_value*threshold)\n", "filtered_fourier_noise = fourier_noise * (abs(fourier_noise) > max_value*threshold)\n", "filtered_fourier_noised_sig = fourier_noised_sig * (abs(fourier_noised_sig) > max_value*threshold)\n", "\n", "ax7.plot(t, abs(filtered_fourier_sig))\n", "ax7.set_title(\"Fourier coefficients\\nafter filtering (signal)\")\n", "\n", "ax8.plot(t, abs(filtered_fourier_noise))\n", "ax8.set_title(\"Fourier coefficients\\nafter filtering (noise)\")\n", "\n", "ax9.plot(t, abs(filtered_fourier_noised_sig))\n", "ax9.set_title(\"Fourier coefficients\\nafter filtering (noised signal)\")\n", "\n", "\n", "# Do the reverse transform #############\n", "\n", "if shift:\n", " filtered_fourier_sig = np.fft.ifftshift(filtered_fourier_sig)\n", " filtered_fourier_noise = np.fft.ifftshift(filtered_fourier_noise)\n", " filtered_fourier_noised_sig = np.fft.ifftshift(filtered_fourier_noised_sig)\n", "\n", "filtered_sig = np.fft.ifft(filtered_fourier_sig)\n", "filtered_noise = np.fft.ifft(filtered_fourier_noise)\n", "filtered_noised_sig = np.fft.ifft(filtered_fourier_noised_sig)\n", "\n", "ax10.plot(t, abs(filtered_sig))\n", "ax10.set_title(\"Filtered signal\")\n", "\n", "ax11.plot(t, abs(filtered_noise))\n", "ax11.set_title(\"Filtered noise\")\n", "\n", "ax12.plot(t, abs(filtered_noised_sig))\n", "ax12.set_title(\"Filtered noised signal\")\n", "\n", "# SAVE FILES ######################\n", "\n", "#ax1.set_axis_off()\n", "#ax2.set_axis_off()\n", "#ax3.set_axis_off()\n", "#ax4.set_axis_off()\n", "#ax5.set_axis_off()\n", "#ax6.set_axis_off()\n", "\n", "#ax7.set_axis_off()\n", "#ax8.set_axis_off()\n", "#ax9.set_axis_off()\n", "\n", "#ax10.set_axis_off()\n", "#ax11.set_axis_off()\n", "#ax12.set_axis_off()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# %load %load https://raw.githubusercontent.com/jeremiedecock/snippets/master/python/numpy/fft_transform/plot.py\n", "\n", "# Fast Fourier Transform snippet\n", "# \n", "# Additional documentation:\n", "# - Numpy implementation: http://docs.scipy.org/doc/numpy/reference/routines.fft.html\n", "# - Scipy implementation: http://docs.scipy.org/doc/scipy/reference/fftpack.html\n", "\n", "# SET OPTIONS ########################################################\n", "\n", "shift = True # Shift the zero to the center\n", "threshold = 0.0001 # The threshold value (between 0 and 1)\n", "noise_coefficient = 0.5 # The noise coefficient (between 0 and 1)\n", "file_path = \"./sample-images/julie_lebrun.jpeg\" # The file image to filter\n", "\n", "# GET DATA ###########################################################\n", "\n", "# Make the signal\n", "t = np.arange(5. * -math.pi, 5. * math.pi, 0.01)\n", "#t = np.arange(0., 10. * math.pi, 0.01)\n", "#sig = np.sin(t) + 2. * np.sin(2. * t) + 4. * np.sin(3. * t)\n", "sig = np.cos(t) + 2. * np.cos(2. * t) + 4. * np.cos(3. * t)\n", "\n", "sig1 = np.cos(t)\n", "sig2 = 2. * np.cos(2. * t)\n", "sig3 = 4. * np.cos(3. * t)\n", "\n", "# Add noise\n", "# TODO: lets the user choose the noise method : uniform, normal, poisson, ...\n", "#noise = np.random.poisson(1. * noise_coefficient, size=sig.shape)\n", "noise = np.random.rand(*sig.shape) * noise_coefficient\n", "\n", "noised_sig = sig + noise\n", "\n", "## Init plots\n", "#fig, (ax1, ax2, ax3, ax4) = plt.subplots(1, 4, figsize=(16, 5))\n", "#\n", "#N = 8\n", "#\n", "## Display images\n", "#ax1.plot(t, sig)\n", "#ax1.set_title(r\"$f(t)$ = \")\n", "#ax1.set_ylim([-N, N])\n", "#\n", "#ax2.plot(t, sig1)\n", "#ax2.set_title(r\"$a_1 \\times \\cos(t)$ +\")\n", "#ax2.set_ylim([-N, N])\n", "#\n", "#ax3.plot(t, sig2)\n", "#ax3.set_title(r\"$a_2 \\times \\cos(2t)$ +\")\n", "#ax3.set_ylim([-N, N])\n", "#\n", "#ax4.plot(t, sig3)\n", "#ax4.set_title(r\"$a_3 \\times \\cos(3t)$\")\n", "#ax4.set_ylim([-N, N])\n", "\n", "# Init plots\n", "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 5))\n", "\n", "# Display images\n", "ax1.plot(t, sig)\n", "ax1.set_title(r\"Direct space\")\n", "#ax1.set_ylim([-N, N])\n", "\n", "ax2.bar([1, 2, 3], [1, 2, 4], 0.01)\n", "ax2.set_title(r\"Fourier space\")\n", "ax2.set_xlim([0, 4])\n", "ax2.set_ylim([0, 5])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# %load %load https://raw.githubusercontent.com/jeremiedecock/snippets/master/python/numpy/fft_transform/rfft2.py\n", "\n", "# Fast Fourier Transform snippet\n", "# \n", "# Additional documentation:\n", "# - Numpy implementation: http://docs.scipy.org/doc/numpy/reference/routines.fft.html\n", "# - Scipy implementation: http://docs.scipy.org/doc/scipy/reference/fftpack.html\n", "\n", "# SET OPTIONS ########################################################\n", "\n", "shift = True # Shift the zero to the center\n", "threshold = 0.0001 # The threshold value (between 0 and 1)\n", "file_path = \"./sample-images/julie_lebrun.jpeg\" # The file image to filter\n", "\n", "# GET DATA ###########################################################\n", "\n", "# Open the image and convert it to grayscale\n", "signal = np.array(pil_img.open(file_path).convert('L'))\n", "\n", "# Init plot\n", "fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(14, 8))\n", "\n", "# Plot\n", "ax1.imshow(signal, interpolation='nearest', cmap=cm.gray)\n", "ax1.set_title(\"Original image\")\n", "\n", "\n", "# FOURIER TRANSFORM WITH NUMPY #######################################\n", "\n", "# Do the fourier transform #############\n", "\n", "transformed_signal = np.fft.rfft2(signal)\n", "\n", "if shift:\n", " transformed_signal = np.fft.fftshift(transformed_signal)\n", "\n", "ax2.imshow(np.log10(abs(transformed_signal)),\n", " interpolation='nearest',\n", " cmap=cm.gray)\n", "ax2.set_title(\"Fourier coefficients before filtering\")\n", "\n", "\n", "# Filter ###############################\n", "\n", "max_value = np.max(abs(transformed_signal))\n", "filtered_transformed_signal = transformed_signal * (abs(transformed_signal) > max_value*threshold)\n", "\n", "ax3.imshow(np.log10(abs(filtered_transformed_signal)),\n", " interpolation='nearest',\n", " cmap=cm.gray)\n", "ax3.set_title(\"Fourier coefficients after filtering\")\n", "\n", "\n", "# Do the reverse transform #############\n", "\n", "if shift:\n", " filtered_transformed_signal = np.fft.ifftshift(filtered_transformed_signal)\n", "\n", "filtered_signal = np.fft.irfft2(filtered_transformed_signal)\n", "\n", "ax4.imshow(abs(filtered_signal), interpolation='nearest', cmap=cm.gray)\n", "ax4.set_title(\"Filtered image\")\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python [default]", "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.5.3" } }, "nbformat": 4, "nbformat_minor": 1 }