Model_testing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Note: \n",
    " In this example code, Python is used for image processing (registration, skull stripping and N4 bias correction via Nipype), applying pretrained PyTorch models (CNN for whole tumor segmentation, CNN classifier for IDH status prediction), and extraction of loci and shape features. \n",
    "  R is used for application of the pretrained radiomics classifier to shpae and loci features, and application of the pretrained logistic model that connects the predicted probability from CNN classifier and radiomics classfier and yield the final probability for IDH status, which is the final step of hybrid model. \n",
    "  \n",
    "  \n",
    "### Python module requirements \n",
    "#### 1. Nipype \n",
    "#### 2. FSL\n",
    "#### 3. ANTs\n",
    "#### 4. PyRadiomics\n",
    "#### 5. PyTorch\n",
    "#### 6. rpy2 (with R)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from img_processing import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Image processing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Resample T1C to 1mm isovoxel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "T1C_original_file = 'example_t1c.nii.gz'            # file path of original T1C \n",
    "T2_original_file = 'example_t2.nii.gz'\n",
    "FLAIR_original_file = 'example_flair.nii.gz'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Isovoxel resampling, registration, skull stripping, bias correction, normalization, and resizing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "resampling T1C - completed\n",
      "register T2 to T1C_isovoxel - completed\n",
      "register FLAIR to T1C_isovoxel - completed\n",
      "191129-10:52:31,681 interface INFO:\n",
      "\t stdout 2019-11-29T10:52:31.681127:Could not find a supported file with prefix \"T1C_bet0_temp_tmp_unbiased_forskull_skull\"\n",
      "Acquired BET mask...\n",
      "Skull stripping of T1C, T2, FLAIR... - done\n",
      "T1C bias correction done...\n",
      "T2 bias correction done...\n",
      "FLAIR bias correction done...\n",
      "Image SI normalization & resizing :  done...\n"
     ]
    }
   ],
   "source": [
    "# filename for isovoxel registered images\n",
    "T1C_iso_file = 't1c_isovoxel.nii.gz'\n",
    "T2_iso_file = 't2_isovoxel.nii.gz'\n",
    "FLAIR_iso_file = 'flair_isovoxel.nii.gz'\n",
    "# filename for skull-stripped images   \n",
    "T1C_bet_file = 't1c_bet.nii.gz'\n",
    "T2_bet_file = 't2_bet.nii.gz'\n",
    "FLAIR_bet_file = 'flair_bet.nii.gz'\n",
    "# filename for bias-corrected images   \n",
    "T1C_corrected_file = 't1c_corrected.nii.gz'\n",
    "T2_corrected_file = 't2_corrected.nii.gz'\n",
    "FLAIR_corrected_file = 'flair_corrected.nii.gz'\n",
    "## filename for preliminary skull-stripped T1C                                       \n",
    "T1C_bet_temp_file = 'T1C_bet0_temp.nii.gz'\n",
    "\n",
    "(t1c_unet_arr, flair_unet_arr, cropdown_info) = func_img_proc(T1C_original_file, T2_original_file, FLAIR_original_file, \n",
    "                                                              T1C_iso_file, T2_iso_file, FLAIR_iso_file,\n",
    "                                                              T1C_bet_file, T2_bet_file, FLAIR_bet_file,\n",
    "                                                              T1C_corrected_file, T2_corrected_file, FLAIR_corrected_file,\n",
    "                                                              T1C_bet_temp_file)\n",
    "\n",
    "brain_mask_file = T1C_bet_temp_file[:len(T1C_bet_temp_file)-len('.nii.gz')] + '_mask.nii.gz'\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####  Model 1 : Automatic tumor segmentation  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calling pretrained Model 1...\n",
      "Acquiring predicted tumor segmentation mask...\n"
     ]
    }
   ],
   "source": [
    "predmask_arr = func_get_predmask(t1c_unet_arr, flair_unet_arr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>res_orig</th>\n",
       "      <th>x_min2</th>\n",
       "      <th>x_max2</th>\n",
       "      <th>x_len</th>\n",
       "      <th>y_min2</th>\n",
       "      <th>y_max2</th>\n",
       "      <th>y_len</th>\n",
       "      <th>z_min2</th>\n",
       "      <th>z_max2</th>\n",
       "      <th>z_len</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>(172, 172, 140)</td>\n",
       "      <td>31</td>\n",
       "      <td>203</td>\n",
       "      <td>240</td>\n",
       "      <td>29</td>\n",
       "      <td>201</td>\n",
       "      <td>240</td>\n",
       "      <td>16</td>\n",
       "      <td>156</td>\n",
       "      <td>170</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          res_orig  x_min2  x_max2  x_len  y_min2  y_max2  y_len  z_min2  \\\n",
       "0  (172, 172, 140)      31     203    240      29     201    240      16   \n",
       "\n",
       "   z_max2  z_len  \n",
       "0     156    170  "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cropdown_info"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "predmask_isovoxel_arr = func_mask_back2iso(predmask_arr, cropdown_info)\n",
    "predmask_isovoxel_arr_sitk = np.transpose(predmask_isovoxel_arr, (2,1,0))\n",
    "predmask_isovoxel_img = sitk.GetImageFromArray(predmask_isovoxel_arr_sitk)\n",
    "\n",
    "predmask_isovoxel_file = 'predmask_isovoxel.nii.gz' #filename for predicted mask of isovoxel resolution\n",
    "sitk.WriteImage(predmask_isovoxel_img, predmask_isovoxel_file)   # save the automatic segmentation of isovoxel resolution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    " ### normalization for ResNet  -> resize to 128 x 128 x 128 \n",
    "t1c_corrected_img = nb.load(T1C_corrected_file)\n",
    "t1c_corrected_arr = t1c_corrected_img.get_data()\n",
    "t2_corrected_img = nb.load(T2_corrected_file)\n",
    "t2_corrected_arr = t2_corrected_img.get_data()\n",
    "\n",
    "brain_mask = nb.load(brain_mask_file)\n",
    "brain_mask_arr = brain_mask.get_data()\n",
    "\n",
    "t1c_resnet_arr = func_norm_resnet(t1c_corrected_arr, predmask_isovoxel_arr, brain_mask_arr, cropdown_info)\n",
    "t2_resnet_arr = func_norm_resnet(t2_corrected_arr, predmask_isovoxel_arr, brain_mask_arr, cropdown_info)\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Model 2 : CNN classifier for IDH status prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "from resnet_model import *    \n",
    "def get_IDH_pred(t1c_resnet_arr, t2_resnet_arr, mask_arr):\n",
    "    \n",
    "    print(\"Calling pretrained Model 2...\")\n",
    "    model_resnet = ResNet(3, BasicBlock, [3,4,6,3])  \n",
    "    model_filename = 'MODEL2_ResNet_CNNclassifier.pth'\n",
    "    checkpoint = torch.load(model_filename)\n",
    "    model_resnet.load_state_dict(checkpoint['model_state_dict'])\n",
    "    model_resnet.eval()\n",
    "    model_resnet.cuda()\n",
    "    \n",
    "    print(\"Calculating a predicted probabilitiy...\")\n",
    "    x_arr = get_maxROI(t1c_resnet_arr, t2_resnet_arr, mask_arr)\n",
    "    \n",
    "    x_arr = torch.from_numpy(x_arr).float()\n",
    "    x_arr = x_arr.cuda()\n",
    "    \n",
    "    with torch.no_grad():\n",
    "        output = model_resnet(x_arr)\n",
    "    \n",
    "    output = nn.Softmax(dim=1)(output)\n",
    "    \n",
    "    output_mean = torch.sum(output[:,1])/5\n",
    "    \n",
    "    return output, output_mean\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calling pretrained Model 2...\n",
      "Calculating a predicted probabilitiy...\n"
     ]
    }
   ],
   "source": [
    "resnet_prob, resnet_prob_mean = get_IDH_pred(t1c_resnet_arr, t2_resnet_arr, predmask_arr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model 3 : radiomics classifier "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Shape and location feature extraction\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Shape_3d_Elongation</th>\n",
       "      <th>Shape_3d_Flatness</th>\n",
       "      <th>Shape_3d_LeastAxisLength</th>\n",
       "      <th>Shape_3d_MajorAxisLength</th>\n",
       "      <th>Shape_3d_Maximum2DDiameterColumn</th>\n",
       "      <th>Shape_3d_Maximum2DDiameterRow</th>\n",
       "      <th>Shape_3d_Maximum2DDiameterSlice</th>\n",
       "      <th>Shape_3d_Maximum3DDiameter</th>\n",
       "      <th>Shape_3d_MeshVolume</th>\n",
       "      <th>Shape_3d_MinorAxisLength</th>\n",
       "      <th>Shape_3d_Sphericity</th>\n",
       "      <th>Shape_3d_SurfaceArea</th>\n",
       "      <th>Shape_3d_SurfaceVolumeRatio</th>\n",
       "      <th>Shape_3d_VoxelVolume</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.923991</td>\n",
       "      <td>0.635154</td>\n",
       "      <td>21.481904</td>\n",
       "      <td>33.821559</td>\n",
       "      <td>38.832976</td>\n",
       "      <td>35.510562</td>\n",
       "      <td>37.64306</td>\n",
       "      <td>41.737274</td>\n",
       "      <td>15695.625</td>\n",
       "      <td>31.250816</td>\n",
       "      <td>0.830737</td>\n",
       "      <td>3649.274091</td>\n",
       "      <td>0.232503</td>\n",
       "      <td>15730.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Shape_3d_Elongation  Shape_3d_Flatness  Shape_3d_LeastAxisLength  \\\n",
       "0             0.923991           0.635154                 21.481904   \n",
       "\n",
       "   Shape_3d_MajorAxisLength  Shape_3d_Maximum2DDiameterColumn  \\\n",
       "0                 33.821559                         38.832976   \n",
       "\n",
       "   Shape_3d_Maximum2DDiameterRow  Shape_3d_Maximum2DDiameterSlice  \\\n",
       "0                      35.510562                         37.64306   \n",
       "\n",
       "   Shape_3d_Maximum3DDiameter  Shape_3d_MeshVolume  Shape_3d_MinorAxisLength  \\\n",
       "0                   41.737274            15695.625                 31.250816   \n",
       "\n",
       "   Shape_3d_Sphericity  Shape_3d_SurfaceArea  Shape_3d_SurfaceVolumeRatio  \\\n",
       "0             0.830737           3649.274091                     0.232503   \n",
       "\n",
       "   Shape_3d_VoxelVolume  \n",
       "0               15730.0  "
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_shape = func_shape(predmask_arr)\n",
    "df_shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Acquiring tumor loci information...\n",
      "True\n",
      "/home/baon/PycharmProjects/0_DeepLearning_glioma/hybridIDH_github_publication/mask_2mni.nii.gz\n",
      "atlasquery -a \"MNI Structural Atlas\" -m /home/baon/PycharmProjects/0_DeepLearning_glioma/hybridIDH_github_publication/mask_2mni.nii.gz\n",
      "Cerebellum:8.0845\n",
      "Occipital Lobe:0.2472\n",
      "Temporal Lobe:51.2400\n",
      "\n",
      "3\n",
      "Cerebellum:8.0845\n",
      "10\n",
      "Occipital Lobe:0.2472\n",
      "14\n",
      "Temporal Lobe:51.2400\n",
      "13\n",
      "{'mni_Temporal Lobe': 51.24, 'mni_Occipital Lobe': 0.2472, 'mni_Cerebellum': 8.0845}\n"
     ]
    }
   ],
   "source": [
    "df_loci = func_loci(T1C_iso_file, predmask_isovoxel_file)   # filepaths to T1C_isovoxel and predicted mask from Model 1 that is reconstructed back to isovoxel space\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "sla_features = pd.concat([df_shape, df_loci], axis =1)\n",
    "\n",
    "# Add patient's age\n",
    "sla_features['age'] = pd.Series(50)  ## change 50 to the patient's age.\n",
    "\n",
    "# Concatenate the predicted probability from CNN classifier (Model 2)\n",
    "sla_features['resnet_prob'] = pd.Series(resnet_prob_mean.cpu().numpy())\n",
    "\n",
    "\n",
    "sla_features.to_csv('./sla_features.csv')  ## file path to shape_loci_age feature table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Shape_3d_Elongation</th>\n",
       "      <th>Shape_3d_Flatness</th>\n",
       "      <th>Shape_3d_LeastAxisLength</th>\n",
       "      <th>Shape_3d_MajorAxisLength</th>\n",
       "      <th>Shape_3d_Maximum2DDiameterColumn</th>\n",
       "      <th>Shape_3d_Maximum2DDiameterRow</th>\n",
       "      <th>Shape_3d_Maximum2DDiameterSlice</th>\n",
       "      <th>Shape_3d_Maximum3DDiameter</th>\n",
       "      <th>Shape_3d_MeshVolume</th>\n",
       "      <th>Shape_3d_MinorAxisLength</th>\n",
       "      <th>...</th>\n",
       "      <th>mni_Cerebellum</th>\n",
       "      <th>mni_Frontal Lobe</th>\n",
       "      <th>mni_Insula</th>\n",
       "      <th>mni_Occipital Lobe</th>\n",
       "      <th>mni_Parietal Lobe</th>\n",
       "      <th>mni_Putamen</th>\n",
       "      <th>mni_Temporal Lobe</th>\n",
       "      <th>mni_Thalamus</th>\n",
       "      <th>age</th>\n",
       "      <th>resnet_prob</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.923991</td>\n",
       "      <td>0.635154</td>\n",
       "      <td>21.481904</td>\n",
       "      <td>33.821559</td>\n",
       "      <td>38.832976</td>\n",
       "      <td>35.510562</td>\n",
       "      <td>37.64306</td>\n",
       "      <td>41.737274</td>\n",
       "      <td>15695.625</td>\n",
       "      <td>31.250816</td>\n",
       "      <td>...</td>\n",
       "      <td>8.0845</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0.2472</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>51.24</td>\n",
       "      <td>0</td>\n",
       "      <td>50</td>\n",
       "      <td>0.9999998</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>1 rows × 25 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   Shape_3d_Elongation  Shape_3d_Flatness  Shape_3d_LeastAxisLength  \\\n",
       "0             0.923991           0.635154                 21.481904   \n",
       "\n",
       "   Shape_3d_MajorAxisLength  Shape_3d_Maximum2DDiameterColumn  \\\n",
       "0                 33.821559                         38.832976   \n",
       "\n",
       "   Shape_3d_Maximum2DDiameterRow  Shape_3d_Maximum2DDiameterSlice  \\\n",
       "0                      35.510562                         37.64306   \n",
       "\n",
       "   Shape_3d_Maximum3DDiameter  Shape_3d_MeshVolume  Shape_3d_MinorAxisLength  \\\n",
       "0                   41.737274            15695.625                 31.250816   \n",
       "\n",
       "   ...  mni_Cerebellum  mni_Frontal Lobe  mni_Insula  mni_Occipital Lobe  \\\n",
       "0  ...          8.0845                 0           0              0.2472   \n",
       "\n",
       "   mni_Parietal Lobe  mni_Putamen  mni_Temporal Lobe  mni_Thalamus  age  \\\n",
       "0                  0            0              51.24             0   50   \n",
       "\n",
       "   resnet_prob  \n",
       "0    0.9999998  \n",
       "\n",
       "[1 rows x 25 columns]"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sla_features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### R for radiomics classfier and logit model -> final predicted probability of hybrid model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%load_ext rpy2.ipython\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R -i df -w 5 -h 5 --units in -r 200"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sla_features <- read.csv(\"./sla_features.csv\", header=T, sep=\",\") # file path to the table with shape+loci+age+ prob from Model2 (CNN classifier)\n",
    "radiomics_classifier <- readRDS(\"./radiomics_classifier.rds\")\n",
    "mylogit <- readRDS(\"./hybrid_logit.rds\")\n",
    "\n",
    "\n",
    "sla_prob <- predict(radiomics_classifier, newdata = sla_features, type = \"prob\")[2]\n",
    "sla_prob <- as.numeric(sla_prob)\n",
    "resnet_prob <- sla_features$resnet_prob\n",
    "table_test <- as.data.frame(cbind(sla_prob, resnet_prob))\n",
    "final_prob <- predict(mylogit, newdata=table_test, type=\"response\")[1]\n",
    "print(paste(\"final probability for IDH mutation is :\", round(final_prob,3)))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}