random_forest_classification/GIF_feature_extraction_by_path_9class.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "e9da3965-aa65-4f64-9acd-dfd7f3a0cac3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import nibabel as nib\n",
    "import os\n",
    "import numpy as np\n",
    "from glob import glob\n",
    "from natsort import natsorted\n",
    "import skimage\n",
    "from skimage import segmentation\n",
    "from skimage import measure\n",
    "from scipy.spatial import distance\n",
    "import csv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "72fc72a5-4c34-4a44-8e33-98564b78190e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# define paths\n",
    "\n",
    "# paths to original GIF output\n",
    "parc_paths = natsorted(glob(os.path.join(\"parcellation_data\", \"GIF_with_tumour_9classes\", \"vs_gk*\", \"vs_gk_seg_GIF_multiclass_refT1.nii.gz\")))\n",
    "output_foldername = os.path.join(\"random_forest_classification\", \"features\", \"extracted_from_GIF_output\")\n",
    "\n",
    "# # paths to nnU-Net output trained with T1\n",
    "# parc_paths = natsorted(glob(os.path.join(\"parcellation_data\", \"GIF_UNet_prediction_9_classes_T1\", \"vs_gk*\", \"vs_gk_seg_GIF_multiclass_mergedCerebellum_refT1.nii.gz\")))\n",
    "# output_foldername = os.path.join(\"random_forest_classification\", \"features\", \"extracted_from_nnUNet_prediction\", \"T1\")\n",
    "\n",
    "# # paths to nnU-Net output trained with T2\n",
    "# parc_paths = natsorted(glob(os.path.join(\"parcellation_data\", \"GIF_UNet_prediction_9_classes_T2\", \"vs_gk*\", \"vs_gk_seg_GIF_multiclass_mergedCerebellum_refT1.nii.gz\")))\n",
    "# output_foldername = os.path.join(\"random_forest_classification\", \"features\", \"extracted_from_nnUNet_prediction\", \"T2\")\n",
    "\n",
    "# # paths to nnU-Net output trained with T1T2\n",
    "# parc_paths = natsorted(glob(os.path.join(\"parcellation_data\", \"GIF_UNet_prediction_9_classes_T1T2\", \"vs_gk*\", \"vs_gk_seg_GIF_multiclass_mergedCerebellum_refT1.nii.gz\")))\n",
    "# output_foldername = os.path.join(\"random_forest_classification\", \"features\", \"extracted_from_nnUNet_prediction\", \"T1T2\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "bed3f109-d399-44cd-8641-23e94d4ca48b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# name the labels that are present in the segmentations\n",
    "label_dict = {\n",
    " '0': 'Background',\n",
    " '1':'Tumour',\n",
    " '2':'Pons',\n",
    " '3':'Brainstem',\n",
    " '4':'Cerebellar Vermal Lobules I-V',\n",
    " '5':'Cerebellar Vermal Lobules VI-VII',\n",
    " '6':'Cerebellar Vermal Lobules VIII-X',\n",
    " '7':'Right Cerebellum',\n",
    " '8':'Left Cerebellum',\n",
    " }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "1d8903f3-8db5-4abc-93b3-07e72795e1ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "# function to get the neighbouring voxels of voxel with coordinates x,y,z\n",
    "def get_neighbors(x, y, z):\n",
    "    neighbors = []\n",
    "    for i,j,k in [(1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1)]: \n",
    "        neighbors.append([x-i, y-j, z-k])\n",
    "        \n",
    "    return neighbors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "37dc943e-488a-49e9-98e6-5c36349b2e91",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calculate smallest distance btween tumour and parcellations\n",
      "Iterate over parcels and calculate minimum distance to VS for each parcel...\n",
      "current label: 0\n",
      "current label: 1\n",
      "current label: 2\n",
      "current label: 3\n",
      "current label: 4\n",
      "current label: 5\n",
      "current label: 6\n",
      "current label: 7\n",
      "current label: 8\n",
      "Features case 8:\n",
      "label\tname                                               vol                                                surfVS                                             distVS                                             \n",
      "0\tBackground                                         3378828.567123413                                  646.3290405273438                                  0.41015625                                         \n",
      "1\tTumour                                             483.7400436401367                                  0                                                  0.0                                                \n",
      "2\tPons                                               8614.96353149414                                   0                                                  5.292426748079662                                  \n",
      "3\tBrainstem                                          7608.622741699219                                  0                                                  4.676500767008183                                  \n",
      "4\tCerebellar Vermal Lobules I-V                      3517.398262023926                                  0                                                  18.953334722278466                                 \n",
      "5\tCerebellar Vermal Lobules VI-VII                   1521.8759536743164                                 0                                                  30.116084002157386                                 \n",
      "6\tCerebellar Vermal Lobules VIII-X                   2709.6508026123047                                 0                                                  19.850777314386246                                 \n",
      "7\tRight Cerebellum                                   60328.71723175049                                  0                                                  21.847798377532516                                 \n",
      "8\tLeft Cerebellum                                    60113.469314575195                                 142.70883178710938                                 0.41015625                                         \n"
     ]
    }
   ],
   "source": [
    "for parc_path in parc_paths:\n",
    "\n",
    "    parc_nii = nib.load(parc_path)\n",
    "    parc_data = np.array(parc_nii.get_fdata(), dtype=int)\n",
    "    spacing = parc_nii.header['pixdim'][1:4]\n",
    "\n",
    "    # get unique values\n",
    "    label_number_list = np.unique(parc_data)\n",
    "\n",
    "    tumour_label_number = 1\n",
    "\n",
    "    ## get boundary pixels of tumour and set the value to the same label as neighbor\n",
    "    border_pixels_img = np.ones_like(parc_data)*(-1)\n",
    "\n",
    "    tumour_voxel_inds = np.where(parc_data == tumour_label_number)\n",
    "    for x, y, z in zip(*tumour_voxel_inds):\n",
    "\n",
    "        if parc_data[x,y,z] == tumour_label_number:\n",
    "            neighbors = get_neighbors(x,y,z)\n",
    "            for nei in neighbors:\n",
    "                n = int(parc_data[tuple([[c] for c in nei])])\n",
    "                if n !=tumour_label_number:\n",
    "                    border_pixels_img[x,y,z] = n\n",
    "\n",
    "\n",
    "    ## calculate surface between tumour and each other tissue\n",
    "    tumour_surface_voxel_inds = np.where(border_pixels_img != -1)\n",
    "    present_neighbour_labels = np.unique(border_pixels_img)\n",
    "\n",
    "    print(\"Calculate smallest distance btween tumour and parcellations\")\n",
    "    surface_areas = {n:0 for n in present_neighbour_labels}\n",
    "    for n in present_neighbour_labels:\n",
    "        binary_3d_surf_data = np.zeros_like(border_pixels_img)\n",
    "        binary_3d_surf_data[border_pixels_img==n] = 1\n",
    "\n",
    "        # Use marching cubes to obtain the surface mesh\n",
    "        verts, faces, normals, values = measure.marching_cubes(binary_3d_surf_data, level=0, spacing=spacing)\n",
    "\n",
    "        surface_areas[n] = skimage.measure.mesh_surface_area(verts, faces)\n",
    "\n",
    "\n",
    "    ## calculate smallest distance btween tumour and parcellations\n",
    "    dists = {}\n",
    "\n",
    "    tumour_coords = list(zip(*tumour_voxel_inds))\n",
    "\n",
    "    print(\"Iterate over parcels and calculate minimum distance to VS for each parcel...\")\n",
    "    for n in label_number_list:\n",
    "        print(\"current label:\", n)\n",
    "        # extract only the voxels on the surface of the nth parcellation\n",
    "        single_parc_binary_img = parc_data==n\n",
    "        boundaries = skimage.segmentation.find_boundaries(single_parc_binary_img, connectivity=1, mode='inner', background=0)\n",
    "\n",
    "        single_parc_surf_voxel_inds = list(zip(*np.where(boundaries)))\n",
    "\n",
    "        # for parcels with many voxels, the resulting distance matrix of scipy.spatial.distance.cdist does not fit in memory\n",
    "        # therefore, define the maximum matrix size, split up larger matrices and iterate over them\n",
    "        max_matrix_elem_in_memory = 776421036  # maximum matrix size\n",
    "        steps=int(np.ceil(len(tumour_coords)*len(single_parc_surf_voxel_inds)/max_matrix_elem_in_memory))  # number of parts into which parcel is split\n",
    "        max_elem_in_parc_vector = int(np.ceil(len(single_parc_surf_voxel_inds)/(steps)))  # length of one sub-parcel\n",
    "\n",
    "        dist = np.Inf\n",
    "        for s in range(steps): # iterate over sub-parcels and save only the minimum distance\n",
    "            part_of_single_parc_surf_voxel_inds = single_parc_surf_voxel_inds[s*max_elem_in_parc_vector:(s+1)*max_elem_in_parc_vector]\n",
    "            min_dist = np.min(distance.cdist(np.multiply(tumour_coords, spacing), np.multiply(part_of_single_parc_surf_voxel_inds, spacing), 'euclidean'))\n",
    "            dist = min_dist if min_dist < dist else dist\n",
    "        dists[n] = dist\n",
    "\n",
    "\n",
    "    ## calculate volumes of all parcels\n",
    "    vols = {}\n",
    "    for n in label_number_list:\n",
    "        nb_voxels = np.count_nonzero(parc_data==n)\n",
    "        vols[n] = nb_voxels*np.prod(spacing)\n",
    "\n",
    "    ## collect all features\n",
    "    features = {n:{'name':'0', 'vol':0, 'surfVS':0, 'distVS':np.nan} for n in label_number_list}\n",
    "    for n in label_number_list:\n",
    "\n",
    "        features[n]['vol'] = vols[n]\n",
    "        features[n]['name'] = label_dict[str(n)]\n",
    "        if n in surface_areas:\n",
    "            features[n]['surfVS'] = surface_areas[n]\n",
    "        if n in dists:\n",
    "            features[n]['distVS'] = dists[n]\n",
    "\n",
    "    ## print all features\n",
    "    print(f\"Features case {n}:\")\n",
    "    col_width = 50\n",
    "    # Print the names of the columns.\n",
    "    print (\"label\\t\"+(\"{:<{col_width}} \"*len(list(features.values())[0].keys())).format(*list(features.values())[0].keys(), col_width=col_width))\n",
    "\n",
    "    # print each data item.\n",
    "    for key, value in features.items():\n",
    "        #name, vol, surfVS = *value.values()\n",
    "        print (str(key)+\"\\t\"+(\"{:<{col_width}} \"*len(value.keys())).format(*value.values(), col_width=col_width))\n",
    "\n",
    "    ##save features as csv\n",
    "    for key in features:\n",
    "        features[key][\"labelnumber\"] = key\n",
    "\n",
    "    csv_columns = list(list(features.values())[0].keys())\n",
    "    case_name = parc_path.split(os.sep)[-2]\n",
    "    csv_file = os.path.join(output_foldername, case_name+\".csv\")\n",
    "\n",
    "    if not os.path.isdir(os.path.dirname(csv_file)):\n",
    "        os.makedirs(os.path.dirname(csv_file), exist_ok=True)\n",
    "\n",
    "    try:\n",
    "        with open(csv_file, 'w') as csvfile:\n",
    "            writer = csv.DictWriter(csvfile, fieldnames=[\"labelnumber\", \"name\", \"vol\", \"surfVS\", \"distVS\"])\n",
    "            writer.writeheader()\n",
    "            for data in list(features.values()):\n",
    "                writer.writerow(data)\n",
    "    except IOError:\n",
    "        print(\"I/O error\")"
   ]
  }
 ],
 "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}