Conv3D-HSI.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "cd68aa96",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-11-19T12:26:27.655503Z",
     "start_time": "2021-11-19T12:26:24.091203Z"
    }
   },
   "outputs": [],
   "source": [
    "############ IMPORTS ####################\n",
    "import sys\n",
    "sys.path.append(\"./\")\n",
    "import os\n",
    "from os import path\n",
    "import numpy as np\n",
    "import random\n",
    "import torch\n",
    "import torch.utils.data as dataf\n",
    "import torch.nn as nn\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import io as sio\n",
    "from sklearn.decomposition import PCA\n",
    "from torch.nn.parameter import Parameter\n",
    "import torchvision.transforms.functional as TF\n",
    "import torch.nn.functional as F\n",
    "import time\n",
    "from PIL import Image\n",
    "import math\n",
    "from sklearn.model_selection import train_test_split\n",
    "from operator import truediv\n",
    "from sklearn.metrics import confusion_matrix, accuracy_score, cohen_kappa_score\n",
    "from torchsummary import summary\n",
    "import record"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "79f77acb",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-11-19T12:26:27.682324Z",
     "start_time": "2021-11-19T12:26:27.657748Z"
    }
   },
   "outputs": [],
   "source": [
    "def loadData(name):\n",
    "    \n",
    "    if name == \"Houston\":\n",
    "        dataHSI = sio.loadmat('./Houston/houston.mat')['houston']\n",
    "\n",
    "        \n",
    "\n",
    "        \n",
    "        labels = sio.loadmat('./Houston/houston_gt.mat')['houston_gt_te']\n",
    "        \n",
    "        labels += sio.loadmat('./Houston/houston_gt.mat')['houston_gt_tr']\n",
    "\n",
    "        \n",
    "    elif(name == \"Trento\"):\n",
    "        \n",
    "        dataHSI = sio.loadmat('./Trento/HSI.mat')['HSI']\n",
    "\n",
    "        \n",
    "        \n",
    "        labels = sio.loadmat('./Trento/TSLabel.mat')['TSLabel']\n",
    "        labels += sio.loadmat('./Trento/TRLabel.mat')['TRLabel']\n",
    "        \n",
    "    elif(name == \"MUUFL\"):\n",
    "        \n",
    "        dataHSI = sio.loadmat('./MUUFL/muufl_share.mat')['hsi_img']\n",
    "\n",
    "        \n",
    "        labels = sio.loadmat('./MUUFL/muufl_share.mat')['labels']\n",
    "    return dataHSI,labels\n",
    "\n",
    "\n",
    "def splitTrainTestSet(X, y, testRatio, randomState=345):\n",
    "    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=testRatio, random_state=randomState,\n",
    "                                                        stratify=y)\n",
    "    return X_train, X_test, y_train, y_test\n",
    "\n",
    "def applyPCA(X, numComponents=75):\n",
    "    newX = np.reshape(X, (-1, X.shape[2]))\n",
    "    pca = PCA(n_components=numComponents, whiten=True)\n",
    "    newX = pca.fit_transform(newX)\n",
    "    newX = np.reshape(newX, (X.shape[0],X.shape[1], numComponents))\n",
    "    return newX, pca\n",
    "\n",
    "def normalizeHSI(X):\n",
    "    for i in range(X.shape[2]):\n",
    "        minimal = X[:, :, i].min()\n",
    "        maximal = X[:, :, i].max()\n",
    "        X[:, :, i] = (X[:, :, i] - minimal)/(maximal - minimal)\n",
    "    return X\n",
    "        \n",
    "def normalizeLIDAR(X):\n",
    "    minimal = X.min()\n",
    "    maximal = X.max()\n",
    "    X = (X - minimal)/(maximal - minimal)\n",
    "    return X\n",
    "\n",
    "def padWithZeros(X, margin=2):\n",
    "    newX = np.zeros((X.shape[0] + 2 * margin, X.shape[1] + 2* margin, X.shape[2]))\n",
    "    x_offset = margin\n",
    "    y_offset = margin\n",
    "    newX[x_offset:X.shape[0] + x_offset, y_offset:X.shape[1] + y_offset, :] = X\n",
    "    return newX\n",
    "\n",
    "def createImageCubes(X, y, windowSize=5, removeZeroLabels = True):\n",
    "    margin = int((windowSize - 1) / 2)\n",
    "    zeroPaddedX = padWithZeros(X, margin=margin)\n",
    "    print(zeroPaddedX.shape)\n",
    "    # split patches\n",
    "    patchesData = np.zeros((X.shape[0] * X.shape[1], windowSize, windowSize, X.shape[2]), dtype=np.float32)\n",
    "    patchesLabels = np.zeros((X.shape[0] * X.shape[1]), dtype=int)\n",
    "    patchIndex = 0\n",
    "    for r in range(margin, zeroPaddedX.shape[0] - margin):\n",
    "        for c in range(margin, zeroPaddedX.shape[1] - margin):\n",
    "            patch = zeroPaddedX[r - margin:r + margin + 1, c - margin:c + margin + 1]   \n",
    "            patchesData[patchIndex, :, :, :] = patch\n",
    "            patchesLabels[patchIndex] = y[r-margin, c-margin]\n",
    "            patchIndex = patchIndex + 1\n",
    "    if removeZeroLabels:\n",
    "        patchesData = patchesData[patchesLabels>0,:,:,:]\n",
    "        patchesLabels = patchesLabels[patchesLabels>0]\n",
    "        \n",
    "    return patchesData, patchesLabels\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "84f9808a",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-11-19T12:26:27.697940Z",
     "start_time": "2021-11-19T12:26:27.683918Z"
    }
   },
   "outputs": [],
   "source": [
    "class CNN(nn.Module):\n",
    "    def __init__(self, FM, Classes, patchsize, NC):\n",
    "        super(CNN, self).__init__()\n",
    "        self.conv1 = nn.Sequential(\n",
    "            nn.Conv3d(\n",
    "                in_channels = 1,\n",
    "                out_channels = FM,\n",
    "                kernel_size = (3, 3, 7),\n",
    "                stride = 1,\n",
    "                padding = (0,0,0)\n",
    "            ),\n",
    "            nn.BatchNorm3d(FM),\n",
    "            nn.ReLU(),\n",
    "            nn.MaxPool3d(kernel_size=(1,1,2)),\n",
    "#             nn.Dropout(0.5),\n",
    "        )\n",
    "        self.final_bands = (NC - 6) // 2\n",
    "        \n",
    "        self.conv2 = nn.Sequential(\n",
    "            nn.Conv3d(FM, FM*2, (3, 3, 7 ), 1, (0,0,0)),\n",
    "            nn.BatchNorm3d(FM*2),\n",
    "            nn.ReLU(),\n",
    "            nn.MaxPool3d(kernel_size=(1,1,2))\n",
    "\n",
    "        )\n",
    "        self.final_bands = (self.final_bands - 6) // 2\n",
    "        self.conv3 = nn.Sequential(\n",
    "            nn.Conv3d(FM*2, FM*4, (3, 3, 7), 1, (0,0,0)),\n",
    "            nn.BatchNorm3d(FM*4),\n",
    "            nn.ReLU(),\n",
    "            nn.MaxPool3d(kernel_size=(1,1,2))\n",
    "\n",
    "        )\n",
    "        self.final_bands = (self.final_bands - 6) // 2\n",
    "        \n",
    "        self.final_patch_size = patchsize - 6\n",
    "        \n",
    "        \n",
    "        \n",
    "        self.out1 =  nn.Linear(self.final_patch_size * self.final_patch_size * FM * 4 * self.final_bands, Classes)\n",
    "#         self.out1 =  nn.Linear(19200, Classes)\n",
    "        \n",
    "    def forward(self, x1):\n",
    "        x1 = x1.unsqueeze(1)\n",
    "        x1 = self.conv1(x1)\n",
    "        x1 = self.conv2(x1)\n",
    "        x1 = self.conv3(x1)\n",
    "        x1 = x1.reshape(x1.shape[0], -1)  # flatten the output of conv2 to (batch_size, 32 * 7 * 7)\n",
    "        out1 = self.out1(x1)\n",
    "        return out1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "32081e72",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-11-19T12:26:49.949167Z",
     "start_time": "2021-11-19T12:26:27.700623Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Current patch size =  11\n",
      "Current training ratio =  0.1\n",
      "Current dataset =  Houston\n",
      "(359, 1915, 144)\n",
      "Train data shape =  torch.Size([1502, 11, 11, 144])\n",
      "Train label shape =  torch.Size([1502])\n",
      "Test data shape =  torch.Size([13527, 11, 11, 144])\n",
      "Test label shape =  torch.Size([13527])\n",
      "------------------------------------------------------------------------------------------\n",
      "Layer (type:depth-idx)                   Output Shape              Param #\n",
      "==========================================================================================\n",
      "├─Sequential: 1-1                        [-1, 16, 9, 9, 69]        --\n",
      "|    └─Conv3d: 2-1                       [-1, 16, 9, 9, 138]       1,024\n",
      "|    └─BatchNorm3d: 2-2                  [-1, 16, 9, 9, 138]       32\n",
      "|    └─ReLU: 2-3                         [-1, 16, 9, 9, 138]       --\n",
      "|    └─MaxPool3d: 2-4                    [-1, 16, 9, 9, 69]        --\n",
      "├─Sequential: 1-2                        [-1, 32, 7, 7, 31]        --\n",
      "|    └─Conv3d: 2-5                       [-1, 32, 7, 7, 63]        32,288\n",
      "|    └─BatchNorm3d: 2-6                  [-1, 32, 7, 7, 63]        64\n",
      "|    └─ReLU: 2-7                         [-1, 32, 7, 7, 63]        --\n",
      "|    └─MaxPool3d: 2-8                    [-1, 32, 7, 7, 31]        --\n",
      "├─Sequential: 1-3                        [-1, 64, 5, 5, 12]        --\n",
      "|    └─Conv3d: 2-9                       [-1, 64, 5, 5, 25]        129,088\n",
      "|    └─BatchNorm3d: 2-10                 [-1, 64, 5, 5, 25]        128\n",
      "|    └─ReLU: 2-11                        [-1, 64, 5, 5, 25]        --\n",
      "|    └─MaxPool3d: 2-12                   [-1, 64, 5, 5, 12]        --\n",
      "├─Linear: 1-4                            [-1, 15]                  288,015\n",
      "==========================================================================================\n",
      "Total params: 450,639\n",
      "Trainable params: 450,639\n",
      "Non-trainable params: 0\n",
      "------------------------------------------------------------------------------------------\n",
      "Input size (MB): 0.07\n",
      "Forward/backward pass size (MB): 4.85\n",
      "Params size (MB): 1.72\n",
      "Estimated Total Size (MB): 6.63\n",
      "------------------------------------------------------------------------------------------\n",
      "\n",
      "Epoch:  0 | train loss: 1.8393 | test accuracy: 0.10\n",
      "Epoch:  1 | train loss: 0.1246 | test accuracy: 0.08\n",
      "Time taken to train =  3.272874355316162 s\n",
      "--------Houston Training Finished-----------\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/users/jgecvision/.conda/envs/purb37/lib/python3.7/site-packages/torch/serialization.py:402: UserWarning: Couldn't retrieve source code for container of type CNN. It won't be checked for correctness upon loading.\n",
      "  \"type \" + obj.__name__ + \". It won't be checked \"\n"
     ]
    }
   ],
   "source": [
    "# CONFIGS\n",
    "os.environ[\"CUDA_VISIBLE_DEVICES\"]=\"1\"\n",
    "datasetNames = [\"Houston\"]\n",
    "testSizeNumber = 2500\n",
    "batchsize = 64\n",
    "EPOCH = 2\n",
    "LR = 0.001\n",
    "trainRatios = [0.1]\n",
    "patchSizes = [11]\n",
    "NUM_ITERATIONS = 1\n",
    "FM = 16\n",
    "\n",
    "\n",
    "def AA_andEachClassAccuracy(confusion_matrix):\n",
    "    counter = confusion_matrix.shape[0]\n",
    "    list_diag = np.diag(confusion_matrix)\n",
    "    list_raw_sum = np.sum(confusion_matrix, axis=1)\n",
    "    each_acc = np.nan_to_num(truediv(list_diag, list_raw_sum))\n",
    "    average_acc = np.mean(each_acc)\n",
    "    return each_acc, average_acc\n",
    "\n",
    "def reports (xtest,ytest,name):\n",
    "    pred_y = np.empty((len(ytest)), dtype=np.float32)\n",
    "    number = len(ytest) // testSizeNumber\n",
    "    for i in range(number):\n",
    "        temp = xtest[i * testSizeNumber:(i + 1) * testSizeNumber, :, :, :]\n",
    "        temp = temp.cuda()\n",
    "\n",
    "        temp2 = cnn(temp)\n",
    "        temp3 = torch.max(temp2, 1)[1].squeeze()\n",
    "        pred_y[i * testSizeNumber:(i + 1) * testSizeNumber] = temp3.cpu()\n",
    "        del temp, temp2, temp3\n",
    "\n",
    "    if (i + 1) * testSizeNumber < len(ytest):\n",
    "        temp = xtest[(i + 1) * testSizeNumber:len(ytest), :, :, :]\n",
    "        temp = temp.cuda()\n",
    "\n",
    "\n",
    "        temp2 = cnn(temp)\n",
    "        temp3 = torch.max(temp2, 1)[1].squeeze()\n",
    "        pred_y[(i + 1) * testSizeNumber:len(ytest)] = temp3.cpu()\n",
    "        del temp, temp2, temp3\n",
    "\n",
    "    pred_y = torch.from_numpy(pred_y).long()\n",
    "    \n",
    "    if name == 'Houston':\n",
    "        target_names = ['Healthy grass', 'Stressed grass', 'Synthetic grass'\n",
    "                        ,'Trees', 'Soil', 'Water', \n",
    "                        'Residential', 'Commercial', 'Road', 'Highway',\n",
    "                        'Railway', 'Parking Lot 1', 'Parking Lot 2', 'Tennis Court',\n",
    "                        'Running Track']\n",
    "    elif name == 'Trento':\n",
    "        target_names = ['Apples','Buildings','Ground','Woods','Vineyard',\n",
    "                        'Roads']\n",
    "    elif name == 'MUUFL':\n",
    "        target_names = ['Trees','Grass_Pure','Grass_Groundsurface','Dirt_And_Sand', 'Road_Materials','Water',\"Buildings'_Shadow\",\n",
    "                    'Buildings','Sidewalk','Yellow_Curb','ClothPanels']\n",
    "    elif name == 'IP':\n",
    "        target_names = ['Alfalfa', 'Corn-notill', 'Corn-mintill', 'Corn'\n",
    "                ,'Grass-pasture', 'Grass-trees', 'Grass-pasture-mowed', \n",
    "                'Hay-windrowed', 'Oats', 'Soybean-notill', 'Soybean-mintill',\n",
    "                'Soybean-clean', 'Wheat', 'Woods', 'Buildings-Grass-Trees-Drives',\n",
    "                'Stone-Steel-Towers']\n",
    "    elif name == 'SA':\n",
    "        target_names = ['Brocoli_green_weeds_1','Brocoli_green_weeds_2','Fallow','Fallow_rough_plow','Fallow_smooth',\n",
    "                        'Stubble','Celery','Grapes_untrained','Soil_vinyard_develop','Corn_senesced_green_weeds',\n",
    "                        'Lettuce_romaine_4wk','Lettuce_romaine_5wk','Lettuce_romaine_6wk','Lettuce_romaine_7wk',\n",
    "                        'Vinyard_untrained','Vinyard_vertical_trellis']\n",
    "    elif name == 'UP':\n",
    "        target_names = ['Asphalt','Meadows','Gravel','Trees', 'Painted metal sheets','Bare Soil','Bitumen',\n",
    "                        'Self-Blocking Bricks','Shadows']\n",
    "\n",
    "    \n",
    "    oa = accuracy_score(ytest, pred_y)\n",
    "    confusion = confusion_matrix(ytest, pred_y)\n",
    "    each_acc, aa = AA_andEachClassAccuracy(confusion)\n",
    "    kappa = cohen_kappa_score(ytest, pred_y)\n",
    "\n",
    "    return confusion, oa*100, each_acc*100, aa*100, kappa*100\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "for patchsize in patchSizes:\n",
    "    for trainRatio in trainRatios:\n",
    "        for datasetName in datasetNames:\n",
    "            print(\"Current patch size = \",patchsize)\n",
    "            print(\"Current training ratio = \",trainRatio)\n",
    "            print(\"Current dataset = \",datasetName)\n",
    "            \n",
    "            try:\n",
    "                os.makedirs(datasetName)\n",
    "            except FileExistsError:\n",
    "                pass\n",
    "            \n",
    "            X1, _, y = loadData(datasetName)\n",
    "            X1 = normalizeHSI(X1)\n",
    "\n",
    "            X1, yPatch = createImageCubes(X1, y, windowSize=patchsize)\n",
    "            TrainPatch, TestPatch, TrainLabel, TestLabel = splitTrainTestSet(X1, yPatch, 1. - trainRatio, randomState= 42)\n",
    "            TrainPatch = TrainPatch.astype(np.float32)\n",
    "            TestPatch = TestPatch.astype(np.float32)\n",
    "            NC = TrainPatch.shape[3]\n",
    "\n",
    "            TrainPatch = torch.from_numpy(TrainPatch)\n",
    "            TrainLabel = torch.from_numpy(TrainLabel)-1\n",
    "            TrainLabel = TrainLabel.long()\n",
    "            TestPatch = torch.from_numpy(TestPatch)\n",
    "            TestLabel = torch.from_numpy(TestLabel)-1\n",
    "            TestLabel = TestLabel.long()\n",
    "\n",
    "            Classes = len(np.unique(TrainLabel))\n",
    "\n",
    "            dataset = dataf.TensorDataset(TrainPatch, TrainLabel)\n",
    "            train_loader = dataf.DataLoader(dataset, batch_size=batchsize, shuffle=True)\n",
    "\n",
    "            print(\"Train data shape = \", TrainPatch.shape)\n",
    "            print(\"Train label shape = \", TrainLabel.shape)\n",
    "            print(\"Test data shape = \", TestPatch.shape)\n",
    "            print(\"Test label shape = \", TestLabel.shape)\n",
    "\n",
    "            KAPPA = []\n",
    "            OA = []\n",
    "            AA = []\n",
    "            ELEMENT_ACC = np.zeros((NUM_ITERATIONS, Classes))\n",
    "\n",
    "            for iterNum in range(NUM_ITERATIONS):    \n",
    "                cnn = CNN(FM, Classes, patchsize,NC)\n",
    "                cnn = cnn.cuda()\n",
    "                summary(cnn, (patchsize, patchsize, NC))\n",
    "                optimizer = torch.optim.Adam(cnn.parameters(), lr=LR)\n",
    "                loss_func = nn.CrossEntropyLoss()  # the target label is not one-hotted\n",
    "\n",
    "                BestAcc = 0\n",
    "                torch.cuda.synchronize()\n",
    "                start = time.time()\n",
    "\n",
    "                # train and test the designed model\n",
    "                for epoch in range(EPOCH):\n",
    "                    for step, (b_x1,b_y) in enumerate(train_loader):\n",
    "                        # move train data to GPU\n",
    "                        b_x1 = b_x1.cuda()\n",
    "                        b_y = b_y.cuda()\n",
    "\n",
    "                        out1 = cnn(b_x1)\n",
    "                        loss = loss_func(out1, b_y)\n",
    "\n",
    "                        optimizer.zero_grad()  # clear gradients for this training step\n",
    "                        loss.backward()  # backpropagation, compute gradients\n",
    "                        optimizer.step()  # apply gradients\n",
    "\n",
    "                        if step == len(train_loader) - 1:\n",
    "                            cnn.eval()\n",
    "                            pred_y = np.empty((len(TestLabel)), dtype='float32')\n",
    "                            number = len(TestLabel) // testSizeNumber\n",
    "                            for i in range(number):\n",
    "                                temp = TestPatch[i * testSizeNumber:(i + 1) * testSizeNumber, :, :, :]\n",
    "                                temp = temp.cuda()\n",
    "\n",
    "                                temp2 = cnn(temp)\n",
    "                                temp3 = torch.max(temp2, 1)[1].squeeze()\n",
    "                                pred_y[i * testSizeNumber:(i + 1) * testSizeNumber] = temp3.cpu()\n",
    "                                del temp, temp2, temp3\n",
    "\n",
    "\n",
    "                            if (i + 1) * testSizeNumber < len(TestLabel):\n",
    "                                temp = TestPatch[(i + 1) * testSizeNumber:len(TestLabel), :, :, :]\n",
    "                                temp = temp.cuda()\n",
    "\n",
    "                                temp2 = cnn(temp)\n",
    "                                temp3 = torch.max(temp2, 1)[1].squeeze()\n",
    "                                pred_y[(i + 1) * testSizeNumber:len(TestLabel)] = temp3.cpu()\n",
    "                                del temp, temp2, temp3\n",
    "\n",
    "                            pred_y = torch.from_numpy(pred_y).long()\n",
    "                            accuracy = torch.sum(pred_y == TestLabel).type(torch.FloatTensor) / TestLabel.size(0)\n",
    "\n",
    "                            print('Epoch: ', epoch, '| train loss: %.4f' % loss.data.cpu().numpy(), '| test accuracy: %.2f' % accuracy)\n",
    "\n",
    "                            # save the parameters in network\n",
    "                            if accuracy > BestAcc:\n",
    "                                BestAcc = accuracy\n",
    "                                torch.save(cnn.state_dict(), datasetName+'/net_params_checkpoint.pkl')\n",
    "                            cnn.train()\n",
    "\n",
    "                torch.cuda.synchronize()\n",
    "                end = time.time()\n",
    "                print(\"Time taken to train = \",end - start, \"s\")\n",
    "                Train_time = end - start\n",
    "\n",
    "                cnn.load_state_dict(torch.load(datasetName+'/net_params_checkpoint.pkl'))\n",
    "                cnn.eval()\n",
    "\n",
    "\n",
    "                confusion, oa, each_acc, aa, kappa = reports(TestPatch,TestLabel,datasetName)\n",
    "                KAPPA.append(kappa)\n",
    "                OA.append(oa)\n",
    "                AA.append(aa)\n",
    "                ELEMENT_ACC[iterNum, :] = each_acc\n",
    "                torch.save(cnn, datasetName+'/best_model_Conv3D-HSI_Iter'+str(iterNum)+'.pt')\n",
    "            print(\"--------\" + datasetName + \" Training Finished-----------\")\n",
    "            record.record_output(OA, AA, KAPPA, ELEMENT_ACC,'./' + datasetName +'/Conv3D-HSI_Report_' + datasetName +'.txt')\n",
    "                        \n",
    "                        \n",
    "        "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}