VGG16 as Feature Extractor for Population Screening.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This code is meant to use pretrained models to extract the features from the optimal layer for the normal and TB-infected chest radiographs to aid in improved population screening. However, you can use these codes as the skeleton to make use of pretrained models as feature extractors for your task of interest.Simply use this skeleton and extract the features from the most optimal layer from the model of your interest for the underlying data. You shall optimize the model hyperparameters to suit your data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To begin with, let us define a few functions to load the data and convert them to Keras compatible targets. We will load the libraries to begin with."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# load libraries\n",
    "import cv2\n",
    "import numpy as np\n",
    "import os\n",
    "from keras.utils import np_utils\n",
    "import matplotlib.pyplot as plt\n",
    "import itertools\n",
    "import time\n",
    "from keras.models import Sequential\n",
    "from keras.models import Model\n",
    "from keras import applications\n",
    "from keras.callbacks import ModelCheckpoint, TensorBoard\n",
    "from keras.layers import Conv2D, Activation, Dense, MaxPooling2D, Flatten, Dropout, GlobalAveragePooling2D\n",
    "from sklearn.metrics import log_loss\n",
    "from sklearn.utils import class_weight\n",
    "from keras.optimizers import SGD\n",
    "from sklearn.metrics import roc_curve, auc\n",
    "from sklearn.metrics import classification_report,confusion_matrix, accuracy_score\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.metrics import average_precision_score\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We performed 5-fold cross validation at the patient level. we had train and test splits for each fold to ensure that none of the patienet information in the training data leaks into the test data. We randomly split 10% of the training data for validation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#define data directories\n",
    "train_data_dir = 'f1_cxr/train' #path to your data\n",
    "valid_data_dir = 'f1_cxr/valid'\n",
    "test_data_dir = 'f1_cxr/test'\n",
    "\n",
    "# declare the number of samples in each category\n",
    "nb_train_samples = 475 #  modify for your dataset\n",
    "nb_valid_samples = 53 #  modify for your dataset\n",
    "nb_test_samples = 134 # modify for your dataset\n",
    "num_classes = 2 # binary classification \n",
    "\n",
    "img_rows_orig = 224 # modify these values for your requirements \n",
    "img_cols_orig = 224"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets define functions to load and resize the training, validation and test data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def load_training_data():\n",
    "    labels = os.listdir(train_data_dir)\n",
    "    total = len(labels)\n",
    "    X_train = np.ndarray((nb_train_samples, img_rows_orig, img_cols_orig, 3), dtype=np.uint8)\n",
    "    Y_train = np.zeros((nb_train_samples,), dtype='uint8')\n",
    "    i = 0\n",
    "    print('-'*30)\n",
    "    print('Creating training images...')\n",
    "    print('-'*30)\n",
    "    j = 0\n",
    "    for label in labels:\n",
    "        image_names_train = os.listdir(os.path.join(train_data_dir, label))\n",
    "        total = len(image_names_train)\n",
    "        print(label, total)\n",
    "        for image_name in image_names_train:\n",
    "            img = cv2.imread(os.path.join(train_data_dir, label, image_name), cv2.IMREAD_COLOR)\n",
    "            img = np.array([img])\n",
    "            X_train[i] = img\n",
    "            Y_train[i] = j\n",
    "            if i % 100 == 0:\n",
    "                print('Done: {0}/{1} images'.format(i, total))\n",
    "            i += 1\n",
    "        j += 1    \n",
    "    print(i)                \n",
    "    print('Loading done.')\n",
    "    print('Transform targets to keras compatible format.')\n",
    "    Y_train = np_utils.to_categorical(Y_train[:nb_train_samples], num_classes)\n",
    "    np.save('imgs_train.npy', X_train, Y_train) #save as numpy files\n",
    "    return X_train, Y_train\n",
    "    \n",
    "def load_validation_data():\n",
    "    # Load validation images\n",
    "    labels = os.listdir(valid_data_dir)\n",
    "    X_valid = np.ndarray((nb_valid_samples, img_rows_orig, img_cols_orig, 3), dtype=np.uint8)\n",
    "    Y_valid = np.zeros((nb_valid_samples,), dtype='uint8')\n",
    "    i = 0\n",
    "    print('-'*30)\n",
    "    print('Creating validation images...')\n",
    "    print('-'*30)\n",
    "    j = 0\n",
    "    for label in labels:\n",
    "        image_names_valid = os.listdir(os.path.join(valid_data_dir, label))\n",
    "        total = len(image_names_valid)\n",
    "        print(label, total)\n",
    "        for image_name in image_names_valid:\n",
    "            img = cv2.imread(os.path.join(valid_data_dir, label, image_name), cv2.IMREAD_COLOR)\n",
    "            img = np.array([img])\n",
    "            X_valid[i] = img\n",
    "            Y_valid[i] = j\n",
    "            if i % 100 == 0:\n",
    "                print('Done: {0}/{1} images'.format(i, total))\n",
    "            i += 1\n",
    "        j += 1\n",
    "    print(i)            \n",
    "    print('Loading done.')\n",
    "    print('Transform targets to keras compatible format.');\n",
    "    Y_valid = np_utils.to_categorical(Y_valid[:nb_valid_samples], num_classes)\n",
    "    np.save('imgs_valid.npy', X_valid, Y_valid) #save as numpy files\n",
    "    return X_valid, Y_valid\n",
    "\n",
    "def load_test_data():\n",
    "    labels = os.listdir(test_data_dir)\n",
    "    X_test = np.ndarray((nb_test_samples, img_rows_orig, img_cols_orig, 3), dtype=np.uint8)\n",
    "    Y_test = np.zeros((nb_test_samples,), dtype='uint8')\n",
    "    i = 0\n",
    "    print('-'*30)\n",
    "    print('Creating test images...')\n",
    "    print('-'*30)\n",
    "    j = 0\n",
    "    for label in labels:\n",
    "        image_names_test = os.listdir(os.path.join(test_data_dir, label))\n",
    "        total = len(image_names_test)\n",
    "        print(label, total)\n",
    "        for image_name in image_names_test:\n",
    "            img = cv2.imread(os.path.join(test_data_dir, label, image_name), cv2.IMREAD_COLOR)\n",
    "            img = np.array([img])\n",
    "            X_test[i] = img\n",
    "            Y_test[i] = j\n",
    "            if i % 100 == 0:\n",
    "                print('Done: {0}/{1} images'.format(i, total))\n",
    "            i += 1\n",
    "        j += 1\n",
    "    print(i)            \n",
    "    print('Loading done.')\n",
    "    print('Transform targets to keras compatible format.');\n",
    "    Y_test = np_utils.to_categorical(Y_test[:nb_test_samples], num_classes)\n",
    "    np.save('imgs_test.npy', X_test, Y_test) #save as numpy files\n",
    "    return X_test, Y_test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will define functions to resize the original images to that dimensions required for the pretrained models using the functions defined below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def load_resized_training_data(img_rows, img_cols):\n",
    "\n",
    "    X_train, Y_train = load_training_data()\n",
    "    X_train = np.array([cv2.resize(img, (img_rows,img_cols)) for img in X_train[:nb_train_samples,:,:,:]])\n",
    "    \n",
    "    return X_train, Y_train\n",
    "    \n",
    "def load_resized_validation_data(img_rows, img_cols):\n",
    "\n",
    "    X_valid, Y_valid = load_validation_data()\n",
    "    X_valid = np.array([cv2.resize(img, (img_rows,img_cols)) for img in X_valid[:nb_valid_samples,:,:,:]])\n",
    "        \n",
    "    return X_valid, Y_valid   \n",
    "\n",
    "def load_resized_test_data(img_rows, img_cols):\n",
    "\n",
    "    X_test, Y_test = load_test_data()\n",
    "    X_test = np.array([cv2.resize(img, (img_rows,img_cols)) for img in X_test[:nb_test_samples,:,:,:]])\n",
    "    \n",
    "    return X_test, Y_test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An evaluation script has been written to compute the confusion matrix for the performance of the trained model. This function prints and plots the confusion matrix. Normalization can be applied by setting 'normalize=True'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def plot_confusion_matrix(cm, classes,\n",
    "                          normalize=False, #if true all values in confusion matrix is between 0 and 1\n",
    "                          title='Confusion matrix',\n",
    "                          cmap=plt.cm.Blues):\n",
    "    plt.imshow(cm, interpolation='nearest', cmap=cmap)\n",
    "    plt.title(title)\n",
    "    plt.colorbar()\n",
    "    tick_marks = np.arange(len(classes))\n",
    "    plt.xticks(tick_marks, classes, rotation=45)\n",
    "    plt.yticks(tick_marks, classes)\n",
    "\n",
    "    if normalize:\n",
    "        cm = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis]\n",
    "        print(\"Normalized confusion matrix\")\n",
    "    else:\n",
    "        print('Confusion matrix, without normalization')\n",
    "\n",
    "    print(cm)\n",
    "    thresh = cm.max() / 2.\n",
    "    for i, j in itertools.product(range(cm.shape[0]), range(cm.shape[1])):\n",
    "        plt.text(j, i, cm[i, j],\n",
    "                 horizontalalignment=\"center\",\n",
    "                 color=\"white\" if cm[i, j] > thresh else \"black\")\n",
    "    plt.tight_layout()\n",
    "    plt.ylabel('True label')\n",
    "    plt.xlabel('Predicted label')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now proceed to extract the features from our dataset using the VGG16 pretrained model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "img_rows=224 \n",
    "img_cols=224\n",
    "channel = 3 #RGB\n",
    "num_classes = 2 #binary classification\n",
    "batch_size = 16 # modify based on the GPUs in your system\n",
    "num_epoch = 100 # modify depending on the model's convergence with your data\n",
    "\n",
    "#load data\n",
    "X_train, Y_train = load_resized_training_data(img_rows, img_cols)\n",
    "X_valid, Y_valid = load_resized_validation_data(img_rows, img_cols)\n",
    "X_test, Y_test = load_resized_test_data(img_rows, img_cols)\n",
    "\n",
    "\n",
    "#print the shape of the data\n",
    "print(X_train.shape, Y_train.shape)\n",
    "print(X_valid.shape, Y_valid.shape)\n",
    "print(X_test.shape, Y_test.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "base_model = applications.VGG16(weights='imagenet', include_top=False, input_shape=(img_rows, img_cols, channel))\n",
    "\n",
    "#extract feature from the optimal layer for your data\n",
    "base_model = Model(input=base_model.input, output=base_model.get_layer('block5_conv1').output) #depends on your data \n",
    "\n",
    "#get the model summary\n",
    "base_model.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets modify the architecture by adding a global spatial average pooling layer and a fully-connected layer with a dropout ratio of 0.5 to prevent overfitting and help model generalization. We will train only the top layers which are randomly initialized, freeze all the convolutional layers to prevent large gradient updates wrecking the learned weights. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x = base_model.output\n",
    "x = GlobalAveragePooling2D()(x)\n",
    "x = Dense(1024, activation='relu')(x)\n",
    "x = Dropout(0.5)(x)\n",
    "predictions = Dense(num_classes, activation='softmax', name='predictions')(x)\n",
    "\n",
    "# this is the model we will train\n",
    "model = Model(inputs=base_model.input, outputs=predictions)\n",
    "\n",
    "# Train only the top layers (which were randomly initialized)\n",
    "# i.e. freeze all convolutional layers to prevent large gradient updates wrecking the learned weights\n",
    "for layer in base_model.layers:\n",
    "    layer.trainable = False\n",
    "\n",
    "#fix the optimizer\n",
    "sgd = SGD(lr=0.0001, decay=1e-6, momentum=0.9, nesterov=True) \n",
    "\n",
    "#compile the model\n",
    "model.compile(optimizer=sgd,\n",
    "              loss='categorical_crossentropy',\n",
    "              metrics=['accuracy'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Its time to train the model. We will store only the best model weights by initializing callbacks. Also we can view the performance of our model during run-time by visualizing the performance graphs with Tensorboard. Create a log directory named 'logs' to store the training logs and a separate folder named 'weights' to store the model weights. You can visualize tensorboard graphs simply by navigating to your working directory and do:\n",
    "\n",
    "$tensorboard --logdir=path/to/log-directory/ --port 6006\n",
    "\n",
    "Then open localhost:6006 in your browser to view the performance graphs, model architecture and other parameters of your interest."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "filepath = 'weights/' + model.name + '.{epoch:02d}-{val_acc:.4f}.hdf5'\n",
    "checkpoint = ModelCheckpoint(filepath, monitor='val_acc', verbose=1, \n",
    "                             save_weights_only=True, save_best_only=True, mode='max', period=1)\n",
    "tensor_board = TensorBoard(log_dir='logs/', histogram_freq=0, batch_size=batch_size)\n",
    "callbacks_list = [checkpoint, tensor_board]\n",
    "\n",
    "#compute training time\n",
    "t=time.time()\n",
    "hist = model.fit(X_train, Y_train, batch_size=batch_size, \n",
    "                 callbacks=callbacks_list,\n",
    "                 epochs=num_epoch, verbose=1, \n",
    "                 shuffle=True, validation_data=[X_valid, Y_valid])\n",
    "\n",
    "#compute the training time\n",
    "print('Training time: %s' % (time.time()-t))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you want to visualize the performance of the model in the console other than with Tensorboard, you can use the following script."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "train_loss=hist.history['loss']\n",
    "val_loss=hist.history['val_loss']\n",
    "train_acc=hist.history['acc']\n",
    "val_acc=hist.history['val_acc']\n",
    "xc=range(num_epoch)\n",
    "\n",
    "plt.figure(1,figsize=(20,10), dpi=100)\n",
    "plt.plot(xc,train_loss)\n",
    "plt.plot(xc,val_loss)\n",
    "plt.xlabel('num of Epochs')\n",
    "plt.ylabel('loss')\n",
    "plt.title('train_loss vs val_loss')\n",
    "plt.grid(True)\n",
    "plt.legend(['train','val'])\n",
    "plt.style.use(['classic'])\n",
    "\n",
    "plt.figure(2,figsize=(20,10), dpi=100)\n",
    "plt.plot(xc,train_acc)\n",
    "plt.plot(xc,val_acc)\n",
    "plt.xlabel('num of Epochs')\n",
    "plt.ylabel('accuracy')\n",
    "plt.title('train_acc vs val_acc')\n",
    "plt.grid(True)\n",
    "plt.legend(['train','val'])\n",
    "plt.style.use(['classic'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once the model is trained, load the best model weights to predict on the test data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model.load_weights('weights/model_2.01-0.8546.hdf5') #path to your model weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#predict on the test data\n",
    "X_test, Y_test = load_resized_test_data(img_rows, img_cols)\n",
    "print(X_test.shape, Y_test.shape)\n",
    "print('-'*30)\n",
    "print('Predicting on the test data...')\n",
    "print('-'*30)\n",
    "y_pred = model.predict(X_test, batch_size=batch_size, verbose=1)\n",
    "\n",
    "# compute the accuracy\n",
    "Test_accuracy = accuracy_score(Y_test.argmax(axis=-1),y_pred.argmax(axis=-1))\n",
    "print(\"Test_Accuracy = \",Test_accuracy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us now compute the performance metrics for the pretrained VGG16 model with the test data. The performance metrics involve computing the ROC-AUC values, cross-entropy loss score, average precision score, prediction probabilities and storing these values and plotting the ROC curves."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#compute the ROC-AUC values\n",
    "fpr = dict()\n",
    "tpr = dict()\n",
    "roc_auc = dict()\n",
    "for i in range(num_classes):\n",
    "    fpr[i], tpr[i], _ = roc_curve(Y_test[:, i], y_pred[:, i])\n",
    "    roc_auc[i] = auc(fpr[i], tpr[i])\n",
    "    \n",
    "# Compute micro-average ROC curve and ROC area\n",
    "fpr[\"micro\"], tpr[\"micro\"], _ = roc_curve(Y_test.ravel(), y_pred.ravel())\n",
    "roc_auc[\"micro\"] = auc(fpr[\"micro\"], tpr[\"micro\"])\n",
    "\n",
    "#Plot ROC curves\n",
    "plt.figure(figsize=(20,10), dpi=100)\n",
    "lw = 1\n",
    "plt.plot(fpr[1], tpr[1], color='red',\n",
    "         lw=lw, label='ROC curve (area = %0.4f)' % roc_auc[1])\n",
    "plt.plot([0, 1], [0, 1], color='black', lw=lw, linestyle='--')\n",
    "plt.xlim([0.0, 1.0])\n",
    "plt.ylim([0.0, 1.05])\n",
    "plt.xlabel('False Positive Rate')\n",
    "plt.ylabel('True Positive Rate')\n",
    "plt.title('Receiver operating characteristics')\n",
    "plt.legend(loc=\"lower right\")\n",
    "plt.show()\n",
    "\n",
    "# computhe the cross-entropy loss score\n",
    "score = log_loss(Y_test,y_pred)\n",
    "print(score)\n",
    "\n",
    "# compute the average precision score\n",
    "prec_score = average_precision_score(Y_test,y_pred)  \n",
    "print(prec_score)\n",
    "\n",
    "# transfer it back\n",
    "y_pred = np.argmax(y_pred, axis=1)\n",
    "Y_test = np.argmax(Y_test, axis=1)\n",
    "print(y_pred)\n",
    "print(Y_test)\n",
    "\n",
    "#save the predictions as a CSV file for further analysis\n",
    "np.savetxt('vgg16_model_y_pred.csv',y_pred,fmt='%i',delimiter = \",\")\n",
    "np.savetxt('vgg16_model_Y_test.csv',Y_test,fmt='%i',delimiter = \",\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let us plot the confusion matrix of the model's performance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "target_names = ['class 0(abnormal)', 'class 1(normal)'] #decide the labels for your own data\n",
    "print(classification_report(Y_test,y_pred,target_names=target_names))\n",
    "print(confusion_matrix(Y_test,y_pred))\n",
    "cnf_matrix = (confusion_matrix(Y_test,y_pred))\n",
    "np.set_printoptions(precision=4)\n",
    "\n",
    "# Plot non-normalized confusion matrix\n",
    "plt.figure(figsize=(20,10), dpi=100)\n",
    "plot_confusion_matrix(cnf_matrix, classes=target_names,\n",
    "                  title='Confusion matrix')\n",
    "plt.show()"
   ]
  }
 ],
 "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}