Test_Dataset.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# matplotlib plots within notebook\n",
    "%matplotlib inline\n",
    "\n",
    "import platform\n",
    "print(\"python: \"+platform.python_version())\n",
    "\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import time\n",
    "\n",
    "import os, shutil, sys\n",
    "\n",
    "import pydicom as dicom\n",
    "\n",
    "from scipy import ndimage\n",
    "from skimage import morphology\n",
    "from skimage import measure\n",
    "\n",
    "import pandas as pd\n",
    "\n",
    "sys.path.insert(0, '.')\n",
    "from DeepAttCorr_lib import GAN_3D_lib as GAN\n",
    "from DeepAttCorr_lib import data_handling as DH\n",
    "from DeepAttCorr_lib import file_manage_utils as File_mng\n",
    "\n",
    "import tensorflow as tf\n",
    "import tensorflow.keras as keras\n",
    "import tensorflow.keras.backend as K\n",
    "print('Using TensorFlow version: '+tf.__version__)\n",
    "\n",
    "from tensorflow.python.client import device_lib\n",
    "print(device_lib.list_local_devices())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Paths and Definitions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Output path of generated dataframe\n",
    "OUTPUT_PATH = './test_outputs'\n",
    "\n",
    "# --------------------------------------------------------------\n",
    "# -------------------- Choose the model ------------------------\n",
    "# --------------------------------------------------------------\n",
    "\n",
    "# -------------- Sucessfull topologies\n",
    "\n",
    "#-- 3D U-Net (Sup. Pre-Train No GAN.Loss.)\n",
    "TEST_NAME = '3D_UNet'\n",
    "CHECKPOINT_PATH = './trained_models/Sup_Pre_Train_No_GAN_Loss/'\n",
    "CHECKPOINT_NAME = 'Sup_Pre_Train_No_GAN_Loss'\n",
    "\n",
    "# #-- Sup. Pre-Train RestrictedGrad.\n",
    "# TEST_NAME = 'Sup_Pre_Train_RestrictedGrad'\n",
    "# CHECKPOINT_PATH = './trained_models/Sup_Pre_Train_RestrictedGrad'\n",
    "# CHECKPOINT_NAME = 'Sup_Pre_Train_RestrictedGrad_generator_model'\n",
    "\n",
    "# -------------- Other tests\n",
    "\n",
    "# #-- Baseline\n",
    "# TEST_NAME = 'Baseline'\n",
    "# CHECKPOINT_PATH = './trained_models/Baseline'\n",
    "# CHECKPOINT_NAME = 'Baseline_generator_model'\n",
    "\n",
    "# #-- Sup. Pre-Train\n",
    "# TEST_NAME = 'Sup_Pre_Train'\n",
    "# CHECKPOINT_PATH = './trained_models/Sup_Pre_Train'\n",
    "# CHECKPOINT_NAME = 'Sup_Pre_Train_generator_model'\n",
    "\n",
    "# #-- Sup. No Pre-Train\n",
    "# TEST_NAME = 'Sup_Not_Pre_Train'\n",
    "# CHECKPOINT_PATH = './trained_models/Sup_Not_Pre_Train'\n",
    "# CHECKPOINT_NAME = 'Sup_Not_Pre_Train_generator_model'\n",
    "\n",
    "# #-- Progressive Growing Pix2Pix GAN\n",
    "# TEST_NAME = 'ProGAN_Pix2Pix_3D'\n",
    "# CHECKPOINT_PATH = './trained_models/ProGAN_Pix2Pix_3D'\n",
    "# CHECKPOINT_NAME = 'ProGAN_Pix2Pix_3D_generator_model'\n",
    "\n",
    "\n",
    "# --------------------------------------------------------------\n",
    "# -------------------- Choose the dataset ----------------------\n",
    "# --------------------------------------------------------------\n",
    "\n",
    "# DATASET_NAME = 'HNSCC_Validation_Samples'\n",
    "DATASET_NAME = 'NSCLC Radiogenomics'\n",
    "# DATASET_NAME = 'TCGA-HNSC'\n",
    "# DATASET_NAME = 'TCGA-LUAD'\n",
    "# DATASET_NAME = 'TCGA-THCA'\n",
    "# DATASET_NAME = 'CPTAC-LUAD'\n",
    "# DATASET_NAME = 'CPTAC-PDA'\n",
    "# DATASET_NAME = 'CPTAC-UCEC'\n",
    "# DATASET_NAME = 'CPTAC-LSCC'\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Input data info"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This data is from the \"Generate_train_and_validation_dataset\" notebook:\n",
    "# Objective size of the volume\n",
    "Objective_size = [256, 256, 512]\n",
    "# Voxel size used in the dataset (and hence the NN)\n",
    "Voxel_size_NORM = [2.734375, 2.734375, 2.5473046004770703]\n",
    "\n",
    "# Set different thresholds for the segmentation of the CT images\n",
    "# Thresholds are in Hounsfield Units\n",
    "HU_Air_limit = -1000\n",
    "HU_Lung_upper_limit = -500\n",
    "HU_Fluids_Fat_lower_limit = -125\n",
    "HU_Fluids_Fat_upper_limit = 10.0\n",
    "HU_Parenchyma_upper_limit = 90.0\n",
    "HU_Parenchyma_lower_limit = 10.0\n",
    "HU_Bone_low_limit = 90.0\n",
    "HU_Bone_upper_limit = 1300\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Set-up"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Build paths\n",
    "DATASET_PATH = os.path.join('./datasets/DICOM/',DATASET_NAME)\n",
    "OUTPUT_PATH = os.path.join(OUTPUT_PATH, TEST_NAME)\n",
    "OUTPUT_PATH = os.path.join(OUTPUT_PATH, DATASET_NAME)\n",
    "\n",
    "File_mng.check_create_path('OUTPUT_PATH', OUTPUT_PATH, clear_folder=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "custom_layers_dict = {'HeScale': GAN.layers.HeScale,\n",
    "                      'BiasLayer': GAN.layers.BiasLayer,\n",
    "                      'PixelNormalization': GAN.layers.PixelNormalization}\n",
    "synth_model = GAN.train_support.load_model(CHECKPOINT_PATH, \n",
    "                                           CHECKPOINT_NAME, \n",
    "                                           custom_obj_dict = custom_layers_dict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set input volume size\n",
    "input_shape_list = synth_model.input.shape.as_list()\n",
    "\n",
    "voxels_X = input_shape_list[1]\n",
    "voxels_Y = input_shape_list[2]\n",
    "voxels_Z = input_shape_list[3]\n",
    "input_size = (voxels_X,voxels_Y,voxels_Z)\n",
    "\n",
    "# Check if it is a segmentator model also\n",
    "output_shape_list = synth_model.output.shape.as_list()\n",
    "SEGM_NET = False\n",
    "if output_shape_list[-1] > 1:\n",
    "    print('Segmentation output found.')\n",
    "    SEGM_NET = True\n",
    "    \n",
    "    \n",
    "# Data full size\n",
    "voxels_X_Full = int(input_shape_list[1])\n",
    "voxels_Y_Full = int(input_shape_list[2])\n",
    "voxels_Z_Full = int(Objective_size[2]*(voxels_X_Full/Objective_size[0]))\n",
    "input_size_Full = (voxels_X_Full,voxels_Y_Full,voxels_Z_Full)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "synth_model.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "tf.keras.utils.plot_model(synth_model, \n",
    "                          to_file=os.path.join(OUTPUT_PATH,TEST_NAME+'.png'), \n",
    "                          show_shapes=True, \n",
    "                          show_layer_names=False,\n",
    "                          rankdir='TB', \n",
    "                          expand_nested=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load and Preprocess Samples\n",
    "\n",
    "Here all the DICOM files in the dataset will be loaded and preprocessed to be tested with the model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Build Dicom Paths"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "PathDicom_list = list()\n",
    "DicomNameList = list()\n",
    "\n",
    "list_paths = sorted(os.listdir(DATASET_PATH))\n",
    "\n",
    "# Explore each series\n",
    "for idx, series_name in enumerate(list_paths):\n",
    "    \n",
    "    print(series_name)\n",
    "    DicomNameList.append(series_name)\n",
    "    \n",
    "    path_series = os.path.join(DATASET_PATH,series_name)\n",
    "    \n",
    "    # Explore each study\n",
    "    for study_name in os.listdir(path_series):\n",
    "        print('\\t%s'%study_name)\n",
    "        path_study = os.path.join(path_series,study_name)\n",
    "        PathDicom_list.append(path_study)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reserve memory\n",
    "num_test_samples = len(PathDicom_list)\n",
    "\n",
    "NAC_PET_images = np.zeros((num_test_samples,voxels_X_Full,voxels_Y_Full,voxels_Z_Full))\n",
    "CT_OBJ_images = np.zeros((num_test_samples,voxels_X_Full,voxels_Y_Full,voxels_Z_Full))\n",
    "test_scale_limits_list = list()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load DICOM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Limits margin\n",
    "HU_Lung_upper_limit = HU_Lung_upper_limit\n",
    "HU_Bone_upper_limit = HU_Bone_upper_limit\n",
    "HU_fat_lower_limit = HU_Fluids_Fat_lower_limit * 1.2\n",
    "\n",
    "img_count = 0\n",
    "for PathDicom_test in PathDicom_list:\n",
    "\n",
    "    print('(%d/%d) Checking path: %s'%(img_count+1,num_test_samples,PathDicom_test))\n",
    "\n",
    "    COPY_NAC_OK = False\n",
    "    COPY_CT_OK = False\n",
    "\n",
    "\n",
    "    for dirName, subdirList, fileList in os.walk(PathDicom_test):\n",
    "\n",
    "        # Scan for all image modalities and get limits\n",
    "        image_path_list = list()\n",
    "        nac_pet_limits = [0,0]\n",
    "        ct_limits = [0,0]\n",
    "        for subDirName in subdirList:\n",
    "            # Build path\n",
    "            PathDicom = os.path.join(dirName, subDirName)\n",
    "            list_paths = sorted(os.listdir(PathDicom))\n",
    "\n",
    "            test_slices = len(list_paths)\n",
    "\n",
    "            # Read all slice locations\n",
    "            test_locations = np.zeros((test_slices))\n",
    "            for idx, filename in enumerate(list_paths):\n",
    "                if \".dcm\" in filename.lower():  # check whether the file's DICOM\n",
    "                    RefDs = dicom.read_file(os.path.join(PathDicom,filename))\n",
    "                    test_locations[idx] = np.array(RefDs.ImagePositionPatient)[2]\n",
    "\n",
    "            # Save extrema\n",
    "            if RefDs.Modality == 'PT':\n",
    "                if (not 'ATTN' in RefDs.CorrectedImage) or 'NAC' in PathDicom or 'uncorrected' in PathDicom or 'NOAC' in PathDicom:\n",
    "                    image_path_list.append(PathDicom)\n",
    "                    nac_pet_limits[0] = test_locations.min()\n",
    "                    nac_pet_limits[1] = test_locations.max()\n",
    "\n",
    "                else:\n",
    "                    image_path_list.append(PathDicom)\n",
    "                    continue\n",
    "\n",
    "            elif RefDs.Modality == 'CT':\n",
    "                image_path_list.append(PathDicom)\n",
    "                ct_limits[0] = test_locations.min()\n",
    "                ct_limits[1] = test_locations.max()\n",
    "\n",
    "        # Get absolute limits\n",
    "        loc_min = np.max([ct_limits[0],nac_pet_limits[0]])\n",
    "        loc_max = np.min([ct_limits[1],nac_pet_limits[1]])\n",
    "\n",
    "        for PathDicom in image_path_list:\n",
    "\n",
    "            list_paths = sorted(os.listdir(PathDicom))\n",
    "            RefDs = dicom.read_file(os.path.join(PathDicom,list_paths[0]))\n",
    "\n",
    "            IS_CT = False\n",
    "            IS_PT = False\n",
    "            IS_NAC_PT = False\n",
    "\n",
    "            # Check image type\n",
    "            if RefDs.Modality == 'PT':\n",
    "                if (not 'ATTN' in RefDs.CorrectedImage) or 'NAC' in PathDicom or 'uncorrected' in PathDicom or 'NOAC' in PathDicom:\n",
    "                    IS_NAC_PT = True\n",
    "                else:\n",
    "                    IS_PT = True\n",
    "            elif RefDs.Modality == 'CT':\n",
    "                IS_CT = True                \n",
    "\n",
    "            if IS_PT:\n",
    "                print('PET image: %d files'%(len(list_paths)))\n",
    "                continue\n",
    "            elif IS_NAC_PT:\n",
    "                print('NAC PET image: %d files'%(len(list_paths)))\n",
    "            elif IS_CT:\n",
    "                print('CT image: %d files'%(len(list_paths)))\n",
    "\n",
    "\n",
    "\n",
    "            test_cols = RefDs.Columns\n",
    "            test_rows = RefDs.Rows\n",
    "            test_slices = len(list_paths)\n",
    "            test_x_size = RefDs.PixelSpacing[0]\n",
    "            test_y_size = RefDs.PixelSpacing[1]\n",
    "            test_FOV_X = test_cols*test_x_size\n",
    "            test_FOV_Y = test_rows*test_y_size\n",
    "            test_FOV_X_min = -test_FOV_X/2.0\n",
    "            test_FOV_X_max = test_FOV_X/2.0\n",
    "            test_FOV_Y_min = -test_FOV_Y/2.0\n",
    "            test_FOV_Y_max = test_FOV_Y/2.0\n",
    "\n",
    "\n",
    "            test_input = np.zeros((test_cols,test_rows,test_slices), dtype=np.float32)\n",
    "            test_locations = np.zeros((test_slices))\n",
    "\n",
    "\n",
    "            for idx, filename in enumerate(list_paths):\n",
    "                if \".dcm\" in filename.lower():  # check whether the file's DICOM\n",
    "                    RefDs = dicom.read_file(os.path.join(PathDicom,filename))\n",
    "                    test_input_aux = ((RefDs.pixel_array*RefDs.RescaleSlope)+RefDs.RescaleIntercept) \n",
    "                    test_input[:,:,idx] = np.reshape(test_input_aux,(test_cols,test_rows))\n",
    "                    test_locations[idx] = np.array(RefDs.ImagePositionPatient)[2]\n",
    "\n",
    "            # re-order\n",
    "            order_idx = np.argsort(test_locations)\n",
    "            test_input = test_input[:,:,order_idx]\n",
    "\n",
    "            test_z_size = np.abs(loc_max-loc_min)/float(len(list_paths))\n",
    "            test_FOV_Z = test_slices*test_z_size\n",
    "\n",
    "\n",
    "\n",
    "            Volume_size_INPUT = [test_cols, test_rows, test_slices]\n",
    "            Voxel_size_INPUT = [test_x_size, test_y_size, test_z_size]\n",
    "\n",
    "            image_position_patient = np.array(RefDs.ImagePositionPatient)\n",
    "\n",
    "\n",
    "            if IS_CT:\n",
    "                # Remove couch\n",
    "                HU_Air_limit = -1000\n",
    "                B_n_d = 4\n",
    "                B_n_e = 2\n",
    "                B_d_d = 15\n",
    "                B_d_e = 13\n",
    "                limits_ct_use = [0,0]\n",
    "                corr_disk = 6\n",
    "\n",
    "                couch_mask = DH.CT_couch_removal_mask_multicore(test_input, \n",
    "                                                                   samples_use = 20, \n",
    "                                                                   n_corr_disk = corr_disk, \n",
    "                                                                   limits_CT = limits_ct_use,\n",
    "                                                                   B_n_dilatation_disk_rad = B_n_d,\n",
    "                                                                   B_n_erosion_disk_rad = B_n_e,\n",
    "                                                                   B_d_dilatation_disk_rad = B_d_d,\n",
    "                                                                   B_d_erosion_disk_rad = B_d_e) \n",
    "\n",
    "\n",
    "                for z_idx in range(0,test_input.shape[2]):\n",
    "                    test_input[:,:,z_idx][couch_mask] = HU_Air_limit               \n",
    "\n",
    "                # Change units from Housefield to mu\n",
    "                test_input = ((test_input/1000.0)+1)*DH.mu_agua_120kev/10.0\n",
    "\n",
    "\n",
    "            # Normalize Size\n",
    "            test_input_norm, test_limits = DH.normalize_volume_single(test_input,\n",
    "                                                                         Voxel_size_INPUT, \n",
    "                                                                         Objective_size, \n",
    "                                                                         Voxel_size_NORM,\n",
    "                                                                        order_interpol=1)\n",
    "            post_scaling = 2.0\n",
    "            test_input_norm_scale = np.zeros((int(test_input_norm.shape[0]/int(post_scaling)),\n",
    "                                            int(test_input_norm.shape[1]/int(post_scaling)),\n",
    "                                                int(test_input_norm.shape[2]/int(post_scaling))))\n",
    "\n",
    "            test_input_norm_scale = ndimage.interpolation.zoom(input=test_input_norm, zoom=[1/post_scaling,1/post_scaling,1/post_scaling], order=2)\n",
    "\n",
    "            input_net_size = [int(Objective_size[0]/post_scaling),int(Objective_size[1]/post_scaling),int(Objective_size[2]/post_scaling)]\n",
    "\n",
    "            test_scale_limits = [int(test_limits[0]/2), int(test_limits[1]/2)]\n",
    "\n",
    "            if IS_CT:\n",
    "                LIM_SUP_CT = (((HU_Bone_upper_limit/1000.0)+1)*DH.mu_agua_120kev/10.0) \n",
    "                LIM_INF_CT = (((HU_fat_lower_limit/1000.0)+1)*DH.mu_agua_120kev/10.0) \n",
    "\n",
    "                test_input_norm_scale[test_input_norm_scale<LIM_INF_CT]=LIM_INF_CT\n",
    "                test_input_norm_scale[test_input_norm_scale>LIM_SUP_CT]=LIM_SUP_CT\n",
    "\n",
    "            test_input_norm_scale = (test_input_norm_scale - test_input_norm_scale.min()) / (test_input_norm_scale.max() - test_input_norm_scale.min())\n",
    "\n",
    "            # Save\n",
    "            if IS_NAC_PT:\n",
    "                NAC_PET_images[img_count,:,:,:] = np.copy(test_input_norm_scale)\n",
    "                test_scale_limits_list.append(test_scale_limits)\n",
    "                COPY_NAC_OK = True\n",
    "\n",
    "            if IS_CT:\n",
    "\n",
    "                CT_OBJ_images[img_count,:,:,:] = np.copy(test_input_norm_scale)\n",
    "                COPY_CT_OK = True\n",
    "\n",
    "    if not COPY_NAC_OK:\n",
    "        print('Missing NAC PET scan. Aborting')\n",
    "    if not COPY_CT_OK:\n",
    "        print('Missing CT scan. Aborting')\n",
    "\n",
    "    if COPY_NAC_OK and COPY_CT_OK:\n",
    "\n",
    "        print('Image limits [ %d ; %d ]'%(test_scale_limits_list[img_count][0],\n",
    "                                         test_scale_limits_list[img_count][1]))\n",
    "\n",
    "        slice_X = 64\n",
    "        slice_Y = 64\n",
    "        slice_Z = 128\n",
    "\n",
    "        fig, ax = plt.subplots(nrows=2, ncols=3, dpi=50, sharex=True, sharey=True)\n",
    "        ax[0][0].imshow(NAC_PET_images[img_count,slice_X,:,:], aspect='equal', vmin=0.0, vmax=1.0)\n",
    "        ax[0][0].axis('off')\n",
    "        ax[0][1].imshow(NAC_PET_images[img_count,:,slice_Y,:], aspect='equal', vmin=0.0, vmax=1.0)\n",
    "        ax[0][1].axis('off')\n",
    "        ax[0][2].imshow(NAC_PET_images[img_count,:,:,slice_Z], aspect='equal', vmin=0.0, vmax=1.0)\n",
    "        ax[0][2].axis('off')\n",
    "        ax[1][0].imshow(CT_OBJ_images[img_count,slice_X,:,:], aspect='equal', vmin=0.0, vmax=1.0)\n",
    "        ax[1][0].axis('off')\n",
    "        ax[1][1].imshow(CT_OBJ_images[img_count,:,slice_Y,:], aspect='equal', vmin=0.0, vmax=1.0)\n",
    "        ax[1][1].axis('off')\n",
    "        ax[1][2].imshow(CT_OBJ_images[img_count,:,:,slice_Z], aspect='equal', vmin=0.0, vmax=1.0)\n",
    "        ax[1][2].axis('off')\n",
    "        plt.tight_layout()\n",
    "        plt.plot()\n",
    "        plt.show()\n",
    "\n",
    "        img_count += 1\n",
    "\n",
    "# hold\n",
    "num_test_samples = img_count\n",
    "NAC_PET_images = NAC_PET_images[:num_test_samples,:,:,:]\n",
    "CT_OBJ_images = CT_OBJ_images[:num_test_samples,:,:,:]\n",
    "test_scale_limits_list = test_scale_limits_list[:num_test_samples]\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "slice_X = 64\n",
    "slice_Y = 64\n",
    "slice_Z = 128\n",
    "img_count = 0\n",
    "fig, ax = plt.subplots(nrows=2, ncols=3, dpi=150, sharex=True, sharey=True)\n",
    "ax[0][0].imshow(NAC_PET_images[img_count,slice_X,:,:], aspect='equal', vmin=0.0, vmax=1.0)\n",
    "ax[0][0].axis('off')\n",
    "ax[0][1].imshow(NAC_PET_images[img_count,:,slice_Y,:], aspect='equal', vmin=0.0, vmax=1.0)\n",
    "ax[0][1].axis('off')\n",
    "ax[0][2].imshow(NAC_PET_images[img_count,:,:,slice_Z], aspect='equal', vmin=0.0, vmax=1.0)\n",
    "ax[0][2].axis('off')\n",
    "ax[1][0].imshow(CT_OBJ_images[img_count,slice_X,:,:], aspect='equal', vmin=0.0, vmax=1.0)\n",
    "ax[1][0].axis('off')\n",
    "ax[1][1].imshow(CT_OBJ_images[img_count,:,slice_Y,:], aspect='equal', vmin=0.0, vmax=1.0)\n",
    "ax[1][1].axis('off')\n",
    "ax[1][2].imshow(CT_OBJ_images[img_count,:,:,slice_Z], aspect='equal', vmin=0.0, vmax=1.0)\n",
    "ax[1][2].axis('off')\n",
    "plt.tight_layout()\n",
    "plt.plot()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compute Test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute whole volumes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "CT_SYNTH_images = np.zeros((num_test_samples,voxels_X_Full,voxels_Y_Full,voxels_Z_Full))\n",
    "\n",
    "SEGMENTED_images = np.zeros((num_test_samples,voxels_X_Full,voxels_Y_Full,voxels_Z_Full))\n",
    "\n",
    "MARGIN = 3\n",
    "data_size = [128,128,256]\n",
    "\n",
    "for idx_img in range(num_test_samples):\n",
    "\n",
    "\n",
    "    out_list = GAN.train_support.compute_whole_volume(NAC_PET_images[idx_img,:,:,:].astype(np.float32), \n",
    "                                                          input_size, \n",
    "                                                          synth_model, \n",
    "                                                          np.copy(test_scale_limits_list[idx_img]), \n",
    "                                                          MARGIN, \n",
    "                                                          sample_num=idx_img, \n",
    "                                                          segm_net=SEGM_NET, \n",
    "                                                          s_model=synth_model, \n",
    "                                                          criticize=False, \n",
    "                                                          d_model='', \n",
    "                                                          slice_val_cap_use=1.0)\n",
    "\n",
    "    if SEGM_NET:\n",
    "        CT_synth_aux = out_list[0]\n",
    "        segmentation_aux = out_list[1]\n",
    "    else:\n",
    "        CT_synth_aux = out_list\n",
    "    # Save\n",
    "    CT_SYNTH_images[idx_img,:,:,:] = np.copy(CT_synth_aux)\n",
    "    if SEGM_NET:\n",
    "        aux_segm = np.argmax(np.copy(segmentation_aux), axis = -1)\n",
    "        aux_segm[:,:,:test_scale_limits_list[idx_img][0]]=3\n",
    "        aux_segm[:,:,test_scale_limits_list[idx_img][1]:]=3\n",
    "        SEGMENTED_images[idx_img,:,:,:] = aux_segm\n",
    "\n",
    "\n",
    "# Save synthetic CT samples\n",
    "np.save(os.path.join(OUTPUT_PATH,'CT_SYNTH_images'), CT_SYNTH_images)\n",
    "np.save(os.path.join(OUTPUT_PATH,'SEGMENTED_images'), SEGMENTED_images)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_out = int(data_size[0]/2)\n",
    "Y_out = int(data_size[1]/2)\n",
    "Z_out = int(data_size[2]/2)\n",
    "IMG_SHOW =  np.random.randint(NAC_PET_images.shape[0])\n",
    "print(IMG_SHOW)\n",
    "\n",
    "plt.figure(dpi=300)\n",
    "ax = plt.subplot(1,3,1)\n",
    "ax.imshow(NAC_PET_images[IMG_SHOW, :,:,Z_out])\n",
    "plt.axis('off')\n",
    "ax = plt.subplot(1,3,2)\n",
    "ax.imshow(NAC_PET_images[IMG_SHOW, :,Y_out,:])\n",
    "plt.axis('off')\n",
    "ax = plt.subplot(1,3,3)\n",
    "ax.imshow(NAC_PET_images[IMG_SHOW, X_out,:,:])\n",
    "plt.tight_layout()\n",
    "plt.axis('off')\n",
    "plt.draw()\n",
    "plt.savefig(os.path.join(OUTPUT_PATH, 'NAC_PET_input'))\n",
    "plt.plot()\n",
    "\n",
    "ct_vmax=CT_OBJ_images.max()\n",
    "ct_vmin=CT_OBJ_images.min()\n",
    "\n",
    "plt.figure(dpi=300)\n",
    "ax = plt.subplot(1,3,1)\n",
    "ax.imshow(CT_SYNTH_images[IMG_SHOW, :,:,Z_out], vmax=ct_vmax, vmin=ct_vmin)\n",
    "plt.axis('off')\n",
    "ax = plt.subplot(1,3,2)\n",
    "ax.imshow(CT_SYNTH_images[IMG_SHOW,:,Y_out,:], vmax=ct_vmax, vmin=ct_vmin)\n",
    "plt.axis('off')\n",
    "ax = plt.subplot(1,3,3)\n",
    "ax.imshow(CT_SYNTH_images[IMG_SHOW,X_out,:,:], vmax=ct_vmax, vmin=ct_vmin)\n",
    "plt.tight_layout()\n",
    "plt.axis('off')\n",
    "plt.draw()\n",
    "plt.savefig(os.path.join(OUTPUT_PATH, 'sCT'))\n",
    "plt.plot()\n",
    "\n",
    "\n",
    "plt.figure(dpi=300)\n",
    "ax = plt.subplot(1,3,1)\n",
    "ax.imshow(CT_OBJ_images[IMG_SHOW,:,:,Z_out], vmax=ct_vmax, vmin=ct_vmin)\n",
    "plt.axis('off')\n",
    "ax = plt.subplot(1,3,2)\n",
    "ax.imshow(CT_OBJ_images[IMG_SHOW,:,Y_out,:], vmax=ct_vmax, vmin=ct_vmin)\n",
    "plt.axis('off')\n",
    "ax = plt.subplot(1,3,3)\n",
    "ax.imshow(CT_OBJ_images[IMG_SHOW,X_out,:,:], vmax=ct_vmax, vmin=ct_vmin)\n",
    "plt.tight_layout()\n",
    "plt.axis('off')\n",
    "plt.draw()\n",
    "plt.savefig(os.path.join(OUTPUT_PATH, 'objective_CT'))\n",
    "plt.plot()\n",
    "\n",
    "\n",
    "if SEGM_NET:\n",
    "    plt.figure(dpi=300)\n",
    "    ax = plt.subplot(1,3,1)\n",
    "    ax.imshow(SEGMENTED_images[IMG_SHOW,:,:,Z_out])\n",
    "    plt.axis('off')\n",
    "    ax = plt.subplot(1,3,2)\n",
    "    ax.imshow(SEGMENTED_images[IMG_SHOW,:,Y_out,:])\n",
    "    plt.axis('off')\n",
    "    ax = plt.subplot(1,3,3)\n",
    "    ax.imshow(SEGMENTED_images[IMG_SHOW,X_out,:,:])\n",
    "    plt.tight_layout()\n",
    "    plt.axis('off')\n",
    "    plt.draw()\n",
    "    plt.savefig(os.path.join(OUTPUT_PATH, 'Tissue_Clases'))\n",
    "    plt.plot()\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compute Metrics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Number of slices skipped at the begining and end of the sample\n",
    "MARGEN_CALC = 16\n",
    "\n",
    "# Create dataframe\n",
    "col_names =  ['SSIN', 'MSE_Whole', 'MSE_RoI', 'PSNR_Whole', 'PSNR_RoI', 'MAE', 'MAE_HU', 'NMSE_Whole', 'NMSE_RoI', 'NCC_Whole', 'NCC_RoI']\n",
    "df_meassurements = pd.DataFrame(columns=col_names) \n",
    "\n",
    "\n",
    "# Compute mertic for each sample\n",
    "for idx_img in range(num_test_samples):\n",
    "\n",
    "    limits_act = test_scale_limits_list[idx_img]\n",
    "\n",
    "    \n",
    "    X_in = np.squeeze(CT_OBJ_images[idx_img,:,:,:])\n",
    "    X_in = np.expand_dims(X_in, axis=-1)\n",
    "    X_in = X_in[:,:,limits_act[0]+MARGEN_CALC:limits_act[1]-MARGEN_CALC,:]\n",
    "    X_in_nonZero = X_in[X_in>0.1]\n",
    "    \n",
    "    Y_in = np.squeeze(CT_SYNTH_images[idx_img,:,:,:])\n",
    "    Y_in = np.expand_dims(Y_in, axis=-1)\n",
    "    Y_in = Y_in[:,:,limits_act[0]+MARGEN_CALC:limits_act[1]-MARGEN_CALC,:]\n",
    "    Y_in_nonZero = Y_in[X_in>0.1]\n",
    "\n",
    "    # SSIM\n",
    "    ssim_aux = measure.compare_ssim(X_in, Y_in, sigma=5, multichannel=True)\n",
    "    # MSE\n",
    "    mse_whole_aux = measure.compare_mse(X_in, Y_in)\n",
    "    mse_RoI_aux = measure.compare_mse(X_in_nonZero, Y_in_nonZero)\n",
    "    # PSNR\n",
    "    psnr_whole_aux = measure.compare_psnr(X_in, Y_in)\n",
    "    psnr_RoI_aux = measure.compare_psnr(X_in_nonZero, Y_in_nonZero)\n",
    "    # MAE\n",
    "    mae_aux = np.sum(np.divide(np.abs(X_in_nonZero-Y_in_nonZero),X_in_nonZero))/(X_in_nonZero.size)\n",
    "    # NMSE\n",
    "    nmse_whole_aux = np.sqrt(np.sum((X_in-Y_in)**2))/(np.sqrt(np.sum(X_in**2)))\n",
    "    nmse_RoI_aux = np.sqrt(np.sum((X_in_nonZero-Y_in_nonZero)**2))/(np.sqrt(np.sum(X_in_nonZero**2)))\n",
    "    # NCC\n",
    "    X_in_zm = X_in-np.mean(X_in)\n",
    "    Y_in_zm = Y_in-np.mean(Y_in)\n",
    "    ncc_whole_aux = (np.sum(X_in_zm*Y_in_zm)/X_in.size) / (np.std(X_in)*np.std(Y_in))\n",
    "    \n",
    "    X_in_zm = X_in_nonZero-np.mean(X_in_nonZero)\n",
    "    Y_in_zm = Y_in_nonZero-np.mean(Y_in_nonZero)\n",
    "    ncc_RoI_aux = (np.sum(X_in_zm*Y_in_zm)/X_in_nonZero.size) / (np.std(X_in_nonZero)*np.std(Y_in_nonZero))\n",
    "    \n",
    "    \n",
    "    # Convert to HU\n",
    "    ct_out_scale = (X_in_nonZero*(LIM_SUP_CT-LIM_INF_CT))+LIM_INF_CT\n",
    "    X_in_nonZero_HU = (((ct_out_scale*10.0)/DH.mu_agua_120kev)-1)*1000.0\n",
    "\n",
    "    ct_out_scale = (Y_in_nonZero*(LIM_SUP_CT-LIM_INF_CT))+LIM_INF_CT\n",
    "    Y_in_nonZero_HU = (((ct_out_scale*10.0)/DH.mu_agua_120kev)-1)*1000.0\n",
    "\n",
    "    mae_HU_aux = np.sum((np.abs(X_in_nonZero_HU-Y_in_nonZero_HU)))/(X_in_nonZero_HU.size)\n",
    "    \n",
    "    \n",
    "    df_meassurements.loc[len(df_meassurements)] = [ssim_aux, \n",
    "                                                   mse_whole_aux, \n",
    "                                                   mse_RoI_aux, \n",
    "                                                   psnr_whole_aux, \n",
    "                                                   psnr_RoI_aux, \n",
    "                                                   mae_aux, \n",
    "                                                   mae_HU_aux,\n",
    "                                                   nmse_whole_aux, \n",
    "                                                   nmse_RoI_aux, \n",
    "                                                   ncc_whole_aux, \n",
    "                                                   ncc_RoI_aux]\n",
    "\n",
    "# Save dataframe\n",
    "df_meassurements.to_pickle(os.path.join(OUTPUT_PATH,'meassurements_dataframe.pkl'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Mean metric values:')\n",
    "print(df_meassurements.mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(dpi=200)\n",
    "\n",
    "ax = plt.subplot(1,3,1)\n",
    "df_meassurements[['PSNR_RoI']].boxplot()\n",
    "plt.ylabel('PSNR (db)')\n",
    "plt.ylim([14.0, 24.0])\n",
    "\n",
    "ax = plt.subplot(1,3,2)\n",
    "df_meassurements[['MAE_HU']].boxplot()\n",
    "plt.ylabel('MAE (HU)')\n",
    "plt.ylim([60.0, 160.0])\n",
    "\n",
    "\n",
    "ax = plt.subplot(1,3,3)\n",
    "df_meassurements[['NCC_Whole']].boxplot()\n",
    "plt.ylabel('NCC (-)')\n",
    "plt.ylim([0.45, 1.0])\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Save as DICOM images\n",
    "\n",
    "Note that the produced DICOM format is only to facilitate viewing the result with a DICOM reader program such as \"Slicer\". It should not be regarded as a valid medical image of any nature."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convert CT to HU units"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "CT_SYNTH_images_HU = np.zeros_like(CT_SYNTH_images)\n",
    "\n",
    "for idx_img in range(num_test_samples):\n",
    "    \n",
    "    # Get current volume\n",
    "    CT_synth_in = np.squeeze(CT_SYNTH_images[idx_img,:,:,:])\n",
    "    \n",
    "    # Voxel dynamic range limits\n",
    "    LIM_INF_CT = (((HU_fat_lower_limit/1000.0)+1)*DH.mu_agua_120kev/10.0)\n",
    "    LIM_SUP_CT = (((HU_Bone_upper_limit/1000.0)+1)*DH.mu_agua_120kev/10.0)\n",
    "    \n",
    "    # Change voxel scale\n",
    "    ct_out_scale = (CT_synth_in*(LIM_SUP_CT-LIM_INF_CT))+LIM_INF_CT\n",
    "    ct_HU = (((ct_out_scale*10.0)/DH.mu_agua_120kev)-1)*1000.0\n",
    "\n",
    "    CT_SYNTH_images_HU[idx_img,:,:,:] = np.copy(ct_HU)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Save as DICOM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for idx_img in range(num_test_samples):\n",
    "    \n",
    "    # Get current volume\n",
    "    CT_synth_in = np.squeeze(CT_SYNTH_images_HU[idx_img,:,:,:])\n",
    "    # Get current sample name\n",
    "    SAMPLE_NAME = DicomNameList[idx_img]\n",
    "    \n",
    "    # Set image to standard CT size\n",
    "\n",
    "    # Netwoerk voxel size\n",
    "    X_VOX_SIZE = Voxel_size_NORM[0]*2\n",
    "    Y_VOX_SIZE = Voxel_size_NORM[1]*2\n",
    "    Z_VOX_SIZE = Voxel_size_NORM[2]*2\n",
    "    # Objective CT voxel size\n",
    "    X_CT_VOX_SIZE = 1.367\n",
    "    Y_CT_VOX_SIZE = 1.367\n",
    "    Z_CT_VOX_SIZE = 3.27\n",
    "\n",
    "    X_zoom = X_VOX_SIZE/X_CT_VOX_SIZE\n",
    "    Y_zoom = Y_VOX_SIZE/Y_CT_VOX_SIZE\n",
    "    Z_zoom = Z_VOX_SIZE/Z_CT_VOX_SIZE\n",
    "\n",
    "    # Resample\n",
    "    ct_HU_resampled = np.copy(CT_synth_in)\n",
    "    ct_HU_resampled = ndimage.interpolation.zoom(input=ct_HU_resampled, zoom=[X_zoom,Y_zoom,Z_zoom], order=1)\n",
    "    \n",
    "    # SAve\n",
    "    DH.saveDICOM_sliced_CT(ct_HU_resampled, \n",
    "                           [X_CT_VOX_SIZE, Y_CT_VOX_SIZE, Z_CT_VOX_SIZE], \n",
    "                           image_position_patient,\n",
    "                           OUTPUT_PATH, \n",
    "                           SAMPLE_NAME,\n",
    "                           patient_ID_in = SAMPLE_NAME,\n",
    "                           StudyInstanceUID = ('1.3.6.1.4.1.14519.5.2.1.3320.3273.2314458152346821195388%d'%np.random.randint(low=63169983, high=93169983)),\n",
    "                           SeriesInstanceUID = ('1.3.6.1.4.1.14519.5.2.1.3320.3273.330516103388011054891%d'%np.random.randint(low=344582212, high=944582212)),\n",
    "                           TYPE_INT = np.int16, \n",
    "                           vervose = True)"
   ]
  }
 ],
 "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}