{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Computation and comparision of the bispectrum and the rotational bispectrum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We show how to compute the bispectrum and the rotational bispectrum, as presented in the paper \n",
    "\n",
    "- _Image processing in the semidiscrete group of rototranslations_ by D. Prandi, U. Boscain and J.-P. Gauthier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy import fft\n",
    "from numpy import linalg as LA\n",
    "from scipy import ndimage\n",
    "from scipy import signal\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.cm as cm\n",
    "import os\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Auxiliary functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def int2intvec(a):\n",
    "    \"\"\"\n",
    "    Auxiliary function to recover a vector with the digits of a \n",
    "    given integer (in inverse order)\n",
    "    \n",
    "    `a` : integer\n",
    "    \"\"\"\n",
    "    digit = a%10\n",
    "    vec = np.array([digit],dtype=int)\n",
    "    a = (a-digit)/10\n",
    "    while a!=0:\n",
    "        digit = a%10\n",
    "        vec = np.append(vec,int(digit))\n",
    "        a = (a-digit)/10\n",
    "    return vec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ALPHABET7 = \"0123456\"\n",
    "ALPHABET10 = \"0123456789\"\n",
    "\n",
    "def base_encode(num, alphabet):\n",
    "    \"\"\"\n",
    "    Encode a number in Base X\n",
    "\n",
    "    `num`: The number to encode\n",
    "    \"\"\"\n",
    "    if (str(num) == alphabet[0]):\n",
    "        return int(0)\n",
    "    arr = []\n",
    "    base = len(alphabet)\n",
    "    while num:\n",
    "        rem = num % base\n",
    "        num = num // base\n",
    "        arr.append(alphabet[rem])\n",
    "    arr.reverse()\n",
    "    return int(''.join(arr))\n",
    "\n",
    "def base7to10(num):\n",
    "    \"\"\"\n",
    "    Convert a number from base 10 to base 7\n",
    "    \n",
    "    `num`: The number to convert    \n",
    "    \"\"\"\n",
    "    arr = int2intvec(num)\n",
    "    num = 0\n",
    "    for i in range(len(arr)):\n",
    "        num += arr[i]*(7**(i))\n",
    "    return num\n",
    "    \n",
    "def base10to7(num):\n",
    "    \"\"\"\n",
    "    Convert a number from base 7 to base 10\n",
    "    \n",
    "    `num`: The number to convert    \n",
    "    \"\"\"\n",
    "    \n",
    "    return base_encode(num, ALPHABET7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def rgb2gray(rgb):\n",
    "    \"\"\"\n",
    "    Convert an image from RGB to grayscale\n",
    "    \n",
    "    `rgb`: The image to convert    \n",
    "    \"\"\"\n",
    "    r, g, b = rgb[:,:,0], rgb[:,:,1], rgb[:,:,2]\n",
    "    gray = 0.2989 * r + 0.5870 * g + 0.1140 * b\n",
    "    return gray"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def oversampling(image, factor = 7):\n",
    "    \"\"\"\n",
    "    Oversample a grayscale image by a certain factor, dividing each\n",
    "    pixel in factor*factor subpixels with the same intensity.\n",
    "    \n",
    "    `image`:  The image to oversample\n",
    "    `factor`: The oversampling factor\n",
    "    \"\"\"\n",
    "    old_shape = image.shape\n",
    "    new_shape = (factor*old_shape[0], factor*old_shape[1])\n",
    "    new_image = np.zeros(new_shape, dtype = image.dtype)\n",
    "    for i in range(old_shape[0]):\n",
    "        for j in range(old_shape[1]):\n",
    "            new_image[factor*i:factor*i+factor,factor*j:factor*j+factor] = image[i,j]*np.ones((factor,factor))\n",
    "    return new_image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Spiral architecture implementation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Spiral architecture has been introduced by Sheridan in \n",
    "\n",
    "- [_Spiral Architecture for Machine Vision_](https://opus.lib.uts.edu.au/research/handle/2100/280), PhD thesis\n",
    "- [_Pseudo-invariant image transformations on a hexagonal lattice_](http://www.sciencedirect.com/science/article/pii/S0262885600000366), P. Sheridan, T. Hintz, and D. Alexander, Image Vis. Comput. 18, 907 (2000).\n",
    "\n",
    "The implementation with hyperpels that we use in the following is presented in \n",
    "\n",
    "- [_A New Simulation of Spiral Architecture_](http://ww1.ucmss.com/books/LFS/CSREA2006/IPC8173.pdf), X. He, T. Hintz, Q. Wu, H. Wang, and W. Jia, Proceedings of International Conference on Image Processing, Computer Vision, and Pattern Recognition (2006).\n",
    "- [_Hexagonal structure for intelligent vision_](http://ieeexplore.ieee.org/xpl/login.jsp?tp=&arnumber=1598543&url=http%3A%2F%2Fieeexplore.ieee.org%2Fiel5%2F10652%2F33619%2F01598543), X. He and W. Jia, in Proc. 1st Int. Conf. Inf. Commun. Technol. ICICT 2005 (2005), pp. 52–64.\n",
    "\n",
    "For a more detailed implementation, see the notebook [Hexagonal grid](http://nbviewer.ipython.org/github/dprn/dprn.github.io/blob/master/docs/notebooks/Hexagonal%20Grid.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start by defining the centered hyperpel, which is defined on a 9x9 grid and is composed of 56 pixels. It has the shape\n",
    "    # o o x x x x x o o\n",
    "    # o x x x x x x x o\n",
    "    # o x x x x x x x o\n",
    "    # x x x x x x x x x\n",
    "    # x x x C x x x x x \n",
    "    # o x x x x x x x o\n",
    "    # o x x x x x x x o\n",
    "    # o o x x x x x o o"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# The centered hyperpel\n",
    "hyperpel = np.array([\\\n",
    "                [-1,4],[0,4],[1,4],[2,4],[3,4],\\\n",
    "                [-2,3],[-1,3], [0,3], [1,3], [2,3], [3,3], [4,3],\\\n",
    "                [-2,2],[-1,2], [0,2], [1,2], [2,2], [3,2], [4,2],\\\n",
    "                [-3,1],[-2,1],[-1,1], [0,1], [1,1], [2,1], [3,1], [4,1],[5,1],\\\n",
    "                [-3,0],[-2,0],[-1,0], [0,0], [1,0], [2,0], [3,0], [4,0],[5,0],\\\n",
    "                [-2,-1],[-1,-1], [0,-1], [1,-1], [2,-1], [3,-1], [4,-1],\\\n",
    "                [-2,-2],[-1,-2], [0,-2], [1,-2], [2,-2], [3,-2], [4,-2],\\\n",
    "                [-1,-3], [0,-3], [1,-3], [2,-3], [3,-3]])\n",
    "\n",
    "hyperpel_sa = hyperpel - np.array([1,1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now compute, in `sa2hex`, the address of the center of the hyperpel corresponding to a certain spiral address."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def sa2hex(spiral_address):\n",
    "    # Split the number in basic unit and call the auxiliary function\n",
    "    # Here we reverse the order, so that the index corresponds to the \n",
    "    # decimal position\n",
    "    digits = str(spiral_address)[::-1] \n",
    "    \n",
    "    hex_address = np.array([0,0])\n",
    "    \n",
    "    for i in range(len(digits)):\n",
    "        if int(digits[i])<0 or int(digits[i])>6:\n",
    "            print(\"Invalid spiral address!\")\n",
    "            return \n",
    "        elif digits[i]!= '0':\n",
    "            hex_address += sa2hex_aux(int(digits[i]),i)\n",
    "    return hex_address\n",
    "        \n",
    "# This computes the row/column positions of the base cases,\n",
    "# that is, in the form a*10^(zeros).\n",
    "def sa2hex_aux(a, zeros):\n",
    "    # Base cases\n",
    "    if zeros == 0:\n",
    "        if a == 0:\n",
    "            return np.array([0,0])\n",
    "        elif a == 1:\n",
    "            return np.array([0,8])\n",
    "        elif a == 2:\n",
    "            return np.array([-7,4])\n",
    "        elif a == 3:\n",
    "            return np.array([-7,-4])\n",
    "        elif a == 4:\n",
    "            return np.array([0,-8])\n",
    "        elif a == 5:\n",
    "            return np.array([7,-4])\n",
    "        elif a == 6:\n",
    "            return np.array([7,4])\n",
    "    \n",
    "    return sa2hex_aux(a,zeros-1)+ 2*sa2hex_aux(a%6 +1,zeros-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we compute the value of the hyperpel corresponding to the spiral address, by averaging the values on the subpixels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def sa_value(oversampled_image,spiral_address):\n",
    "    \"\"\"\n",
    "    Computes the value of the hyperpel corresponding to the given\n",
    "    spiral coordinate.\n",
    "    \"\"\"\n",
    "    hp = hyperpel_sa + sa2hex(spiral_address)\n",
    "    val = 0.\n",
    "    for i in range(56):\n",
    "        val += oversampled_image[hp[i,0],hp[i,1]]\n",
    "    \n",
    "    return val/56"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Spiral addition and multiplication"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def spiral_add(a,b,mod=0):\n",
    "    addition_table = [\n",
    "    [0,1,2,3,4,5,6],\n",
    "    [1,63,15,2,0,6,64],\n",
    "    [2,15,14,26,3,0,1],\n",
    "    [3,2,26,25,31,4,0],\n",
    "    [4,0,3,31,36,42,5],\n",
    "    [5,6,0,4,42,41,53],\n",
    "    [6,64,1,0,5,53,52]\n",
    "    ]\n",
    "    \n",
    "    dig_a = int2intvec(a)\n",
    "    dig_b = int2intvec(b) \n",
    "    \n",
    "    if (dig_a<0).any() or (dig_a>7).any() \\\n",
    "      or (dig_b<0).any() or (dig_b>7).any():\n",
    "        print(\"Invalid spiral address!\")\n",
    "        return\n",
    "    \n",
    "    if len(dig_a) == 1 and len(dig_b)==1:\n",
    "        return addition_table[a][b]\n",
    "    \n",
    "    if len(dig_a) < len(dig_b):\n",
    "        dig_a.resize(len(dig_b))\n",
    "    elif len(dig_b) < len(dig_a):\n",
    "        dig_b.resize(len(dig_a))\n",
    "        \n",
    "    res = 0\n",
    "    \n",
    "    for i in range(len(dig_a)):\n",
    "        \n",
    "        if i == len(dig_a)-1:\n",
    "            res += spiral_add(dig_a[i],dig_b[i])*(10**i)\n",
    "        else:\n",
    "            temp = spiral_add(dig_a[i],dig_b[i])\n",
    "            res += (temp%10)*(10**i)\n",
    "        \n",
    "            carry_on = spiral_add(dig_a[i+1],(temp - temp%10)/10)\n",
    "            dig_a[i+1] = str(carry_on)\n",
    "    \n",
    "    if mod!=0:\n",
    "        return res%mod\n",
    "    \n",
    "    return res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def spiral_mult(a,b, mod=0):\n",
    "    multiplication_table = [\n",
    "    [0,0,0,0,0,0,0],\n",
    "    [0,1,2,3,4,5,6],\n",
    "    [0,2,3,4,5,6,1],\n",
    "    [0,3,4,5,6,1,2],\n",
    "    [0,4,5,6,1,2,3],\n",
    "    [0,5,6,1,2,3,4],\n",
    "    [0,6,1,2,3,4,5],\n",
    "    ]\n",
    "    \n",
    "    dig_a = int2intvec(a)\n",
    "    dig_b = int2intvec(b) \n",
    "    \n",
    "    if (dig_a<0).any() or (dig_a>7).any() \\\n",
    "      or (dig_b<0).any() or (dig_b>7).any():\n",
    "        print(\"Invalid spiral address!\")\n",
    "        return\n",
    "    \n",
    "    sa_mult = int(0)\n",
    "    \n",
    "    for i in range(len(dig_b)):\n",
    "        for j in range(len(dig_a)):\n",
    "            temp = multiplication_table[dig_a[j]][dig_b[i]]*(10**(i+j))\n",
    "            sa_mult=spiral_add(sa_mult,temp)\n",
    "    \n",
    "    if mod!=0:\n",
    "        return sa_mult%mod\n",
    "    \n",
    "    return sa_mult"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Computation of the bispectrum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start by computing the vector $\\omega_f(\\lambda)$, where $\\lambda$ is a certain spiral address."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def omegaf(fft_oversampled, sa):\n",
    "    \"\"\"\n",
    "    Evaluates the vector omegaf corresponding to the given \n",
    "    spiral address sa.\n",
    "    \n",
    "    `fft_oversampled`: the oversampled FFT of the image\n",
    "    `sa`: the spiral address where to compute the vector\n",
    "    \"\"\"\n",
    "    \n",
    "    omegaf = np.zeros(6, dtype=fft_oversampled.dtype)\n",
    "    \n",
    "    for i in range(1,7):\n",
    "        omegaf[i-1] = sa_value(fft_oversampled,spiral_mult(sa,i))\n",
    "    \n",
    "    return omegaf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we can compute the \"generalized invariant\" corresponding to $\\lambda_1$, $\\lambda_2$ and $\\lambda_3$, starting from the FFT of the image.\n",
    "That is\n",
    "\n",
    "$$\n",
    "    I^3_f(\\lambda_1,\\lambda_2,\\lambda_3) = \\langle\\omega_f(\\lambda_1)\\odot\\omega_f(\\lambda_2),\\omega_f(\\lambda_3)\\rangle.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def invariant(fft_oversampled, sa1,sa2,sa3):\n",
    "    \"\"\"\n",
    "    Evaluates the generalized invariant of f on sa1, sa2 and sa3\n",
    "    \n",
    "    `fft_oversampled`: the oversampled FFT of the image\n",
    "    `sa1`, `sa2`, `sa3`: the spiral addresses where to compute the invariant\n",
    "    \"\"\"\n",
    "    \n",
    "    omega1 = omegaf(fft_oversampled,sa1)\n",
    "    omega2 = omegaf(fft_oversampled,sa2)\n",
    "    omega3 = omegaf(fft_oversampled,sa3)\n",
    "    \n",
    "    # Attention: np.vdot uses the scalar product with the complex \n",
    "    # conjugation at the first place!\n",
    "    return np.vdot(omega1*omega2,omega3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Finally, this function computes the bispectrum (or the rotational bispectrum) corresponding to the spiral addresses in the following picture.\n",
    "<img src=\"./pixels.png\" alt=\"Hexagonal pixels\" style=\"width: 200px;\"/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def bispectral_inv(fft_oversampled_example, rotational = False):\n",
    "    \"\"\"\n",
    "    Computes the (rotational) bispectral invariants for any sa1 \n",
    "    and any sa2 in the above picture.\n",
    "    \n",
    "    `fft_oversampled_example`: oversampled FFT of the image\n",
    "    `rotational`: if True, we compute the rotational bispectrum\n",
    "    \"\"\"\n",
    "    \n",
    "    if rotational == True:\n",
    "        bispectrum = np.zeros(9**2*6,dtype = fft_oversampled_example.dtype)\n",
    "    else:\n",
    "        bispectrum = np.zeros(9**2,dtype = fft_oversampled_example.dtype)\n",
    "        \n",
    "    indexes = [0,1,10,11,12,13,14,15,16]\n",
    "    \n",
    "    count = 0\n",
    "    for i in range(9):\n",
    "        sa1 = indexes[i]\n",
    "        sa1_base10 = base7to10(sa1)\n",
    "        for k in range(9):\n",
    "            sa2 = indexes[k]\n",
    "            if rotational == True:\n",
    "                for r in range(6):\n",
    "                    sa2_rot = spiral_mult(sa2,r)\n",
    "                    sa2_rot_base10 = base7to10(sa2_rot)\n",
    "                    sa3 = base10to7(sa1_base10+sa2_rot_base10)\n",
    "                    bispectrum[count]=invariant(fft_oversampled_example,sa1,sa2,sa3)\n",
    "                    count += 1\n",
    "            else:\n",
    "                sa2_base10 = base7to10(sa2)\n",
    "                sa3 = base10to7(sa1_base10+sa2_base10)\n",
    "                bispectrum[count]=invariant(fft_oversampled_example,sa1,sa2,sa3)\n",
    "                count += 1\n",
    "    \n",
    "    return bispectrum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some timing tests."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "example =  1 - rgb2gray(plt.imread('./test-images/butterfly.png'))\n",
    "fft_example = np.fft.fftshift(np.fft.fft2(example))\n",
    "fft_oversampled_example = oversampling(fft_example)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loops, best of 3: 372 ms per loop\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "bispectral_inv(fft_oversampled_example)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loops, best of 3: 2.35 s per loop\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "bispectral_inv(fft_oversampled_example, rotational=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tests"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we define various functions to batch test the images in the test folder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "folder = './test-images'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def evaluate_invariants(image, rot = False):\n",
    "    \"\"\"\n",
    "    Evaluates the invariants of the given image.\n",
    "    \n",
    "    `image`: the matrix representing the image (not oversampled)\n",
    "    `rot`: if True we compute the rotational bispectrum\n",
    "    \"\"\"\n",
    "    \n",
    "    # compute the normalized FFT\n",
    "    fft = np.fft.fftshift(np.fft.fft2(image))\n",
    "    fft /= fft / LA.norm(fft)\n",
    "    \n",
    "    # oversample it\n",
    "    fft_oversampled = oversampling(fft)\n",
    "    \n",
    "    return bispectral_inv(fft_oversampled, rotational = rot)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some timing tests."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loops, best of 3: 1.07 s per loop\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "evaluate_invariants(example)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loops, best of 3: 3.09 s per loop\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "evaluate_invariants(example, rot = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def bispectral_folder(folder_name = folder, rot = False): \n",
    "    \"\"\"\n",
    "    Evaluates all the invariants of the images in the selected folder, \n",
    "    storing them in a dictionary with their names as keys.\n",
    "    \n",
    "    `folder_name`: path to the folder\n",
    "    `rot`: if True we compute the rotational bispectrum\n",
    "    \"\"\"\n",
    "    \n",
    "    # we store the results in a dictionary\n",
    "    results = {}\n",
    "    \n",
    "    for filename in os.listdir(folder_name):\n",
    "        infilename = os.path.join(folder_name, filename)\n",
    "        if not os.path.isfile(infilename): \n",
    "            continue\n",
    "\n",
    "        base, extension = os.path.splitext(infilename)\n",
    "        if extension == '.png':\n",
    "            test_img = 1 - rgb2gray(plt.imread(infilename))\n",
    "            bispectrum = evaluate_invariants(test_img, rot = rot)\n",
    "            \n",
    "            results[os.path.splitext(filename)[0]] = bispectrum\n",
    "            \n",
    "    return results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def bispectral_comparison(bispectrums, comparison = 'triangle', plot = True, log_scale = True):\n",
    "    \"\"\"\n",
    "    Returns the difference of the norms of the given invariants w.r.t. the \n",
    "    comparison element.\n",
    "    \n",
    "    `bispectrums`: a dictionary with as keys the names of the images and \n",
    "                    as values their invariants\n",
    "    `comparison`:  the element to use as comparison\n",
    "    \"\"\"\n",
    "    \n",
    "    if comparison not in bispectrums:\n",
    "        print(\"The requested comparison is not in the folder\")    \n",
    "        return\n",
    "    \n",
    "    \n",
    "    bispectrum_diff = {}\n",
    "    for elem in bispectrums:\n",
    "        diff = LA.norm(bispectrums[elem]-bispectrums[comparison])\n",
    "        # we remove nan results\n",
    "        if not np.isnan(diff):\n",
    "            bispectrum_diff[elem] = diff\n",
    "        \n",
    "    return bispectrum_diff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def bispectral_plot(bispectrums, comparison = 'triangle', log_scale = True):\n",
    "    \"\"\"\n",
    "    Plots the difference of the norms of the given invariants w.r.t. the \n",
    "    comparison element (by default in logarithmic scale).\n",
    "    \n",
    "    `bispectrums`: a dictionary with as keys the names of the images and \n",
    "                    as values their invariants\n",
    "    `comparison`:  the element to use as comparison\n",
    "    `log_scale`:   wheter the plot should be in log_scale\n",
    "    \"\"\"\n",
    "    \n",
    "    bispectrum_diff = bispectral_comparison(bispectrums, comparison = comparison)\n",
    "\n",
    "    plt.plot(bispectrum_diff.values(),'ro')\n",
    "    if log_scale == True:\n",
    "        plt.yscale('log')\n",
    "    for i in range(len(bispectrum_diff.values())):\n",
    "        # if we plot in log scale, we do not put labels on items that are\n",
    "        # too small, otherwise they exit the plot area.\n",
    "        if log_scale and bispectrum_diff.values()[i] < 10**(-3):\n",
    "            continue\n",
    "        plt.text(i,bispectrum_diff.values()[i],bispectrum_diff.keys()[i][:3])\n",
    "        plt.title(\"Comparison with as reference '\"+ comparison +\"'\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Construction of the table for the paper"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "comparisons_paper = ['triangle', 'rectangle', 'ellipse', 'etoile', 'diamond']\n",
    "\n",
    "def extract_table_values(bispectrums, comparisons = comparisons_paper):\n",
    "    \"\"\"\n",
    "    Extract the values for the table of the paper.\n",
    "    \n",
    "    `bispectrums`: a dictionary with as keys the names of the images and \n",
    "                    as values their invariants\n",
    "    `comparison`:  list of elements to use as comparison\n",
    "    \n",
    "    Returns a list of tuples. Each tuple contains the name of the comparison \n",
    "    element, the maximal value of the difference of the norm of the invariants \n",
    "    with its rotated and the minimal values of the same difference with the \n",
    "    other images.\n",
    "    \"\"\"\n",
    "    table_values = []\n",
    "    for elem in comparisons:\n",
    "        diff = bispectral_comparison(bispectrums, comparison= elem, plot=False)\n",
    "\n",
    "        l = len(elem)\n",
    "        match = [x for x in diff.keys() if x[:l]==elem]\n",
    "        not_match = [x for x in diff.keys() if x[:l]!=elem]\n",
    "\n",
    "        max_match = max([ diff[k] for k in match ])\n",
    "        min_not_match = min([ diff[k] for k in not_match ])\n",
    "        \n",
    "        table_values.append((elem,'%.2E' % (max_match),'%.2E' % min_not_match))\n",
    "        \n",
    "    return table_values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python2.7/site-packages/IPython/kernel/__main__.py:6: RuntimeWarning: invalid value encountered in divide\n"
     ]
    }
   ],
   "source": [
    "bispectrums = bispectral_folder()\n",
    "bispectrums_rotational = bispectral_folder(rot=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('triangle', '9.23E+10', '7.03E+12'),\n",
       " ('rectangle', '7.93E+10', '8.22E+12'),\n",
       " ('ellipse', '7.42E+10', '7.11E+12'),\n",
       " ('etoile', '7.27E+10', '5.54E+12'),\n",
       " ('diamond', '3.78E+10', '5.47E+12')]"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "extract_table_values(bispectrums)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('triangle', '2.26E+11', '1.72E+13'),\n",
       " ('rectangle', '1.94E+11', '2.01E+13'),\n",
       " ('ellipse', '1.82E+11', '1.74E+13'),\n",
       " ('etoile', '1.78E+11', '1.36E+13'),\n",
       " ('diamond', '9.27E+10', '1.34E+13')]"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "extract_table_values(bispectrums_rotational)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}