{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# MAC0317/MAC5920\n",
    "## Introdução ao Processamento de Sinais Digitais\n",
    "### Seção 2.7: A DFT bidimensional"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "import math as m\n",
    "from urllib.request import urlopen\n",
    "import numpy as np\n",
    "from imageio import imread\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import gridspec"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Definição da DFT bidimensional\n",
    "\n",
    "A transformada de Fourier de uma imagem $A\\in\\mathcal{M}_{m,n}(\\mathbb{C})$ segue bem de perto o exemplo 1.24, e corresponde aos coeficientes não-normalizados da mudança de base da representação usual (base canônica) para a base das formas básicas de onda bidimensionais dadas pelas matrizes $\\mathcal{E}_{k,l}$, definidas na seção 1.7.2 como\n",
    "\n",
    "$$\\left(\\mathcal{E}_{k,l}\\right)_{r,s} = e^{i2\\pi(kr/m+ls/n)},\\quad r=0,\\ldots,m-1,\\quad s=0,\\ldots,n-1,$$\n",
    "\n",
    "onde escrevemos $A=\\displaystyle\\sum_{k=0}^{m-1}\\sum_{l=0}^{n-1}c_{k,l}\\mathcal{E}_{k,l}$ e obtemos os coeficientes $c_{k,l}$ pelo teorema 1.8.3 como $c_{k,l}=\\frac{\\left(A,\\mathcal{E}_{k,l}\\right)}{mn}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Definição 2.7.1 (DFT bidimensional)\n",
    "\n",
    "Dada uma matriz $A\\in\\mathcal{M}_{m,n}(\\mathbb{C})$, definimos sua DFT $\\hat{A}=DFT(A)\\in\\mathcal{M}_{m,n}(\\mathbb{C})$ como a matriz cujas componentes são\n",
    "\n",
    "$$\\hat{A}_{k,l} = \\left(A,\\mathcal{E}_{k,l}\\right) = \\sum_{r=0}^{m-1}\\sum_{s=0}^{n-1}A_{r,s}e^{-i2\\pi(kr/m+ls/n)},\\quad \\genfrac{}{}{0pt}{}{k=0,\\ldots,m-1,}{l=0,\\ldots,n-1}.$$\n",
    "\n",
    "**Observe** que nessa definição $k$ e $l$ são índices associados às *frequências de varredura por linhas e por colunas*, respectivamente, enquanto $r$ e $s$ são os índices de linhas e colunas dentro das matrizes $A$ e $\\mathcal{E}_{k,l}$.\n",
    "\n",
    "**Observe** também que a única diferença entre os coeficientes da DFT $\\hat{A}_{k,l}$ e os coeficientes $c_{k,l}$ da combinação linear que expressa $A$ em função das formas básicas de onda é a ausência da constante $\\frac{1}{mn}$ nos coeficientes da DFT. Isso é análogo à definição usada no caso unidimensional, sendo que necessariamente essa constante aparecerá na equação da IDFT."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Observação 2.5\n",
    "\n",
    "Note que a expressão\n",
    "$$\\hat{A}_{k,l} = \\left(A,\\mathcal{E}_{k,l}\\right) = \\sum_{r=0}^{m-1}\\sum_{s=0}^{n-1}A_{r,s}e^{-i2\\pi(kr/m+ls/n)}$$\n",
    "é periódica em relação a $k$ e $l$, com períodos $m$ e $n$ respectivamente, ou seja\n",
    "$$\\hat{A}_{k+pm,l+qn} = \\hat{A}_{k,l},\\ \\forall p,q\\in\\mathbb{Z}.$$\n",
    "\n",
    "Em função disso, é possível visualizar a DFT bidimensional em outras faixas de frequência além de $k=0,\\ldots,m-1,\\quad l=0,\\ldots,n-1$. Uma escolha comum, análoga ao caso unidimensional, consiste em utilizar a faixa\n",
    "$$-\\frac{m}{2}<k\\le\\frac{m}{2},\\quad\\quad-\\frac{n}{2}<l\\le\\frac{n}{2}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Definição 2.7.2 (IDFT bidimensional)\n",
    "\n",
    "Dada uma matriz $\\hat{A}\\in\\mathcal{M}_{m,n}(\\mathbb{C})$ definimos sua IDFT $A=IDFT(\\hat{A})\\in\\mathcal{M}_{m,n}(\\mathbb{C})$ pela expressão\n",
    "\n",
    "$$A=\\frac{1}{mn}\\sum_{k=0}^{m-1}\\sum_{l=0}^{n-1}\\hat{A}_{k,l}\\mathcal{E}_{k,l},$$\n",
    "\n",
    "sendo que as componentes de $A$ podem ser expressas em função de $\\hat{A}$ como\n",
    "\n",
    "$$A_{r,s} = \\frac{1}{mn}\\sum_{k=0}^{m-1}\\sum_{l=0}^{n-1}\\hat{A}_{k,l}e^{i2\\pi(kr/m+ls/n)},\\quad r=0,\\ldots,m-1,\\quad s=0,\\ldots,n-1.$$\n",
    "\n",
    "**Observe** que mais uma vez nessa definição $k$ e $l$ são índices associados às *frequências de varredura por linhas e por colunas*, respectivamente, enquanto $r$ e $s$ são os índices de linhas e colunas dentro das matrizes $A$ e $\\mathcal{E}_{k,l}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Formulação matricial da DFT 2D\n",
    "\n",
    "Uma relação muito importante entre a DFT unidimensional e a DFT bidimensional é a seguinte: \n",
    "\n",
    "### Proposição 2.7.1\n",
    "se $\\hat{A} = DFT(A)$, onde $A,\\hat{A}\\in\\mathcal{M}_{m,n}(\\mathbb{C})$, então\n",
    "\n",
    "$$\\hat{A} = F_mAF_n^T,$$\n",
    "\n",
    "onde $F_m\\in\\mathcal{M}_{m,m}(\\mathbb{C})$ e $F_n\\in\\mathcal{M}_{n,n}(\\mathbb{C})$ são as matrizes que representam as DFTs unidimensionais nos espaços $\\mathbb{C}^m$ e $\\mathbb{C}^n$, respectivamente."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Para ver que isso é verdade, basta calcular o produto de matrizes acima. Vamos começar calculando a matriz $F_mA$; usando a definição do produto de matrizes:\n",
    "\n",
    "$$\\begin{array}{ll}\n",
    "\\left(F_mA\\right)_{k,s} & = (\\mbox{\"linha k de $F_m$\" $\\quad,\\quad$ \"coluna s de $A$\")}\\\\\n",
    "& = \\displaystyle\\sum_{r=0}^{m-1}A_{r,s}e^{-i2\\pi kr/m}.\n",
    "\\end{array}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Usando agora essa expressão no produto $F_mAF_n^T$, temos:\n",
    "$$\\begin{array}{ll}\n",
    "\\left(F_mAF_n^T\\right)_{k,l} = ((F_mA)F_n^T)_{k,l}\n",
    "& = \\mbox{\"linha k de $F_mA$\" $\\quad\\times\\quad$ \"coluna l de $F_n^T$\"}\\\\\n",
    "& = \\mbox{\"linha k de $F_mA$\" $\\quad\\times\\quad$ \"linha l de $F_n$\"}\\\\\n",
    "& = \\displaystyle\\sum_{s=0}^{n-1}(FmA)_{k,s}e^{-i2\\pi ls/n}\\\\\n",
    "& = \\displaystyle\\sum_{s=0}^{n-1}\\left(\\sum_{r=0}^{m-1}A_{r,s}e^{-i2\\pi kr/m}\\right)e^{-i2\\pi ls/n}\\\\\n",
    "& = \\displaystyle\\sum_{r=0}^{m-1}\\sum_{s=0}^{n-1}A_{r,s}e^{-i2\\pi(kr/m+ls/n)}\\\\\n",
    "& = \\hat{A}_{k,l}\n",
    "\\end{array}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "A interpretação algorítmica da identidade $DFT(A) = F_mAF_n^T$ é a seguinte:\n",
    "\n",
    "- A parte $F_mA$ pode ser vista como uma aplicação da DFT nas **colunas** de A:\n",
    "\n",
    "$$F_mA = \\left[\\begin{array}{l|l|l|l}\n",
    "\\rule{0mm}{0mm}\\\\\n",
    "F_mA^0&F_mA^1&\\ldots&F_mA^{n-1}\\\\\n",
    "\\rule{0mm}{0mm}\n",
    "\\end{array}\\right]\n",
    " = \\left[\\begin{array}{l|l|l|l}\n",
    "\\rule{0mm}{0mm}\\\\\n",
    "DFT(A^0)&DFT(A^1)&\\ldots&DFT(A^{n-1})\\\\\n",
    "\\rule{0mm}{0mm}\n",
    "\\end{array}\\right]$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "- O produto subsequente por $F_n^T$ pode ser visto como uma aplicação da DFT nas **linhas** da matriz resultante $X=F_mA$, pois\n",
    "\n",
    "$$XF_n^T = \\left[\\begin{array}{c}\n",
    "\\rule{2cm}{0mm}X_0F_n^T\\rule{2cm}{0mm}\\\\\\hline\n",
    "X_1F_n^T\\\\\\hline\n",
    "\\vdots\\\\\\hline\n",
    "X_{m-1}F_n^T\n",
    "\\end{array}\\right]\n",
    " = \\left[\\begin{array}{c}\n",
    "\\rule{2cm}{0mm}DFT(X_0^T)^T\\rule{2cm}{0mm}\\\\\\hline\n",
    "DFT(X_1^T)^T\\\\\\hline\n",
    "\\vdots\\\\\\hline\n",
    "DFT(X_{m-1}^T)^T\n",
    "\\end{array}\\right]$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "- Finalmente, a associatividade do produto de matrizes $$F_mAF_n^T=(F_mA)F_n^T=F_m(AF_n^T)$$ permite ver que o mesmo resultado poderia ser atingido calculando-se primeiro a DFT das linhas de $A$ e depois a DFT das colunas resultantes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Seção 2.7.1: Interpretação e exemplos da DFT 2D\n",
    "\n",
    "- Abaixo veremos uma imagem \"natural\" e uma imagem sintética, e observaremos como os coeficientes de Fourier se distribuem no âmbito de frequências\n",
    "$$k=0,\\ldots,m-1,\\quad l=0,\\ldots,n-1$$\n",
    "e também\n",
    "$$-\\frac{m}{2}<k\\le\\frac{m}{2},\\quad\\quad-\\frac{n}{2}<l\\le\\frac{n}{2}.$$\n",
    "\n",
    "- Como a magnitude dos coeficientes possui grande variação, visualizaremos também as imagens correspondentes ao mapeamento $\\log\\left(1+|\\hat{A}|\\right)$.\n",
    "\n",
    "- Exploraremos também a ressíntese da imagem aplicando um esquema simples de compressão similar ao que usamos no caso unidimensional (som do sino): qualquer coeficiente com amplitude inferior a $\\frac{\\max|\\hat{A}_{k,l}|}{1000}$ será eliminado. O objetivo é comparar as taxas de compressão atingidas e o resultado perceptual desse esquema simples de compressão."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "url = \"http://sutherncharm.files.wordpress.com/2009/09/double-ferris.jpg\"\n",
    "M = imread(urlopen(url).read())\n",
    "# FFT da imagem\n",
    "N = np.fft.fft2(M,axes=(0,1)); absN = np.abs(N)\n",
    "# Centralizando fft naorigem\n",
    "halfheight, halfwidth = int(absN.shape[0]/2), int(absN.shape[1]/2)\n",
    "shiftN = np.roll(absN, (halfheight,halfwidth), axis=(0,1))\n",
    "# Normalizando pelom ln\n",
    "Nlog = np.log(1 + shiftN)\n",
    "maxi, mini = Nlog.max(), Nlog.min(); Nshow = (Nlog - mini) / (maxi - mini)\n",
    "# comprimindo a imagem\n",
    "C = 0.0001 * absN.max()\n",
    "R = 0\n",
    "for val in np.nditer(N, op_flags=['readwrite']):\n",
    "    if abs(val) < C:\n",
    "        val[...] = 0\n",
    "        R += 1\n",
    "# calcula compressão em relação a cada componente R, G ou B anulada\n",
    "Compression = R/(3*N.shape[0]*N.shape[1])\n",
    "reM = np.round(np.real(np.fft.ifft2(N, axes=(0,1))))\n",
    "reMin, reMax = reM.min(), reM.max()\n",
    "reMshow = (reM - reM.min())/ (reM.max() - reM.min())\n",
    "Distortion = np.linalg.norm(reM-M)**2/np.linalg.norm(M)**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Comandos para plotar imagens\n",
    "f, ax = plt.subplots(1,5,figsize=(15,4.5))\n",
    "f.suptitle(\"Figura 2.13 - Imagem e várias apresentações da DFT\")\n",
    "ax[0].set_title(\"Imagem original\"); ax[0].imshow(M)\n",
    "ax[1].set_title(\"Espectro\"); ax[1].imshow(absN / absN.max())\n",
    "ax[2].set_title(\"Espectro centralizado em (0,0)\"); ax[2].imshow(shiftN / shiftN.max(), extent=[ -halfwidth, halfwidth,-halfheight, halfheight])\n",
    "ax[3].set_title(\"Log-Espectro normalizado\"); ax[3].imshow(Nshow, extent=[ -halfwidth, halfwidth,-halfheight, halfheight])\n",
    "ax[4].set_title(\"Ressíntese\"); ax[4].imshow(reMshow)\n",
    "plt.show()\n",
    "print(\"Compressão na ressíntese = {0:.3f}%\".format(100*Compression))\n",
    "print(\"Distorção na ressíntese = {0:.3f}%\".format(100*Distortion))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "O = np.zeros((350,450,3)); O[50:90,40:60,0:3] += 1\n",
    "#FFT da imagem\n",
    "N = np.fft.fft2(O,axes=(0,1)); absN = abs(N)\n",
    "# Centralizando o espectro na origem\n",
    "halfheight, halfwidth = int(absN.shape[0]/2), int(absN.shape[1]/2)\n",
    "shiftN = np.roll(absN, (halfheight, halfwidth), axis=(0,1))\n",
    "# Normalização com ln\n",
    "Nlog = np.log(1 + shiftN)\n",
    "maxi, mini = Nlog.max(), Nlog.min()\n",
    "Nshow = (Nlog - mini) / (maxi - mini)\n",
    "# comprimindo a imagem\n",
    "C = 0.0001 * absN.max()\n",
    "R = 0\n",
    "for val in np.nditer(N, op_flags=['readwrite']):\n",
    "    if abs(val) < C:\n",
    "        val[...] = 0\n",
    "        R += 1\n",
    "# calcula compressão\n",
    "Compression = R/(N.shape[0]*N.shape[1])\n",
    "reM = np.round(np.real(np.fft.ifft2(N, axes=(0,1))))\n",
    "reMin, reMax = reM.min(), reM.max()\n",
    "reMshow = (reM - reM.min())/ (reM.max() - reM.min())\n",
    "Distortion = np.linalg.norm(reM-O)**2/np.linalg.norm(O)**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# comandos para plotar as imagens\n",
    "f, ax = plt.subplots(1,5,figsize=(15,5))\n",
    "f.suptitle(\"Figura 2.14: Imagem artificial e visualizações da DFT\")\n",
    "ax[0].set_title(\"Imagem original\"); ax[0].imshow(O)\n",
    "ax[1].set_title(\"Espectro\"); ax[1].imshow(absN/ absN.max())\n",
    "ax[2].set_title(\"Espectro centralizado em (0,0)\"); ax[2].imshow(shiftN / shiftN.max(), extent=[-halfwidth, halfwidth,-halfheight,halfheight])\n",
    "ax[3].set_title(\"Log-Espectro normalizado\"); ax[3].imshow(Nshow)\n",
    "ax[4].set_title(\"Ressíntese\"); ax[4].imshow(reMshow)\n",
    "plt.show()\n",
    "print(\"Compressão na ressíntese = {0:.3f}%\".format(100*Compression))\n",
    "print(\"Distorção na ressíntese = {0:.3f}%\".format(100*Distortion))"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}