Random Forest Super-Resolution (RFSR).ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!/usr/bin/env python2 or python3\n",
    "# -*- coding: utf-8 -*-\n",
    "\"\"\"\n",
    "Created on Mon Feb 12 2018\n",
    "\n",
    "@author: jshermeyer\n",
    "\n",
    "Random Forest Super-Resolution (RFSR)\n",
    "Super resolves imagery by analyzing the relationships between adjacent LR pixels and the corresponding \n",
    "HR image pixel.  Uses a random forest regressor and a residual training schema.  Operates on the luminance\n",
    "component of a YCbCr converted image.\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "#imports\n",
    "import sys\n",
    "import os\n",
    "import glob\n",
    "import gdal\n",
    "import numpy as np\n",
    "import scipy\n",
    "from scipy import misc\n",
    "import sklearn\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "from sklearn.metrics import mean_squared_error\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.externals import joblib\n",
    "import datetime\n",
    "import cv2\n",
    "from tqdm import tqdm_notebook as tqdm\n",
    "\n",
    "#Various functions called throughout\n",
    "\n",
    "#Crop out the center of an array.\n",
    "def crop_center(array,cropx,cropy):\n",
    "    y,x = array.shape\n",
    "    startx = x//2-(cropx//2)\n",
    "    starty = y//2-(cropy//2)    \n",
    "    return array[starty:starty+cropy,startx:startx+cropx]\n",
    "\n",
    "#Convert an RGB image to an YCbCr array\n",
    "def Calc_YCbCr(RGB_Image):\n",
    "    srcRaster=gdal.Open(RGB_Image)\n",
    "    Blue = srcRaster.GetRasterBand(3).ReadAsArray()\n",
    "    Green = srcRaster.GetRasterBand(2).ReadAsArray()\n",
    "    Red = srcRaster.GetRasterBand(1).ReadAsArray()\n",
    "    shape=Blue.shape\n",
    "    Y=0+(0.2999 * Red)+(0.587 * Green)+(0.114 * Blue)\n",
    "    CB=128-(0.168736 * Red)-(0.331264 * Green)+(0.5 * Blue)\n",
    "    CR=128+(0.5 * Red)-(0.418688 * Green)-(0.081312 * Blue)\n",
    "    outArray=np.asarray([CR,Y,CB])    \n",
    "    return outArray\n",
    "\n",
    "    \n",
    "#Convert an RGB array to an YCbCr array    \n",
    "def Calc_YCbCr_Array(RGB_Array):\n",
    "    srcRaster=RGB_Array\n",
    "    #print(srcRaster.shape)\n",
    "    Blue = srcRaster[2,:,:]\n",
    "    Green = srcRaster[1,:,:]\n",
    "    Red = srcRaster[0,:,:]\n",
    "    shape=Blue.shape\n",
    "    #print(shape)\n",
    "    Y=0+(0.2999 * Red)+(0.587 * Green)+(0.114 * Blue)\n",
    "    CB=128-(0.168736 * Red)-(0.331264 * Green)+(0.5 * Blue)\n",
    "    CR=128+(0.5 * Red)-(0.418688 * Green)-(0.081312 * Blue)\n",
    "    outArray=np.asarray([CR,Y,CB])    \n",
    "    return outArray\n",
    "\n",
    "#Convert an YCbCr array to a RGB array  \n",
    "def Calc_RGB_Array(YcBcR_Array):\n",
    "    srcRaster=YcBcR_Array\n",
    "    #print(srcRaster.shape)\n",
    "    Cb = srcRaster[2,:,:]\n",
    "    Y = srcRaster[1,:,:]\n",
    "    Cr = srcRaster[0,:,:]\n",
    "    shape=Y.shape\n",
    "    #print(shape)\n",
    "    R  = Y + (Cr - 128) *  1.40200\n",
    "    G  = Y + (Cb - 128) * -0.34414 + (Cr - 128) * -0.71414\n",
    "    B  = Y + (Cb - 128) *  1.77200\n",
    "    outArray=np.asarray([R,G,B])    \n",
    "    return outArray\n",
    "     \n",
    "\n",
    "#Upscale a lower resolution array to match an HR arrays shape using bicubic interpolation\n",
    "def Upscale_Array(LR_Array,HR_Shape):\n",
    "    RasHolder=[]\n",
    "    for i in range(LR_Array.shape[0]):\n",
    "        LR_band=LR_Array[i]\n",
    "        LR_Resample=scipy.misc.imresize(LR_band,HR_Shape,interp='bicubic')\n",
    "        RasHolder.append(LR_Resample)\n",
    "    RasHolder=np.asarray(RasHolder)\n",
    "    return RasHolder\n",
    "\n",
    "#Shift an array around to create a 3D cube. \n",
    "#A shift_dimension of 1 takes into account the neighboring 1 pixels to the central pixel \n",
    "#https://i.stack.imgur.com/CWIHi.jpg (queen neighborhood)\n",
    "#A shift_dimension of 2 takes into account the neighboring 2 pixels to the central pixel\n",
    "#https://upload.wikimedia.org/wikipedia/commons/thumb/e/e1/Ising_model_5x5_0.svg/2000px-Ising_model_5x5_0.svg.png\n",
    "#(5x5 grid)\n",
    "#shift_dimension of 3==7x7 grid, 4===9x9 grid\n",
    "#The larger the grid the more memory intensive and more difficult training becomes\n",
    "#Recommend a dimension of 2 or 3, but typically 2 is sufficient.\n",
    "\n",
    "def Shift_Array(input_array,shift_dimension=2):\n",
    "    shift_holder=[]\n",
    "    current_shift_dimension_x=-shift_dimension\n",
    "    current_shift_dimension_y=-shift_dimension    \n",
    "    while current_shift_dimension_y<=shift_dimension:\n",
    "        if current_shift_dimension_x < shift_dimension:\n",
    "            Extent=np.roll(input_array,current_shift_dimension_y,axis=0)\n",
    "            Extent=np.roll(Extent,current_shift_dimension_x,axis=1)\n",
    "            shift_holder.append(Extent)           \n",
    "            current_shift_dimension_x+=1\n",
    "        elif current_shift_dimension_x == shift_dimension:\n",
    "            Extent=np.roll(input_array,current_shift_dimension_y,axis=0)\n",
    "            Extent=np.roll(Extent,current_shift_dimension_x,axis=1)\n",
    "            shift_holder.append(Extent)\n",
    "            current_shift_dimension_x=-shift_dimension\n",
    "            current_shift_dimension_y+=1\n",
    "            \n",
    "    shift_holder=np.asarray(shift_holder).astype(int)\n",
    "    return shift_holder\n",
    "\n",
    "#Downscale an array by a set factor.\n",
    "#If an image is 100x100 and the factor is 2, it will be downscaled to 50x50\n",
    "#Set blur>0 to mimic the PSF of a camera as it moves further away from an object or the surface\n",
    "#blur= the sigma of a gaussian blur you intend to mimic, we use 1.\n",
    "#set inter_area==1 to use the more robust and accurate inter area decimation when degrading an image\n",
    "#Otherwise, uses a bicubic decimation.\n",
    "def Downscale_Array(HR_Array,factor=2,blur=0,inter_area=0):\n",
    "    RasHolder=[]\n",
    "    #print(HR_Array.shape)\n",
    "    if blur > 0:\n",
    "        HR_Array=np.swapaxes(HR_Array,0,2)\n",
    "        blur_level=(factor/2)*blur\n",
    "        #print(HR_Array.shape)\n",
    "        HR_Array = cv2.GaussianBlur(HR_Array, (0, 0), blur_level, blur_level, 0)\n",
    "        HR_Array=np.swapaxes(HR_Array,2,0)\n",
    "        #print(HR_Array.shape)\n",
    "    if inter_area > 0:\n",
    "        del RasHolder\n",
    "        #print(HR_Array.shape)\n",
    "        HR_Array=np.swapaxes(HR_Array,0,2)\n",
    "        #print(HR_Array.shape)\n",
    "        #print(factor)\n",
    "        x=(1/float(factor))\n",
    "        #print(x)\n",
    "        RasHolder=cv2.resize(HR_Array, (0,0), fx=x,fy=x, interpolation=cv2.INTER_AREA)\n",
    "        RasHolder=np.swapaxes(RasHolder,0,2)\n",
    "    else:\n",
    "        for i in range(HR_Array.shape[0]):\n",
    "            HR_band=HR_Array[i]\n",
    "            HR_Shape=(int(HR_band.shape[0]/factor), int(HR_band.shape[1]/factor))\n",
    "            HR_Resample=scipy.misc.imresize(HR_band,HR_Shape,interp='bicubic')\n",
    "            RasHolder.append(HR_Resample)\n",
    "        RasHolder=np.asarray(RasHolder)\n",
    "    return RasHolder\n",
    "\n",
    "\n",
    "#Output multi or single band geotiffs\n",
    "def CreateMultiBandGeoTiff(Array, Name):\n",
    "    driver=gdal.GetDriverByName('GTiff')\n",
    "    DataSet = driver.Create(Name, Array.shape[2], Array.shape[1], Array.shape[0], gdal.GDT_Float32)\n",
    "    for i, image in enumerate(Array, 1):\n",
    "        DataSet.GetRasterBand(i).WriteArray( image )\n",
    "    del DataSet\n",
    "    return Name\n",
    "\n",
    "def CreateSingleBandGeoTiff(Array, Name):\n",
    "    driver=gdal.GetDriverByName('GTiff')\n",
    "    DataSet = driver.Create(Name, Array.shape[1], Array.shape[0], 1, gdal.GDT_Float32)\n",
    "    DataSet.GetRasterBand(1).WriteArray(Array)\n",
    "    del DataSet\n",
    "    return Name   \n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Create training data as an input to train RFSR\n",
    "#Index is the index of the luminance channel in a YCbCr converted image\n",
    "#Set =1 if using my functions for conversion, 0 for most others.\n",
    "def Create_SR_TrainingData(LR_Array,HR_Array,shift_dimension=2,index=1):\n",
    "    shift_window=((shift_dimension * 2) + 1) ** 2\n",
    "\n",
    "    #Upscale the raster to match HR_Raster\n",
    "    size1=HR_Array.shape[1]\n",
    "    size2=HR_Array.shape[2]\n",
    "    Bicube_LR=Upscale_Array(LR_Array,(size1,size2))\n",
    "    \n",
    "    HR_Shape=(LR_Array.shape[0],size1,size2)\n",
    "\n",
    "    #Convert to YCbCr\n",
    "    Bicube_LR=(Calc_YCbCr_Array(Bicube_LR)[index,:,:])\n",
    "    HR_Array=(Calc_YCbCr_Array(HR_Array)[index,:,:]) \n",
    "    \n",
    "    \n",
    "    #If imagery is irregularly shaped, you may want to uncomment this\n",
    "    #cropx=(HR_Shape[2])\n",
    "    #cropy=(HR_Shape[1])\n",
    "    #HR_Array=crop_center(HR_Array,cropx,cropy)\n",
    "    \n",
    "    #Optionally Pad\n",
    "    Bicube_LR = np.pad(Bicube_LR, pad_width=shift_dimension, mode='constant', constant_values=0)\n",
    "    HR_Array = np.pad(HR_Array, pad_width=shift_dimension, mode='constant', constant_values=0)\n",
    "\n",
    "    #Reshape the HR to 3D\n",
    "    HR_Array=HR_Array.reshape(1,HR_Shape[1]+(shift_dimension*2),HR_Shape[2]+(shift_dimension*2))\n",
    "\n",
    "    #Create patches\n",
    "    BC_LR_Shift=Shift_Array(Bicube_LR,shift_dimension=shift_dimension)\n",
    "    \n",
    "    #Subtract the Bicube_LR upscaled image\n",
    "    BC_LR_Shift=BC_LR_Shift-Bicube_LR\n",
    "    HR_Array=HR_Array-Bicube_LR    \n",
    "    \n",
    "    #Append the Arrays to ensure they are consistently located in 2d space before converting to 1d\n",
    "    Append_Array=np.append(BC_LR_Shift,HR_Array, axis=0)\n",
    "    nfeatures, nx, ny = Append_Array.shape\n",
    "    Append_Array=Append_Array.reshape(nfeatures,(nx*ny))\n",
    "    Append_Array=np.swapaxes(Append_Array,0,1)\n",
    "    Train=Append_Array[:,0:shift_window]\n",
    "    Target=Append_Array[:,shift_window]\n",
    "    \n",
    "    return Train, Target"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "## After importing all functions, start here\n",
    "\n",
    "#Path to your high res imagery.\n",
    "HR_Path=\"/HR_Images\"\n",
    "os.chdir(HR_Path)\n",
    "\n",
    "#The scale of enhancement\n",
    "ScalingFactor=2\n",
    "\n",
    "#The sigma of your gaussian blur,if blur = 0, no blurring happens, if > 0 sigma=blur\n",
    "blur=1\n",
    "\n",
    "## If -inter_area=1, cv2 interarea, otherwise bicubic decimation\n",
    "inter_area=1\n",
    "\n",
    "#Change this to whatever your extension is.\n",
    "HR=glob.glob('*.tif')\n",
    "\n",
    "#shift_dimension- see description above in the functions, I'd use 2, and perhaps 3 if you have time to kill.\n",
    "shift_dimension=2  \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#Create a training and testing dataset for scoring\n",
    "HR.sort()\n",
    "HR_Rasters=[]\n",
    "for raster in HR:\n",
    "    raster=gdal.Open(raster).ReadAsArray()\n",
    "    HR_Rasters.append(raster)\n",
    "    \n",
    "LR_Rasters=[]\n",
    "for raster in HR_Rasters:\n",
    "    LR_raster=Downscale_Array(raster,factor=ScalingFactor,blur=blur,inter_area=inter_area)\n",
    "    LR_Rasters.append(LR_raster)\n",
    "\n",
    "Train_TrainSet=[]\n",
    "Train_TestSet=[]\n",
    "Target_TrainSet=[]\n",
    "Target_TestSet=[]\n",
    "count=1\n",
    "t1=datetime.datetime.now()\n",
    "print(datetime.datetime.now())\n",
    "for LR, HR in (zip(LR_Rasters,HR_Rasters)):\n",
    "    train, target = Create_SR_TrainingData(LR,HR,shift_dimension=shift_dimension)\n",
    "    Train_TrainSet_temp, Train_TestSet_temp, Target_TrainSet_temp, Target_TestSet_temp = train_test_split(train, target, test_size=0.1,train_size=0.1)\n",
    "    \n",
    "    Train_TrainSet.append(Train_TrainSet_temp)\n",
    "    Train_TestSet.append(Train_TestSet_temp)\n",
    "    Target_TrainSet.append(Target_TrainSet_temp)\n",
    "    Target_TestSet.append(Target_TestSet_temp)\n",
    "    \n",
    "    Train_TrainSet=[np.concatenate(Train_TrainSet,axis=0)]\n",
    "    Train_TestSet=[np.concatenate(Train_TestSet,axis=0)]  \n",
    "    Target_TrainSet=[np.concatenate(Target_TrainSet,axis=0)]\n",
    "    Target_TestSet=[np.concatenate(Target_TestSet,axis=0)]\n",
    "    if count % 100 == 0:\n",
    "        print(count)\n",
    "    count+=1\n",
    "\n",
    "print(datetime.datetime.now())\n",
    "t2=datetime.datetime.now()\n",
    "print(t2-t1)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(20330310, 25)\n",
      "(20330310,)\n"
     ]
    }
   ],
   "source": [
    "#Check to ensure the first dimension of the shapes match.\n",
    "print(Train_TrainSet[0].shape)\n",
    "print(Target_TrainSet[0].shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Sample the training and testing set at 10% rate.  May want to increase this depending on your data.\n",
    "t1=datetime.datetime.now()\n",
    "print(datetime.datetime.now())\n",
    "#Default settings listed, adjust as you see fit.\n",
    "rf = RandomForestRegressor(n_estimators=100, max_depth=12,min_samples_split=200, n_jobs=-1,verbose=1,oob_score=True)\n",
    "print(rf)\n",
    "OutputModel=rf.fit(Train_TrainSet[0],Target_TrainSet[0])\n",
    "print(datetime.datetime.now())\n",
    "t2=datetime.datetime.now()\n",
    "print(t2-t1)\n",
    "mse=mean_squared_error( Target_TestSet[0],rf.predict(Train_TestSet[0]))\n",
    "#Training and testing scores, these will be close to what you will expect to see for a larger dataset\n",
    "print(\"MSE:\",mse)\n",
    "print(\"PSNR:\",20 * np.log10(255 / np.sqrt((mse)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Save your model somewhere\n",
    "joblib.dump(rf, '/models/rfmodel1.pkl')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#Inference code\n",
    "#shift_dimension again, must be identical to how you trained your model\n",
    "shift_dimension=2\n",
    "#Scaling factor\n",
    "SF=2\n",
    "#Load a model\n",
    "rf = joblib.load(\"/RF_Models/RF_2x_1Sigma_30cm_10percent_5x5PatchSize_Pad_100Depth_IA.pkl\")\n",
    "#These should be your low resolution images.  This code does not degrade then reupsample imagery.\n",
    "LR_Path=\"/images_to_be_enhanced\"\n",
    "#output directory\n",
    "output_dir=\"/output_path\"\n",
    "os.chdir(LR_Path)\n",
    "#Again change this to your extension.\n",
    "LR=glob.glob('*.tif')\n",
    "\n",
    "\n",
    "LR.sort()\n",
    "LR_Rasters=[]\n",
    "driver = gdal.GetDriverByName(\"GTiff\")\n",
    "for image in tqdm(LR):\n",
    "    raster=gdal.Open(image)\n",
    "    geo = raster.GetGeoTransform()\n",
    "    pixW=float(geo[1])/SF\n",
    "    pixH=float(geo[5])/SF\n",
    "    geo=[geo[0],pixW,geo[2],geo[3],geo[4],pixH]\n",
    "    proj = raster.GetProjection()\n",
    "    raster=raster.ReadAsArray()\n",
    "    #print(raster.shape)\n",
    "    size1=raster.shape[1]*SF\n",
    "    size2=raster.shape[2]*SF\n",
    "    raster=Upscale_Array(raster,(size1,size2))\n",
    "    \n",
    "    #Convert to YCbCr\n",
    "    YCbCr=(Calc_YCbCr_Array(raster))\n",
    "    Y=YCbCr[1,:,:]\n",
    "    Cr=YCbCr[0,:,:]\n",
    "    Cb=YCbCr[2,:,:]\n",
    "    \n",
    "    # Pad\n",
    "    Y = np.pad(Y, pad_width=shift_dimension, mode='constant', constant_values=0)\n",
    "    Cr = np.pad(Cr, pad_width=shift_dimension, mode='constant', constant_values=0)\n",
    "    Cb = np.pad(Cb, pad_width=shift_dimension, mode='constant', constant_values=0)\n",
    "    \n",
    "    #print(Y.shape,Cb.shape,Cr.shape)\n",
    "\n",
    "    #Create patches\n",
    "    Y_Shift=Shift_Array(Y,shift_dimension=shift_dimension)\n",
    "    \n",
    "    #Subtract the Bicube_LR upscaled image\n",
    "    Y_Shift=Y_Shift-Y \n",
    "    \n",
    "    #Create model ready input\n",
    "    nfeatures, nx, ny = Y_Shift.shape\n",
    "    Input=Y_Shift.reshape(nfeatures,(nx*ny))\n",
    "    Input=np.swapaxes(Input,0,1)\n",
    "    \n",
    "    #Infer\n",
    "    Output=rf.predict(Input)\n",
    "    Output=Output.reshape(Output.shape[0]//ny,-1)\n",
    "    #print(Output.shape)\n",
    "    Output=Output+Y\n",
    "    Stack=np.array([Cr,Output,Cb])\n",
    "    Stack=Stack[:, shift_dimension:-shift_dimension, shift_dimension:-shift_dimension]\n",
    "    Stack=Calc_RGB_Array(Stack)\n",
    "    \n",
    "    #Save\n",
    "    out=output_dir+str(image)\n",
    "    #print(out)\n",
    "    DataSet = driver.Create(out, Stack.shape[2], Stack.shape[1], Stack.shape[0], gdal.GDT_Byte)\n",
    "    for i, B in enumerate(Stack, 1):\n",
    "        DataSet.GetRasterBand(i).WriteArray( B )\n",
    "    DataSet.SetProjection(proj)\n",
    "    DataSet.SetGeoTransform(geo)\n",
    "    #DataSet.SetNoDataValue(0)\n",
    "    del DataSet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Code for simply testing your performance across multiple scaling factors for a different or whole dataset.\n",
    "\n",
    "from skimage.measure import compare_ssim\n",
    "\n",
    "shift_dimension=2\n",
    "SFs=[2,4,8]\n",
    "rfs=[]\n",
    "#Load models.\n",
    "rf = joblib.load(\"/RF_Models/RF_2x_1Sigma_120cm_10percent_5x5PatchSize_Pad_100Depth_IA.pkl\")\n",
    "rfs.append(rf)\n",
    "rf = joblib.load(\"/RF_Models/RF_4x_1Sigma_120cm_10percent_5x5PatchSize_Pad_100Depth_IA.pkl\")\n",
    "rfs.append(rf)\n",
    "rf = joblib.load(\"/RF_Models/RF_8x_1Sigma_120cm_10percent_5x5PatchSize_Pad_100Depth_IA.pkl\")\n",
    "rfs.append(rf)\n",
    "#This code DOES degrade and then upscale.  If you want to test performance for 30cm SR output feed in 30cm data\n",
    "LR_Path=\"/ImageryPath/\"\n",
    "os.chdir(LR_Path)\n",
    "#extension.\n",
    "LR=glob.glob('*.tif')\n",
    "LR.sort()\n",
    "\n",
    "for SF,rf in zip(SFs,rfs):\n",
    "    LR_Rasters=[]\n",
    "    mse=[]\n",
    "    ssim=[]\n",
    "    driver = gdal.GetDriverByName(\"GTiff\")\n",
    "    count=0\n",
    "    for image in tqdm(LR):\n",
    "        count+=1\n",
    "        raster=gdal.Open(image)\n",
    "        raster=raster.ReadAsArray()\n",
    "        YCbCr=(Calc_YCbCr_Array(raster))\n",
    "        HR_Y=YCbCr[1,:,:]\n",
    "        raster=Downscale_Array(raster,factor=SF,blur=1,inter_area=1)\n",
    "        raster=Upscale_Array(raster,HR_Y.shape)\n",
    "\n",
    "\n",
    "        #Convert to YCbCr\n",
    "        YCbCr=(Calc_YCbCr_Array(raster))\n",
    "        Y=YCbCr[1,:,:]\n",
    "\n",
    "        # Pad\n",
    "        Y = np.pad(Y, pad_width=shift_dimension, mode='constant', constant_values=0)\n",
    "\n",
    "        #Create patches\n",
    "        Y_Shift=Shift_Array(Y,shift_dimension=shift_dimension)\n",
    "\n",
    "        #Subtract the Bicube_LR upscaled image\n",
    "        Y_Shift=Y_Shift-Y \n",
    "\n",
    "        #Create model ready input\n",
    "        nfeatures, nx, ny = Y_Shift.shape\n",
    "        Input=Y_Shift.reshape(nfeatures,(nx*ny))\n",
    "        Input=np.swapaxes(Input,0,1)\n",
    "\n",
    "        #Infer\n",
    "        Output=rf.predict(Input)\n",
    "        Output=Output.reshape(Output.shape[0]//ny,-1)\n",
    "        Output=Output+Y\n",
    "        Output=Output[shift_dimension:-shift_dimension, shift_dimension:-shift_dimension]\n",
    "\n",
    "\n",
    "        #Test\n",
    "        mse.append(mean_squared_error(HR_Y,Output))\n",
    "        ssim.append(compare_ssim(HR_Y, Output, data_range=Output.max() - Output.min()))\n",
    "        #if count % 50 == 0:\n",
    "            #print(\"MSE:\",np.mean(mse))\n",
    "            #print(\"SSIM:\",np.nanmean(ssim))\n",
    "\n",
    "\n",
    "    print(\"SF:\", SF)\n",
    "    print(\"MSE:\",np.mean(mse))\n",
    "    print(\"PSNR:\",20 * np.log10(255 / np.sqrt(np.mean(mse))))\n",
    "    print(\"SSIM:\",np.nanmean(ssim))"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}