3.Methods/Metrics.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Metrics.ipynb",
      "provenance": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "code",
      "metadata": {
        "id": "tLn9kBueMbxJ"
      },
      "source": [
        "import numpy as np\r\n",
        "import math\r\n",
        "import scipy\r\n",
        "import pandas as pd\r\n",
        "import PIL\r\n",
        "import gdal\r\n",
        "import matplotlib.pyplot as plt\r\n",
        "plt.style.use('ggplot')\r\n",
        "import sys, os\r\n",
        "from pathlib import Path\r\n",
        "import time\r\n",
        "import random\r\n",
        "import collections, functools, operator\r\n",
        "import csv\r\n",
        "import subprocess\r\n",
        "import datetime\r\n",
        "\r\n",
        "from osgeo import gdal,osr\r\n",
        "from gdalconst import *\r\n",
        "import subprocess\r\n",
        "from osgeo.gdalconst import GA_Update\r\n",
        "\r\n",
        "import torch\r\n",
        "import torch.nn as nn\r\n",
        "\r\n",
        "import skimage\r\n",
        "from skimage import io, transform\r\n",
        "import sklearn\r\n",
        "import sklearn.metrics\r\n",
        "from sklearn.feature_extraction import image\r\n",
        "from sklearn import svm"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zWPCTZBRNQh8"
      },
      "source": [
        "# Metrics and utils for model training and evaluation"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "pGxuT8PXMfqm"
      },
      "source": [
        "def pcc(img1, img2): # Pearson's Correlation Coefficient\r\n",
        "    num = np.sum( (img1-np.mean(img1)) * (img2-np.mean(img2)) )\r\n",
        "    denom = np.sqrt(np.sum((img1-np.mean(img1))**2)) * np.sqrt(np.sum((img2-np.mean(img2))**2))\r\n",
        "    if denom==0.:\r\n",
        "        return 0\r\n",
        "    PCC = np.divide(num, denom)\r\n",
        "    return PCC\r\n",
        "def rmse(img1, img2): # Relative Mean Squared Error\r\n",
        "    RMSE = np.sqrt(np.mean((img1 - img2) ** 2))\r\n",
        "    return RMSE\r\n",
        "def psnr(img1, img2): # Peak Signal-to-Noise Ratio\r\n",
        "    PSNR = skimage.metrics.peak_signal_noise_ratio(img1, img2, data_range=img1.max() - img2.min())\r\n",
        "    return PSNR\r\n",
        "def ssim(img1, img2): # Structural Similarity\r\n",
        "    SSIM = skimage.metrics.structural_similarity(img1, img2, data_range=img1.max() - img1.min())\r\n",
        "    return SSIM\r\n",
        "def sam(px1, px2): # Spectral Angle Mapper\r\n",
        "    num = np.dot(px1, px2)\r\n",
        "    denom = np.sqrt(np.dot(px1, px1)) * np.sqrt(np.dot(px2, px2))\r\n",
        "    SAM = math.acos(num / denom)\r\n",
        "    return SAM\r\n",
        "def sid(px1, px2): # Spectral Information Divergence\r\n",
        "    if all(px2) == 0.:\r\n",
        "        return 1.\r\n",
        "    p = np.divide(px1, np.sum(px1))\r\n",
        "    q = np.divide(px2, np.sum(px2))\r\n",
        "    Dxy = np.sum(p * np.log(np.divide(p,q)))\r\n",
        "    Dyx = np.sum(q * np.log(np.divide(q,p)))\r\n",
        "    SID = Dxy + Dyx\r\n",
        "    return SID\r\n",
        "\r\n",
        "def bandWise_metrics(batch_real, batch_predicted):\r\n",
        "    #BAND-WISE EVALUATION\r\n",
        "    PCC = []\r\n",
        "    RMSE = []\r\n",
        "    PSNR = []\r\n",
        "    SSIM = []\r\n",
        "    for band in range(batch_real.shape[1]): # Iterate through bands\r\n",
        "        PCCband = 0\r\n",
        "        RMSEband = 0\r\n",
        "        PSNRband = 0\r\n",
        "        SSIMband = 0\r\n",
        "        for real_band, predicted_band in zip(batch_real[:,band], batch_predicted[:,band]): # Iterate through all patches in batch\r\n",
        "            PCCband += pcc(real_band, predicted_band)\r\n",
        "            RMSEband += rmse(real_band, predicted_band)\r\n",
        "            PSNRband += psnr(real_band, predicted_band)\r\n",
        "            SSIMband += ssim(real_band, predicted_band)\r\n",
        "        PCC.append(PCCband / batch_real.shape[0])\r\n",
        "        RMSE.append(RMSEband / batch_real.shape[0])\r\n",
        "        PSNR.append(PSNRband / batch_real.shape[0])\r\n",
        "        SSIM.append(SSIMband / batch_real.shape[0])\r\n",
        "    bwm = {'PCC': np.array(PCC), 'RMSE': np.array(RMSE), 'PSNR': np.array(PSNR), 'SSIM': np.array(SSIM)}\r\n",
        "    return bwm\r\n",
        "\r\n",
        "def pixelWise_metrics(batch_real, batch_predicted):\r\n",
        "    #PIXEL-WISE EVALUATION\r\n",
        "    SAM = []\r\n",
        "    SID = []\r\n",
        "    for i in range(batch_real.shape[2]):\r\n",
        "        for j in range(batch_real.shape[3]):\r\n",
        "            SAMpix = 0\r\n",
        "            SIDpix = []\r\n",
        "            for real_pix, predicted_pix in zip(batch_real[:,:,i,j], batch_predicted[:,:,i,j]):\r\n",
        "                SAMpix += sam(real_pix, predicted_pix)\r\n",
        "                SIDpix.append(sid(real_pix, predicted_pix))\r\n",
        "            SAM.append(SAMpix / batch_real.shape[0])\r\n",
        "            SIDpix = np.array(SIDpix)\r\n",
        "            SIDpix[np.isnan(SIDpix)] = np.mean(SIDpix[~np.isnan(SIDpix)])\r\n",
        "            SIDpix = np.sum(SIDpix)\r\n",
        "            SID.append(SIDpix / batch_real.shape[0])\r\n",
        "    pwm = {'SAM': np.array(SAM), 'SID': np.array(SID)}\r\n",
        "    return pwm\r\n",
        "\r\n",
        "\r\n",
        "def calc_metrics(real, predicted, verbose=False):\r\n",
        "    #BAND-WISE EVALUATION\r\n",
        "    bwm = bandWise_metrics(real, predicted)\r\n",
        "    PCC = bwm['PCC']\r\n",
        "    RMSE = bwm['RMSE']\r\n",
        "    PSNR = bwm['PSNR']\r\n",
        "    SSIM = bwm['SSIM']\r\n",
        "    \r\n",
        "    #PIXEL-WISE EVALUATION\r\n",
        "    pwm = pixelWise_metrics(real, predicted)\r\n",
        "    SAM = pwm['SAM']\r\n",
        "    SID = pwm['SID']\r\n",
        "\r\n",
        "    metrics = {'PCC': PCC,\r\n",
        "               'RMSE': RMSE,\r\n",
        "               'PSNR': PSNR,\r\n",
        "               'SSIM': SSIM,\r\n",
        "               'SAM': SAM,\r\n",
        "               'SID': SID}\r\n",
        "    if verbose:\r\n",
        "        show_metrics(metrics)\r\n",
        "    return metrics\r\n",
        "\r\n",
        "def show_metrics(metrics):\r\n",
        "    print('PCC:',np.mean(metrics['PCC']))\r\n",
        "    print('RMSE:',np.mean(metrics['RMSE']))\r\n",
        "    print('PSNR:',np.mean(metrics['PSNR']))\r\n",
        "    print('SSIM:',np.mean(metrics['SSIM']))\r\n",
        "    print('SAM:',np.mean(metrics['SAM']))\r\n",
        "    print('SID:',np.mean(metrics['SID']))\r\n",
        "\r\n",
        "    fig, axs = plt.subplots(2,3)\r\n",
        "    axs[0, 0].boxplot(metrics['PCC'], showmeans=True)\r\n",
        "    axs[0, 0].set_title('PCC')\r\n",
        "    axs[0, 1].boxplot(metrics['RMSE'], showmeans=True)\r\n",
        "    axs[0, 1].set_title('RMSE')\r\n",
        "    axs[1, 0].boxplot(metrics['PSNR'], showmeans=True)\r\n",
        "    axs[1, 0].set_title('PSNR')\r\n",
        "    axs[1, 1].boxplot(metrics['SSIM'], showmeans=True)\r\n",
        "    axs[1, 1].set_title('SSIM')\r\n",
        "    axs[0, 2].boxplot(metrics['SAM'], showmeans=True)\r\n",
        "    axs[0, 2].set_title('SAM')\r\n",
        "    axs[1, 2].boxplot(metrics['SID'], showmeans=True)\r\n",
        "    axs[1, 2].set_title('SID')\r\n",
        "    plt.show()\r\n",
        "\r\n",
        "\r\n",
        "\r\n",
        "def show_patches(input, prediction, target, saveSinglePatches=False): # Shows a few 64x64 patches from a batch\r\n",
        "    fig, axs = plt.subplots(3,10)\r\n",
        "    fig.set_figwidth(15)\r\n",
        "    for patch in range(10):\r\n",
        "        i = [ (input[patch][4]*65535)/4000, (input[patch][3]*65535)/4000, (input[patch][2]*65535)/4000 ]\r\n",
        "        p = [ (prediction[patch][160]*65535)/4000, (prediction[patch][45]*65535)/4000, (prediction[patch][9]*65535)/4000 ]\r\n",
        "        t = [ (target[patch][160]*65535)/4000,(target[patch][45]*65535)/4000,(target[patch][9]*65535)/4000 ]\r\n",
        "\r\n",
        "        for y in (i, p, t):\r\n",
        "            for x in y:\r\n",
        "                x[x<0] = 0\r\n",
        "                x[x>1] = 1\r\n",
        "\r\n",
        "        i = np.dstack((i[0],i[1],i[2]))\r\n",
        "        p = np.dstack((p[0],p[1],p[2]))\r\n",
        "        t = np.dstack((t[0],t[1],t[2]))\r\n",
        "        axs[0, patch].imshow(i)\r\n",
        "        axs[0, patch].set_title('Input')\r\n",
        "        axs[1, patch].imshow(p)\r\n",
        "        axs[1, patch].set_title('Prediction')\r\n",
        "        axs[2, patch].imshow(t)\r\n",
        "        axs[2, patch].set_title('Target')\r\n",
        "    plt.show()\r\n",
        "\r\n",
        "    if saveSinglePatches:\r\n",
        "        for i, patch in enumerate([0, 2, 4, 6]):\r\n",
        "            p = [ (prediction[patch][160]*65535)/4000, (prediction[patch][45]*65535)/4000, (prediction[patch][9]*65535)/4000 ] # Bands 45 is 37, 21 is 13, 14 is 6\r\n",
        "\r\n",
        "            for y in (p, t):\r\n",
        "                for x in y:\r\n",
        "                    x[x<0] = 0\r\n",
        "                    x[x>1] = 1\r\n",
        "                    \r\n",
        "            p = np.dstack((p[0],p[1],p[2]))\r\n",
        "            print('predicted')\r\n",
        "            img = plt.imshow(p)\r\n",
        "            plt.axis('off')\r\n",
        "            plt.savefig(os.getcwd() + '/drive/My Drive/TFG/PatchesViz/' + f\"SAMGAN_epoch400_p{i}.png\", bbox_inches='tight')\r\n",
        "            plt.show()\r\n",
        "\r\n",
        "            \r\n",
        "            print('real')\r\n",
        "            img = plt.imshow(t)\r\n",
        "            plt.axis('off')\r\n",
        "            plt.savefig(os.getcwd() + '/drive/My Drive/TFG/PatchesViz/' + f\"GT_p{i}.png\", bbox_inches='tight')\r\n",
        "            plt.show()\r\n",
        "    \r\n",
        "\r\n",
        "\r\n",
        "\r\n",
        "def weights_init(m):\r\n",
        "    if isinstance(m, nn.Conv2d) or isinstance(m, nn.ConvTranspose2d):\r\n",
        "        torch.nn.init.normal_(m.weight, 0.0, 0.02)\r\n",
        "    if isinstance(m, nn.BatchNorm2d):\r\n",
        "        torch.nn.init.normal_(m.weight, 0.0, 0.02)\r\n",
        "        torch.nn.init.constant_(m.bias, 0)\r\n",
        "\r\n",
        "\r\n",
        "\r\n",
        "class svmClassifier(object):\r\n",
        "    def __init__(self):\r\n",
        "        self.clf = svm.SVC(kernel='poly', degree=2, gamma='auto', C=1.0, decision_function_shape = \"ovr\")\r\n",
        "        self.PCs = [2, 6, 14, 23, 36, 43, 81, 110, 153]\r\n",
        "\r\n",
        "    def train(self, train_data):\r\n",
        "        spectral_data=[]\r\n",
        "        crop_class=[]\r\n",
        "        for sample in train_data:\r\n",
        "            spectral_data.append(net(sample['input'].resize(64*64,9,5,5).to(device)).detach().cpu().numpy())\r\n",
        "            crop_class.append(sample['crop'].numpy())\r\n",
        "\r\n",
        "        \r\n",
        "\r\n",
        "        spectral_data = np.array(spectral_data).reshape(len(spectral_data), 64, 64, 170).transpose(0,3,1,2)\r\n",
        "        crop_class = np.array(crop_class)\r\n",
        "        spectral_data = spectral_data.transpose(0,2,3,1).reshape(spectral_data.shape[0]*spectral_data.shape[2]*spectral_data.shape[3], spectral_data.shape[1])\r\n",
        "        crop_class = crop_class.transpose(0,2,3,1).reshape(crop_class.shape[0]*crop_class.shape[2]*crop_class.shape[3], crop_class.shape[1])\r\n",
        "\r\n",
        "        \r\n",
        "        spectral_PCs = []\r\n",
        "        for i in self.PCs: \r\n",
        "            spectral_PCs.append(spectral_data[:,i])\r\n",
        "        spectra_data = np.array(spectral_PCs)\r\n",
        "        \r\n",
        "\r\n",
        "\r\n",
        "        self.clf.fit(spectral_data, crop_class)\r\n",
        "    \r\n",
        "    def test(self, crop_class, spectral_data):\r\n",
        "        spectral_data = spectral_data.transpose(0,2,3,1).reshape(spectral_data.shape[0]*spectral_data.shape[2]*spectral_data.shape[3], spectral_data.shape[1])\r\n",
        "        crop_class = crop_class.transpose(0,2,3,1).reshape(crop_class.shape[0]*crop_class.shape[2]*crop_class.shape[3], crop_class.shape[1])\r\n",
        "\r\n",
        "        \r\n",
        "        PCs = [2, 6, 14, 23, 36, 43, 81, 110, 153]\r\n",
        "        spectral_PCs = []\r\n",
        "        for i in self.PCs: \r\n",
        "            spectral_PCs.append(spectral_data[:,i])\r\n",
        "        spectra_data = np.array(spectral_PCs)\r\n",
        "\r\n",
        "        return crop_class, self.clf.predict(spectral_data)"
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}