brisque.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "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.6"
    },
    "pycharm": {
      "stem_cell": {
        "cell_type": "raw",
        "metadata": {
          "collapsed": false
        },
        "source": []
      }
    },
    "colab": {
      "name": "Cópia de brisque.ipynb",
      "provenance": [],
      "collapsed_sections": [
        "JyT9GuV4lGTx",
        "C1s1dnFflGUR"
      ]
    }
  },
  "cells": [
    {
      "cell_type": "code",
      "metadata": {
        "pycharm": {
          "is_executing": false
        },
        "id": "P3VyW-0ylGTU"
      },
      "source": [
        "!pip install opencv-python>=3.4.2.17\n",
        "!pip install libsvm>=3.23.0"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "pycharm": {
          "is_executing": false
        },
        "id": "7vpEahFjlGTa"
      },
      "source": [
        "import collections\n",
        "from itertools import chain\n",
        "import urllib.request as request\n",
        "import pickle \n",
        "\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "\n",
        "import scipy.signal as signal\n",
        "import scipy.special as special\n",
        "import scipy.optimize as optimize\n",
        "\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "import skimage.io\n",
        "import skimage.transform\n",
        "\n",
        "import cv2\n",
        "import glob\n",
        "import os\n",
        "\n",
        "\n",
        "from libsvm import svmutil\n",
        "import pathlib\n",
        "import IPython.display as display\n",
        "from PIL import Image\n",
        "\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "49YurxrZlGTd"
      },
      "source": [
        "# Image Quality Assessment\n",
        "https://www.sciencedirect.com/science/article/pii/S0923596514001490\n",
        "\n",
        "Image quality is a notion that highly depends on observers. Generally, \n",
        "it is linked to the conditions in which it is viewed; therefore, a highly subjective topic. Image quality assessment aims to quantitatively represent the human perception of quality. These metrics commonly are used to analyze the performance of algorithms in different fields of computer vision like image compression, image transmission, and image processing [1].\n",
        "\n",
        "Image quality assessment (IQA) is mainly divided into two areas of research (1) reference-based evaluation and (2)  no-reference evaluation. The main difference is that reference-based methods depend on a high-quality image as a source to evaluate the difference between images. An example of reference-based evaluations is the Structural Similarity Index (SSIM) [2].\n",
        "\n",
        "## No-reference Image Quality Assessment\n",
        "\n",
        "No-reference image quality assessment does not require a base image to evaluate an image quality, the only information that the algorithm receives is the distorted image whose quality is being assessed.\n",
        "\n",
        "Blind methods are mostly comprised of two steps. The first step calculates features that describe the image's structure and the second step relates the features with the human opinion of the image quality. TID2008 is a famous database created following a methodology that describes how to measure human opinion scores from referenced images [3], it is widely used to compare the performance of IQA algorithms.\n",
        "\n",
        "## Blind/referenceless image spatial quality evaluator (BRISQUE)\n",
        "\n",
        "BRISQUE [4] is a model that only used the image pixels to calculate features (other methods are based on image transformation to other spaces like wavelet or DCT). It is demonstrated to be highly efficient as it does not need any transformation to calculate its features.\n",
        "\n",
        "It relies on spatial Natural Scene Statistics (NSS) model of locally normalized luminance coefficients in the spatial domain, as well as the model for pairwise products of these coefficients. \n",
        "\n",
        "## Methodology\n",
        "### Natural Scene Statistics in the Spatial Domain\n",
        "Given an image $I(i, j)$, first, compute the locally normalized luminances $\\hat{I}(i,j)$ via local mean subtraction $\\mu(i,j)$ and divide it by the local deviation $\\sigma(i, j)$. $C$ is added to avoid zero divisions. \n",
        "\n",
        "$$\\hat{I}(i,j) = \\frac{I(i,j) - \\mu(i,j)}{\\sigma(i,j) + C}$$\n",
        "\n",
        "*Hint: If $I(i,j)$'s domain is [0,255] then $C=1$ if the domain is [0,1] then $C=1/255$.* \n",
        "\n",
        "To calculate the locally normalized luminance, also known as mean substracted contrast normalized (MSCN) coefficients, first, we need to calculate the local mean\n",
        "\n",
        "$$\\mu(i,j) = \\sum_{k=-K}^{K}\\sum_{l=-L}^{L}w_{k,l}I_{k,l}(i,j)$$\n",
        "\n",
        "where $w$ is a Gaussian kernel of size (K, L).\n",
        "\n",
        "The way that the author displays the local mean could be a little bit confusing but it is just applying a Gaussian filter to the image."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "pycharm": {
          "is_executing": false
        },
        "id": "wveqfVj2lGTd"
      },
      "source": [
        "def normalize_kernel(kernel):\n",
        "    return kernel / np.sum(kernel)\n",
        "\n",
        "def gaussian_kernel2d(n, sigma):\n",
        "    Y, X = np.indices((n, n)) - int(n/2)\n",
        "    gaussian_kernel = 1 / (2 * np.pi * sigma ** 2) * np.exp(-(X ** 2 + Y ** 2) / (2 * sigma ** 2)) \n",
        "    return normalize_kernel(gaussian_kernel)\n",
        "\n",
        "def local_mean(image, kernel):\n",
        "    return signal.convolve2d(image, kernel, 'same')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "x4bCzXholGTg"
      },
      "source": [
        "Then, we calculate the local deviation\n",
        "\n",
        "$$ \\sigma(i,j) = \\sqrt{\\sum_{k=-K}^{K}\\sum_{l=-L}^{L}w_{k,l}(I_{k,l}(i, j) - \\mu(i, j))^2 } $$"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "pycharm": {
          "is_executing": false
        },
        "id": "vIK80XPvlGTg"
      },
      "source": [
        "def local_deviation(image, local_mean, kernel):\n",
        "    \"Vectorized approximation of local deviation\"\n",
        "    sigma = image ** 2\n",
        "    sigma = signal.convolve2d(sigma, kernel, 'same')\n",
        "    return np.sqrt(np.abs(local_mean ** 2 - sigma))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QldI_4RSlGTk"
      },
      "source": [
        "Finally, we calculate the MSCN coefficients\n",
        "\n",
        "$$\\hat{I}(i,j) = \\frac{I(i,j) - \\mu(i,j)}{\\sigma(i,j) + C}$$"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "pycharm": {
          "is_executing": false
        },
        "id": "BRgVOrqYlGTl"
      },
      "source": [
        "def calculate_mscn_coefficients(image, kernel_size=6, sigma=7/6):\n",
        "    C = 1/255\n",
        "    kernel = gaussian_kernel2d(kernel_size, sigma=sigma)\n",
        "    local_mean = signal.convolve2d(image, kernel, 'same')\n",
        "    local_var = local_deviation(image, local_mean, kernel)\n",
        "    \n",
        "    return (image - local_mean) / (local_var + C)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "uM_xiGoLlGTo"
      },
      "source": [
        "The author found that the MSCN coefficients are distributed as a Generalized Gaussian Distribution (GGD) for a broader spectrum of distorted image.\n",
        "\n",
        "$$f(x; \\alpha, \\sigma^2) = \\frac{\\alpha}{2\\beta\\Gamma(1/\\alpha)}e^{-\\big(\\frac{|x|}{\\beta}\\big)^\\alpha}$$\n",
        "\n",
        "where\n",
        "\n",
        "$$\\beta = \\sigma \\sqrt{\\frac{\\Gamma\\big(\\frac{1}{\\alpha}\\big)}{\\Gamma\\big(\\frac{3}{\\alpha}\\big)}}$$\n",
        "\n",
        "and $\\Gamma$ is the gamma function.\n",
        "\n",
        "The shape $\\alpha$ controls the shape and $\\sigma^2$ th variance."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "pycharm": {
          "is_executing": false
        },
        "id": "PQLXfUqclGTo"
      },
      "source": [
        "def generalized_gaussian_dist(x, alpha, sigma):\n",
        "    beta = sigma * np.sqrt(special.gamma(1 / alpha) / special.gamma(3 / alpha))\n",
        "    \n",
        "    coefficient = alpha / (2 * beta() * special.gamma(1 / alpha))\n",
        "    return coefficient * np.exp(-(np.abs(x) / beta) ** alpha)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "nrsyo-y-lGTs"
      },
      "source": [
        "*Pairwise* products of neighboring MSCN coefficients\n",
        "\n",
        "The signs of adjacent coefficients also exhibit a regular structure, which gets disturbed in the presence of distortion. The author proposes the model of pairwise products of neighboring MSCN coefficients along four directions (1) horizontal $H$, (2) vertical $V$, (3) main-diagonal $D1$ and (4) secondary-diagonal $D2$.\n",
        "\n",
        "$$H(i,j) = \\hat{I}(i,j) \\hat{I}(i, j + 1)$$\n",
        "$$V(i,j) = \\hat{I}(i,j) \\hat{I}(i + 1, j)$$\n",
        "$$D1(i,j) = \\hat{I}(i,j) \\hat{I}(i + 1, j + 1)$$\n",
        "$$D2(i,j) = \\hat{I}(i,j) \\hat{I}(i + 1, j - 1)$$"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "pycharm": {
          "is_executing": false
        },
        "id": "5qBXkox8lGTt"
      },
      "source": [
        "def calculate_pair_product_coefficients(mscn_coefficients):\n",
        "    return collections.OrderedDict({\n",
        "        'mscn': mscn_coefficients,\n",
        "        'horizontal': mscn_coefficients[:, :-1] * mscn_coefficients[:, 1:],\n",
        "        'vertical': mscn_coefficients[:-1, :] * mscn_coefficients[1:, :],\n",
        "        'main_diagonal': mscn_coefficients[:-1, :-1] * mscn_coefficients[1:, 1:],\n",
        "        'secondary_diagonal': mscn_coefficients[1:, :-1] * mscn_coefficients[:-1, 1:]\n",
        "    })"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pu8xtiHUlGTv"
      },
      "source": [
        "The author mentions that the Generalized Gaussian Distribution does not provide good fit to the empirical histograms of coefficient producs. Thus, they propose the Asymmetric Generalized Gaussian Distribution (AGGD) model [5].\n",
        "\n",
        "$$\n",
        "f(x; \\nu, \\sigma_l^2, \\sigma_r^2) =  \n",
        "   \\begin{cases} \n",
        "      \\frac{\\nu}{(\\beta_l + \\beta_r)\\Gamma\\big(\\frac{1}{\\nu}\\big)}e^{\\big(-\\big(\\frac{-x}{\\beta_l}\\big)^\\nu\\big)} & x < 0 \\\\\n",
        "        \\frac{\\nu}{(\\beta_l + \\beta_r)\\Gamma\\big(\\frac{1}{\\nu}\\big)}e^{\\big(-\\big(\\frac{x}{\\beta_r}\\big)^\\nu\\big)} & x >= 0\n",
        "\\end{cases}\n",
        "$$\n",
        "\n",
        "where\n",
        "\n",
        "$$\\beta_{side} = \\sigma_{side} \\sqrt{\\frac{\\Gamma\\big(\\frac{1}{\\nu}\\big)}{\\Gamma\\big(\\frac{3}{\\nu}\\big)}}$$\n",
        "\n",
        "and $side$ can be either $r$ or $l$.\n",
        "\n",
        "Another parameter that is not reflected in the previous formula is the mean\n",
        "\n",
        "$$\\eta = (\\beta_r - beta_l) \\frac{\\Gamma\\big(\\frac{2}{\\nu}\\big)}{\\Gamma\\big(\\frac{1}{\\nu}\\big)}$$"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "pycharm": {
          "is_executing": false
        },
        "id": "96oyUY_-lGTv"
      },
      "source": [
        "def asymmetric_generalized_gaussian(x, nu, sigma_l, sigma_r):\n",
        "    def beta(sigma):\n",
        "        return sigma * np.sqrt(special.gamma(1 / nu) / special.gamma(3 / nu))\n",
        "    \n",
        "    coefficient = nu / ((beta(sigma_l) + beta(sigma_r)) * special.gamma(1 / nu))\n",
        "    f = lambda x, sigma: coefficient * np.exp(-(x / beta(sigma)) ** nu)\n",
        "        \n",
        "    return np.where(x < 0, f(-x, sigma_l), f(x, sigma_r))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "JyT9GuV4lGTx"
      },
      "source": [
        "### Fitting Asymmetric Generalized Gaussian Distribution\n",
        "\n",
        "The methodology to fit an Asymmetric Generalized Gaussian Distribution is described in [5].\n",
        "\n",
        "1. Calculate $\\hat{\\gamma}$ where $N_l$ is the number of negative samples and $N_r$ is the number of positive samples. \n",
        "\n",
        "$$\n",
        "\\hat{\\gamma} = \\frac{\\sqrt{\\frac{1}{N_l - 1}\\sum_{k=1, x_k < 0}^{N_l} x_k^2}\n",
        "}{\\sqrt{\\frac{1}{N_r - 1}\\sum_{k=1, x_k >= 0}^{N_r} x_k^2}\n",
        "}\n",
        "$$\n",
        "\n",
        "2. Calculate $\\hat{r}$.\n",
        "\n",
        "$$\\hat{r} = \\frac{\\big(\\frac{\\sum|x_k|}{N_l + N_r}\\big)^2}{\\frac{\\sum{x_k ^ 2}}{N_l + N_r}} $$\n",
        "\n",
        "3. Calculate $\\hat{R}$ using $\\hat{\\gamma}$ and $\\hat{r}$ estimations.\n",
        "\n",
        "$$\\hat{R} = \\hat{r} \\frac{(\\hat{\\gamma}^3 + 1)(\\hat{\\gamma} + 1)}{(\\hat{\\gamma}^2 + 1)^2}$$\n",
        "\n",
        "4. Estimate $\\alpha$ using the approximation of the inverse generalized Gaussian ratio.\n",
        "\n",
        "$$\\hat{\\alpha} = \\hat{\\rho} ^ {-1}(\\hat{R})$$\n",
        "\n",
        "$$\\rho(\\alpha) = \\frac{\\Gamma(2 / \\alpha) ^ 2}{\\Gamma(1 / \\alpha) \\Gamma(3 / \\alpha)}$$\n",
        "\n",
        "5. Estimate left and right scale parameters.\n",
        "$$\\sigma_l = \\sqrt{\\frac{1}{N_l - 1}\\sum_{k=1, x_k < 0}^{N_l} x_k^2}$$\n",
        "$$\\sigma_r = \\sqrt{\\frac{1}{N_r - 1}\\sum_{k=1, x_k >= 0}^{N_r} x_k^2}$$"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "pycharm": {
          "is_executing": false
        },
        "id": "4QUT1BiglGTy"
      },
      "source": [
        "def asymmetric_generalized_gaussian_fit(x):\n",
        "    def estimate_phi(alpha):\n",
        "        numerator = special.gamma(2 / alpha) ** 2\n",
        "        denominator = special.gamma(1 / alpha) * special.gamma(3 / alpha)\n",
        "        return numerator / denominator\n",
        "\n",
        "    def estimate_r_hat(x):\n",
        "        size = np.prod(x.shape)\n",
        "        return (np.sum(np.abs(x)) / size) ** 2 / (np.sum(x ** 2) / size)\n",
        "\n",
        "    def estimate_R_hat(r_hat, gamma):\n",
        "        numerator = (gamma ** 3 + 1) * (gamma + 1)\n",
        "        denominator = (gamma ** 2 + 1) ** 2\n",
        "        return r_hat * numerator / denominator\n",
        "\n",
        "    def mean_squares_sum(x, filter = lambda z: z == z):\n",
        "        filtered_values = x[filter(x)]\n",
        "        squares_sum = np.sum(filtered_values ** 2)\n",
        "        return squares_sum / ((filtered_values.shape))\n",
        "\n",
        "    def estimate_gamma(x):\n",
        "        left_squares = mean_squares_sum(x, lambda z: z < 0)\n",
        "        right_squares = mean_squares_sum(x, lambda z: z >= 0)\n",
        "\n",
        "        return np.sqrt(left_squares) / np.sqrt(right_squares)\n",
        "\n",
        "    def estimate_alpha(x):\n",
        "        r_hat = estimate_r_hat(x)\n",
        "        gamma = estimate_gamma(x)\n",
        "        R_hat = estimate_R_hat(r_hat, gamma)\n",
        "\n",
        "        solution = optimize.root(lambda z: estimate_phi(z) - R_hat, [0.2]).x\n",
        "\n",
        "        return solution[0]\n",
        "\n",
        "    def estimate_sigma(x, alpha, filter = lambda z: z < 0):\n",
        "        return np.sqrt(mean_squares_sum(x, filter))\n",
        "    \n",
        "    def estimate_mean(alpha, sigma_l, sigma_r):\n",
        "        return (sigma_r - sigma_l) * constant * (special.gamma(2 / alpha) / special.gamma(1 / alpha))\n",
        "    \n",
        "    alpha = estimate_alpha(x)\n",
        "    sigma_l = estimate_sigma(x, alpha, lambda z: z < 0)\n",
        "    sigma_r = estimate_sigma(x, alpha, lambda z: z >= 0)\n",
        "    \n",
        "    constant = np.sqrt(special.gamma(1 / alpha) / special.gamma(3 / alpha))\n",
        "    mean = estimate_mean(alpha, sigma_l, sigma_r)\n",
        "    \n",
        "    return alpha, mean, sigma_l, sigma_r"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8rI7OenblGT0"
      },
      "source": [
        "### Calculate BRISQUE features\n",
        "\n",
        "The features needed to calculate the image quality are the result of fitting the MSCN coefficients and shifted products to the Generalized Gaussian Distributions. First, we need to fit the MSCN coefficients to the GDD, then the pairwise products to the AGGD. A summary of the features is the following:\n",
        "\n",
        "| Feature ID      | Feature Description                            | Computation Procedure            |\n",
        "|-----------------|------------------------------------------------|----------------------------------|\n",
        "| $f_1-f_2$       | Shape and variance                             | Fit GGD to MSCN coefficients     |\n",
        "| $f_3-f_6$       | Shape, mean, left variance, right variance     | Fit AGGD to H pairwise products  |\n",
        "| $f_7-f_{10}$    | Shape, mean, left variance, right variance     | Fit AGGD to V pairwise products  |\n",
        "| $f_{11}-f_{14}$ | Shape, mean, left variance, right variance     | Fit AGGD to D1 pairwise products |\n",
        "| $f_{15}-f_{18}$ | Shape, mean, left variance, right variance     | Fit AGGD to D2 pairwise products |"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "pycharm": {
          "is_executing": false
        },
        "id": "bUCTOXx6lGT1"
      },
      "source": [
        "def calculate_brisque_features(image, kernel_size=7, sigma=7/6):\n",
        "    def calculate_features(coefficients_name, coefficients, accum=np.array([])):\n",
        "        alpha, mean, sigma_l, sigma_r = asymmetric_generalized_gaussian_fit(coefficients)\n",
        "\n",
        "        if coefficients_name == 'mscn':\n",
        "            var = (sigma_l ** 2 + sigma_r ** 2) / 2\n",
        "            return [alpha, var]\n",
        "        \n",
        "        return [alpha, mean, sigma_l ** 2, sigma_r ** 2]\n",
        "    \n",
        "    mscn_coefficients = calculate_mscn_coefficients(image, kernel_size, sigma)\n",
        "    coefficients = calculate_pair_product_coefficients(mscn_coefficients)\n",
        "    \n",
        "    features = [calculate_features(name, coeff) for name, coeff in coefficients.items()]\n",
        "    flatten_features = list(chain.from_iterable(features))\n",
        "    return np.array(flatten_features)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tcTKWbK9lGT3"
      },
      "source": [
        "# Hands-on\n",
        "\n",
        "After creating all the functions needed to calculate the brisque features, we can estimate the image quality for a given image. In [4], they use an image that comes from the Kodak dataset [6]."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "NqrxLHKBlGT4"
      },
      "source": [
        "## Auxiliary Functions"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "pycharm": {
          "is_executing": false
        },
        "id": "AMt89i2ilGT4"
      },
      "source": [
        "def load_image(url):\n",
        "    image_stream = request.urlopen(url)\n",
        "    return skimage.io.imread(image_stream, plugin='pil')\n",
        "\n",
        "def plot_histogram(x, label):\n",
        "    n, bins = np.histogram(x.ravel(), bins=50)\n",
        "    n = n / np.max(n)\n",
        "    plt.plot(bins[:-1], n, label=label, marker='o')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "xLElglG2op_o",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "86dfd678-5b9c-4068-8b57-1a88817c123a"
      },
      "source": [
        "#IMPORT DATABASE FROM DRIVE\n",
        "from google.colab import drive\n",
        "drive.mount('/content/drive')\n",
        "\n",
        "#UNZIPE FILES\n",
        "def unzip_file(file_name, local_path):\n",
        "  from zipfile import ZipFile\n",
        "  zip_ref = ZipFile(file_name, 'r')\n",
        "  zip_ref.extractall(local_path)\n",
        "  zip_ref.close()\n",
        "  \n",
        "file_name = '/content/drive/My Drive/TREINA_DFT/data_base.zip'\n",
        "local_path='.'\n",
        "\n",
        "unzip_file(file_name, local_path)\n",
        "\n"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Go to this URL in a browser: https://accounts.google.com/o/oauth2/auth?client_id=947318989803-6bn6qk8qdgf4n4g3pfee6491hc0brc4i.apps.googleusercontent.com&redirect_uri=urn%3aietf%3awg%3aoauth%3a2.0%3aoob&scope=email%20https%3a%2f%2fwww.googleapis.com%2fauth%2fdocs.test%20https%3a%2f%2fwww.googleapis.com%2fauth%2fdrive%20https%3a%2f%2fwww.googleapis.com%2fauth%2fdrive.photos.readonly%20https%3a%2f%2fwww.googleapis.com%2fauth%2fpeopleapi.readonly%20https%3a%2f%2fwww.googleapis.com%2fauth%2fdrive.activity.readonly%20https%3a%2f%2fwww.googleapis.com%2fauth%2fexperimentsandconfigs%20https%3a%2f%2fwww.googleapis.com%2fauth%2fphotos.native&response_type=code\n",
            "\n",
            "Enter your authorization code:\n",
            "4/1AX4XfWjmyEDwyAAIgAlfjijBhMpbIOFwPCSbLPij7seeN3nbUuS2VAk8uHg\n",
            "Mounted at /content/drive\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "T71e2efxUK-W"
      },
      "source": [
        "# EXTRACT IMAGE NAME"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Ws7FrSpWUJZD"
      },
      "source": [
        "IMG_DIR = '/content/data_base/train'\n",
        "\n",
        "def trataNome(filename):\n",
        "    filename = filename[1].split('.')\n",
        "    filename = filename[0]\n",
        "    filename = filename.replace('_1','')\n",
        "    return filename\n",
        "\n",
        "def extract_img_name(dir_entrada):\n",
        "    listaResultado = []\n",
        "    for dir in glob.glob(dir_entrada+'/*'):\n",
        "      #print(dir)\n",
        "      for categoria_dir in glob.glob(dir+'/*'):\n",
        "        head_tail = os.path.split(categoria_dir) \n",
        "        head_tail = trataNome(head_tail)\n",
        "        listaResultado.append(head_tail)\n",
        "    return listaResultado\n",
        "\n",
        "def extract_img_name_path(dir_entrada):\n",
        "    listaResultado = []\n",
        "    for dir in glob.glob(dir_entrada+'/*'):\n",
        "      #print(dir)\n",
        "      for categoria_dir in glob.glob(dir+'/*'):\n",
        "        head_tail = os.path.split(categoria_dir) \n",
        "        head_tail = trataNome(head_tail)\n",
        "        img_name_path = ([categoria_dir,head_tail])\n",
        "        listaResultado.append(img_name_path)\n",
        "        #listaResultado = pd.DataFrame(listaResultado, columns=['path','img_name']).transpose()\n",
        "    return listaResultado\n",
        "\n",
        "def numberOfNonNans(data):\n",
        "    count = 0\n",
        "    for i in data:\n",
        "        if not np.isnan(i):\n",
        "            count += 1\n",
        "    return count \n",
        "\n",
        "def extract_unique_img_name(image_path):\n",
        "  return trataNome(os.path.split(image_path)) \n",
        "\n",
        "\n",
        "#EXPORT DATA\n",
        "#lista  = pd.DataFrame(extract_img_name(IMG_DIR),columns=['train'])\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "gzysgYualGT8"
      },
      "source": [
        "## 1. Load image"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "pycharm": {
          "is_executing": false
        },
        "id": "FHWKCk4vlGT8"
      },
      "source": [
        "%matplotlib inline\n",
        "from skimage import img_as_float\n",
        "from skimage import io\n",
        "import pathlib\n",
        "plt.rcParams[\"figure.figsize\"] = 12, 9\n",
        "\n",
        "list_of_img=[]\n",
        "category_name = 'andador'\n",
        "v_prod = list(pathlib.Path(IMG_DIR).glob(category_name+'/*'))\n",
        "for image_path in v_prod[:6]:\n",
        "  img_name = extract_unique_img_name(image_path)\n",
        "  image =  skimage.io.imread(image_path, plugin='pil')\n",
        "  gray_image = skimage.color.rgb2gray(image)\n",
        "  img_with_name = [img_name,gray_image]\n",
        "  print(img_with_name)\n",
        "  list_of_img.append(img_with_name)\n",
        "  #display.display(Image.open(str(image_path)))\n",
        "  #skimage.io.imshow(image)\n",
        "list_of_img_to_df = pd.DataFrame(list_of_img, columns=['img_name','img_gray'])#.transpose()\n",
        "\n",
        "\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "IKUt3c9Uj5eA"
      },
      "source": [
        ""
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "hbDY3PjPkiyK",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 115
        },
        "outputId": "dfb6da8a-5ca3-47fa-b582-146cd2c73584"
      },
      "source": [
        "#list_of_img_to_df.head()\n",
        "list_of_img_to_df['img_gray'].head()"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "0    [[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,...\n",
              "1    [[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,...\n",
              "2    [[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,...\n",
              "3    [[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,...\n",
              "4    [[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,...\n",
              "Name: img_gray, dtype: object"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 66
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "V_nI4XhqlGT_"
      },
      "source": [
        "## 2. Calculate Coefficients\n",
        "\n",
        "```\n",
        "# Isto está formatado como código\n",
        "```\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "pycharm": {
          "is_executing": false
        },
        "id": "0uAVTXRXlGT_",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 50
        },
        "outputId": "17f7c627-8dcf-43be-b1b2-699fc0b6b5ce"
      },
      "source": [
        "%%time \n",
        "\n",
        "mscn_coefficients = calculate_mscn_coefficients(gray_image, 7, 7/6)\n",
        "coefficients = calculate_pair_product_coefficients(mscn_coefficients)\n",
        "#print(coefficients)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "CPU times: user 21.1 ms, sys: 0 ns, total: 21.1 ms\n",
            "Wall time: 22.2 ms\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ar6sSlR_lGUC"
      },
      "source": [
        "After calculating the MSCN coefficients and the pairwise products, we can \n",
        "\n",
        "---\n",
        "\n",
        "verify that the distributions are in fact different."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "pycharm": {
          "is_executing": false
        },
        "id": "BxsljqwilGUC"
      },
      "source": [
        "%matplotlib inline\n",
        "plt.rcParams[\"figure.figsize\"] = 12, 11\n",
        "\n",
        "for name, coeff in coefficients.items():\n",
        "    plot_histogram(coeff.ravel(), name)\n",
        "\n",
        "plt.axis([-2.5, 2.5, 0, 1.05])\n",
        "plt.legend()\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "e9j5ttaJlGUF"
      },
      "source": [
        "## 3. Fit Coefficients to Generalized Gaussian Distributions "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "pycharm": {
          "is_executing": false
        },
        "id": "5w9GdL5FlGUF"
      },
      "source": [
        "%%time \n",
        "\n",
        "brisque_features = calculate_brisque_features(gray_image, kernel_size=7, sigma=7/6)\n",
        "#print(brisque_features)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "HuhHWVImlGUI"
      },
      "source": [
        "## 4. Resize Image and Calculate BRISQUE Features"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "pycharm": {
          "is_executing": false
        },
        "id": "2_67VM21lGUJ"
      },
      "source": [
        "%%time\n",
        "\n",
        "downscaled_image = cv2.resize(gray_image, None, fx=1/2, fy=1/2, interpolation = cv2.INTER_CUBIC)\n",
        "downscale_brisque_features = calculate_brisque_features(downscaled_image, kernel_size=7, sigma=7/6)\n",
        "\n",
        "brisque_features = np.concatenate((brisque_features, downscale_brisque_features))\n",
        "#print(brisque_features)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RV1mKO4ilGUM"
      },
      "source": [
        "## 5. Scale Features and Feed the SVR\n",
        "The author provides a pretrained SVR model to calculate the quality assessment. However, in order to have good results, we need to scale the features to [-1, 1]. For the latter, we need the same parameters the author used to scale the features vector."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "pycharm": {
          "is_executing": false
        },
        "id": "ZRXWeMBelGUN"
      },
      "source": [
        "def scale_features(features):\n",
        "    with open('/content/sample_data/normalize.pickle', 'rb') as handle:\n",
        "        scale_params = pickle.load(handle)\n",
        "    \n",
        "    min_ = np.array(scale_params['min_'])\n",
        "    max_ = np.array(scale_params['max_'])\n",
        "    \n",
        "    return -1 + (2.0 / (max_ - min_) * (features - min_))\n",
        "\n",
        "def calculate_image_quality_score(brisque_features):\n",
        "    model = svmutil.svm_load_model('/content/sample_data/brisque_svm.txt')\n",
        "    scaled_brisque_features = scale_features(brisque_features)\n",
        "    \n",
        "    x, idx = svmutil.gen_svm_nodearray(\n",
        "        scaled_brisque_features,\n",
        "        isKernel=(model.param.kernel_type == svmutil.PRECOMPUTED))\n",
        "    \n",
        "    nr_classifier = 1\n",
        "    prob_estimates = (svmutil.c_double * nr_classifier)()\n",
        "    \n",
        "    return svmutil.libsvm.svm_predict_probability(model, x, prob_estimates)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ulijQYtzlGUP"
      },
      "source": [
        "The scaled used to represent image quality goes from 0 to 100. An image quality of 100 means that the image's quality is very bad. In the case of the analyzed image, we get that it is a good quality image. It makes sense because we are using the reference image."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "pycharm": {
          "is_executing": false
        },
        "id": "KPBiQRDRlGUP",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 66
        },
        "outputId": "8d2a1deb-610d-4344-f87a-f3bfd762aa0c"
      },
      "source": [
        "%%time\n",
        "\n",
        "calculate_image_quality_score(brisque_features)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "CPU times: user 11.6 ms, sys: 456 µs, total: 12.1 ms\n",
            "Wall time: 16.2 ms\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "43.9675731706412"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 136
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "C1s1dnFflGUR"
      },
      "source": [
        "## Conclusion\n",
        "\n",
        "This method was tested with the TID2008 database and performs well; even compared with referenced IQA methods. Although, the future work is to check the performance of other machine learning algorithms like XGBoost, LightGBM, for the pattern recognition step.\n",
        "\n",
        "## References\n",
        "\n",
        "[1] Maître, H. (2017). From Photon to pixel: the digital camera handbook. John Wiley & Sons.\n",
        "\n",
        "[2] Wang, Z., Bovik, A. C., Sheikh, H. R., & Simoncelli, E. P. (2004). Image quality assessment: from error visibility to structural similarity. IEEE transactions on image processing, 13(4), 600-612.\n",
        "\n",
        "[3] Ponomarenko, N., Lukin, V., Zelensky, A., Egiazarian, K., Carli, M., & Battisti, F. (2009). TID2008-a database for evaluation of full-reference visual quality assessment metrics. Advances of Modern Radioelectronics, 10(4), 30-45.\n",
        "\n",
        "[4] Mittal, A., Moorthy, A. K., & Bovik, A. C. (2012). No-reference image quality assessment in the spatial domain. IEEE Transactions on Image Processing, 21(12), 4695-4708.\n",
        "\n",
        "[5] Lasmar, N. E., Stitou, Y., & Berthoumieu, Y. (2009). Multiscale skewed heavy-tailed model for texture analysis. Proceedings - International Conference on Image Processing, ICIP, (1), 2281–2284."
      ]
    }
  ]
}