TP4 - Descripteurs d'image.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Séance 4 : descripteurs d'image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Cette dernière séance s'intéresse aux possibilités de reconnaître des objets dans une image ou de trouver des images similaires à une image requête. Comme pour l'alignement d'images à partir de points d'intérêt, nous allons pour ce faire utiliser des **descripteurs**. Ces descripteurs ne porteront toutefois pas sur un point d'intérêt précis mais sur l'image entière."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Template matching"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Le *template matching* consiste à localiser un objet, décrit par un modèle (*template*) dans une image. En pratique, on aura donc l'image d'un modèle $T$ de dimensions $w \\times h$ et une image $I$ de dimensions $W \\times H$ dans laquelle effectuer la recherche.\n",
    "\n",
    "La méthode classique consiste à chercher la position $(x^*,y^*)$ idéale dans l'image qui minimise la distance entre la fenêtre de taille $w \\times h$ centrée en $(x,y)$ et le *template*:\n",
    "\n",
    "$$x^*, y^* = \\arg\\min_{x,y} E(x,y) = \\arg\\min \\sum_{i=1}^w \\sum_{j=1}^h \\left( I[x+i,y+j] - T[i, j] \\right)^2$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Cette recherche est équivalente à chercher la position $(x^*,y^*)$ qui maximise la **corrélation croisée** entre la fenêtre considérée et le *template*:\n",
    "\n",
    "$$x^*, y^* = \\arg\\max_{x,y} \\sum_{i=1}^w \\sum_{j=1}^h I[x+i,y+j] \\cdot T[i,j] = \\arg\\max_{x,y} F^{(x,y)} \\star T$$\n",
    "\n",
    "où $F^{(x,y)}$ est la fenêtre centrée en $(x,y)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Pour trouver la position la plus probable du *template* dans l'image, on peut donc procéder à l'aide de l'algorithme suivant :\n",
    "1. calculer la corrélation croisée pour toutes les positions $(x,y)$\n",
    "2. renvoyer les valeurs $(x^*, y^*)$ pour lesquelles la corrélation est maximale"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## À la main"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Mettons cette logique en application sur un exemple. Pour simplifier les choses, nous allons travailler en niveaux de gris mais on aurait pu tout aussi bien définir la corrélation croisée en sommant la corrélation pour les canaux rouge, bleu et vert."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Les images avec lesquelles nous allons faire nos tests sont des captures d'écran de Pokémon Rouge, sorti en France en 1999 sur Game Boy ."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "import skimage\n",
    "import matplotlib.pyplot as plt\n",
    "from skimage import io\n",
    "\n",
    "template = io.imread(\"images/pokeball.png\", as_gray=True)\n",
    "image = io.imread(\"images/pokemon_safari.jpg\", as_gray=True)\n",
    "\n",
    "fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(12, 6), gridspec_kw={'width_ratios': [2.5, 1]})\n",
    "axes[0].imshow(image, cmap=\"gray\") and axes[0].set_title(\"Image où chercher l'objet\")\n",
    "axes[1].imshow(template, cmap=\"gray\") and axes[1].set_title(\"Template\") and plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Remarquons que les dimensions du *template* correspondent à la dimension de l'objet dans l'image."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Nous pourrions calculer la corrélation manuellement en utilisant la formule précédemment détaillée mais ce n'est pas particulièrement intéressant. Utilisons plutôt SciPy pour le faire grâce à la fonction `scipy.signal.correlate2d`. Celle-ci calcule la corrélation croisée entre deux matrices 2D (par exemple, des images en niveaux de gris)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from scipy.signal import correlate2d\n",
    "\n",
    "correlation_map = correlate2d(image, template, mode=\"same\")\n",
    "# Affichons la carte de corrélation\n",
    "fig = plt.figure(figsize=(10, 8))\n",
    "plt.imshow(correlation_map, cmap=\"magma\")\n",
    "plt.colorbar(shrink=0.5)\n",
    "plt.title(\"Carte de corrélation croisée\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Oups. Le résultat obtenu n'est pas tout à fait ce que l'on escomptait. En effet, dans la corrélation croisée, les pixels brillants sont très favorisés. Toutes les zones de terrain blanc apparaissent donc avec une forte corrélation, plus forte même que la position réelle de l'objet.\n",
    "\n",
    "L'astuce consiste donc à ne pas considérer le template tel quel mais à lui retirer sa valeur moyenne, de sorte à débiaiser les résultats :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "template -= template.mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "correlation_map = correlate2d(image, template, mode=\"same\")\n",
    "# Affichons la carte de corrélation\n",
    "fig = plt.figure(figsize=(10, 8))\n",
    "plt.imshow(correlation_map, cmap=\"magma\")\n",
    "plt.colorbar(shrink=0.5)\n",
    "plt.title(\"Carte de corrélation croisée\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "C'est nettement mieux ! On peut maintenant utiliser un peu de magie avec NumPy pour extraire les coordonnées du maximum (le pixel le plus brillant dans la carte de corrélation croisée) :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "# y, x sont inversés entre NumPy et Matplotlib\n",
    "y_max, x_max = np.argwhere(correlation_map == correlation_map.max())[0]\n",
    "# On affiche la carte avec le point correspondant au maximum\n",
    "fig, axes = plt.subplots(ncols=2, nrows=1, figsize=(16, 8))\n",
    "axes[0].imshow(correlation_map, cmap=\"magma\")\n",
    "axes[0].plot(x_max, y_max, \"w*\", markersize=10)\n",
    "axes[0].set_title(\"Carte de corrélation croisée\")\n",
    "\n",
    "axes[1].imshow(image, cmap=\"gray\")\n",
    "w, h = template.shape\n",
    "rect = plt.Rectangle((x_max - w//2, y_max - h//2), w, h, edgecolor=\"r\", facecolor=\"none\")\n",
    "axes[1].add_patch(rect)\n",
    "axes[1].set_title(\"Détection dans l'image\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "À noter que cette technique est plutôt robuste au bruit et aux imperfections. On peut aussi bien l'utiliser sur une image dégradée qu'avec un template approximatif. Appliquons par exemple un bruit gaussien :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from skimage.util import random_noise\n",
    "noisy = random_noise(image, mode=\"gaussian\", var=0.5)\n",
    "correlation_map = correlate2d(noisy, template, mode=\"same\")\n",
    "y_max, x_max = np.argwhere(correlation_map == correlation_map.max())[0]\n",
    "\n",
    "fig, axes = plt.subplots(ncols=2, nrows=1, figsize=(16, 8))\n",
    "axes[0].imshow(correlation_map, cmap=\"magma\")\n",
    "axes[0].plot(x_max, y_max, \"w*\", markersize=10)\n",
    "axes[0].set_title(\"Carte de corrélation croisée\")\n",
    "axes[1].imshow(noisy, cmap=\"gray\")\n",
    "w, h = template.shape\n",
    "rect = plt.Rectangle((x_max - w//2, y_max - h//2), w, h, edgecolor=\"r\", facecolor=\"none\")\n",
    "axes[1].add_patch(rect)\n",
    "axes[1].set_title(\"Détection dans l'image\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Détections multiples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "L'algorithme précédent fonctionne bien mais le fait de ne considérer que la valeur maximale fait que l'on peut potentiellement rater la présence de plusieurs objets correspondants au *template*. Prenons l'image suivante :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "image2 = io.imread(\"images/pokemon_choice.png\", as_gray=True)\n",
    "plt.imshow(image2, cmap=\"gray\") and plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Que faire dans cette situation ? L'objet correspondant au *template* s'y trouve à trois reprises. On aimerait donc pouvoir détecter ces trois occurrences mais rien ne dit que la corrélation sera strictement égale au maximum pour les trois en simultanée. Une solution est donc de ne pas considérer la corrélation maximale mais de considérer les corrélations « presque maximales », par exemple entre 95% et 100% de la valeur maximale."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "On remarque par ailleurs que les objets sont plus petits dans cette image, il faut donc redimensionner le *template* (on verra plus tard comment s'en tirer quand on ne connaît pas la taille des objets que l'on cherche) :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from skimage.transform import rescale\n",
    "small_template = rescale(template, scale=0.5)\n",
    "w, h = small_template.shape\n",
    "\n",
    "correlation_map = correlate2d(image2, small_template, mode=\"same\")\n",
    "tolerance = 0.95 # on cherche les maxima au-delà d'un seuil de 95%\n",
    "\n",
    "fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(16, 8))\n",
    "axes[0].imshow(correlation_map) and axes[0].set_title(\"Corrélations croisées\")\n",
    "axes[1].imshow(image2, cmap=\"gray\") and axes[1].set_title(\"Détection dans l'image\")\n",
    "for y_max, x_max in np.argwhere(correlation_map > tolerance * correlation_map.max()):\n",
    "    axes[0].plot(x_max, y_max, \"w*\", markersize=12)\n",
    "    rect = plt.Rectangle((x_max - w//2, y_max - h//2), w, h, edgecolor=\"r\", facecolor=\"none\")\n",
    "    axes[1].add_patch(rect)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Et le tour est joué ! Essayez d'abaisser ou d'augmenter le seuil de tolérance pour voir l'impact que cela peut avoir sur les détections dans l'image."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Détections à plusieurs échelles"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Comme nous venons de le voir, une limite du *template matching* est que l'on a besoin de connaître au préalable la taille de l'objet dans l'image pour savoir à quelle échelle redimensionner notre *template*. C'est assez gênant car en général, on ne peut pas vraiment faire d'hypothèse : si l'objet est plus loin de l'appareil photo, il apparaîtra de la même façon mais plus petit."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Une solution de contournement consiste à tester le *template* à plusieurs échelles. Pour ce faire, une possibilité est de redimensionner le *template* pour réduire sa taille :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# les échelles de \"zoom\" à tester : 25%, 50%, 75% et 100%\n",
    "scales = [0.25, 0.5, 0.75, 1.0]\n",
    "tolerance = 0.95\n",
    "\n",
    "fig, ax = plt.subplots(ncols=len(scales), nrows=1, figsize=(16, 6))\n",
    "for idx, scale in enumerate(scales):\n",
    "    scaled_template = rescale(template, scale)\n",
    "    correlation_map = correlate2d(image2, scaled_template, mode=\"same\")\n",
    "    ax[idx].imshow(correlation_map)\n",
    "    for y_max, x_max in np.argwhere(correlation_map > tolerance * correlation_map.max()):\n",
    "        ax[idx].plot(x_max, y_max, \"w*\", markersize=12)\n",
    "    ax[idx].set_title(f\"Échelle = {scale},\\n corrélation max. = {correlation_map.max():.1f}\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "*Note* : Pour agrandir le *template*, on va plutôt faire l'opposé et réduire la taille de l'image (la corrélation croisée est plus rapide dans ce sens)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Maintenant que l'on a effectué le *template matching* à plusieurs échelles, on se retrouve avec plein de détections sur les bras. Mais comment peut-on décider quelle échelle est la plus pertinente pour nous ?\n",
    "\n",
    "**Question** : comment évoluent les corrélations maximales quand on augmente les dimensions du *template* ? Pourquoi ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Une solution à ce problème consiste à *normaliser* la corrélation croisée en divisant par la taille du *template*:\n",
    "    \n",
    "$$x^*, y^* = \\arg\\max_{x,y} \\frac{1}{wh} \\sum_{i=1}^w \\sum_{j=1}^h I[x+i,y+j] \\cdot T[i,j] = \\arg\\max_{x,y} \\frac{1}{wh} F^{(x,y)} \\star T$$\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "En Python, cela donne :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# les échelles de \"zoom\" à tester : 25%, 50%, 75% et 100%\n",
    "scales = [0.25, 0.5, 0.75, 1.0]\n",
    "tolerance = 0.95\n",
    "\n",
    "fig, ax = plt.subplots(ncols=len(scales), nrows=1, figsize=(16, 6))\n",
    "for idx, scale in enumerate(scales):\n",
    "    scaled_template = rescale(template, scale)\n",
    "    correlation_map = correlate2d(image2, scaled_template, mode=\"same\")\n",
    "    # Normalisation ->\n",
    "    correlation_map /= scaled_template.size\n",
    "    ax[idx].imshow(correlation_map)\n",
    "    for y_max, x_max in np.argwhere(correlation_map > tolerance * correlation_map.max()):\n",
    "        ax[idx].plot(x_max, y_max, \"w*\", markersize=12)\n",
    "    ax[idx].set_title(f\"Échelle = {scale},\\n corrélation max. = {correlation_map.max():.1f}\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "On peut alors considérer les détections avec les corrélation croisées les plus hautes pour toutes les échelles. Cela permet même de détecter des objets de différentes taille dans une seule image."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Avec scikit-image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Vous avez l'habitude, `scikit-image` fait tout le boulot pour nous. Il suffit d'aller chercher la fonction appropriée dans le bon sous-module. En l'occurrence, il s'agit de `skimage.feature.match_template`, qui implémente une version efficace de l'algorithme présenté ci-dessus."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from skimage.feature import match_template\n",
    "\n",
    "result = match_template(image, template, pad_input=True)\n",
    "y, x = np.unravel_index(np.argmax(result), result.shape)\n",
    "\n",
    "fig, ax = plt.subplots(ncols=3, nrows=1, figsize=(16, 6), gridspec_kw={\"width_ratios\": [1, 3, 3]})\n",
    "ax[0].imshow(template, cmap=plt.cm.gray) and ax[0].set_title('Template')\n",
    "ax[0].set_axis_off()\n",
    "\n",
    "ax[1].imshow(image, cmap=plt.cm.gray) and ax[1].set_title('Image')\n",
    "ax[1].set_axis_off()\n",
    "h, w = template.shape\n",
    "rect = plt.Rectangle((x-w//2, y-h//2), w, h, edgecolor='r', facecolor='none')\n",
    "ax[1].add_patch(rect)\n",
    "\n",
    "ax[2].imshow(result) and ax[2].set_title('Résultat de match_template')\n",
    "ax[2].set_axis_off()\n",
    "ax[2].autoscale(False)\n",
    "ax[2].plot(x, y, 'w*', markersize=15)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Cet algorithme fonctionne bien entendu de la même façon sur des images en couleur :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "airport = io.imread(\"images/orly_gmaps.jpg\")\n",
    "plane = io.imread(\"images/avion.jpg\")\n",
    "result = match_template(airport, plane, pad_input=True)\n",
    "h, w, _ = plane.shape\n",
    "tolerance = 0.7\n",
    "\n",
    "fig, ax = plt.subplots(ncols=3, nrows=1, figsize=(16, 6), gridspec_kw={\"width_ratios\": [1, 3, 3]})\n",
    "ax[0].imshow(plane, cmap=plt.cm.gray) and ax[0].set_title('Template')\n",
    "ax[0].set_axis_off()\n",
    "ax[1].imshow(airport, cmap=plt.cm.gray) and ax[1].set_title('Image')\n",
    "ax[1].set_axis_off()\n",
    "ax[2].imshow(skimage.exposure.rescale_intensity(result, out_range=np.uint8))\n",
    "ax[2].set_axis_off()\n",
    "ax[2].set_title('Résultat de match_template')\n",
    "ax[2].autoscale(False)\n",
    "for y, x, _ in np.argwhere(result > tolerance * result.max()):\n",
    "    rect = plt.Rectangle((x-w//2, y-h//2), w, h, edgecolor='r', facecolor='none')\n",
    "    ax[1].add_patch(rect)\n",
    "    ax[2].plot(x, y, 'w*', markersize=15)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "On voit tout de même apparaître quelques limitations du *template matching*:\n",
    "- pas d'invariance à la rotation (il faudrait tester le *template* selon différentes orientations et pour espérer détecter des avions positionnés différemment),\n",
    "- détections multiples d'un même objet (un algorithme de *non-maximal suppression* (NMS) peut résoudre ce problème),\n",
    "- sensible à des variations légères d'apparence."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Descripteurs d'image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Histogrammes de couleur\n",
    "\n",
    "Du point de vue perceptuel, une image est avant tout caractérisée par les couleurs qui la composent. Une façon de « résumer » rapidement le contenu d'une image est donc de lister les couleurs qui y apparaissent et le nombre de fois que chaque couleur apparaît dans un pixel. Ces valeurs sont représentées par les **histogrammes de couleur**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import skimage\n",
    "from skimage.exposure import histogram\n",
    "\n",
    "cat = skimage.data.cat()\n",
    "\n",
    "fig = plt.figure(figsize=(12, 6))\n",
    "fig.add_subplot(121)\n",
    "plt.imshow(cat)\n",
    "fig.add_subplot(122)\n",
    "for idx, col in enumerate([\"Rouge\", \"Vert\", \"Bleu\"]):\n",
    "    hist, hist_centers = histogram(cat[:,:,idx], normalize=True)\n",
    "    plt.bar(hist_centers, hist, alpha=0.5, label=col)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "L'image est ainsi décrite par la combinaison des 3 histogrammes :\n",
    "$$\\mathbf{v} = h_\\text{rouge}::h_\\text{vert}::h_\\text{bleu}$$\n",
    "\n",
    "où $::$ désigne la concaténation de vecteurs. Pour des histogrammes complets à 256 couleurs, cela réduit donc une image de dimensions $W \\times H$ à un vecteur de $3 \\times 256 = 768$ valeurs. Cela permet ainsi de faire une recherche d'images : pour une image requête $I_r$, je peux trouver l'image de ma base de données $\\{I_1, I_2, \\ldots, I_n\\}$ qui a le profil colorimétrique le plus proche."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Voyons un exemple sur une recherche de logo :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "database = [io.imread(f\"images/wow{idx}.jpg\") for idx in range (5)]\n",
    "print(\"Base de données :\")\n",
    "fig, axes = plt.subplots(nrows=1, ncols=len(database), figsize=(16, 6))\n",
    "for img, ax in zip(database, axes):\n",
    "    ax.imshow(img) and ax.axis(\"off\")\n",
    "plt.show()\n",
    "\n",
    "query = io.imread(\"images/wow_petit.jpg\")\n",
    "plt.imshow(query) and plt.axis(\"off\")\n",
    "plt.title(\"Image requête\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "À l'aide de la fonction suivante, on peut calculer les histogrammes pour toutes les images de la base de données ainsi que pour l'image requête : "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def color_histogram(image):\n",
    "    hists = [histogram(image[:,:,idx], normalize=True)[0][1:-1] for idx in range(3)]\n",
    "    return np.concatenate(hists)\n",
    "                       \n",
    "db_hists = [color_histogram(img) for img in database]\n",
    "query_hist = color_histogram(query)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "On calcule les scores (ici on utilise la distance euclidienne mais pour les histogrammes, il serait préférable d'utiliser la [divergence de Kullback-Leibler](https://fr.wikipedia.org/wiki/Divergence_de_Kullback-Leibler) qui est nettement plus adaptée) :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "scores = [np.linalg.norm(db_hist - query_hist) for db_hist in db_hists]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Il ne reste plus qu'à trouver et afficher l'image pour laquelle le score est le plus faible :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# on trouve l'image pour laquelle le score est le minimum\n",
    "closest_idx = np.argmin(scores)\n",
    "\n",
    "fig, axes = plt.subplots(nrows=1, ncols=len(database), figsize=(16, 6))\n",
    "for img, ax, score in zip(database, axes, scores):\n",
    "    ax.imshow(img)\n",
    "    ax.set_title(f\"Distance: {score:.3f}\")\n",
    "    ax.axis(\"off\")\n",
    "plt.show()\n",
    "\n",
    "plt.imshow(database[closest_idx])\n",
    "plt.title(\"Image la plus proche de la requête\")\n",
    "plt.axis(\"off\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "L'inconvénient de ce descripteur, c'est qu'il ignore complètement la structure de l'image. On peut mélanger les pixels et obtenir exactement le même histogramme. Pour l'image de Chelsea le chat, cela donne :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "shuffle = np.reshape(cat, (-1, cat.shape[2]))\n",
    "np.random.shuffle(shuffle)\n",
    "shuffle = np.reshape(shuffle, cat.shape)\n",
    "\n",
    "fig = plt.figure(figsize=(12, 6))\n",
    "fig.add_subplot(121)\n",
    "plt.imshow(shuffle) and plt.title(\"Image mélangée\")\n",
    "fig.add_subplot(122)\n",
    "for idx, col in enumerate([\"Rouge\", \"Vert\", \"Bleu\"]):\n",
    "    hist, hist_centers = histogram(shuffle[:,:,idx])\n",
    "    plt.bar(hist_centers, hist, alpha=0.5, label=col)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Ce qui est exactement le même histogramme que précédemment. Autrement dit, ce descripteur ne comporte que de l'information concernant couleur mais pas du tout concernant la structure spatiale de l'image. On va donc s'intéresser à deux autres types de descripteurs qui prennent en compte l'organisation d'une image, c'est-à-dire :\n",
    "- la texture (les structures haute fréquence),\n",
    "- et la forme (les structures basse fréquence)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "*Note* : les histogrammes de couleur que l'on a calculé ne sont pas tout fait des histogrammes de couleur 3D. Ce sont plutôt des triples histogrammes d'intensité RVB. Le problème fondamental ne change pas."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Descripteurs de forme : histogrammes de gradients orientés\n",
    "\n",
    "Les histogrammes de gradients orientés (HOG) sont des descripteurs d'image particulièrement populaires pour la détection d'objets. Ils ont été proposés par [Dalal et Triggs en 2005](http://ieeexplore.ieee.org/document/1467360/) notamment pour détecter des piétons dans des images."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "L'algorithme fonctionne en quatre temps :\n",
    "1. calculer les gradients de l'image selon $x$ et $y$\n",
    "2. calculer les histogrammes des intensités des gradients pour chaque orientation\n",
    "3. normaliser les histogrammes dans un bloc\n",
    "4. concaténer les histogrammes en un vecteur de caractéristique"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "L'étape 1 se fait de la façon habituelle (filtre de Prewitt) et permet de calculer la norme et l'orientation des gradients. Pour une image couleur, on utilise le gradient de la couleur qui est le plus intense localement (ce qui permet d'être invariant à un changement de couleur)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "L'étape 2 divise l'image en cellules (*cells*) sur lesquelles calculer l'histogramme. Pour chaque cellule on calcule l'histogramme des intensités des gradients selon $k$ orientations (par exemple, $k=4$ pour considérer les orientations à 90° haut/bas/gauche/droite, $k=8$ pour aussi considérer les diagonales). Cela produit pour chaque cellules autant d'histogrammes que d'orientation. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "L'étape 3 considère un bloc, c'est-à-dire un ensemble de $(c_x, c_y)$ cellules. Pour un bloc, on normalise les histogrammes des cellules de ce bloc en divisant par la norme euclidienne de l'histogramme global du bloc."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "L'étape 4 collecte les histogrammes des cellules de tous les blocs et les concatène. Une cellule peut appartenir à plusieurs blocs, auquel cas son histogramme apparaît donc plusieurs fois dans la caractéristique finale (mais avec des normalisations différentes)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Voyons ce que cela donne sur un exemple :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "from skimage.feature import hog\n",
    "from skimage import data, exposure\n",
    "\n",
    "image = data.astronaut()\n",
    "\n",
    "descriptor, hog_image = hog(image, orientations=8, pixels_per_cell=(16, 16),\n",
    "                    cells_per_block=(1, 1), visualize=True, channel_axis=2)\n",
    "\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4), sharex=True, sharey=True)\n",
    "ax1.axis('off')\n",
    "ax1.imshow(image, cmap=\"gray\") and ax1.set_title(\"Image\")\n",
    "\n",
    "# Amélioration du contraste\n",
    "hog_image = exposure.rescale_intensity(hog_image, in_range=(0, 10))\n",
    "ax2.axis('off')\n",
    "ax2.imshow(hog_image, cmap=\"gray\") and ax2.set_title(\"HOG\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Ici, le paramètre `visualize=True` permet d'obtenir la visualisation à droite. En général, on se contente du descripteur :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "descriptor.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "L'image $512\\times512$ devient donc décrite par un vecteur de taille $8192$ qui contient (on l'espère) toutes les informations utiles concernant la forme des objets contenus dans l'image. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Applications des descripteurs : reconnaissance d'image\n",
    "\n",
    "Une fois les descripteurs d'image obtenus, leur principal intérêt est de pouvoir reconnaître automatiquement le contenu d'une image ou de retrouver une image similaire (comme nous l'avons fait avec les histogrammes de couleurs). La reconnaissance d'objet à partir des descripteurs ou des caractéristiques extraites de l'image passe généralement par des algorithmes d'apprentissage automatique (*machine learning*), ce qui sort du cadre de ce cours.\n",
    "\n",
    "Néanmoins, on peut évoquer l'algorithme le plus simple : la classification par $k$ plus proches voisins."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## k-nearest neighbors classification\n",
    "\n",
    "Imaginons une collection de $n$ images $I_1, I_2, \\dots, I_n$. Ces images sont réparties dans $p$ groupes (par exemple, les images de chat, de chien, de raton-laveur, etc.) que l'on représente par une étiquette $y_i$ (un entier).\n",
    "\n",
    "Grâce aux descripteurs $v_1, v_2, \\dots, v_n$ calculés sur les images, on peut calculer une distance entre chaque paire d'images, définie par :\n",
    "$$d(I_i, I_j) = \\| v_i - v_j \\|$$\n",
    "\n",
    "À partir de là, admettons que nous recevons une nouvelle image $I_{n+1}$ dont **on ne connaît pas le groupe**. On aimerait savoir quel animal est représenté dans cette image. Pour ce faire, on peut aller rechercher dans la banque d'exemples les $k$ images les plus proches (au sens de la distance $d$) de $I_{n+1}$. On rassemble ensuite les $k$ étiquettes correspondantes et, par un vote majoritaire, on décide que l'étiquette $y_{n+1}$ doit être l'étiquette majoritaire parmi les $k$ voisins trouvés."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Commençons par charger nos collections d'images. On les redimensionne en (512, 512) pour plus de simplicité. Nous avons ici deux groupes : les images de voiture et les images de vélo."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from skimage.transform import resize\n",
    "\n",
    "size = 256\n",
    "\n",
    "bikes = [resize(io.imread(f\"images/velo{idx}.jpg\"), (size, size)) for idx in range(1,6)]\n",
    "cars = [resize(io.imread(f\"images/voiture{idx}.jpg\"), (size, size)) for idx in range(1,6)]\n",
    "images = bikes + cars"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(12, 6))\n",
    "axes[0].imshow(bikes[0]) and axes[0].axis('off')\n",
    "axes[1].imshow(cars[0]) and axes[1].axis('off')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Calculons les descripteurs HOG pour toutes ces images :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "bikes_hog = [hog(img, orientations=8, pixels_per_cell=(16, 16),\n",
    "                    cells_per_block=(1, 1), channel_axis=2) for img in bikes]\n",
    "cars_hog = [hog(img, orientations=8, pixels_per_cell=(16, 16),\n",
    "                    cells_per_block=(1, 1), channel_axis=2) for img in cars]\n",
    "\n",
    "# étiquettes (1 pour les vélos, 2 pour les voitures)\n",
    "bikes_labels = 1 * np.ones((len(bikes_hog), ))\n",
    "cars_labels = 2 * np.ones((len(cars_hog), ))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "On créé la base de données :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "dataset = np.concatenate((bikes_hog, cars_hog))\n",
    "labels = np.concatenate((bikes_labels, cars_labels)).astype(\"int\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Chargeons une nouvelle image. On la redimensionne en (256, 256), on calcule ses HOG puis on calcule les scores pour chaque image de la base de données avant de trier les images par distance croissante :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "query = resize(io.imread(\"images/velo6.jpg\"), (size, size))\n",
    "plt.imshow(query) and plt.show()\n",
    "\n",
    "hog_query = hog(query, orientations=8, pixels_per_cell=(16, 16),\n",
    "                    cells_per_block=(1, 1), channel_axis=2)\n",
    "scores = np.array([np.linalg.norm(hog_query - hog) for hog in dataset])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Les images les plus proches sont les suivantes : "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "closest = np.argsort(scores)\n",
    "print(closest)\n",
    "fig, axes = plt.subplots(nrows=1, ncols=len(closest), figsize=(20, 8))\n",
    "for idx, img_idx in enumerate(closest):\n",
    "    axes[idx].imshow(images[img_idx]) and axes[idx].axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "On peut conserver seulement les $k$ plus proches, par exemple $k=3$. On cherche alors l'étiquette majoritaire des les étiquettes des images les plus similaires à la requête :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "k = 3\n",
    "closest_labels = labels[closest[:k]]\n",
    "print(f\"Étiquettes des {k} images les plus proches : {closest_labels}\")\n",
    "print(f\"Étiquette majoritaire : {np.argmax(np.bincount(closest_labels))}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Vous pouvez tester avec l'image `\"voiture6.jpg\"` pour un autre exemple."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Détection de visage : méthode de Viola et Jones\n",
    "\n",
    "Pour terminer, parlons d'une application très courante lorsque l'on cherche à construire un outil interactif : la détection de visage. La méthode la plus populaire est celle introduite en 2004 par [Viola et Jones](https://www.merl.com/publications/docs/TR2004-043.pdf).\n",
    "\n",
    "Pour détecter rapidement des visages, Viola et Jones proposent des caractéristiques rapides à calculer : les caractéristiques pseudo-Haar."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Dans une fenêtre rectangulaire, ces caractéristiques vont calculer la différences des sommes des pixels contenus dans les régions blanches et dans les régions noires. Les régions sont rectangulaires et peuvent être définies de différentes façons :\n",
    "![image.png](attachment:image.png)\n",
    "\n",
    "Par exemple pour (1), la caractéristique est obtenue en calculant :\n",
    "- la somme des pixels dans la partie blanche à gauche,\n",
    "- la somme des pixels dans la partie noire à droite,\n",
    "- la différence entre la somme « noire » et la somme « blanche ».\n",
    "\n",
    "Viola et Jones utilisent toutes les combinaisons pour construire un ensemble de caractéristiques."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Pour pouvoir calculer plus rapidement ces caractéristiques, Viola et Jones proposent d'utiliser l'image intégrale, définie par :\n",
    "\n",
    "$$II[x,y] = \\sum_{x'\\leq x, y' \\leq y} I[x', y']$$\n",
    "\n",
    "c'est-à-dire que le pixel $(x,y)$ de l'image intégrale est obtenu en sommant tous les pixels au-dessus de lui et à sa gauche dans l'image initiale. En utilisant l'image intégrale, le calcul des caractéristiques de pseudo-Haar se fait en temps constant."
   ]
  },
  {
   "attachments": {
    "image.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbgAAABwCAYAAABsBmXeAAAABHNCSVQICAgIfAhkiAAAIABJREFUeF7tXQeYVEXWrdcTJSNgAIQBEUwzjAQxIEExL2YwobCiqJgwBxCGIArGVde0BnYx6xrXXzASFBNJskhGlJwkTeiu/5yafm1PT+fp6TB97/ed6X7vVXrnvalb99ataqVEhAFhQBgQBoQBYUAYEAaEAWFAGBAGhAFhQBgQBoQBYUAYEAaEAWFAGBAGhAFhIJYMWH4Kq9ugQYNb69at29TPNTkVQwZKIevXr5+IIn+MYbFSlDAgDCQHA+2aNm36QkZGRsPkaE7NbYXL5Srbt2/fY1u2bHnV+y4rKTg8kFcHDhzYNz8/P6vm0pEcd+Z0OtUdd9yxYd26dYehRX8mR6ukFcKAMBALBtCX/vTqq692KiwsjEVxUkYQBnbt2qW6du26Zc2aNa2RbKedNNNPnk5XX311VosWLfxcklOxZmDixIkOKLg8lDs/1mVLecKAMJA4BuCgadqhQwdVv379xDUiTWpu2LChat68uRMKrgFu2aPgHGly/0l7m5ZVyYhO2rZKw4QBYUAYSCUGRMGl0tOStgoDwoAwIAyEzYAouLCpkoTCgDAgDAgDqcSAKLhUelrSVmFAGEg7Bj777DO1//77+wWC1MLm491331WXXHKJSf/jjz8aJFLi0QZ/QSaJvGepO/4M1EaVZwJnAIcCBwI5cWiGE3VsBNYBXwIfA3/Eod5oquBAsDNwHpAPNAPqRVNQFHm2Iw95mQ186P7UUZQTryyNUVFvoBfASLUDgHj0M6WoZwOwCpgE/B+wA0h56dmzp1q+fLm5j5dffllR4b355pvmOCcn/H9VRBmqNm3amHzvv/+++Tz22GPNZyIkHm2Ix4uXCO6kztAMHIQkQ4GBwH7u5Lvx+RuwLXT2KqfIRQntgBOBvsCzwCfA/cDPVS49NgXw/2MAQJ7y3EXaHSmVc3UrGipWDjioVM8GyM2vwGjgNcAFJIu0dbfrQnxmuBvF94jKOR5LYGqhHsbjdwOuBIoBcjQSWAOkrGRlZSlGCVJq1aqlvI95bvDgwcYye/jhh9VLL72k9uzZo8aMGaOWLl2qGA1/9913Kyz7UogwNFYbQ+o/+OADUx7PX3bZZer7779XTz75pMI6MnXppZeqAQMGmOu0ELt3766w3MGUO3bsWPXVV1+pSZMmKUaIjho1SmGdn8IaNHPthx9+MHUOGzZMtWzZUn3zzTdq9uzZKjMzU/3vf/8z55jnl19+qdQGU2F1C9ZuLFm9erUWiQ8DZ599NjtKdmDxlOtRGTsddtCfA+wQOPJOhFC5svOeAFB50LL7JxD+0FSp15HeezEtFdPLAJV4tFKAjAsBcrQauA84OtrCqpiPobbtgRHA7wDbRIuOSiUSmYDE3pxQgb4AVGVNEJXZeMB+dnwW5wP0DCRCGCZOPxw9AuRpD3A3EG64Mjk+zqfhHOB09zkX8rBJkybrtm/fHtOO5JlnntHoMyqUCYWiu3TpoqHc9O7duzUUkn7ggQf0ypUr9dNPP62hxEz61157zeSFEtNQYAa//fab/vnnnzVC7PXHH3+sp0+fro8++mj9yiuvmDzM26lTJz116lR9++236/3220/fe++9esqUKbp9+/aedBdeeKG+9tpr9dy5czWUqz700EM11vmaOpln+PDhGkpUn3rqqXrIkCGV2lDhhqI8OOGEE2jBB3+XRcFFyW6U2eKs4NgZPQfwH38ZQDdSMgktui/c7fsGn+EqXVo1Hb1uhAq83IcT3d3RxcYBQAnAztG2cKMrLba56qC4IqAM2ApE8gxpGdPCsWUAvnzkdRzpV7pp6Q7k+zQDoBJOJjkejVngbh8Vb3YYjXsFaW7xSkercCkQyYDLZI+ngps8ebKnB5o2bZouLi7Wmzdv1m+99ZbGOjxzzVZw/H7PPfcYULCxhx43bpz5Tnn99dd1r169zHcqOMzdme8zZ87U2OVKY9cQcwzrTt9333161apV5jzW/ZnzFLhC9bfffmvqpCK0hUrYVtDebfAkqMIXfwpOXJRhvPE1KMk43Mu1wFTgImBzkt3bL2gP5wIfBW4G6Ec5GaCiCSbM1wqYBdQF7gR6BMsQ5Bo7xXcBKjgquilB0ibi0i5UWgQwQoCdNuflTgDCcesuQbo8YBpAFzEtk7OAaISDJdZ/OvAfYBBAt2AyyXdoDJ8n23mpu31/D9FA+12yk/F/5i533hBZE3f54IMP9lQOa8u4HWGVqcMO4yZJwQWWnnE7vv32256E9lwdTxx4IL3kGB1kZ6t69eope+0u3Y5QdmZ+EApVHXfcX4YvF7eXlJT/27IdtuTm8rWLn4iCix/Xia7pMjTgduAHgEpkX6IbFKB+Wib2CJpK7glgcIC09ml2SnnuA97jRCCaeRcGRFBhsA3kaKa7zGT8YBDFOQBdzHTH0X3q2cEhQINp6ea5r5Fj3ivPRSMjkYmuZc5zDQBoxSWjcKDC4CDyNQDgM6ULPJDwXbJ76lPxnX1k+YRVoBxJdH7RokXqscceU7/++qtq1KiRgstQffRRcCOdyvHMM880822UZcuWGYQrzE/FxXk8Kj0Kg2A4v/f777+HW0y1pEuogoN/WM2fP7/SjcFfq5o1a1bpfLgnbrzxRkMszGMF/6/JxolVjkCINBSSwJEoJ/z5zx5IuZ3r5oYdX6LlNjSAc5O0DJ4B6GoKJHQhMe1BAANWOgVKGOL8MFxvAnAwEEi5cahMK3g0wPnCRAotMSoqBujQamUQSjChBdcDaATQkjk2WOIg11riGntDcnQ1EEi5cWDykxtBiqv2S3xOfC9o4Y8CaNExAMaf8F3KAxzAA8BV/hIl6zkqmrKyMjVnzhzT9z300EMKrsNKiobBKt99952i9QYXpQkqYURl7dq11fXXX6+uu+66sG/x8MMPV5i3U/fff7/q37+/grtUPfvsswrzckHLYFt/+ukn04ZWreiAib3wISZMGImzceNG1bp16wqoihL64osvTOjs+eefr6hAbWGUEdeBpKlwTop+gvuA9QE4qIPzVHCfBrge79PslNhBsvMcE6JydtzsdIvcafeESO/vMien+V9Ni+gNfwnc527E53+BRCs3u4nP4wvdlRwQUDkHE3be7EnomnwECGXxBSqrCBc4UiQXgQZLbXCtNUAFlwyyA42gdb8/wMFAIKFFy9G1PcgJNrAKVEa1nT/yyCPVGWfQufCX0B3JdXIU9qUIRFHjx483n4888oi64YYb1HvvvWeWCNh5qYgOOOAAhXk1xWUIdlpGXA4aNEhddVW5Xr/ooovUQQdx3AjiUAfrsoUKkaDLklYiIykRaKLmzZtnljIw2tO7TubzPr7iiitU48aNTRviJvEMMuFkI6OCAgnCWvW///1v/c4772iQZ5IxSgjhpiYKhxOYWEuh9+7d6yli27ZtZmIV60XMJCsFhGuEtOqbb75ZM0J01qxZGia4/vzzzzUsSJMGGx5rbHysYVqbOuIlcQoyWYwXiD4Hzpv4Cq07dpK+YAeVBdTyzRDj4+ND1PFvXKcyKZ8I8F85raoNwFTA8p8k5FlaP1SmbI8/YWfuy5H9317fX4YYnaNCuhzgXFmgASnnwdh27+AIf9WznVuAGUHK8pfP+xwHQhxAfBIgId8bX554zPeMEzDVNQnDdnHO9GKgvLf330D2phzkBeKSuZYDVGy0dKOW6ggyiVe/lIr1+AsyCfaQo36wsciIcFWz/mLhwoXG1cgJTCgyBaWlLr/8ctW7d29jhtOdec0115gq165dqxAuqzjJyvUYRx11lPElQ3mpHTt2GFN469at6sUXX1T9+vVTiOIxpjvN5OOPP96MPDjS4foOpq8u4UiHbSRwP1QidKl1B2iFRCqHIwNnhw8JkLEdzjPN+4A/q4Mj8X7AFz7YjuPTAI70q0PY4fUA6A5t7K+CWqM6HeMoaETrzJF5aP3bc8ce073e/Ud19pOWrtd3AXbwft1l2UXt/1NrZGd2gIGErtvfgO8DJOCwmVGL3jyxbZR5AfJU9TRdiKyPCrw/QNdaRRlzcrPa9xWUKYe106qd1Z8cEeqfPdjh+wpf6v8AQwCX70UeZ41s/0jOqMJgc558J/js+M75Ew73Ganp+z6V4RzrJWIt7MemACcDRwNcQtHAt5K6Yzp2zWhei27KA3O6NbqePNUZ2+VI33Q4puunCOBgwL+M7do2p6jgHTX0mGj+Z/2XKWern4F4W3C4I3ZIHsD/awYPxxxzjP700089A4lu3brpN954w1hgTE8LjDJjxgyNNpvvUHQaiwg9eW699VYTykqhtUiLj8I6zjnnHE+4K8um1WcL0z733HOe41h/ge+7wj173f+kKJ8wXU8s8xfgCp8y6J7ktYp+jb8ScVRPZeZP6K5c6M5LS2488D/gUcC2Wtri+7+AfwCXAlSmlB7Aq8ATQNPyUxX+0vphJ8m6/a5dyS7K35F939GGK8dRDXT2iAKdXVSw109ZIU9lDy+4tzx/4baskflP+iiAhiiAHT47/0DC+2ZQha8MxAkGMgx3XyAHnNVnoIs9x5Xjvk4+zgcucaelS5H8MC0Vh688hRMsn0Lrmy5FKhePQHG/wftytKun4SvS2cPyy3ka1uEC73Thft9vTMFxhqeR+fuyR7Z/Uw3t5jtwYufPZ+J73q6iK75M91Mf3X4/uJGPTx5zbvVjgC5T27twEr7Tai8COHK1FRV5fRMYCfiusSvAOXJuC63LE72OzdfskQWurGvamPcpo+dBhqecosKVvunCPc4ZXrCkvIz2v2QNL+zvm08suFj3nMHL82fBZfo+lHgf00ds+3tZN+fPGHKKxYdmopJb01Dwy9dq8eLFCjdhfl8JCwnNeQajwKVovtNqW7JkiSdwhSv3bd+0SeAltNjscFf6gCdMmOCZo6Olx3JSSC5CW+cCVDbspBl0QOVxF2Arl3VB7ocd8N1e14vxnR0vraKeAK0azl+Q6AEA5zBuAJiGHRQtp10A07NDzgWoDOlaY5uoKAsBdi620Bohgls/WRic75eh9J+lXlkj/1oyat6DOcMLb9LKdbClrZuyNm29QY3M/yGjRN2z74H5HKnTtRmMI1ZKnmkh2ML75f2PAviiHgbcCtBS5Av6EsDOnFYwuRsDjAPI1TsAleZogC8budgGeM9Z8RqVJ4XKcjPgV8Fb9eEI0FrpXWXKfI9S9g6b933WyILplsuCotEXZ2dv72uhI3c4M0bvfWAO5yb5PnEwQKs5kDTHBe/3iZYx75cuZMpy4BXgvwCVIRX5eQCtqxcBDtJYD9/hyQBdsH2AIQAjR6kA+Sxs4TvE8xxcnQlQKfL/obLUyzbn9M6qvU8so7g067zsrNJFYL0txhYTMNh4Rln6nRKN/7uieRsrVy5n4s1AwhUco3nsbWjsm2cUELd/oQvRXoPBbWKorHgtkHBS88orr1Qnn0xPBRZPYR0GywklDG297bbbzAQtha5QrvdIgHRCnZ9HWS/dT/zHpnCE+3dgALC2/JSZowokVDzsXG0pcX9Zj09aDdsBds50PbFjOQZgZ9QZYOdsW57s8Cm0UH4GqNQodQB2+MvcxxF9WLXxmqLjdkt2VlFBVBxBuW2CiqW7D/0Q5mC0dbwrU03NvCJvS9nEVTwdjCNep2Ly5om9JDsyNo6dOIWWxqkAFZv9EtGaoquY+an4BwNHAXUBWmTkczHAjtxbwbFsXr8fONcNfPiR2m6lthtNooLLLBsBnq73kzL0KVcG3gfeEnS+C1RZ6ghnpvN1WEDPl4xfvFbtLuWAIPA/Yvk1b56opPlO0cVL4RzelQAHT1RcHBgcAJwNUBH+CFBmuD/74ZO8dAE2Ab3c530/+J7tD7Asvm+VBk9WbXd/sKtcweGdOCja94mbt2ht7bCULv+/07oW3q/+2Q51JQYFk9SzwcYAvk2X4+pgIOEKzt9NUeH06NHDhLFS8ezcudPMkTEq0rbc/OXD6nv19ddfm7BXLkDs27evyWev76Di8ifYPkZxzo/5qUjhsjR7qR1ySCAvjL9SYnKOo86+UZREi8nX1cLRPkfDvMaRL12KgUaV7HxeCFHvv3CdvShdQZ8BjYBaAOuxxTYd2NF4D5HvwTE7puikGAZDXferaiHqea+OhiOs1bE+oOnhEe0otRx6ctnCPbY1YbtdA7WTipzur0BCBTUWoAX9FmC7ItmrUrlRAnFEa87XbcCO+v+At4ETAHvgYQqqIMXO8sOc8g5cK+ejpfsctC4jlqxc53Aofw5mjEDBUc0ZaxfK7T6covIIJhwYhXqfvkaaacB0wFZkgd4nDtjse+cjvMancrooydUXwEyACvRCoJKCU3yXKLnlPEGFbyzZF937pBy1EPC9d5W3XwJOtC0ZOvPFfa3aF+FVoaWaMGEkIz1RCK7ztGH06NFm2QD3g4yH0IvGKMxgkfFMw+UMjOSMtSRUwdWpUyfgbtjPP/+8YogqNw9lwAfdmAxnxd5pCtvCeHhwOBye46FDhyqGnjIUFfugmcAU++HSqsOcnAmNpdXovaL+H//4h+rTp49CtKb6888/jWI877zzYs11OOVxVOw98g0nD9Nw1Es3I/97+Q9+F2C7g25zF0LL5Vf3d98P+AFNSLe3UBlSSdHKoLCDJSmrgFcBKlQqZM69NQGYrjdABTgFYOf+X4D1fgh8AEQubreb42AaMhD4g9RD8yPmiHNLztLybargmv5du/TTJaPmjkOJ5Kypu2HGugsiB+GaN0/sSNmZkz/yfyzAe+XAgvfP85Q5AK04umqvBVjnfIAd9zcAuabbbbv7PD6MDAPIGy3j5u5zq/BZQU/zvO3CteqW60+rLGeXemh2xDypm9rUs/arfQPdncpybNeWc2Lp5r3D1FPLbKfeAFTHe+UAh5acP+Ggyvd9ooVrv0/0E7YDaInxO3tbKqPvAA4Qngc6AB0ByjSAyo/vOTke6f5uLkL4cjwCHO8+hnvVuIwriS9PeEQu9dDcyHlCyZz/xNuYC+Xv0i41KzvDunf3iPlflo/s8Mib8N8iccIdRjhtY8sff/yhEMNQLYok0F2edNJJZpNlLvoOJFyvF+x6oHzhnE+ogrN/8sFfQ7nwD+H8atOmTYqK0F6wzW1fsPeZJ4v3cd26dc0O1VSInF/zdjPSEiQojMD0Fs7jIVjFRGhS+RHVKbQMuXMABRGcJaj3R3ylG+e9KOqle6cr8E+Ao+udPmWwc6Wws2BH4StUqvwvYKfiLQ/igG4hKst+AMumVbIaeA24EeD/8hCAnTMtOSpXdr5UbEcA3wIsn8R7W3Q49Mj3+LbP+4T9HbrsNbVyb3t0tidkNMz+FW7FZdrhoEKJWMpK9JMZjowpDqfzjj2j583yKeAPHP8OsFMNJFTmVFq2VcZ05I1K6D/ASwC5IjfdAFpPUwBaErQ4HgCuAhYCVKR83gOBNwF20FSMnwO+0h0nCFs4yNhtH2SUOj5yZTjr67W7u1kNspSV5ZhGi6LYaS3yLSic4+xGtfEeuBZlOnLv3zPiJ96Dr8zGib8DfJ8+872IYyppwvd9GoRzkwBadvwkH3wflwFPAOSGli8HT1TqvwB02+JujFv3MYCWHi1hvnvewudA3qkg+f7x/4jvYAWxtP5QL99JxVvgqJc1B+/TevQTvu+Cbzb/x4ietJSzh8vKeL50y593eQ0A/Kf3Ocvobm6NxSmUc889V3GxNOMMuPs+B/YUpqF1w8E5z3OX/i+//NKkYxpuw0WPE71WXB/HaHME3Jkoc38/o7NhwwaTB0uvzHo2LupmVDm36SooKDAD/LZt26oLLrjAtI0R5lxPbE/dMPobwXgmWp39V8eO5eMPrmc+4ogjTH/NNrO92KDZ/OIAN9jg+mP2sZyKwobOZt0bAnDUxRdfbAwOWne8Nwp/7YDHjIPgTiz0pvH7hx/y3wP+f/zSgT1tZU5EKvGMogweE5MeV2OwDo4j4GBCXwxH2rZFFyxtpNfqIgMVITw9JgqOipodSCxlFApjJ9e+SoX26ROKJ3au7BxtS6lK1flkvhXHbdzn7sTn0FgWjrIOB8jRuCqXG5on+u3JU+z9SeVW3S3ue6iPT3ocsqp8TxULoOLjgKpOlcp9ZQCt1KASKIoSHbbmLwBgMbbGT9yYnfy5Hpc7+ds7/7P3Y9Q3FJHpCBnZzQ2MGSXOTYqhLPSCBQtMVDkG9hpKRePnbvQpp5xiduyncI2xnZ9R4lBuGkF9GkukNBSa2ZCZa4Sx2Fqffvrp+sEHH9RQOLpz5876lltuMZswYzCuEctg1iFDoZlIdbab5/ETOp62QcHpm266ybSNvx7Ae+T6ZW7CzF8PwIYeJh+mjMwaZUa346dzzC8OsA3kg8L8bBumqEwUPTd45jX+QgEMFI1twcwvIPgTf1GUlR6QKDh/1FXfuRgouErP0M+J53COnVJnP9eqemoMCuDImaPo26tamE9+BmlsAGj1VLfQSqKSoKUQa+mJAml9kCe6G9l5x1JoPbLthbEsNEhZdKvSU3BAkDTRXOIghNYwLTsOlujWjaW0Q2H0JLwdy0IDlRVIwcHLpLFGV3NTCgosKI09JEMqOC57soWdPXYc8SybwnphcwkWnUa9Gp6hCgqOS6m4wYUt/HkdbpJB5QKLydMWbNHlUYpMS+WJ327TiGjXsOY8+T/55BPzMzoUKl/+hI4tVGKwAM0hlRProFCxIRLefOdmGpgm0rAGKyk4Kltb+B0R7p7jq6++usJSMM8FfPGn4BLqogz0Ysj5mDPwEEocAHCe4mTAHZEQk3rooiOqQ+gWZSd6U3UU7lMmLdwpwGCAA4KlMazza5TFOczqEEZg9gfow6EbNR7COTC6J0cD18awwhKUdWUMy/Muil4Gvv8OgIOyhAmD2ui2y8vLM24+uhsZL8AfG/UW34hx2yXINLCy1Asv0NsLUxRTOHRRUuz5rBUrVphjCrdDpLuRG2LYS6M4F0d3IoP26PKz4xr46RvjwGVb3FSDQX90PVIYoc4ybLHr5zGj3RnI4ivMy2ASWJ4mRoJgMKCvYA205xQ3i6ab9qmnnjLneC+MxQhXRMGFy1Rqp1uF5nOOjnNhdClyrijZ5Vx3O2kZMrggHnI3KqGl9T5AhbQjHpVWoQ5GMdDlxgHLvVUoJ9KsnCucDHBejVbphEgLSED6e1Dn39xtnZeA+j1Vci6LHTajBjl/Bteg2W2J+0TCEPGk84414Enms4VxBvbyKioTKkN7J3/OqzH+wFaYdkwB4xOoDClUFFREnM+zlZ6ncD9fmI/Bd0VFReYq62TQii0M9gsmVGSnnXaa2XmKQX0sjzEToYRtf/zxx82elxT8gGxYS7887QpVgVyvMQzQyqKy4BwQkcxCK3MisBlggMtf//XV22q6xcgTh8MMdGhQvdVVqfTGyE2rrRXAeSsGZMRT/o7KfgMYTHJ+PCuOoq6rkIfWJgNvhkSRP6ZZaL0xEISdN3ZUUl27djUKgwqLSo3BIFQ8vkF4DA5hUAmXTTHKHO470y5aWIw2p3CzeZbBDTFsoTLhtocMHGEd3BgDc3WmjnCFQSWMhqSiYhlPPPGE2dA5lFDxMTKdFh8DBqnkqNiwU5RZbxxK4P40gS20KqnAuWk+LclwRSy4cJlK/XR8mxjRR3cZ30wGJrBjtNdnRXSHXzQ7qmdWVkYlRVla6ny417qFrCMa4TCQEXKPAmzvxQCj7OIpdOceClwNfA+Qs6jnAKfnFb7swmJi7xvA4uJVPVYtoCs0WmH4Gufy8oDHAbpU4y0cvtPKpquSbRkBkLuyaBoyvVXn3i5dfL1vXp3tHNFj6UJG80YjOcj0IMAgHypjtjfhVjk3o+BPyjBykAqAkYtUZowspOLjHrqMLuTSKFostjASkVGQVAyMQOR6X1p+VGDs9LlMipGR2NLQWHAsz/4ZGu7UxF8QoGKi0BJjaD7dlN4h+r4/U1ZYWGjcnlwjPHjwYJOWViZ+pVthc3pTFpdlee8Y5X1M9yujMrnvLt2T3OeXCpK/WMBITrZ5xIgRnjb41o/gEvMLBTzPevmzPrz3qEWCTLynLav/e5yCTLzfh4Y4YKdEq+h34DogYktlal6Hoql5BboyOhRF8fIxuKA3wBB0tosLZI+OopxYZrkDhbGz5mQCIwbbRVM4+Nnry9H0VvnRdrLkhEN1uiQ5X8Vnl2jhYGABwOfGgUBfgMseIpLpLds/7csTj79qE3Tj50B10PfVH+AkEdtFt3OFQUagjLE8HyjIxO5V4FY0EYrhiPdeuow8tIVRlNi60BzCugpZFCzFkGlCJahKGd5tD1WP9/Vw6pQgk1i+valb1jY0/QzgcoCT7QyPfxJgJ7AK4Ej3TyCoPLh59Ykts6iXKsrq0pITcYZzWaEkFwmaAs2BrgAjJmlN0hJgJGNUliXyxUoeQUGca6JVQsuC4NqsOcB6gNGd7DyDyus7NsJLwviGvwQDUVoW4XBEi/ZA4GCAVhuVCYVzhEOBxe7jRH5wMNIZoEVK9+5bwB5gOrAWWAfQGg8q47et7dA0g7RUlMXFu0/DmdCTNeWL5vkucUEV30G+X7QyOQjgoCAqyxL5qk24Bi6crQR9GxBoviucuTR7ns63zEiOq1JGoLaHqr8qdVYoWyy4SMYVVU+bAAvO+3mzR+Gq0tcB+kLYYccbxaiTFiUVSCMgGYUd+DhgKRBvfuz6aB2NBdonI0HuNnGQMhD4GKBSSwRXHBgxKKkfELEl6b6PmHyEsuAi6T24Ro7Wmq8w3J7h/iKyTCAmL20NK4TKhfMnBIUKhtZC5aG0O4H9Mapxq0EHZWcO8jmt1peUvTB888oXfM/7OaabjVYQUTlW2E+GBJ7iHBBBq4udZjOgfjjt+VeztjPQ+VQ0dbXaPeiPX7uHkx9paHH/DvwVQhdmxgQk49o4WksEzVZan0TIuf5xTVrd1TArk+7NCrLMuXfc+A2/hRNFS5ctLetNvmXUhGMGovgTBqokaFtBf81JunP7uIb0AAAMuElEQVQhX7yka7E0qDoZ4EKciotxAtTWvU5DuH8qe33aZWf+oTavnBUgW004TctkWbg30jYrt5Lydliayr0mc0R6aL1R4RAh5YRa9Ta5sHeWr7TJrL0GCq6mc+V723IcIwZEwcWIyLQrJqPOS1bpVtPRa8vxBPabNeHXOqvO1LTjIsgNO6ycvtq5t75DZfZ2OfQe5XJ+mZHl2BwkS1pe2s+V+fgeq/h7c/MZjrHK6bqPX5vkujgPKiIMRMWAKLioaJNM3ZdPYwDBq2RiWquCB7qtXGC+i1Rk4KSV5RsWT219THP86syObquEJ3/vSKe1sxisQqiprQru6r5GePLHk5yLjIHgy88jK0tSCwPCgDAgDAgDScNAJQsOoaZbuY2L/fMFSdPSGtgQbq+DyCgOMhhIUNPkRdzQqYD3RF0RjstXh9a0uy3fVaQAt+U958bQeW5pJVKRAS69uAzgGkNbuAL5qYrJUvsoKytrBRZ0H4z9HiuuE0nt20rK1nOnFOxZSZ63ejewkoKDchuM0PV3scM015GIVCMDWLyI8YTFDv+3aqwmUUUzSrAHsDJRDYhjvRykMMytJZDo9XtxvO2oq+qGnFyewq3Raqxgg+NLsev/Y1jrxkhSkeploBhbmFVaP1tJwaHTnQsl16Z62yKl13AG9sf9MYw+HZQbHyV3OWEovyi30C82l0wcBcwNnTTlU/yGrbQqLX1I+btKoRuQObgUelgp1NQuaCsX/dqBA/z0v5AnhW4qSFN5v6197veIIOnT+VI+br4WsNiLr3PSmRC59+pjoJIFV31VSclpxAB3/uCuG6PS5J55v1wm8a80ud+q3CYHA9zb85aqFCJ5hYFwGBALLhyWJE2kDLATK1/TFGnO1EzP+63R80kxfCwcDET76wAxbIYUlQ4MiIJLh6cc/3vshCpnxr/ahNTIYCzOWUf9kzoJaXXiKj0OVYuCSxz/aVWzKLi0etxxuVkGEVwDVAjXjUvNiakkC9VeDFTetywx7Un2Woehgdy0WkQYqHYGZA6u2ilOuwq46e1HaXTX/GkhWesW/gP/b/hJJaUwUDUGxIKrGn+SWxgQBoQBYSBJGRALLkkfTGo1y7HflFb5fVKrzZG1to6V/UOnFbPWRJYr9qmntGzfUzlcjWNfcvKUaCkr5M81JU9rpSXJzIAouGR+OinSNvyawDPolPiL0zVUnHX3KCfD2hO6lm9aiy5HWmrfKJeV+W0NJdrclkPpW2vy/cm9xY8BUXDx47rG1tR9xeyiGntzuLGZrTu22Ossey3R95iZuSej1GWt6b5izj2JbovULwykAgMyB5cKT0naKAwIA8KAMBAxA6LgIqZMMggDwoAwIAykAgOi4FLhKUkbhQFhQBgQBiJmQBRcxJRJBmFAGBAGhIFUYEAUXCo8JWmjMCAMCAPCQMQMiIKLmDLJIAwIA8KAMJAKDIiCS4WnJG0UBoQBYUAYiJgBWQcXMWWSIR0YeFv1yT641dJpZZnZ1ymny7Nx9NetOwzJ0GWF3VbOGxCKh+lt21+hSl29y7p3vFJN+9kk39ijR53Fq7e9nllmvXTi2rkfhirj68LCBo7t+tPc3Lr9nSU7PMmn5nUocmhn45NW/3xjqDLkujCQrgxY6Xrjct/CQCgGprXKn2K5rG5O5frWcjjqW0o30No6RGdbF/RY+vP7ofJrbMoxvVX7fVq5XJa2vtSWKrOU40ylVVm3VXP5q9ZhyfS89gu00ke6lDXZcujaymnh18J141oqu0Pn1TPnhFWIJBIG0pABcVGm4UOXWw6PgTKnut2llGVZjq5QSvlUbg5L7wxHubEGjB5dUHKTlLZy8HkWyjhHa1cW9qJ6L7wWlKeyHBnDkd+Cgj1DudRJytKN0ajfRblFwqKkTUcGRMGl41OXew6LgZPXzJ/lgCLxTuzUOqItu3Jz6t7hnR9KTzesVXZ7WA1wJ+q6YvZ7sNx2VijH4XoykjIkrTCQjgyIgkvHpy73HDYDOlM/bSd2wPzak7kvon0guyz5ZqlDqxV2GdiUelH+ggUbwm6AJ6Pjdc9Xy1HadcX8hyMuQzIIA2nGgCi4NHvgcruRMXDSsnkPOxxWKXNph2PmWcuWVbCkwinNZVnj7XT45YWicPL4pnHWte6lgjXtUHoS3Z++aeRYGBAGKjIgCk7eCGEgCANQJGUul+tzzIGpkpKyiKw3u9juq35+Ht/3WQ71Z7eVc94NUl3ASz3nzt3ucqhZTODr9gyYSS4IA2nOgCi4NH8B5PZDM5CbW/9WRHhs6LVu4dehU/tPoZX1gdPK8LgZ/acKftZp0YpTK+j2DJ5SrgoDwgAZkGUC8h6kGgN8Z88Gmsaz4d3qNDh42q7tf0Rb5xG1a9ffWVZWuq64eE+0ZTBfVdsRRd38cdWFUeSTLMJAwhkQBZfwRyANiISBJk2a3NipU6cHzjvvvHqR5JO0kTNQUlLiGj9+/Ja1a9cWIneFaNLIS5McwoAwIAwIA0EZaNGixeTp06drkfgwcPPNN3P7FFrMIsJAyjEgc3Ap98jSvsHidYjjK+BwSBcRR7qlqhgzIG9vjAmV4oQBYUAYEAaSgwFRcMnxHKQVwoAwIAwIAzFmQBRcjAmV4oQBYUAYEAaSgwFRcMnxHKQVwoAwIAwIAzFmQBRcjAmV4oQBYUAYEAaSgwFRcMnxHKQVwoAwIAwIAzFmQBRcjAmV4oQBYUAYEAaSgwFRcMnxHKQVwoAwIAwIAzFmQBRcjAmV4oQBYUAYEAaSgwFRcMnxHKQVwoAwIAwIAzFmIDPG5UlxwkBaMNC3b1/PfVr4LZ3GjRurs846S519dvBtG/fs2aMGDx6sJkyYEJSn3bt3q8zMTJWTkxM0nVwUBoSBwAyIBReYG7kiDARk4J133lFdunRRffr0URdddJFq3ry5uuSSSxTPBxPu7ZiXlxcsibk2aNAg9fLLL4dMJwmEAWEgMANiwQXmRq4IA0EZOO2001R+fr4nDX5WRk2aNMkoPcqaNWvU5MmTVW5urjr//PNVnTp1VEZGhios5K/PlMvOnTvVJ598orZt26bOPfdc1axZM7VkyRKTd/bs2Wr58uXq0EMP9aSXL8KAMBA+A2LBhc+VpBQGAjKAH69RK1asUK1atTJpZsyYobp3765Wr16tPvvsM9WxY0e1Y8cORRflgAEDTJrNmzer4447TuHnf4wi69Chg5ozZ45av3692rp1q1FymzZtClinXBAGhIHgDIgFF5wfuSoMBGSAbkRaZZSVK1eqRo0aqdtuu80c33nnneqxxx4zlhvlb3/7m3rzzTeNG9OWhx9+WJ1xxhkmHQU/5mq+T5w40Vh5Xbt2NQpQRBgQBqJjQBRcdLxJLmFADRkyRLVp00aVlZWpRYsWqVtuuUV98803qlevXmrmzJnq0UcfVc8884xhihZau3btKrD2448/Ggvt1FNPNee3bNliAktEhAFhIDYMyH9TbHiUUtKQgSOPPNIzB8eAk48++sjMm3FujvNuo0ePVi1btjTM7N27V9WvX78CS7Vq1VL9+vVTdkRmaWmpoqtTRBgQBmLDgMzBxYZHKUUYUPvvv7+ZO6PQKvvqq69U69atjevxsssuU/Pnz6/Akp3mkEMOMZGV48aNM+5JCoNRdu3aJawKA8JAFRgQBVcF8iRr+jJAxZWdnV2BgMMOO8wElVCeeuopEzxCFyZdk71791ZnnnlmhfQ33HCDatGihVGCDE5hJOXQoUNNGlqBY8eOVW+88UaFPHIgDAgDwoAwUEMZgEL4DIoDnrzUEFhhuri42NNYRE5quCorNH7fvn0aC7uT8oYwz7gDr1Lw1es19F2T20p9BmQOLvWfodxBEjNQu3ZtT+ugyNTAgQPV4YcfXqHFsltJEj9AaVpKMyAKLqUfnzQ+lRjgvBoXgZ9yyimp1GxpqzCQsgyIgkvZRycNTzUGsrKy1OWXX55qzZb2CgMpy4AEmaTso5OGCwPCgDAgDARjQBRcMHbkmjAgDAgDwkDKMiAKLmUfnTRcGBAGhAFhIBgDouCCsSPXhAFhQBgQBlKWAVFwKfvopOHCgDAgDAgDwRgQBReMHbkmDAgDwoAwkLIMiIJL2UcnDRcGhAFhQBgIxoAouGDsyDVhQBgQBoSBlGVAFFzKPjppuDAgDAgDwkAwBkTBBWNHrgkDwoAwIAykLAOi4FL20aVnw10u1w75nbT4PXv8+oETte2NX41SkzAQOwas2BUlJQkDcWHgiMaNG39pWdZf2/THpdq0rESD53kbN27sibunohMRBoQBYUAYEAaEAWFAGBAGhAFhQBgQBoQBYUAYEAaEAWFAGBAGhAFhQBgQBoQBYUAYSGMG/h/h85IW49Y8VgAAAABJRU5ErkJggg=="
    }
   },
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "En pratique, les caractéristiques sont calculées à toutes les positions et à toutes les échelles pour une fenêtre de détection (par exemple $24\\times 24$). Cela donne de très nombreuses caractéristiques (plusieurs dizaines de milliers). Celles-ci sont ensuite traitées par *machine learning* pour décider si la fenêtre de détection contient un visage ou pas.\n",
    "\n",
    "On ne détaille pas la partie ML de l'algorithme mais, grossièrement, la méthode de Viola et Jones va passer chaque fenêtre de détection dans une cascade de 3 classifieurs. Le premier classifieur (étage 1) est rapide mais peu précis (il dit souvent oui), le deuxième classifieur (étage 2) est un peu moins rapide mais plus précis et le classifier final (étage 3) est moins rapide mais le plus performant. Cela permet de rapidement éjecter (en peu de calculs) les fenêtres sans visage.\n",
    "\n",
    "![image.png](attachment:image.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Concrètement la méthode de Viola et Jones est une méthode de recherche exhaustive : on va tester toutes les fenêtres possibles de l'image à toutes les échelles (dans une pyramide, comme pour ORB). Cette exhaustivité est possible car le calcul des caractéristiques de pseudo-Haar est très rapide et le rejet des (nombreuses et largement majoritaires) fenêtres négatives (sans visage) se fait de façon efficace."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Une version préentraînée de la méthode de Viola et Jones pour la détection de visages est disponible dans `scikit-image` : "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from skimage.feature import Cascade\n",
    "\n",
    "# Charge le fichier contenant le modèle préentraîné\n",
    "trained = skimage.data.lbp_frontal_face_cascade_filename()\n",
    "\n",
    "# Détecteur\n",
    "detector = Cascade(trained)\n",
    "\n",
    "image = data.astronaut()\n",
    "\n",
    "detected = detector.detect_multi_scale(image,\n",
    "                                       scale_factor=1.2,\n",
    "                                       step_ratio=1,\n",
    "                                       min_size=(60,60),\n",
    "                                       max_size=(120, 120))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Les paramètres importants du détecteur sont :\n",
    "- `scale_factor` : facteur d'échelle de la pyramide d'images,\n",
    "- `step_ratio` : facteur par lequel est multiplié le pas de la recherche exhaustive à chaque échelle. 1 correspond à une recherche exhaustive. En général, entre 1 et 1,5.\n",
    "- `min_size` : la taille minimale de visage à détecter,\n",
    "- `max_size` : la taille maximale de visage à détecter.\n",
    "\n",
    "Les résultats de la détection donnent :\n",
    "- `r` : la ligne du coin en haut à gauche du rectangle qui englobe le visage,\n",
    "- `c` : idem mais pour la colonne,\n",
    "- `width`/`height`: la largeur et la hauteur du rectangle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "print(detected)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "On peut dessiner avec `matplotlib` le rectangle :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from matplotlib import patches\n",
    "fig, ax = plt.subplots()\n",
    "plt.imshow(image, cmap=\"gray\")\n",
    "\n",
    "for patch in detected:\n",
    "    ax.add_patch(\n",
    "        patches.Rectangle(\n",
    "            (patch['c'], patch['r']),\n",
    "            patch['width'],\n",
    "            patch['height'],\n",
    "            fill=False,\n",
    "            color='r',\n",
    "            linewidth=2\n",
    "        )\n",
    "    )\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### À vous de jouer !\n",
    "\n",
    "Essayez ces algorithmes avec vos propres images. Notamment, essayez d'appliquer des déformations géométriques sur ces images. Les algorithmes sont-ils robustes ? Si non, à quelles transformations ne sont-ils pas invariants ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.9.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}