{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# MAC0317/MAC5920\n",
    "## Introdução ao Processamento de Sinais Digitais\n",
    "### Seção 3.8: Compressão JPEG"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.fftpack as spfft\n",
    "from imageio import imread\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.animation as anim\n",
    "from IPython.display import HTML\n",
    "from urllib.request import urlopen"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Seção 3.8.1: Visão geral da compressão JPEG\n",
    "\n",
    "Vamos considerar como entrada uma imagem RGB $M\\times N$ com 24 bits por pixel, sendo 8 bits para cada canal, em valores inteiros com sinal (ou seja, na faixa $-127,\\ldots,+128$). Se os valores originais forem sem sinal (na faixa $0,\\ldots,255$) basta subtrair $127$ de cada pixel."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Passos do método JPEG:\n",
    "1. Separação em cores (RGB $\\rightarrow R\\in M_{M\\times N}(\\mathbb{R}), G\\in M_{M\\times N}(\\mathbb{R}), B\\in M_{M\\times N}(\\mathbb{R})$ )\n",
    "\n",
    "2. Blocagem (8x8) e DCT em blocos.\n",
    "\n",
    "3. Quantização: esse passo usa uma matriz $E \\in M_{8\\times 8}$ onde $e_{k,l}$ determina a quantização de cada coeficiente DCT de frequências $(k,l)$ através da fórmula:\n",
    "\n",
    "\t$$q(x) = round\\left(\\frac{x}{e_{k,l}}\\right).$$\n",
    "    \n",
    "4. Codificação"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "A reconstrução é dada pelos passos\n",
    "\n",
    "4. Decodificação\n",
    "\n",
    "3. Dequantização:\n",
    "\n",
    "$$\\tilde{x} = e_{k,l}q(x)$$\n",
    "    \n",
    "2. DCT inversa nos blocos e remontagem da matriz\n",
    "\n",
    "1. Combinação das cores"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Seção 3.8.2: Detalhes da quantização na DCT em JPEG\n",
    "\n",
    "Um exemplo de matriz usada na expressão $q(x) = round\\left(\\frac{x}{e_{k,l}}\\right)$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "E = np.array([ [ 16,  11,  10,  16,  24,  40,  51,  61 ],\n",
    "               [ 12,  12,  14,  19,  26,  58,  60,  55 ],\n",
    "               [ 14,  13,  16,  24,  40,  57,  69,  56 ],\n",
    "               [ 14,  17,  22,  29,  51,  87,  80,  62 ],\n",
    "               [ 18,  22,  37,  56,  68, 109, 103,  77 ],\n",
    "               [ 24,  35,  55,  64,  81, 104, 113,  92 ],\n",
    "               [ 49,  64,  78,  87, 103, 121, 120, 101 ],\n",
    "               [ 72,  92,  95,  98, 112, 100, 103,  99 ] ] )\n",
    "plt.imshow(E, cmap='gray');plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Exemplo 3.1: exemplo de bloco (parte de uma imagem):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "B = np.array( \n",
    "    [ [ 47,  32,  75, 148, 192, 215, 216, 207 ],\n",
    "      [ 36,  82, 161, 196, 205, 207, 190, 140 ],\n",
    "      [ 86, 154, 200, 203, 213, 181, 143, 82 ],\n",
    "      [ 154, 202, 209, 203, 159, 145, 147, 127 ],\n",
    "      [ 184, 207, 199, 147, 134, 127, 137, 138 ],\n",
    "      [ 205, 203, 125,  72, 123, 129, 150, 115 ],\n",
    "      [ 209, 167, 126, 107, 111,  94, 105, 107 ],\n",
    "      [ 191, 129, 126, 136, 106,  54,  99, 165 ] ] )\n",
    "plt.imshow(B, cmap='gray')\n",
    "plt.title(\"Figura 3.14: Bloco B de exemplo para quantização\");plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Cálculo de $\\hat{B}=\\mbox{DCT}(B)$ (usando $B-127$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# calcula a DCT 2D fazendo a DCT 1D das linhas e das colunas\n",
    "def dct_2d(m):\n",
    "    D1 = spfft.dct(m.T, norm='ortho')\n",
    "    D2 = spfft.dct(D1.T, norm='ortho')\n",
    "    return D2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "Bhat = dct_2d(B-127)\n",
    "print(np.round(Bhat,0))\n",
    "plt.imshow(Bhat, cmap='gray');plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Bloco com coeficientes DCT quantizados $q(\\hat{B})$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "Bquant = np.round(Bhat/(E))\n",
    "print(Bquant)\n",
    "plt.imshow(Bquant, cmap='gray');plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### DCT $\\tilde{B}$ reconstruída a partir dos coeficientes quantizados"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "Btil = Bquant*(E);\n",
    "print(Btil)\n",
    "plt.imshow(Btil, cmap='gray');plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# calcula a IDCT 2D fazendo a IDCT 1D das linhas e das colunas\n",
    "def idct_2d(m):\n",
    "    D1 = spfft.idct(m.T, norm='ortho')\n",
    "    D2 = spfft.idct(D1.T, norm='ortho')\n",
    "    return D2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Bloco JPEG reconstruído a partir da DCT quantizada"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "Bn = np.round(idct_2d(Btil)+127)\n",
    "print(Bn)\n",
    "plt.imshow(Bn, cmap='gray');plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Compara bloco $B$ original e o bloco reconstruído $\\mbox{IDCT}(\\tilde{B})$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,2,figsize=(15,5))\n",
    "ax[0].set_title(\"Imagem original\")\n",
    "ax[0].imshow(B, cmap='gray')\n",
    "ax[1].set_title(\"Imagem reconstruída: Dist = {0:.2f}%, Comp = {1:.2f}%\".format\n",
    "                (100*np.linalg.norm(B-Bn)**2/np.linalg.norm(B)**2,100*np.sum(np.sum(Bquant!=0))/64))\n",
    "ax[1].imshow(Bn, cmap='gray');plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Flexibilizando a quantização\n",
    "\n",
    "É possível aumentar a flexibilidade do esquema de quantização usando um parâmetro adicional $r>0$ junto com a matriz de quantização $E$, adaptando a fórmula da quantização e dequantização para\n",
    "$$q(x) = round\\left(\\frac{x}{r\\cdot e_{k,l}}\\right)$$\n",
    "e\n",
    "$$\\tilde{x} = r\\cdot e_{k,l}q(x).$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "A interpretação dessa flexibilização é:\n",
    "- para $r>1$ teremos uma quantização mais \"grosseira\", reduzindo o número de bits necessários para armazenar cada valor (e piorando a qualidade da imagem reconstruída);\n",
    "- para $r\\in(0,1)$ teremos uma quantização mais \"fina\", preservando mais detalhes da imagem reconstruída (às custas de uma compressão menos econômica)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Repete o exemplo anterior com r=0.1, 0.5, 1 e 2\n",
    "R = [ 0.1, 0.5, 1, 2 ]\n",
    "fig, ax = plt.subplots(1,len(R),figsize=(15,5))\n",
    "for k in range(len(R)):\n",
    "    Bquant = np.round(Bhat/(R[k]*E)) # Calcula valores quantizados\n",
    "    Btil = Bquant*(R[k]*E) # Retorna à escala original (este é o passo da de(s)quantização)\n",
    "    Bn = np.round(idct_2d(Btil)+127) # Aplica a DCT, voltando os valores à faixa 0...255\n",
    "    ax[k].set_title(f\"r={R[k]}\")\n",
    "    ax[k].set_xlabel(f\"Dist={100*np.linalg.norm(B-Bn)**2/np.linalg.norm(B)**2:.2f}%\\nComp={100*np.sum(np.sum(Bquant!=0))/64:.2f}%\")\n",
    "    ax[k].imshow(Bn/np.amax(Bn), cmap='gray', interpolation='None')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Funções auxiliares: DCT e IDCT em blocos m x n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Adaptado de https://www.rose-hulman.edu/mathDFT/matlabbookroutines/blkdct2.m\n",
    "# Uso: Y = blkdct2(y, m, n) e y = blkidct2(Y, m, n).\n",
    "# DCT e IDCT em blocos sobre a matriz y usando blocos de m x n. As dimensões\n",
    "# das matrizes devem ser divisíveis por m e n respectivamente.\n",
    "\n",
    "def blkdct2(y, m, n):\n",
    "    my,ny = y.shape\n",
    "    Y = 0*y\n",
    "    for j in range(0, my, m):\n",
    "        for k in range(0, ny, n):\n",
    "            Z = y[j : j+m, k : k+n]\n",
    "            Y[j : j+m, k : k+n] = dct_2d(Z)\n",
    "    return Y\n",
    "\n",
    "def blkidct2(y, m, n):\n",
    "    my,ny = y.shape\n",
    "    Y = 0*y\n",
    "    for j in range(0, my, m):\n",
    "        for k in range(0, ny, n):\n",
    "            Z = y[j : j+m, k : k+n]\n",
    "            Y[j : j+m, k : k+n] = idct_2d(Z)\n",
    "    return Y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Seção 3.8.3: JPEG aplicado à roda gigante"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# função auxiliar: converte imagem colorida para cinza\n",
    "def rgb2gray(rgb):\n",
    "    fil = [0.299, 0.587, 0.144]\n",
    "    return np.dot(rgb, fil)\n",
    "# Carrega a imagem da internet e converte para nível de cinza\n",
    "url = \"http://sutherncharm.files.wordpress.com/2009/09/double-ferris.jpg\"\n",
    "M = rgb2gray(imread(urlopen(url).read()))\n",
    "My, Mx = M.shape\n",
    "divsize = 8\n",
    "new_My = My - My % divsize\n",
    "new_Mx = Mx - Mx % divsize\n",
    "M=M[:new_My, :new_Mx]\n",
    "My,Mx = M.shape\n",
    "# Calcula a DCT em blocos de 8x8, e mostra o resultado\n",
    "N=blkdct2(M,8,8)\n",
    "Nlog = np.log(1 + abs(N))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,2,figsize=(15,7))\n",
    "ax[0].imshow(M, cmap='gray');ax[0].axis('off');ax[0].set_title(\"Imagem original\")\n",
    "ax[1].imshow(Nlog / np.amax(Nlog), cmap='gray');ax[1].axis('off');ax[1].set_title(\"DCT em blocos\");plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "Nquant = N.copy()\n",
    "Ntil = N.copy()\n",
    "# Parâmetro de granularidade da quantização\n",
    "# (faz pouca diferença nesse exemplo; experimente R=0.1)\n",
    "R=1\n",
    "# Processa todos os blocos\n",
    "for j in range(0, My, divsize):\n",
    "    for k in range(0, Mx, divsize):\n",
    "        Nquant[j:j+divsize,k:k+divsize] = np.round(N[j:j+divsize,k:k+divsize]/(R*E))\n",
    "        Ntil[j:j+divsize,k:k+divsize] = (R*E)*Nquant[j:j+divsize,k:k+divsize]\n",
    "Mn = np.round(blkidct2(Ntil,8,8))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,2,figsize=(15,7))\n",
    "ax[0].imshow(M, cmap='gray');ax[0].axis('off');ax[0].set_title(\"Imagem original\")\n",
    "ax[1].imshow(Mn, cmap='gray');ax[1].axis('off');\n",
    "ax[1].set_title(\"JPEG: Dist = {0:.2f}% e Comp = {1:.2f}%\".format(100*np.linalg.norm(M-Mn)**2/np.linalg.norm(M)**2,100*np.sum(Nquant!=0)/(My*Mx)));\n",
    "fig.suptitle(\"Aplicação do esquema JPEG em blocos 8x8 à imagem da roda gigante\");plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,2,figsize=(15,8))\n",
    "ax[0].imshow(abs(M-Mn), cmap='gray', vmin=0, vmax=255)\n",
    "ax[0].set_title(\"Imagem diferença\");ax[0].axis('off')\n",
    "ax[1].imshow((abs(M-Mn)/np.amax(M-Mn)), cmap='gray')\n",
    "ax[1].set_title(\"Imagem diferença normalizada\");ax[1].axis('off');\n",
    "fig.suptitle(\"Comparação entre as imagens original e JPEG através das diferenças entre as imagens\");plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Codificação sequencial $\\times$ progressiva\n",
    "\n",
    "- **Sequencial:** varredura tradicional (por linhas, por exemplo) da matriz $\\hat{B}$ com blocos $\\hat{B}_{ij}\\in M_{8\\times 8}(\\mathbb{R})$. Pode-se reconstruir a imagem a partir dos blocos $8\\times 8$ do canto superior esquerdo, na varredura por linhas, sendo que a qualidade da imagem parcial reconstruída é a melhor possível.\n",
    "\n",
    "- **Progressiva:** baseada no agrupamento das DCT's dos blocos por pares de frequência ($k,l$) iguais. Transmite-se inicialmente os coeficientes D.C. ($\\hat{C}_{0,0}$) e progressivamente adicionam-se frequências mais altas. A informação chega na forma de uma imagem inteira e grosseira, e progressivamente aumenta-se a definição."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Reconstrução progressiva usando coeficientes DCT\n",
    "\n",
    "Considere uma imagem $A$ de tamanho $m\\times n$, e seja $\\hat{A}$ sua DCT em blocos agrupados por frequência, onde $\\tilde{A}^{(k,l)}$, é um bloco de tamanho $\\frac{m}{8}\\times\\frac{n}{8}$, para $k,l=0,1,\\ldots,7$.\n",
    "\n",
    "Seja agora a matriz $\\widehat{A^{(k,l)}}$ a matriz $m\\times n$ obtida por manter apenas o bloco $\\tilde{A}^{(k,l)}$ e zerar o resto da matriz."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Definiremos a reconstrução de nível $s=0,1,\\ldots,14$ como a matriz $A_s$ que é a transformada inversa da matriz\n",
    "\n",
    "$$\\widehat{A_s} = \\sum_{k+l=s} \\widehat{A^{(k,l)}} = \\sum_{k=0}^s\\sum_{l=0}^{s-k}\\widehat{A^{(k,l)}}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Assim,\n",
    "\n",
    "$$\\begin{array}{l}\n",
    "A_0 = \\mbox{IDCT_BLK}(\\widehat{A_0}) = \\mbox{IDCT_BLK}(\\widehat{A^{(0,0)}})\\\\\n",
    "A_1 = \\mbox{IDCT_BLK}(\\widehat{A_1}) = \\mbox{IDCT_BLK}(\\widehat{A^{(0,1)}}+\\widehat{A^{(1,0)}})\\\\\n",
    "A_2 = \\mbox{IDCT_BLK}(\\widehat{A_2}) = \\mbox{IDCT_BLK}(\\widehat{A^{(0,2)}}+\\widehat{A^{(1,1)}}+\\widehat{A^{(2,0)}})\\\\\n",
    "\\vdots\n",
    "\\end{array}$$\n",
    "\n",
    "correspondem a níveis de detalhes associados a combinações de frequências $k+l=s$ progressivamente mais altas."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Note que\n",
    "$$\\hat{A} = \\sum_{s=0}^{14}\\widehat{A_s},$$\n",
    "\n",
    "assim, pela linearidade da DCT,\n",
    "\n",
    "$$A = \\sum_{s=0}^{14} A_s.$$\n",
    "\n",
    "A *reconstrução até ordem* $l=0,1,\\ldots,14$ é a imagem\n",
    "\n",
    "$$\\sum_{s=0}^{l} A_s.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "%%capture\n",
    "\n",
    "plt.ion()\n",
    "fig = plt.figure(figsize=(15,8))\n",
    "ax = fig.add_subplot(111)\n",
    "plt.gray()\n",
    "plt.axis(\"off\")\n",
    "ims = []\n",
    "Nl=0*Ntil\n",
    "for l in range(15):\n",
    "    for j in range(l+1):\n",
    "        k = l-j\n",
    "        Nl[j:My:divsize,k:Mx:divsize] = Ntil[j:My:divsize,k:Mx:divsize]\n",
    "    Ml = np.round(blkidct2(Nl,divsize,divsize))\n",
    "    im = plt.imshow(Ml, cmap='gray')\n",
    "    title = f\"Reconstrução progressiva até j+k={l}, distorção introduzida={100*np.linalg.norm(Mn-Ml)**2/np.linalg.norm(M)**2:.2f}%\"\n",
    "    ttl = plt.text(0.5, -0.1, title, horizontalalignment='center', verticalalignment='bottom', fontsize='large', transform=ax.transAxes)\n",
    "    ims.append([im, ttl])\n",
    "    \n",
    "imf = anim.ArtistAnimation(fig, ims, interval=1000, blit=True, repeat=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Exemplo 3.2 e figura 3.17"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "HTML(imf.to_jshtml())"
   ]
  }
 ],
 "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
}