3d-resnet.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Common code and data loading"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#ONLY RUN THIS IF THE GPU SERVER IS BUSY\n",
    "#USES THE SLOWASS CPU INSTEAD\n",
    "\n",
    "import os\n",
    "os.environ[\"CUDA_DEVICE_ORDER\"] = \"PCI_BUS_ID\"   # see issue #152\n",
    "os.environ[\"CUDA_VISIBLE_DEVICES\"] = \"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Run this to choose 2nd GPU\n",
    "\n",
    "import os\n",
    "os.environ[\"CUDA_DEVICE_ORDER\"] = \"PCI_BUS_ID\"   # see issue #152\n",
    "os.environ[\"CUDA_VISIBLE_DEVICES\"] = \"1\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Using TensorFlow backend.\n"
     ]
    }
   ],
   "source": [
    "# Run this if the following error is thrown when training/testing models:\n",
    "\n",
    "# Error : Failed to get convolution algorithm. This is probably because cuDNN failed to initialize, \n",
    "# so try looking to see if a warning log message was printed above.\n",
    "\n",
    "import tensorflow as tf\n",
    "from keras.backend.tensorflow_backend import set_session\n",
    "from keras.backend.tensorflow_backend import clear_session\n",
    "from keras.backend.tensorflow_backend import get_session\n",
    "import gc\n",
    "\n",
    "# Reset Keras Session\n",
    "def reset_keras():\n",
    "    sess = get_session()\n",
    "    clear_session()\n",
    "    sess.close()\n",
    "    sess = get_session()\n",
    "\n",
    "    try:\n",
    "        del model # this is from global space - change this as you need\n",
    "    except:\n",
    "        pass\n",
    "\n",
    "    print(gc.collect()) # if it's done something you should see a number being outputted\n",
    "\n",
    "    # use the same config as you used to create the session\n",
    "    config_tf()\n",
    "\n",
    "\n",
    "def config_tf():\n",
    "    config = tf.ConfigProto()\n",
    "    config.gpu_options.allow_growth = True\n",
    "    sess = tf.Session(config=config)\n",
    "    set_session(sess)\n",
    "    \n",
    "config_tf()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "import SimpleITK as sitk\n",
    "import os, sys\n",
    "sys.path.insert(1, './Models/Resnet-3D')\n",
    "from resnet3d3 import Resnet3DBuilder\n",
    "import resnet3d3\n",
    "\n",
    "import numpy as np\n",
    "import random\n",
    "import keras\n",
    "from sklearn.model_selection import train_test_split\n",
    "import math\n",
    "from keras.callbacks import EarlyStopping\n",
    "\n",
    "#these are imported or the transferring models are not deserialised properly\n",
    "import conv4d\n",
    "from keras.regularizers import l2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "#if this is not here, the loaded notebook will not detect ensuing changes in imported numpy scripts "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/home/jzhe0882/numpydata/HeadNeckCancer/PET/HN-CHUS-048.npy /home/jzhe0882/numpydata/HeadNeckCancer/CT/HN-CHUS-048.npy /home/jzhe0882/Radiomics/HeadNeckCancer/PET/HN-CHUS-048.npy /home/jzhe0882/Radiomics/HeadNeckCancer/CT/HN-CHUS-048.npy [ 64  49 105] [9 8 7] /home/jzhe0882/numpydata/HeadNeckCancer/Mask/HN-CHUS-048.npy\n",
      "/home/jzhe0882/numpydata/HeadNeckCancer/PET/HN-HGJ-035.npy /home/jzhe0882/numpydata/HeadNeckCancer/CT/HN-HGJ-035.npy /home/jzhe0882/Radiomics/HeadNeckCancer/PET/HN-HGJ-035.npy /home/jzhe0882/Radiomics/HeadNeckCancer/CT/HN-HGJ-035.npy [70 38 61] [16 14 33] /home/jzhe0882/numpydata/HeadNeckCancer/Mask/HN-HGJ-035.npy\n",
      "/home/jzhe0882/numpydata/BreastCancer/PET/10368550.npy /home/jzhe0882/numpydata/BreastCancer/CT/10368550.npy /home/jzhe0882/Radiomics/BreastCancer/PET/10368550.npy /home/jzhe0882/Radiomics/BreastCancer/CT/10368550.npy [46 58 62] [3 3 4] /home/jzhe0882/numpydata/BreastCancer/Mask/10368550.npy\n",
      "/home/jzhe0882/numpydata/BreastCancer/PET/10359092.npy /home/jzhe0882/numpydata/BreastCancer/CT/10359092.npy /home/jzhe0882/Radiomics/BreastCancer/PET/10359092.npy /home/jzhe0882/Radiomics/BreastCancer/CT/10359092.npy [77 61 79] [4 4 7] /home/jzhe0882/numpydata/BreastCancer/Mask/10359092.npy\n"
     ]
    }
   ],
   "source": [
    "#load training data from disk\n",
    "\n",
    "pet_data, ct_data, centre_data, size_data, mask_data, ct_radiomics_data, pet_radiomics_data = {'train':{}, 'test':{}}, {'train':{}, 'test':{}}, \\\n",
    "                                                        {'train':{}, 'test':{}}, {'train':{}, 'test':{}}, {'train':{}, 'test':{}}, \\\n",
    "                                                        {'train':{}, 'test':{}}, {'train':{}, 'test':{}}\n",
    "\n",
    "for dataset in ['HeadNeckCancer', 'BreastCancer']:\n",
    "    pet_files = []\n",
    "    ct_files = []\n",
    "    centre_files = []\n",
    "    size_files = []\n",
    "    mask_files = []\n",
    "    pet_radiomics_files = []\n",
    "    ct_radiomics_files = []\n",
    "    \n",
    "    for root, dirs, files in os.walk('/home/jzhe0882/numpydata/' + dataset + '/PET'):\n",
    "        for name in files:\n",
    "            file_path = os.path.join(root, name)\n",
    "            pet_files.append(file_path)\n",
    "\n",
    "    for root, dirs, files in os.walk('/home/jzhe0882/numpydata/' + dataset + '/CT'):\n",
    "        for name in files:\n",
    "            file_path = os.path.join(root, name)\n",
    "            ct_files.append(file_path)\n",
    "\n",
    "    for root, dirs, files in os.walk('/home/jzhe0882/numpydata/' + dataset + '/MaskCentres'):\n",
    "        for name in files:\n",
    "            file_path = os.path.join(root, name)\n",
    "            centre_files.append(file_path)\n",
    "\n",
    "    for root, dirs, files in os.walk('/home/jzhe0882/numpydata/' + dataset + '/MaskSizes'):\n",
    "        for name in files:\n",
    "            file_path = os.path.join(root, name)\n",
    "            size_files.append(file_path)\n",
    "\n",
    "    for root, dirs, files in os.walk('/home/jzhe0882/numpydata/' + dataset + '/Mask'):\n",
    "        for name in files:\n",
    "            file_path = os.path.join(root, name)\n",
    "            mask_files.append(file_path)\n",
    "            \n",
    "    for root, dirs, files in os.walk('/home/jzhe0882/Radiomics/' + dataset + '/CT'):\n",
    "        for name in files:\n",
    "            file_path = os.path.join(root, name)\n",
    "            ct_radiomics_files.append(file_path)\n",
    "            \n",
    "    for root, dirs, files in os.walk('/home/jzhe0882/Radiomics/' + dataset + '/PET'):\n",
    "        for name in files:\n",
    "            file_path = os.path.join(root, name)\n",
    "            pet_radiomics_files.append(file_path)\n",
    "\n",
    "    pet_files = sorted(pet_files)\n",
    "    ct_files = sorted(ct_files)\n",
    "    centre_files = sorted(centre_files)\n",
    "    centres = [np.load(c) for c in centre_files] #can load all of these into memory (other volumes are too large)\n",
    "    size_files = sorted(size_files)\n",
    "    bounding_box_sizes = [np.load(s) for s in size_files] #can also load all of these into memory\n",
    "    mask_files = sorted(mask_files)\n",
    "    pet_radiomics_files= sorted(pet_radiomics_files)\n",
    "    ct_radiomics_files = sorted(ct_radiomics_files)\n",
    "\n",
    "    #Inputs are PET/CT data, outputs are centres, sizes, masks\n",
    "    X_train, X_test, y_train, y_test = train_test_split(list(zip(pet_files, ct_files, pet_radiomics_files, ct_radiomics_files)), \n",
    "                                                        list(zip(centres, bounding_box_sizes, mask_files)), \n",
    "                                                        test_size=0.33, shuffle=True, random_state=9)\n",
    "    \n",
    "    pet_data['train'][dataset], ct_data['train'][dataset], pet_radiomics_data['train'][dataset], ct_radiomics_data['train'][dataset] = zip(*X_train)\n",
    "    pet_data['test'][dataset], ct_data['test'][dataset], pet_radiomics_data['test'][dataset], ct_radiomics_data['test'][dataset] = zip(*X_test)\n",
    "    centre_data['train'][dataset], size_data['train'][dataset], mask_data['train'][dataset] = zip(*y_train)\n",
    "    centre_data['test'][dataset], size_data['test'][dataset], mask_data['test'][dataset] = zip(*y_test)\n",
    "\n",
    "    print(pet_data['train'][dataset][0], ct_data['train'][dataset][0], pet_radiomics_data['train'][dataset][0], ct_radiomics_data['train'][dataset][0],\n",
    "          centre_data['train'][dataset][0], size_data['train'][dataset][0], mask_data['train'][dataset][0])\n",
    "    \n",
    "    print(pet_data['test'][dataset][0], ct_data['test'][dataset][0], pet_radiomics_data['test'][dataset][0], ct_radiomics_data['test'][dataset][0],\n",
    "          centre_data['test'][dataset][0], size_data['test'][dataset][0], mask_data['test'][dataset][0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "# (x,y,z, a) + (x,y,z, b) => (x,y,z, a+b)\n",
    "def zip_np_volumes(vol_a, vol_b):\n",
    "    vol_a = np.expand_dims(np.array(vol_a), axis=-1)\n",
    "    vol_b = np.expand_dims(np.array(vol_b), axis=-1)\n",
    "    return np.concatenate((vol_a, vol_b), axis=-1)\n",
    "\n",
    "def combine_radiomics_data(source_vol, radiomics_vol):\n",
    "    source_vol = np.expand_dims(source_vol, axis=-1)\n",
    "    return np.concatenate((source_vol, radiomics_vol), axis=-1)\n",
    "\n",
    "def normalise_volume(vol):\n",
    "    std = np.std(vol)\n",
    "    if np.isclose(std, 0):\n",
    "        std = 1\n",
    "    vol = np.divide(vol - np.mean(vol), std)\n",
    "    return vol\n",
    "\n",
    "volume_size = np.array([128,128,128], dtype=int)\n",
    "def shift_centre(centre, shift_weight, seed=None):\n",
    "    np.random.seed(seed) \n",
    "    \n",
    "    centre = np.random.randn(3) * shift_weight * np.array([19.69, 4.24, 13.83]) + centre #std taken from NumpyAnalysis\n",
    "    centre = np.maximum(centre, [0,0,0]) #keep centre within image bounds\n",
    "    centre = np.minimum(centre, volume_size)\n",
    "    return centre.astype(int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#source_volume is a 3-dim (w,h,d) or 4-dim array (w,h,d,channels)\n",
    "#generates a bounding box volume given the parameters\n",
    "#centre is assumed to be relative to the source volume and rounded down\n",
    "#extents is size of the bounding_box / 2\n",
    "def get_bounding_box(source_volume, centre, extents):\n",
    "    extents_ceil = np.ceil(extents).astype(int)\n",
    "    extents_floor = extents.astype(int)\n",
    "    centre = np.rint(centre).astype(int)\n",
    "    \n",
    "    maxima = centre + extents_ceil\n",
    "    minima = centre - extents_floor\n",
    "    \n",
    "    #keep bounding box dimensions within the mask dimensions\n",
    "    maxima = np.minimum(maxima, np.array(source_volume.shape)[:3]).astype(int)\n",
    "    minima = np.maximum(minima, [0,0,0]).astype(int) \n",
    "        \n",
    "    bounding_box_values = source_volume[minima[0]:maxima[0],\n",
    "                                      minima[1]:maxima[1],\n",
    "                                      minima[2]:maxima[2]]\n",
    "    \n",
    "    relative_centre = extents_floor\n",
    "    relative_maxima = relative_centre + maxima - centre\n",
    "    relative_minima = relative_centre + minima - centre\n",
    "    \n",
    "    bbox_shape = tuple(extents_ceil + extents_floor)\n",
    "    if len(source_volume.shape) == 4:\n",
    "        bbox_shape += (source_volume.shape[3],)\n",
    "    bounding_box = np.zeros(bbox_shape)\n",
    "    \n",
    "    bounding_box[relative_minima[0]:relative_maxima[0],\n",
    "                relative_minima[1]:relative_maxima[1],\n",
    "                relative_minima[2]:relative_maxima[2]] = bounding_box_values\n",
    "    \n",
    "    return bounding_box\n",
    "\n",
    "#translates the bounding box so that its values have a new reference centre\n",
    "def align_bounding_box(bounding_box, box_centre, target_centre):\n",
    "    displacement = (target_centre - box_centre).astype(int)\n",
    "    new_box = np.copy(bounding_box)\n",
    "    \n",
    "    #boxes to be shifted forward have a trails of zeroes at the end of the array\n",
    "    #boxes to be shifted backward have a trails of zeroes at the beginning of the array\n",
    "    if displacement[0] < 0:\n",
    "        new_box[displacement[0]:, :, :] = 0\n",
    "    else:\n",
    "        new_box[:displacement[0], :, :] = 0\n",
    "                      \n",
    "    if displacement[1] < 0:\n",
    "        new_box[:, displacement[1]:, :] = 0\n",
    "    else:\n",
    "        new_box[:, :displacement[1], :] = 0\n",
    "    \n",
    "    if displacement[2] < 0:\n",
    "        new_box[:, :, displacement[2]:] = 0\n",
    "    else:\n",
    "        new_box[:, :, :displacement[2]] = 0\n",
    "    \n",
    "    new_box = np.roll(new_box, -displacement, axis=(0,1,2))\n",
    "    return new_box\n",
    "\n",
    "#resample bounding_box to target_size while keeping aspect ratios\n",
    "def resize_bounding_box(bounding_box, target_size, mode):\n",
    "    size_ratio = np.amin(np.divide(target_size, bounding_box.shape))\n",
    "    size_ratio_index = np.argmin(np.divide(target_size, bounding_box.shape))\n",
    "    \n",
    "    #the resampling here is slightly inaccurate since new_size is rounded to the nearest int\n",
    "    bounding_box = sitk.GetImageFromArray(np.transpose(bounding_box))\n",
    "    new_size = np.rint(size_ratio * np.array(bounding_box.GetSize()))\n",
    "    new_spacing = tuple(np.multiply(np.divide(bounding_box.GetSpacing(), size_ratio), np.divide(new_size, new_size+1)))\n",
    "    #additionally multiply spacing by size/(size+1) to fill out an extra row of values\n",
    "        \n",
    "    assert new_size[size_ratio_index] == target_size[size_ratio_index]\n",
    "    \n",
    "    resampler = sitk.ResampleImageFilter()\n",
    "    resampler.SetReferenceImage(bounding_box)\n",
    "    resampler.SetSize((int(new_size[0]), int(new_size[1]), int(new_size[2])))\n",
    "    resampler.SetOutputSpacing(new_spacing)\n",
    "\n",
    "    if mode == 'mask':\n",
    "        resampler.SetInterpolator(sitk.sitkNearestNeighbor)\n",
    "    else:\n",
    "        resampler.SetInterpolator(sitk.sitkLinear)\n",
    "        \n",
    "    new_box = np.transpose(sitk.GetArrayFromImage(resampler.Execute(bounding_box)))\n",
    "    \n",
    "    #pads the surroundings with 0's\n",
    "    return get_bounding_box(new_box, 0.5 * new_size, 0.5 * target_size)\n",
    "            \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Detection Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "#generator for (ct,pet)->centre prediction models\n",
    "def detection_generator(mode='train', dataset='HeadNeckCancer', shuffle=True, normalise=True, use_radiomics=True, batch_size=4):\n",
    "    \n",
    "    def helper(pet_files, ct_files, pet_radiomics_files, ct_radiomics_files, centres, shuffle, normalise, use_radiomics, batch_size):\n",
    "        while True:\n",
    "\n",
    "            if shuffle:\n",
    "                z = list(zip(pet_files, ct_files, pet_radiomics_files, ct_radiomics_files, centres))\n",
    "                random.shuffle(z)\n",
    "                pet_files, ct_files, pet_radiomics_files, ct_radiomics_files, centres = zip(*z)\n",
    "\n",
    "            ct_pet_batch = []\n",
    "            centre_batch = []\n",
    "            \n",
    "            for i in range(len(pet_files)):\n",
    "                ct, pet = np.load(ct_files[i]), np.load(pet_files[i])\n",
    "                \n",
    "                if use_radiomics:\n",
    "                    ct_radiomics = np.load(ct_radiomics_files[i])\n",
    "                    ct = combine_radiomics_data(ct, ct_radiomics)\n",
    "                    \n",
    "                    pet_radiomics = np.load(pet_radiomics_files[i])\n",
    "                    pet = combine_radiomics_data(pet, pet_radiomics)\n",
    "            \n",
    "                else:\n",
    "                    ct = np.expand_dims(ct, axis=-1)\n",
    "                    pet = np.expand_dims(pet, axis=-1)\n",
    "                    \n",
    "                ct_pet = np.concatenate((ct, pet), axis=-1)\n",
    "                    \n",
    "                if normalise:                    \n",
    "                    for axis in range(ct_pet.shape[-1]):\n",
    "                        ct_pet[:,:,:,axis] = normalise_volume(ct_pet[:,:,:,axis])\n",
    "                \n",
    "                ct_pet_batch.append(ct_pet)\n",
    "                centre_batch.append(centres[i])\n",
    "\n",
    "                if len(ct_pet_batch) == batch_size:\n",
    "                    yield np.array(ct_pet_batch), np.array(centre_batch)\n",
    "\n",
    "                    ct_pet_batch.clear()\n",
    "                    centre_batch.clear()\n",
    "\n",
    "            if len(ct_pet_batch) > 0:\n",
    "                yield np.array(ct_pet_batch), np.array(centre_batch)\n",
    "    \n",
    "    return helper(pet_data[mode][dataset], ct_data[mode][dataset], pet_radiomics_data[mode][dataset], ct_radiomics_data[mode][dataset],\n",
    "                  centre_data[mode][dataset], shuffle, normalise, use_radiomics, batch_size)\n",
    "\n",
    "# generator for (ct,pet,centre)->(centre,size) prediction models\n",
    "def localisation_generator(mode='train', dataset='HeadNeckCancer', shuffle=True, normalise=True, use_radiomics=True, batch_size=4, \n",
    "                           shift_weight=0.7, region_radius=[16,16,16]):\n",
    "\n",
    "    \n",
    "    def helper(pet_files, ct_files, pet_radiomics_files, ct_radiomics_files, centres, sizes, shuffle, normalise, use_radiomics, batch_size):\n",
    "        val_seeds = list(range(len(pet_files)))\n",
    "\n",
    "        while True:\n",
    "            if shuffle:\n",
    "                z = list(zip(pet_files, ct_files, pet_radiomics_files, ct_radiomics_files, centres, sizes, val_seeds))\n",
    "                random.shuffle(z)\n",
    "                pet_files, ct_files, pet_radiomics_files, ct_radiomics_files, centres, sizes, val_seeds = zip(*z)\n",
    "\n",
    "            ct_pet_batch = []\n",
    "            shift_centre_batch = []\n",
    "            true_centre_size_batch = []\n",
    "            \n",
    "            for i in range(len(pet_files)):\n",
    "                ct, pet = np.load(ct_files[i]), np.load(pet_files[i])\n",
    "                \n",
    "                if use_radiomics:\n",
    "                    ct_radiomics = np.load(ct_radiomics_files[i])\n",
    "                    ct = combine_radiomics_data(ct, ct_radiomics)\n",
    "                    \n",
    "                    pet_radiomics = np.load(pet_radiomics_files[i])\n",
    "                    pet = combine_radiomics_data(pet, pet_radiomics)\n",
    "            \n",
    "                else:\n",
    "                    ct = np.expand_dims(ct, axis=-1)\n",
    "                    pet = np.expand_dims(pet, axis=-1)\n",
    "                    \n",
    "                ct_pet = np.concatenate((ct, pet), axis=-1)\n",
    "                    \n",
    "                if normalise:                    \n",
    "                    for axis in range(ct_pet.shape[-1]):\n",
    "                        ct_pet[:,:,:,axis] = normalise_volume(ct_pet[:,:,:,axis])\n",
    "                \n",
    "                if mode == 'test': #make the validation generator generate predictable centres\n",
    "                    seed = val_seeds[i]\n",
    "                elif mode == 'train':\n",
    "                    seed = None\n",
    "                \n",
    "                shifted_centre = shift_centre(centres[i], shift_weight, seed)\n",
    "                ct_pet = get_bounding_box(ct_pet, shifted_centre, np.array(region_radius))\n",
    "                                \n",
    "                ct_pet_batch.append(ct_pet)\n",
    "                shift_centre_batch.append(shifted_centre)\n",
    "                true_centre_size_batch.append([centres[i], sizes[i]])\n",
    "\n",
    "                if len(ct_pet_batch) == batch_size:\n",
    "                    yield [np.array(ct_pet_batch), np.array(shift_centre_batch)], np.array(true_centre_size_batch)\n",
    "\n",
    "                    ct_pet_batch.clear()\n",
    "                    shift_centre_batch.clear()\n",
    "                    true_centre_size_batch.clear()\n",
    "\n",
    "            if len(ct_pet_batch) > 0:\n",
    "                yield [np.array(ct_pet_batch), np.array(shift_centre_batch)], np.array(true_centre_size_batch)\n",
    "    \n",
    "    return helper(pet_data[mode][dataset], ct_data[mode][dataset], pet_radiomics_data[mode][dataset], ct_radiomics_data[mode][dataset],\n",
    "                  centre_data[mode][dataset], size_data[mode][dataset], shuffle, normalise, use_radiomics, batch_size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10 batches of 10 samples taken over 98 total samples\n",
      "0 (10, 32, 32, 32, 2) (10, 3) (10, 2, 3)\n"
     ]
    }
   ],
   "source": [
    "#test to see if generator works\n",
    "\n",
    "batch_size = 10\n",
    "dataset = 'HeadNeckCancer'\n",
    "mode = 'test'\n",
    "test_generator = localisation_generator(mode=mode, dataset=dataset, batch_size=batch_size, shuffle=True, use_radiomics=False)\n",
    "print('{} batches of {} samples taken over {} total samples'.format(\n",
    "    math.ceil(len(ct_data[mode][dataset])/batch_size), batch_size, len(ct_data[mode][dataset])))\n",
    "\n",
    "for i in range(1):#math.ceil(len(ct_data[mode][dataset])/batch_size)):\n",
    "    ctpet, centresize = next(test_generator)\n",
    "    print(i, ctpet[0].shape, ctpet[1].shape, centresize.shape)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(3, 2, 3) (3, 1, 2, 3)\n",
      "[-0.33333334 -0.03333333]\n",
      "[-0.33333334 -0.2       ]\n",
      "[[0 5]\n",
      " [0 1]]\n"
     ]
    }
   ],
   "source": [
    "import tensorflow as tf\n",
    "import keras.backend as K\n",
    "import keras.losses\n",
    "\n",
    "#y_true, y_pred in the format (batch_size, centre/size, image axes)\n",
    "#GENERALISED IoU :o o_0 ermagerd\n",
    "def GIoU_helper(IoU_only=False):\n",
    "    \n",
    "    def GIoU(y_true, y_pred):\n",
    "        num_axes=3\n",
    "\n",
    "        #calculate volume of cuboid given its max and min (x,y,z) coordinates\n",
    "        def calculate_volume(max_coords, min_coords):\n",
    "            #K.ones throws an uninitialised variable exception something something\n",
    "            vols = K.maximum(tf.constant([0.0]), max_coords[:, 0] - min_coords[:, 0]) #max_coords.shape[0] = batch size\n",
    "\n",
    "             # if min_coord >= max_coord for some sample in the batch, then that sample has non-overlapping boxes\n",
    "                # assign 0 volume to such samples\n",
    "            for i in range(1, num_axes):\n",
    "                vols = vols * K.maximum(tf.constant([0.0]), max_coords[:, i] - min_coords[:, i])\n",
    "\n",
    "            return vols\n",
    "\n",
    "        #extract values from tensors\n",
    "        centre_pred = K.cast(y_pred[:, 0, :], dtype='float32')\n",
    "        size_pred = K.cast(y_pred[:, 1, :], dtype='float32')\n",
    "        centre_true = K.cast(y_true[:, 0, :], dtype='float32')\n",
    "        size_true = K.cast(y_true[:, 1, :], dtype='float32')\n",
    "\n",
    "        #obtain the min/max coordinates of the bounding boxes\n",
    "        pred_max = centre_pred + 0.5 * size_pred\n",
    "        pred_min = centre_pred - 0.5 * size_pred\n",
    "        true_max = centre_true + 0.5 * size_true\n",
    "        true_min = centre_true - 0.5 * size_true\n",
    "\n",
    "        intersect_min = K.maximum(pred_min, true_min) # \"small\" corner of overlapping box\n",
    "        intersect_max = K.minimum(pred_max, true_max) # \"large\" corner of overlapping box\n",
    "        intersect_volume = calculate_volume(intersect_max, intersect_min)\n",
    "\n",
    "        union_volume = calculate_volume(pred_max, pred_min) + calculate_volume(true_max, true_min) - intersect_volume\n",
    "\n",
    "        hull_min = K.minimum(pred_min, true_min)\n",
    "        hull_max = K.maximum(pred_max, true_max)\n",
    "        hull_volume = calculate_volume(hull_max, hull_min)\n",
    "\n",
    "        IoU = intersect_volume / union_volume\n",
    "        extra = (hull_volume - union_volume) / hull_volume\n",
    "\n",
    "        if IoU_only:\n",
    "            return -(IoU)\n",
    "        else:\n",
    "            return -(IoU - extra)\n",
    "    \n",
    "    return GIoU\n",
    "\n",
    "def IoU(y_true, y_pred):\n",
    "    return GIoU_helper(IoU_only=True)(y_true, y_pred)\n",
    "\n",
    "def hybrid_GIoU_MSE(y_true, y_pred, delta=50):\n",
    "    giou = GIoU_helper()(y_true, y_pred)\n",
    "    mse = K.mean(keras.losses.mean_squared_error(y_true, y_pred))\n",
    "    return delta * giou + mse\n",
    "    \n",
    "\n",
    "with tf.Session() as sess:\n",
    "\n",
    "    t = tf.constant([[[1, 1, 1], [2, 2, 2]],\n",
    "                 [[3, 3, 3], [4, 4, 4]],\n",
    "                 [[5, 5, 5], [6, 6, 6]]])\n",
    "    \n",
    "    u = tf.constant([[[1, 1, 1], [2, 2, 2]], \n",
    "                     [[1, 1, 1], [2, 2, 2]]])\n",
    "    v = tf.constant([[[1, 1, 1], [2, 2, 6]],\n",
    "                    [[1, 2, 1], [2, 2, 4]]])\n",
    "    \n",
    "    a = tf.slice(t, [1, 0, 0], [1, 1, 3])  # [[[3, 3, 3]]]\n",
    "    b = tf.slice(t, [1, 0, 0], [1, 2, 3])  # [[[3, 3, 3],\n",
    "                                       #       [4, 4, 4]]]\n",
    "    c = tf.slice(t, [1, 0, 0], [2, 1, 3])  # [[[3, 3, 3]],\n",
    "                                           #  [[5, 5, 5]]]\n",
    "        \n",
    "    d = K.expand_dims(t, axis=1)\n",
    "        \n",
    "    booleans = tf.constant([False, False])\n",
    "          \n",
    "    centres = tf.slice(t, [0,0,0], [3,1,3])\n",
    "    sizes = tf.slice(t, [0,1,0], [3,1,3])\n",
    "    \n",
    "    print(sess.run(t).shape, sess.run(d).shape)\n",
    "    \n",
    "    print(sess.run(GIoU_helper()(u, v)))\n",
    "    \n",
    "    print(sess.run(IoU(u, v)))\n",
    "    \n",
    "    print(sess.run(keras.losses.mean_squared_error(u, v)))\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Evaluate detection model\n",
    "\n",
    "from sklearn.metrics import mean_squared_error\n",
    "\n",
    "\n",
    "#print(detection_model.evaluate_generator(validation_generator, steps=math.ceil(len(ct_data[mode][dataset])/batch_size)))\n",
    "\n",
    "def evaluate_detection_model(model, dataset, function='detection', batch_size=4):\n",
    "    mode = 'test'\n",
    "    \n",
    "    if function == 'detection':\n",
    "        test_generator = detection_generator(mode=mode, dataset=dataset, shuffle=False, normalise=True, batch_size=batch_size, use_radiomics=False)    \n",
    "        detected_centres = model.predict_generator(test_generator, steps=math.ceil(len(ct_data[mode][dataset])/batch_size))\n",
    "        \n",
    "        return str(mean_squared_error(detected_centres, centre_data[mode][dataset], multioutput='raw_values'))\n",
    "    \n",
    "    elif function == 'localisation':\n",
    "        test_generator = localisation_generator(mode=mode, dataset=dataset, shuffle=False, normalise=True, batch_size=batch_size, use_radiomics=False)    \n",
    "        detected_centresizes = model.predict_generator(test_generator, steps=math.ceil(len(ct_data[mode][dataset])/batch_size))\n",
    "        detected_centre_offsets = detected_centresizes[:, 0, :]\n",
    "        detected_sizes = detected_centresizes[:, 1, :]\n",
    "\n",
    "        #print(detected_centres.shape)\n",
    "        centre_mse =  mean_squared_error(detected_centre_offsets, centre_data[mode][dataset], multioutput='raw_values')\n",
    "        size_mse = mean_squared_error(detected_sizes, size_data[mode][dataset], multioutput='raw_values')\n",
    "    \n",
    "        return str(centre_mse) + str(size_mse)\n",
    "\n",
    "    \n",
    "#print(evaluate_detection_model(detection_model, 'BreastCancer', function='localisation'))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "956\n"
     ]
    }
   ],
   "source": [
    "reset_keras()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from keras.models import load_model\n",
    "from keras import losses\n",
    "import utils\n",
    "\n",
    "#cycles_per_epoch = how many times the entire training set should be cycled over for each epoch\n",
    "#total_cycles = how many times the entire training set should be cycled in total\n",
    "\n",
    "def train_detection_model(dataset, save_model=True, batch_size=4, total_epochs=100, normalise=True, shuffle=True,\n",
    "                          use_radiomics=True, transferring_model=None, loss=losses.mean_squared_error, es_patience=20,\n",
    "                         function='detection'):\n",
    "\n",
    "    def get_model_name():\n",
    "        return '{} {} batch={} cycles={} shuffle={} normalise={} use_radiomics={} transferred={} {}.h5'.format(\n",
    "            function, dataset, batch_size, total_epochs, shuffle, normalise, use_radiomics, transferring_model!=None, loss.__name__)\n",
    "        \n",
    "    #define input shapes\n",
    "    if use_radiomics:\n",
    "        num_channels = 16\n",
    "    else:\n",
    "        num_channels = 2\n",
    "        \n",
    "    if function == 'detection':\n",
    "        detection_input_shape = (128, 128, 128, num_channels)\n",
    "        generator_funct = detection_generator\n",
    "            \n",
    "    elif function == 'localisation':\n",
    "        detection_input_shape = [(32, 32, 32, num_channels), (3,)]\n",
    "        generator_funct = localisation_generator\n",
    "\n",
    "    train_generator = generator_funct(mode='train', dataset=dataset, batch_size=batch_size, shuffle=shuffle, normalise=normalise,\n",
    "                                          use_radiomics=use_radiomics)\n",
    "    validation_generator = generator_funct(mode='test', dataset=dataset, batch_size=batch_size, shuffle=False, normalise=normalise, \n",
    "                                               use_radiomics=use_radiomics)\n",
    "\n",
    "    #determine whether to use old model or build new one from scratch\n",
    "    if transferring_model is not None:\n",
    "        detection_model = transferring_model\n",
    "        #for source_layer, target_layer in zip(transferring_model.layers, detection_model.layers):\n",
    "         #   target_layer.set_weights(source_layer.get_weights())\n",
    "    else:\n",
    "        if function == 'localisation':\n",
    "            resnet_function = Resnet3DBuilder.build_resnet_34\n",
    "            metrics = [IoU, GIoU_helper(), 'mean_squared_error']\n",
    "            \n",
    "        elif function == 'detection':\n",
    "            resnet_function = Resnet3DBuilder.build_resnet_18\n",
    "            metrics = []\n",
    "\n",
    "        \n",
    "        detection_model = resnet_function(detection_input_shape, [3,3], mode=function, reg_factor=1e-8)\n",
    "        detection_model.compile(optimizer='adam',\n",
    "                          loss=loss, metrics=metrics)\n",
    "        \n",
    "    get_best = utils.GetBest(monitor='val_loss', verbose=0, mode='min')\n",
    "    history = detection_model.fit_generator(train_generator, validation_data=validation_generator, verbose=1,\n",
    "                                  validation_steps=math.ceil(len(ct_data['test'][dataset])/batch_size),\n",
    "                                  steps_per_epoch=math.ceil(len(ct_data['train'][dataset])/batch_size), \n",
    "                                  epochs=total_epochs,\n",
    "                                    callbacks=[get_best])\n",
    "\n",
    "    loss_hist = history.history['loss']\n",
    "    val_loss_hist = history.history['val_loss']\n",
    "    best_loss = np.amin(val_loss_hist)\n",
    "    best_mse = evaluate_detection_model(detection_model, dataset, function=function) #np.amin(history.history['val_mean_squared_error'])\n",
    "    \n",
    "    if function == 'localisation':\n",
    "        best_giou = np.amin(history.history['val_GIoU'])\n",
    "        best_iou = np.amin(history.history['val_IoU'])\n",
    "    else:\n",
    "        best_giou = 'N/A'\n",
    "        best_iou = 'N/A'\n",
    "    \n",
    "    with np.printoptions(precision=2):\n",
    "        print('loss:', np.array(loss_hist))\n",
    "        print('val_loss:', np.array(val_loss_hist))\n",
    "        print('best_loss:', best_loss)\n",
    "        print('best_mse:', best_mse)\n",
    "        print('best_iou:', best_iou)\n",
    "        print('best_giou:', best_giou)\n",
    "        \n",
    "    with open('Models/keras models/training_log.txt', 'a+') as log:\n",
    "        log.write('{} \\nbest loss/mse/giou/iou: {} {} {} {}\\ntraining: {}\\nvalidation: {}\\n\\n'.format(\n",
    "            get_model_name(), best_loss, best_mse, best_giou, best_iou, loss_hist, val_loss_hist))\n",
    "    \n",
    "    if save_model:\n",
    "        model_name = get_model_name()\n",
    "        detection_model.save('Models/keras models/' + model_name)\n",
    "        \n",
    "    return detection_model\n",
    "    \n",
    "'''detection_model = train_detection_model(dataset='HeadNeckCancer', use_radiomics=False, total_epochs=200, save_model=True, function='detection',\n",
    "                                   loss=losses.mean_squared_error)'''\n",
    "\n",
    "'''transferring_model=load_model('Models/keras models/detection HeadNeckCancer batch=4 cycles=200 shuffle=True normalise=True use_radiomics=False transferred=False mean_squared_error.h5',\n",
    "                             custom_objects={'jaccard_distance_loss' : jaccard_distance_loss,\n",
    "                                    'dice_coef_loss': dice_coef_loss,\n",
    "                                          'dice_metric': dice_metric,\n",
    "                                          'dice_coef': dice_coef,\n",
    "                                       'conv4d' : conv4d.conv4d,\n",
    "                                    'l2': l2,\n",
    "                                    'Resnet3DBuilder': Resnet3DBuilder,\n",
    "                                    'tf' : tf,\n",
    "                                    'upsample3d_helper' : resnet3d3.upsample3d_helper})\n",
    "detection_model = train_detection_model(dataset='BreastCancer', use_radiomics=False, total_epochs=200, save_model=True, function='detection',\n",
    "                                   loss=losses.mean_squared_error, transferring_model=transferring_model)'''\n",
    "\n",
    "\n",
    "\n",
    "localisation_model = train_detection_model(dataset='BreastCancer', use_radiomics=False, total_epochs=200, save_model=True, function='localisation',\n",
    "                                       loss=hybrid_GIoU_MSE, es_patience=190)\n",
    "reset_keras()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(detected_centres[92], centre_data['test']['HeadNeckCancer'][92])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Segmentation Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "#values derived from NumpyAnalysis.ipynb\n",
    "#get_best_bounding_box()\n",
    "maximal_bounding_volume = np.array([32,32,32])\n",
    "\n",
    "#generator for (ct bounding box, pet bounding_box)-> mask bounding box prediction models\n",
    "def mask_bounding_box_predictor_generator(mode='train', dataset='HeadNeckCancer', shift_centres=True, shift_sizes=True, shift_weight=0.7, \n",
    "                                          shuffle=False, normalise=True, use_radiomics=True, batch_size=4, region_radius=[16,16,16]):\n",
    "\n",
    "            \n",
    "    def helper(pet_files, ct_files, pet_radiomics_files, ct_radiomics_files, mask_files, centres, sizes, \n",
    "               shift_centres, shift_sizes, shuffle, normalise, use_radiomics, batch_size):\n",
    "        val_seeds = list(range(len(pet_files)))\n",
    "        \n",
    "        while True:\n",
    "\n",
    "            if shuffle:\n",
    "                z = list(zip(pet_files, ct_files, pet_radiomics_files, ct_radiomics_files, centres, sizes, mask_files, val_seeds))\n",
    "                random.shuffle(z)\n",
    "                pet_files, ct_files, pet_radiomics_files, ct_radiomics_files, centres, sizes, mask_files, val_seeds = zip(*z)\n",
    "                                \n",
    "            ct_pet_batch = []\n",
    "            mask_batch = []\n",
    "            \n",
    "            for i in range(len(pet_files)):\n",
    "                \n",
    "                if mode == 'test': #make the validation generator generate predictable centres\n",
    "                    seed = val_seeds[i]\n",
    "                elif mode == 'train':\n",
    "                    seed = None\n",
    "                \n",
    "                #randomly sample centres/sizes around the true centre/sizes to generate bounding boxes\n",
    "                if shift_centres:\n",
    "                    centre = shift_centre(centres[i], shift_weight, seed)\n",
    "                else:\n",
    "                    centre = centres[i]\n",
    "\n",
    "                ct = get_bounding_box(np.load(ct_files[i]), centre, np.array(region_radius))\n",
    "                pet = get_bounding_box(np.load(pet_files[i]), centre, np.array(region_radius))\n",
    "                mask = get_bounding_box(np.load(mask_files[i]), centre, np.array(region_radius))\n",
    "\n",
    "                ct_pet = zip_np_volumes(ct, pet)\n",
    "                \n",
    "                if use_radiomics:\n",
    "                    ct_radiomics = np.load(ct_radiomics_files[i])\n",
    "                    pet_radiomics = np.load(pet_radiomics_files[i])\n",
    "                    \n",
    "                    for axis in range(ct_radiomics.shape[-1]):\n",
    "                        ct_radiomic_feature = get_bounding_box(ct_radiomics[:,:,:,axis], centre, np.array(region_radius))\n",
    "                        pet_radiomic_feature = get_bounding_box(pet_radiomics[:,:,:,axis], centre, np.array(region_radius))\n",
    "                        ct_radiomic_feature = np.expand_dims(ct_radiomic_feature, axis=-1)\n",
    "                        pet_radiomic_feature = np.expand_dims(pet_radiomic_feature, axis=-1)\n",
    "                        \n",
    "                        ct_pet = np.concatenate((ct_pet, ct_radiomic_feature, pet_radiomic_feature), axis=-1)\n",
    "                \n",
    "                if normalise:\n",
    "                    for axis in range(ct_pet.shape[-1]):\n",
    "                        ct_pet[:,:,:,axis] = normalise_volume(ct_pet[:,:,:,axis])\n",
    "                \n",
    "                ct_pet_batch.append(ct_pet)\n",
    "                mask_batch.append(mask)\n",
    "                \n",
    "                #print(centre - centres[i])\n",
    "\n",
    "                if len(ct_pet_batch) == batch_size:\n",
    "                    yield np.array(ct_pet_batch), np.array(mask_batch)\n",
    "                    ct_pet_batch.clear()\n",
    "                    mask_batch.clear()\n",
    "\n",
    "            if len(ct_pet_batch) > 0:\n",
    "                yield np.array(ct_pet_batch), np.array(mask_batch)\n",
    "                \n",
    "    return helper(pet_data[mode][dataset], ct_data[mode][dataset], pet_radiomics_data[mode][dataset], ct_radiomics_data[mode][dataset],\n",
    "                  mask_data[mode][dataset], centre_data[mode][dataset], size_data[mode][dataset],\n",
    "                 shift_centres, shift_sizes, shuffle, normalise, use_radiomics, batch_size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(10, 32, 32, 32, 16)\n",
      "mask proportion: 0.0015282107708295127\n"
     ]
    }
   ],
   "source": [
    "#check to see if generator works\n",
    "batch_size = 10\n",
    "mode = 'test'\n",
    "dataset = 'BreastCancer'\n",
    "test_generator = mask_bounding_box_predictor_generator(mode, dataset, shuffle=True, normalise=False, \n",
    "                                                       shift_centres=True, shift_sizes=False, batch_size=batch_size, \n",
    "                                                       use_radiomics=True)\n",
    "num_bg = 0.0\n",
    "num_mask = 0.0\n",
    "\n",
    "for i in range(1):\n",
    "    ctpet, mask = next(test_generator)\n",
    "       \n",
    "    print(ctpet.shape)\n",
    "        \n",
    "    unique, counts = np.unique(mask, return_counts=True)\n",
    "    counts = dict(zip(unique, counts))\n",
    "    \n",
    "    num_bg += counts[0]\n",
    "    num_mask += counts.get(1,0)\n",
    "                      \n",
    "print('mask proportion:', num_mask / num_bg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[5. 5. 7.]\n",
      "[ 7.  7. 11.]\n",
      "[11. 11. 15.]\n",
      "[15. 15. 19.]\n",
      "[18. 19. 23.]\n"
     ]
    }
   ],
   "source": [
    "# calculates output shape given input parameters\n",
    "\n",
    "def deconv_calculator(input_size, kernels, stride):\n",
    "    input_size = np.array(input_size)\n",
    "    stride = np.array(stride)\n",
    "\n",
    "    for kernel in kernels:\n",
    "        \n",
    "        kernel = np.array(kernel)\n",
    "        input_size = np.ceil((input_size - 1) * stride + kernel)\n",
    "        print(input_size)\n",
    "\n",
    "deconv_calculator([3,3,3], [[3,3,5],\n",
    "                            [3,3,5],\n",
    "                            [5,5,5],\n",
    "                            [5,5,5],\n",
    "                           [4,5,5]], [1,1,1])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.67476383 1.93050193]\n"
     ]
    }
   ],
   "source": [
    "from sklearn.utils import class_weight\n",
    "\n",
    "#computes the weights where if they are applied to each item, taking the mean would yield 1\n",
    "#used for segmentation mask weighting\n",
    "def get_class_weight(mask_proportion=0.33, precision=1000):\n",
    "    \n",
    "    arr = np.empty(precision)\n",
    "    index = int(mask_proportion * precision)\n",
    "    arr[:index] = 1\n",
    "    arr[index:] = 0\n",
    "    weights = class_weight.compute_class_weight('balanced', [0,1], arr)\n",
    "    return weights\n",
    "\n",
    "class_weights = get_class_weight(0.25996801724912405) #proportion taken from above\n",
    "print(class_weights)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TYPE                 |Almost_right |half right |all_wrong\n",
      "mult [[1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]]\n",
      "jaccard_distance_loss 0.0\n",
      "binary_crossentropy [1.0174631e-07]\n",
      "binary_crossentropy_scaled [1.]\n",
      "weighted_binary_crossentropy 8.226345e-08\n",
      "weighted_binary_crossentropy_scaled 1.0\n",
      "upsampled [[[[0. 0. 1. 1.]\n",
      "   [0. 0. 1. 1.]\n",
      "   [2. 2. 3. 3.]\n",
      "   [2. 2. 3. 3.]]\n",
      "\n",
      "  [[0. 0. 1. 1.]\n",
      "   [0. 0. 1. 1.]\n",
      "   [2. 2. 3. 3.]\n",
      "   [2. 2. 3. 3.]]\n",
      "\n",
      "  [[4. 4. 5. 5.]\n",
      "   [4. 4. 5. 5.]\n",
      "   [6. 6. 7. 7.]\n",
      "   [6. 6. 7. 7.]]\n",
      "\n",
      "  [[4. 4. 5. 5.]\n",
      "   [4. 4. 5. 5.]\n",
      "   [6. 6. 7. 7.]\n",
      "   [6. 6. 7. 7.]]]]\n"
     ]
    }
   ],
   "source": [
    "from scipy.spatial import distance\n",
    "import keras.backend as K\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "#Accepts tensors; from https://github.com/keras-team/keras/issues/3611\n",
    "#Higher is better\n",
    "def dice_metric(y_true, y_pred, smooth=1):\n",
    "    y_pred = K.round(y_pred)\n",
    "\n",
    "    intersection = K.sum(y_true * y_pred, axis=[1,2,3])\n",
    "    union = K.sum(y_true, axis=[1,2,3]) + K.sum(y_pred, axis=[1,2,3])\n",
    "    return K.mean( (2. * intersection + smooth) / (union + smooth), axis=0)\n",
    "\n",
    "#From https://gist.github.com/wassname/7793e2058c5c9dacb5212c0ac0b18a8a\n",
    "def dice_coef(y_true, y_pred, smooth=1):\n",
    "    \"\"\"\n",
    "    Dice = (2*|X & Y|)/ (|X|+ |Y|)\n",
    "         =  2*sum(|A*B|)/(sum(A^2)+sum(B^2))\n",
    "    ref: https://arxiv.org/pdf/1606.04797v1.pdf\n",
    "    \"\"\"\n",
    "    y_true = K.batch_flatten(y_true)\n",
    "    y_pred = K.batch_flatten(y_pred)\n",
    "    \n",
    "    intersection = K.sum(K.abs(y_true * y_pred), axis=-1)\n",
    "    return K.mean((2. * intersection + smooth) / (K.sum(K.square(y_true),-1) + K.sum(K.square(y_pred),-1) + smooth))\n",
    "\n",
    "def dice_coef_loss(y_true, y_pred):\n",
    "    return 1-dice_coef(y_true, y_pred)\n",
    "\n",
    "\n",
    "#From https://gist.github.com/wassname/f1452b748efcbeb4cb9b1d059dce6f96\n",
    "def jaccard_distance_loss(y_true, y_pred, smooth=100):\n",
    "    \"\"\"\n",
    "    Jaccard = (|X & Y|)/ (|X|+ |Y| - |X & Y|)\n",
    "            = sum(|A*B|)/(sum(|A|)+sum(|B|)-sum(|A*B|))\n",
    "    \n",
    "    The jaccard distance loss is usefull for unbalanced datasets. This has been\n",
    "    shifted so it converges on 0 and is smoothed to avoid exploding or disapearing\n",
    "    gradient.\n",
    "    \n",
    "    Ref: https://en.wikipedia.org/wiki/Jaccard_index\n",
    "    \n",
    "    @url: https://gist.github.com/wassname/f1452b748efcbeb4cb9b1d059dce6f96\n",
    "    @author: wassname\n",
    "    \"\"\"\n",
    "    y_true = K.batch_flatten(y_true)\n",
    "    y_pred = K.batch_flatten(y_pred)\n",
    "    intersection = K.sum(K.abs(y_true * y_pred), axis=-1)\n",
    "    sum_ = K.sum(K.abs(y_true) + K.abs(y_pred), axis=-1)\n",
    "    jac = K.mean((intersection + smooth) / (sum_ - intersection + smooth))\n",
    "    return (1 - jac) * smooth\n",
    "\n",
    "\n",
    "\n",
    "#from https://stackoverflow.com/questions/46009619/keras-weighted-binary-crossentropy\n",
    "def weighted_binary_crossentropy(y_true, y_pred):\n",
    "    zero_weight = class_weights[0]\n",
    "    one_weight = class_weights[1]\n",
    "    \n",
    "    # Calculate the binary crossentropy\n",
    "    b_ce = K.binary_crossentropy(y_true, y_pred)\n",
    "    \n",
    "    # Apply the weights\n",
    "    weight_vector = y_true * one_weight + (1. - y_true) * zero_weight\n",
    "    weighted_b_ce = weight_vector * b_ce\n",
    "\n",
    "    # Return the mean error\n",
    "    return K.mean(weighted_b_ce)\n",
    "\n",
    "def blyat():\n",
    "\n",
    "    # Test\n",
    "    # Test\n",
    "    print(\"TYPE                 |Almost_right |half right |all_wrong\")\n",
    "    y_true = np.array([[1,0,0,0,0,0,0,0,0,0,0]])\n",
    "    y_pred = np.array([[1,0,0,0,0,0,0,0,0,0,0]])\n",
    "    \n",
    "    r = (\n",
    "        K.variable(y_true) *\n",
    "        K.variable(y_pred)\n",
    "    ).eval(session=K.get_session())\n",
    "    print('mult',r)\n",
    "    \n",
    "    r = jaccard_distance_loss(\n",
    "        K.variable(y_true),\n",
    "        K.variable(y_pred),\n",
    "    ).eval(session=K.get_session())\n",
    "    print('jaccard_distance_loss',r)\n",
    "    #assert r[0]<r[1]\n",
    "    #assert r[1]<r[2]\n",
    "\n",
    "    r = keras.losses.binary_crossentropy(\n",
    "        K.variable(y_true),\n",
    "        K.variable(y_pred),\n",
    "    ).eval(session=K.get_session())\n",
    "    print('binary_crossentropy',r)\n",
    "    print('binary_crossentropy_scaled',r/r.max())\n",
    "    #assert r[0]<r[1]\n",
    "    #assert r[1]<r[2]\n",
    "    \n",
    "    r = weighted_binary_crossentropy(\n",
    "        K.variable(y_true),\n",
    "        K.variable(y_pred),\n",
    "    ).eval(session=K.get_session())\n",
    "    print('weighted_binary_crossentropy',r)\n",
    "    print('weighted_binary_crossentropy_scaled',r/r.max())\n",
    "    \n",
    "    ups = np.arange(8).reshape((1,2,2,2))\n",
    "    r = resnet3d3.upsample3d(2)(\n",
    "        K.variable(ups)\n",
    "    ).eval(session=K.get_session())\n",
    "    print('upsampled',r)\n",
    "    \n",
    "blyat()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0.735979651667195, 0.7131510984944008, '10790677.npy', 0.9461363390050938)\n"
     ]
    }
   ],
   "source": [
    "import utils\n",
    "from scipy.ndimage.morphology import binary_fill_holes\n",
    "\n",
    "#Evaluate segmentation model\n",
    "\n",
    "#validation_generator = mask_bounding_box_predictor_generator('test', dataset, shuffle=False, shift_centres=False, shift_sizes=False, batch_size=batch_size)\n",
    "#print(segmentation_model.evaluate_generator(validation_generator, steps=math.ceil(len(ct_data['train'][dataset])/batch_size)))\n",
    "\n",
    "#Accepts numpy arrays\n",
    "#Higher is better\n",
    "def dice_eval(true_mask, pred_mask):\n",
    "    pred_mask = np.rint(pred_mask)\n",
    "    true_mask = np.rint(true_mask)\n",
    "    \n",
    "    true_values = np.unique(true_mask)\n",
    "    pred_values = np.unique(pred_mask)\n",
    "    \n",
    "    \n",
    "    if np.array_equal(true_values, [0]) and np.array_equal(pred_values, [0]):\n",
    "        return 1\n",
    "    \n",
    "    pred_mask = pred_mask.flatten()\n",
    "    true_mask = true_mask.flatten()\n",
    "    ret = 1 - distance.dice(pred_mask, true_mask)\n",
    "    return ret\n",
    "\n",
    "#returns dice acc, slice-by-slice acc and best image segmentation\n",
    "def evaluate_segmentation_model(dataset, model, mode='test', output_segmentation=False, fill_holes=False):\n",
    "    \n",
    "    validation_generator = mask_bounding_box_predictor_generator(\n",
    "                    mode, dataset, shuffle=False, shift_centres=False, shift_sizes=False, use_radiomics=False, normalise=True, batch_size=1)\n",
    "    \n",
    "    dice_sum = 0\n",
    "    max_acc = 0\n",
    "    max_patient = ''\n",
    "    \n",
    "    dice_slice_sum = 0\n",
    "    \n",
    "    for i in range(len(ct_data[mode][dataset])):\n",
    "        ctpet, true_mask = next(validation_generator)\n",
    "        true_mask = true_mask[0]\n",
    "        pred_mask = model.predict(ctpet)[0]\n",
    "        \n",
    "        if fill_holes:\n",
    "            pred_mask = binary_fill_holes(pred_mask)\n",
    "        \n",
    "        dice_index = dice_eval(true_mask, pred_mask)\n",
    "        dice_sum += dice_index\n",
    "                \n",
    "        if dice_index > max_acc:\n",
    "            max_acc = dice_index\n",
    "            max_patient = ct_data[mode][dataset][i]\n",
    "        \n",
    "        if output_segmentation:\n",
    "            patient_name = os.path.splitext(os.path.basename(pet_data[mode][dataset][i]))[0]\n",
    "            mask_radius = 0.5 * np.array((32, 32, 32))\n",
    "            utils.output_numpy_mask_to_nrrd(patient_name, np.rint(pred_mask), centre_data[mode][dataset][i], mask_radius, dataset)\n",
    "            utils.output_numpy_mask_to_nrrd(patient_name, true_mask, centre_data[mode][dataset][i], mask_radius, dataset, filename_tag='-true')\n",
    "            \n",
    "        temp_slice_sum = 0\n",
    "        for j in range(pred_mask.shape[2]):\n",
    "            temp_slice_sum += dice_eval(true_mask[:,:,j], pred_mask[:,:,j])\n",
    "        \n",
    "        dice_slice_sum += temp_slice_sum / pred_mask.shape[2]\n",
    "        \n",
    "        \n",
    "    return dice_sum / len(ct_data[mode][dataset]), dice_slice_sum / len(ct_data[mode][dataset]), os.path.basename(max_patient), max_acc\n",
    "\n",
    "#print('dice accuracy, dice slice-by-slice accuracy, most accurate:')\n",
    "print(evaluate_segmentation_model('BreastCancer', segmentation_model, 'test', output_segmentation=False\n",
    "                                 ))\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1713576\n"
     ]
    }
   ],
   "source": [
    "reset_keras()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Model number for session: 8\n",
      "[<tf.Tensor 'multiply_1/mul:0' shape=(?, 4, 4, 4, 1024) dtype=float32>, <tf.Tensor 'multiply_2/mul:0' shape=(?, 8, 8, 8, 512) dtype=float32>, <tf.Tensor 'multiply_3/mul:0' shape=(?, 16, 16, 16, 256) dtype=float32>, <tf.Tensor 'multiply_4/mul:0' shape=(?, 32, 32, 32, 128) dtype=float32>]\n",
      "Tensor(\"conv3d_42/add:0\", shape=(?, 8, 8, 8, 256), dtype=float32)\n",
      "Tensor(\"conv3d_44/add:0\", shape=(?, 8, 8, 8, 256), dtype=float32)\n",
      "Tensor(\"conv3d_46/add:0\", shape=(?, 16, 16, 16, 128), dtype=float32)\n",
      "Tensor(\"conv3d_48/add:0\", shape=(?, 16, 16, 16, 128), dtype=float32)\n",
      "Tensor(\"conv3d_50/add:0\", shape=(?, 32, 32, 32, 64), dtype=float32)\n",
      "Tensor(\"conv3d_52/add:0\", shape=(?, 32, 32, 32, 64), dtype=float32)\n",
      "Tensor(\"activation_51/Relu:0\", shape=(?, 32, 32, 32, 192), dtype=float32)\n",
      "Epoch 1/100\n",
      "27/27 [==============================] - 41s 2s/step - loss: 57.6628 - dice_coef: 0.6886 - dice_metric: 0.6558 - val_loss: 54.1623 - val_dice_coef: 0.6982 - val_dice_metric: 0.6714\n",
      "Epoch 2/100\n",
      "27/27 [==============================] - 15s 554ms/step - loss: 45.3613 - dice_coef: 0.7769 - dice_metric: 0.7483 - val_loss: 55.5285 - val_dice_coef: 0.6533 - val_dice_metric: 0.6358\n",
      "Epoch 3/100\n",
      "27/27 [==============================] - 15s 557ms/step - loss: 41.5582 - dice_coef: 0.7943 - dice_metric: 0.7658 - val_loss: 47.4749 - val_dice_coef: 0.7596 - val_dice_metric: 0.7289\n",
      "Epoch 4/100\n",
      "27/27 [==============================] - 15s 559ms/step - loss: 38.8660 - dice_coef: 0.8091 - dice_metric: 0.7841 - val_loss: 48.9250 - val_dice_coef: 0.7236 - val_dice_metric: 0.7048\n",
      "Epoch 5/100\n",
      "27/27 [==============================] - 15s 558ms/step - loss: 38.4359 - dice_coef: 0.8068 - dice_metric: 0.7837 - val_loss: 53.0506 - val_dice_coef: 0.6802 - val_dice_metric: 0.6621\n",
      "Epoch 6/100\n",
      "27/27 [==============================] - 15s 560ms/step - loss: 36.7935 - dice_coef: 0.8177 - dice_metric: 0.7965 - val_loss: 43.1142 - val_dice_coef: 0.7718 - val_dice_metric: 0.7469\n",
      "Epoch 7/100\n",
      "27/27 [==============================] - 15s 560ms/step - loss: 36.3065 - dice_coef: 0.8164 - dice_metric: 0.7971 - val_loss: 44.8022 - val_dice_coef: 0.7504 - val_dice_metric: 0.7333\n",
      "Epoch 8/100\n",
      "27/27 [==============================] - 15s 561ms/step - loss: 35.9922 - dice_coef: 0.8200 - dice_metric: 0.7999 - val_loss: 42.6556 - val_dice_coef: 0.7713 - val_dice_metric: 0.7479\n",
      "Epoch 9/100\n",
      "27/27 [==============================] - 15s 561ms/step - loss: 33.8695 - dice_coef: 0.8324 - dice_metric: 0.8138 - val_loss: 48.6797 - val_dice_coef: 0.7221 - val_dice_metric: 0.6907\n",
      "Epoch 10/100\n",
      "27/27 [==============================] - 15s 561ms/step - loss: 33.1666 - dice_coef: 0.8381 - dice_metric: 0.8192 - val_loss: 41.0125 - val_dice_coef: 0.7794 - val_dice_metric: 0.7582\n",
      "Epoch 11/100\n",
      "27/27 [==============================] - 15s 563ms/step - loss: 32.2937 - dice_coef: 0.8433 - dice_metric: 0.8257 - val_loss: 41.7178 - val_dice_coef: 0.7730 - val_dice_metric: 0.7524\n",
      "Epoch 12/100\n",
      "27/27 [==============================] - 15s 562ms/step - loss: 31.1984 - dice_coef: 0.8502 - dice_metric: 0.8334 - val_loss: 40.6742 - val_dice_coef: 0.7769 - val_dice_metric: 0.7592\n",
      "Epoch 13/100\n",
      "27/27 [==============================] - 15s 564ms/step - loss: 31.2630 - dice_coef: 0.8490 - dice_metric: 0.8326 - val_loss: 40.5395 - val_dice_coef: 0.7762 - val_dice_metric: 0.7587\n",
      "Epoch 14/100\n",
      "27/27 [==============================] - 15s 563ms/step - loss: 29.4407 - dice_coef: 0.8614 - dice_metric: 0.8453 - val_loss: 42.7350 - val_dice_coef: 0.7635 - val_dice_metric: 0.7433\n",
      "Epoch 15/100\n",
      "27/27 [==============================] - 15s 562ms/step - loss: 28.9474 - dice_coef: 0.8641 - dice_metric: 0.8493 - val_loss: 42.2791 - val_dice_coef: 0.7690 - val_dice_metric: 0.7488\n",
      "Epoch 16/100\n",
      "27/27 [==============================] - 15s 563ms/step - loss: 27.5438 - dice_coef: 0.8734 - dice_metric: 0.8588 - val_loss: 41.3277 - val_dice_coef: 0.7710 - val_dice_metric: 0.7535\n",
      "Epoch 17/100\n",
      "27/27 [==============================] - 15s 563ms/step - loss: 26.2453 - dice_coef: 0.8806 - dice_metric: 0.8671 - val_loss: 42.8278 - val_dice_coef: 0.7623 - val_dice_metric: 0.7434\n",
      "Epoch 18/100\n",
      "27/27 [==============================] - 15s 561ms/step - loss: 25.1833 - dice_coef: 0.8881 - dice_metric: 0.8749 - val_loss: 42.3810 - val_dice_coef: 0.7651 - val_dice_metric: 0.7465\n",
      "Epoch 19/100\n",
      "27/27 [==============================] - 15s 562ms/step - loss: 24.7780 - dice_coef: 0.8903 - dice_metric: 0.8777 - val_loss: 41.0853 - val_dice_coef: 0.7714 - val_dice_metric: 0.7556\n",
      "Epoch 20/100\n",
      "27/27 [==============================] - 15s 564ms/step - loss: 24.5486 - dice_coef: 0.8917 - dice_metric: 0.8785 - val_loss: 41.2261 - val_dice_coef: 0.7709 - val_dice_metric: 0.7548\n",
      "Epoch 21/100\n",
      "27/27 [==============================] - 15s 563ms/step - loss: 23.6699 - dice_coef: 0.8976 - dice_metric: 0.8848 - val_loss: 42.4391 - val_dice_coef: 0.7613 - val_dice_metric: 0.7463\n",
      "Epoch 22/100\n",
      "27/27 [==============================] - 15s 561ms/step - loss: 22.8080 - dice_coef: 0.9024 - dice_metric: 0.8903 - val_loss: 43.3914 - val_dice_coef: 0.7566 - val_dice_metric: 0.7396\n",
      "Epoch 23/100\n",
      "27/27 [==============================] - 15s 564ms/step - loss: 23.6406 - dice_coef: 0.8966 - dice_metric: 0.8842 - val_loss: 42.3158 - val_dice_coef: 0.7612 - val_dice_metric: 0.7467\n",
      "Epoch 24/100\n",
      "27/27 [==============================] - 15s 564ms/step - loss: 22.7983 - dice_coef: 0.9035 - dice_metric: 0.8914 - val_loss: 41.6094 - val_dice_coef: 0.7638 - val_dice_metric: 0.7522\n",
      "Epoch 25/100\n",
      "27/27 [==============================] - 15s 563ms/step - loss: 21.2384 - dice_coef: 0.9123 - dice_metric: 0.9015 - val_loss: 42.5716 - val_dice_coef: 0.7616 - val_dice_metric: 0.7468\n",
      "Epoch 26/100\n",
      "27/27 [==============================] - 15s 564ms/step - loss: 19.9771 - dice_coef: 0.9193 - dice_metric: 0.9091 - val_loss: 41.8122 - val_dice_coef: 0.7689 - val_dice_metric: 0.7535\n",
      "Epoch 27/100\n",
      "27/27 [==============================] - 15s 564ms/step - loss: 19.0705 - dice_coef: 0.9238 - dice_metric: 0.9143 - val_loss: 42.9346 - val_dice_coef: 0.7592 - val_dice_metric: 0.7434\n",
      "Epoch 28/100\n",
      "27/27 [==============================] - 15s 562ms/step - loss: 18.3758 - dice_coef: 0.9280 - dice_metric: 0.9188 - val_loss: 41.3919 - val_dice_coef: 0.7697 - val_dice_metric: 0.7554\n",
      "Epoch 29/100\n",
      "27/27 [==============================] - 15s 563ms/step - loss: 17.7751 - dice_coef: 0.9310 - dice_metric: 0.9220 - val_loss: 41.7936 - val_dice_coef: 0.7662 - val_dice_metric: 0.7524\n",
      "Epoch 30/100\n",
      "27/27 [==============================] - 15s 562ms/step - loss: 17.1072 - dice_coef: 0.9346 - dice_metric: 0.9260 - val_loss: 41.5663 - val_dice_coef: 0.7703 - val_dice_metric: 0.7548\n",
      "Epoch 31/100\n",
      "27/27 [==============================] - 15s 563ms/step - loss: 16.4881 - dice_coef: 0.9376 - dice_metric: 0.9296 - val_loss: 42.6112 - val_dice_coef: 0.7603 - val_dice_metric: 0.7463\n",
      "Epoch 32/100\n",
      "27/27 [==============================] - 15s 562ms/step - loss: 15.9056 - dice_coef: 0.9406 - dice_metric: 0.9328 - val_loss: 40.9936 - val_dice_coef: 0.7728 - val_dice_metric: 0.7592\n",
      "Epoch 33/100\n",
      "27/27 [==============================] - 15s 562ms/step - loss: 15.2693 - dice_coef: 0.9444 - dice_metric: 0.9369 - val_loss: 41.3717 - val_dice_coef: 0.7701 - val_dice_metric: 0.7561\n",
      "Epoch 34/100\n",
      "27/27 [==============================] - 15s 563ms/step - loss: 14.8108 - dice_coef: 0.9464 - dice_metric: 0.9392 - val_loss: 42.3298 - val_dice_coef: 0.7623 - val_dice_metric: 0.7489\n",
      "Epoch 35/100\n",
      "27/27 [==============================] - 15s 561ms/step - loss: 14.0425 - dice_coef: 0.9504 - dice_metric: 0.9436 - val_loss: 41.7850 - val_dice_coef: 0.7663 - val_dice_metric: 0.7528\n",
      "Epoch 36/100\n",
      "27/27 [==============================] - 15s 564ms/step - loss: 13.4404 - dice_coef: 0.9532 - dice_metric: 0.9469 - val_loss: 43.7214 - val_dice_coef: 0.7506 - val_dice_metric: 0.7364\n",
      "Epoch 37/100\n",
      "27/27 [==============================] - 15s 563ms/step - loss: 12.6076 - dice_coef: 0.9578 - dice_metric: 0.9518 - val_loss: 42.2625 - val_dice_coef: 0.7634 - val_dice_metric: 0.7494\n",
      "Epoch 38/100\n",
      "27/27 [==============================] - 15s 562ms/step - loss: 11.8375 - dice_coef: 0.9613 - dice_metric: 0.9559 - val_loss: 42.3756 - val_dice_coef: 0.7619 - val_dice_metric: 0.7480\n",
      "Epoch 39/100\n",
      "27/27 [==============================] - 15s 563ms/step - loss: 11.3892 - dice_coef: 0.9634 - dice_metric: 0.9584 - val_loss: 41.4261 - val_dice_coef: 0.7698 - val_dice_metric: 0.7558\n",
      "Epoch 40/100\n",
      "27/27 [==============================] - 15s 564ms/step - loss: 11.2293 - dice_coef: 0.9644 - dice_metric: 0.9594 - val_loss: 41.9045 - val_dice_coef: 0.7659 - val_dice_metric: 0.7521\n",
      "Epoch 41/100\n",
      "27/27 [==============================] - 15s 563ms/step - loss: 10.4142 - dice_coef: 0.9680 - dice_metric: 0.9636 - val_loss: 41.5238 - val_dice_coef: 0.7686 - val_dice_metric: 0.7548\n",
      "Epoch 42/100\n",
      "27/27 [==============================] - 15s 563ms/step - loss: 9.9005 - dice_coef: 0.9705 - dice_metric: 0.9664 - val_loss: 41.3310 - val_dice_coef: 0.7691 - val_dice_metric: 0.7566\n",
      "Epoch 43/100\n",
      "27/27 [==============================] - 15s 563ms/step - loss: 9.1690 - dice_coef: 0.9738 - dice_metric: 0.9704 - val_loss: 42.1125 - val_dice_coef: 0.7629 - val_dice_metric: 0.7494\n",
      "Epoch 44/100\n",
      "27/27 [==============================] - 15s 564ms/step - loss: 8.4174 - dice_coef: 0.9775 - dice_metric: 0.9743 - val_loss: 41.9545 - val_dice_coef: 0.7635 - val_dice_metric: 0.7508\n",
      "Epoch 45/100\n",
      "27/27 [==============================] - 15s 562ms/step - loss: 7.9552 - dice_coef: 0.9794 - dice_metric: 0.9765 - val_loss: 40.9033 - val_dice_coef: 0.7711 - val_dice_metric: 0.7591\n",
      "Epoch 46/100\n",
      "27/27 [==============================] - 15s 564ms/step - loss: 8.1774 - dice_coef: 0.9785 - dice_metric: 0.9755 - val_loss: 41.0866 - val_dice_coef: 0.7680 - val_dice_metric: 0.7574\n",
      "Epoch 47/100\n",
      "27/27 [==============================] - 15s 565ms/step - loss: 7.7492 - dice_coef: 0.9804 - dice_metric: 0.9777 - val_loss: 41.2437 - val_dice_coef: 0.7676 - val_dice_metric: 0.7568\n",
      "Epoch 48/100\n",
      "27/27 [==============================] - 15s 562ms/step - loss: 7.1863 - dice_coef: 0.9831 - dice_metric: 0.9807 - val_loss: 41.4406 - val_dice_coef: 0.7656 - val_dice_metric: 0.7552\n",
      "Epoch 49/100\n",
      "27/27 [==============================] - 15s 562ms/step - loss: 6.7091 - dice_coef: 0.9852 - dice_metric: 0.9831 - val_loss: 41.2103 - val_dice_coef: 0.7666 - val_dice_metric: 0.7562\n",
      "Epoch 50/100\n",
      "27/27 [==============================] - 15s 564ms/step - loss: 6.6055 - dice_coef: 0.9857 - dice_metric: 0.9836 - val_loss: 40.4410 - val_dice_coef: 0.7712 - val_dice_metric: 0.7616\n",
      "Epoch 51/100\n",
      "27/27 [==============================] - 15s 562ms/step - loss: 6.3748 - dice_coef: 0.9865 - dice_metric: 0.9845 - val_loss: 40.5843 - val_dice_coef: 0.7702 - val_dice_metric: 0.7603\n",
      "Epoch 52/100\n",
      "27/27 [==============================] - 15s 562ms/step - loss: 5.8938 - dice_coef: 0.9886 - dice_metric: 0.9871 - val_loss: 41.6082 - val_dice_coef: 0.7616 - val_dice_metric: 0.7526\n",
      "Epoch 53/100\n",
      "27/27 [==============================] - 15s 563ms/step - loss: 5.8151 - dice_coef: 0.9889 - dice_metric: 0.9873 - val_loss: 41.3552 - val_dice_coef: 0.7642 - val_dice_metric: 0.7543\n",
      "Epoch 54/100\n",
      "27/27 [==============================] - 15s 564ms/step - loss: 5.6282 - dice_coef: 0.9896 - dice_metric: 0.9882 - val_loss: 41.0576 - val_dice_coef: 0.7654 - val_dice_metric: 0.7565\n",
      "Epoch 55/100\n",
      "27/27 [==============================] - 15s 563ms/step - loss: 5.3306 - dice_coef: 0.9909 - dice_metric: 0.9896 - val_loss: 41.1294 - val_dice_coef: 0.7649 - val_dice_metric: 0.7553\n",
      "Epoch 56/100\n",
      "27/27 [==============================] - 15s 565ms/step - loss: 5.2089 - dice_coef: 0.9913 - dice_metric: 0.9901 - val_loss: 41.4214 - val_dice_coef: 0.7611 - val_dice_metric: 0.7525\n",
      "Epoch 57/100\n",
      "27/27 [==============================] - 15s 564ms/step - loss: 5.0826 - dice_coef: 0.9916 - dice_metric: 0.9905 - val_loss: 40.3912 - val_dice_coef: 0.7697 - val_dice_metric: 0.7607\n",
      "Epoch 58/100\n",
      "27/27 [==============================] - 15s 561ms/step - loss: 4.9735 - dice_coef: 0.9922 - dice_metric: 0.9910 - val_loss: 41.9803 - val_dice_coef: 0.7575 - val_dice_metric: 0.7484\n",
      "Epoch 59/100\n",
      "27/27 [==============================] - 15s 564ms/step - loss: 4.9690 - dice_coef: 0.9920 - dice_metric: 0.9909 - val_loss: 40.8729 - val_dice_coef: 0.7661 - val_dice_metric: 0.7573\n",
      "Epoch 60/100\n",
      "27/27 [==============================] - 15s 565ms/step - loss: 4.8339 - dice_coef: 0.9925 - dice_metric: 0.9915 - val_loss: 41.5764 - val_dice_coef: 0.7606 - val_dice_metric: 0.7514\n",
      "Epoch 61/100\n",
      "27/27 [==============================] - 15s 563ms/step - loss: 4.6910 - dice_coef: 0.9930 - dice_metric: 0.9920 - val_loss: 40.1696 - val_dice_coef: 0.7690 - val_dice_metric: 0.7617\n",
      "Epoch 62/100\n",
      "27/27 [==============================] - 15s 563ms/step - loss: 4.5825 - dice_coef: 0.9934 - dice_metric: 0.9924 - val_loss: 40.0896 - val_dice_coef: 0.7697 - val_dice_metric: 0.7618\n",
      "Epoch 63/100\n",
      "27/27 [==============================] - 15s 564ms/step - loss: 4.5501 - dice_coef: 0.9934 - dice_metric: 0.9925 - val_loss: 40.6143 - val_dice_coef: 0.7668 - val_dice_metric: 0.7577\n",
      "Epoch 64/100\n",
      "27/27 [==============================] - 15s 564ms/step - loss: 4.3427 - dice_coef: 0.9943 - dice_metric: 0.9935 - val_loss: 41.5058 - val_dice_coef: 0.7589 - val_dice_metric: 0.7507\n",
      "Epoch 65/100\n",
      "27/27 [==============================] - 15s 564ms/step - loss: 4.1825 - dice_coef: 0.9948 - dice_metric: 0.9941 - val_loss: 40.9095 - val_dice_coef: 0.7633 - val_dice_metric: 0.7552\n",
      "Epoch 66/100\n",
      "27/27 [==============================] - 15s 563ms/step - loss: 4.1016 - dice_coef: 0.9950 - dice_metric: 0.9943 - val_loss: 41.0377 - val_dice_coef: 0.7615 - val_dice_metric: 0.7536\n",
      "Epoch 67/100\n",
      "27/27 [==============================] - 15s 563ms/step - loss: 3.9682 - dice_coef: 0.9955 - dice_metric: 0.9948 - val_loss: 40.0283 - val_dice_coef: 0.7689 - val_dice_metric: 0.7610\n",
      "Epoch 68/100\n",
      "27/27 [==============================] - 15s 561ms/step - loss: 3.9563 - dice_coef: 0.9954 - dice_metric: 0.9947 - val_loss: 42.1616 - val_dice_coef: 0.7528 - val_dice_metric: 0.7443\n",
      "Epoch 69/100\n",
      "27/27 [==============================] - 15s 563ms/step - loss: 4.0362 - dice_coef: 0.9950 - dice_metric: 0.9942 - val_loss: 42.2725 - val_dice_coef: 0.7508 - val_dice_metric: 0.7430\n",
      "Epoch 70/100\n",
      "27/27 [==============================] - 15s 564ms/step - loss: 3.8908 - dice_coef: 0.9956 - dice_metric: 0.9949 - val_loss: 39.5731 - val_dice_coef: 0.7723 - val_dice_metric: 0.7644\n",
      "Epoch 71/100\n",
      "27/27 [==============================] - 15s 562ms/step - loss: 3.8964 - dice_coef: 0.9954 - dice_metric: 0.9948 - val_loss: 41.0975 - val_dice_coef: 0.7605 - val_dice_metric: 0.7531\n",
      "Epoch 72/100\n",
      "27/27 [==============================] - 15s 564ms/step - loss: 3.8952 - dice_coef: 0.9954 - dice_metric: 0.9947 - val_loss: 41.6698 - val_dice_coef: 0.7562 - val_dice_metric: 0.7482\n",
      "Epoch 73/100\n",
      "27/27 [==============================] - 15s 564ms/step - loss: 4.0419 - dice_coef: 0.9947 - dice_metric: 0.9939 - val_loss: 41.9456 - val_dice_coef: 0.7540 - val_dice_metric: 0.7463\n",
      "Epoch 74/100\n",
      "27/27 [==============================] - 15s 562ms/step - loss: 4.1612 - dice_coef: 0.9942 - dice_metric: 0.9933 - val_loss: 44.6164 - val_dice_coef: 0.7324 - val_dice_metric: 0.7245\n",
      "Epoch 75/100\n",
      "27/27 [==============================] - 15s 564ms/step - loss: 4.7054 - dice_coef: 0.9923 - dice_metric: 0.9911 - val_loss: 42.9350 - val_dice_coef: 0.7486 - val_dice_metric: 0.7399\n",
      "Epoch 76/100\n",
      "27/27 [==============================] - 15s 564ms/step - loss: 5.2365 - dice_coef: 0.9904 - dice_metric: 0.9890 - val_loss: 43.9175 - val_dice_coef: 0.7424 - val_dice_metric: 0.7330\n",
      "Epoch 77/100\n",
      "27/27 [==============================] - 15s 565ms/step - loss: 5.3356 - dice_coef: 0.9906 - dice_metric: 0.9890 - val_loss: 41.6278 - val_dice_coef: 0.7580 - val_dice_metric: 0.7498\n",
      "Epoch 78/100\n",
      "27/27 [==============================] - 15s 561ms/step - loss: 5.6197 - dice_coef: 0.9897 - dice_metric: 0.9881 - val_loss: 41.7712 - val_dice_coef: 0.7574 - val_dice_metric: 0.7500\n",
      "Epoch 79/100\n",
      "27/27 [==============================] - 15s 565ms/step - loss: 5.1469 - dice_coef: 0.9919 - dice_metric: 0.9907 - val_loss: 41.6477 - val_dice_coef: 0.7588 - val_dice_metric: 0.7509\n",
      "Epoch 80/100\n",
      "27/27 [==============================] - 15s 563ms/step - loss: 4.4245 - dice_coef: 0.9951 - dice_metric: 0.9943 - val_loss: 40.5995 - val_dice_coef: 0.7661 - val_dice_metric: 0.7584\n",
      "Epoch 81/100\n",
      "27/27 [==============================] - 15s 562ms/step - loss: 4.1026 - dice_coef: 0.9963 - dice_metric: 0.9957 - val_loss: 40.7773 - val_dice_coef: 0.7646 - val_dice_metric: 0.7573\n",
      "Epoch 82/100\n",
      "27/27 [==============================] - 15s 562ms/step - loss: 3.9207 - dice_coef: 0.9968 - dice_metric: 0.9963 - val_loss: 40.5889 - val_dice_coef: 0.7651 - val_dice_metric: 0.7580\n",
      "Epoch 83/100\n",
      "27/27 [==============================] - 15s 563ms/step - loss: 3.7838 - dice_coef: 0.9971 - dice_metric: 0.9967 - val_loss: 41.0962 - val_dice_coef: 0.7598 - val_dice_metric: 0.7529\n",
      "Epoch 84/100\n",
      "27/27 [==============================] - 15s 562ms/step - loss: 3.5843 - dice_coef: 0.9977 - dice_metric: 0.9974 - val_loss: 41.1628 - val_dice_coef: 0.7598 - val_dice_metric: 0.7531\n",
      "Epoch 85/100\n",
      "27/27 [==============================] - 15s 562ms/step - loss: 3.4967 - dice_coef: 0.9978 - dice_metric: 0.9975 - val_loss: 41.3292 - val_dice_coef: 0.7584 - val_dice_metric: 0.7511\n",
      "Epoch 86/100\n",
      "27/27 [==============================] - 15s 564ms/step - loss: 3.3646 - dice_coef: 0.9981 - dice_metric: 0.9979 - val_loss: 42.2871 - val_dice_coef: 0.7496 - val_dice_metric: 0.7428\n",
      "Epoch 87/100\n",
      "27/27 [==============================] - 15s 562ms/step - loss: 3.3333 - dice_coef: 0.9980 - dice_metric: 0.9978 - val_loss: 40.7544 - val_dice_coef: 0.7618 - val_dice_metric: 0.7551\n",
      "Epoch 88/100\n",
      "27/27 [==============================] - 15s 562ms/step - loss: 3.2513 - dice_coef: 0.9981 - dice_metric: 0.9979 - val_loss: 41.2388 - val_dice_coef: 0.7571 - val_dice_metric: 0.7507\n",
      "Epoch 89/100\n",
      "27/27 [==============================] - 15s 564ms/step - loss: 3.1234 - dice_coef: 0.9985 - dice_metric: 0.9983 - val_loss: 41.1691 - val_dice_coef: 0.7574 - val_dice_metric: 0.7506\n",
      "Epoch 90/100\n",
      "27/27 [==============================] - 15s 563ms/step - loss: 3.0949 - dice_coef: 0.9984 - dice_metric: 0.9982 - val_loss: 41.1468 - val_dice_coef: 0.7577 - val_dice_metric: 0.7510\n",
      "Epoch 91/100\n",
      "27/27 [==============================] - 15s 562ms/step - loss: 3.0281 - dice_coef: 0.9984 - dice_metric: 0.9982 - val_loss: 40.8475 - val_dice_coef: 0.7594 - val_dice_metric: 0.7525\n",
      "Epoch 92/100\n",
      "27/27 [==============================] - 15s 563ms/step - loss: 2.9545 - dice_coef: 0.9986 - dice_metric: 0.9985 - val_loss: 40.7828 - val_dice_coef: 0.7588 - val_dice_metric: 0.7523\n",
      "Epoch 93/100\n",
      "27/27 [==============================] - 15s 564ms/step - loss: 2.8859 - dice_coef: 0.9987 - dice_metric: 0.9986 - val_loss: 40.6405 - val_dice_coef: 0.7603 - val_dice_metric: 0.7534\n",
      "Epoch 94/100\n",
      "27/27 [==============================] - 15s 562ms/step - loss: 2.8080 - dice_coef: 0.9988 - dice_metric: 0.9987 - val_loss: 40.1830 - val_dice_coef: 0.7628 - val_dice_metric: 0.7562\n",
      "Epoch 95/100\n",
      "27/27 [==============================] - 15s 564ms/step - loss: 2.7313 - dice_coef: 0.9989 - dice_metric: 0.9989 - val_loss: 40.9922 - val_dice_coef: 0.7566 - val_dice_metric: 0.7501\n",
      "Epoch 96/100\n",
      "27/27 [==============================] - 15s 563ms/step - loss: 2.6836 - dice_coef: 0.9990 - dice_metric: 0.9989 - val_loss: 41.7313 - val_dice_coef: 0.7512 - val_dice_metric: 0.7446\n",
      "Epoch 97/100\n",
      "27/27 [==============================] - 15s 561ms/step - loss: 2.6452 - dice_coef: 0.9989 - dice_metric: 0.9988 - val_loss: 41.0338 - val_dice_coef: 0.7554 - val_dice_metric: 0.7490\n",
      "Epoch 98/100\n",
      "27/27 [==============================] - 15s 564ms/step - loss: 2.5877 - dice_coef: 0.9990 - dice_metric: 0.9989 - val_loss: 40.2878 - val_dice_coef: 0.7617 - val_dice_metric: 0.7552\n",
      "Epoch 99/100\n",
      "27/27 [==============================] - 15s 565ms/step - loss: 2.5525 - dice_coef: 0.9990 - dice_metric: 0.9989 - val_loss: 40.4667 - val_dice_coef: 0.7596 - val_dice_metric: 0.7533\n",
      "Epoch 100/100\n",
      "27/27 [==============================] - 15s 563ms/step - loss: 2.5692 - dice_coef: 0.9988 - dice_metric: 0.9986 - val_loss: 42.5265 - val_dice_coef: 0.7423 - val_dice_metric: 0.7360\n",
      "loss: [57.82 45.43 41.57 38.91 38.31 36.78 36.32 36.02 33.91 33.18 32.3  31.21\n",
      " 31.18 29.47 28.96 27.55 26.15 25.2  24.72 24.45 23.7  22.73 23.65 22.8\n",
      " 21.26 19.95 19.1  18.39 17.81 17.11 16.44 15.94 15.31 14.74 14.07 13.46\n",
      " 12.64 11.8  11.39 11.22 10.42  9.9   9.18  8.42  7.99  8.19  7.77  7.22\n",
      "  6.71  6.6   6.37  5.91  5.75  5.64  5.34  5.22  5.08  4.99  4.97  4.83\n",
      "  4.7   4.59  4.56  4.35  4.19  4.11  3.96  3.96  4.04  3.89  3.9   3.89\n",
      "  4.05  4.14  4.7   5.23  5.34  5.62  5.14  4.43  4.09  3.92  3.78  3.59\n",
      "  3.49  3.36  3.33  3.25  3.12  3.1   3.03  2.96  2.89  2.81  2.73  2.68\n",
      "  2.65  2.59  2.55  2.57]\n",
      "val_loss: [54.16 55.53 47.47 48.92 53.05 43.11 44.8  42.66 48.68 41.01 41.72 40.67\n",
      " 40.54 42.74 42.28 41.33 42.83 42.38 41.09 41.23 42.44 43.39 42.32 41.61\n",
      " 42.57 41.81 42.93 41.39 41.79 41.57 42.61 40.99 41.37 42.33 41.78 43.72\n",
      " 42.26 42.38 41.43 41.9  41.52 41.33 42.11 41.95 40.9  41.09 41.24 41.44\n",
      " 41.21 40.44 40.58 41.61 41.36 41.06 41.13 41.42 40.39 41.98 40.87 41.58\n",
      " 40.17 40.09 40.61 41.51 40.91 41.04 40.03 42.16 42.27 39.57 41.1  41.67\n",
      " 41.95 44.62 42.93 43.92 41.63 41.77 41.65 40.6  40.78 40.59 41.1  41.16\n",
      " 41.33 42.29 40.75 41.24 41.17 41.15 40.85 40.78 40.64 40.18 40.99 41.73\n",
      " 41.03 40.29 40.47 42.53]\n",
      "dice_index: 0.04875987902350325\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "\"segmentation_model = train_segmentation_model('BreastCancer', batch_size=4, total_epochs=100,\\n                                              shift_centres=False, shift_sizes=False, normalise=True, shuffle=True, use_radiomics=False,\\n    transferring_model=load_model('/home/jzhe0882/Models/keras models/segmentation HeadNeckCancer batch=4 cycles=100 shift_centre=False shift_size=False shuffle=True normalise=True use_radiomics=False transferred=False jaccard_distance_loss.h5', \\n                                custom_objects={'jaccard_distance_loss' : jaccard_distance_loss,\\n                                    'dice_coef_loss': dice_coef_loss,\\n                                          'dice_metric': dice_metric,\\n                                          'dice_coef': dice_coef,\\n                                       'conv4d' : conv4d.conv4d,\\n                                    'l2': l2,\\n                                    'Resnet3DBuilder': Resnet3DBuilder,\\n                                    'tf' : tf,\\n                                    'upsample3d_helper' : resnet3d3.upsample3d_helper},\\n                                 compile=True))\""
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from keras.models import load_model\n",
    "\n",
    "\n",
    "#cycles_per_epoch = how many times the entire training set should be cycled over for each epoch\n",
    "#total_cycles = how many times the entire training set should be cycled in total\n",
    "def train_segmentation_model(dataset, batch_size=4, total_epochs=50, shift_centres=False, shift_sizes=False, shuffle=True, \n",
    "                             normalise=True, save_model=True, use_radiomics=True, transferring_model=None, loss=jaccard_distance_loss,\n",
    "                            es_patience=200):\n",
    "    \n",
    "    def get_model_name():\n",
    "        return 'segmentation {} batch={} cycles={} shift_centre={} shift_size={} shuffle={} normalise={} use_radiomics={} transferred={} {}.h5'.format(\n",
    "            dataset, batch_size, total_epochs, shift_centres, shift_sizes, shuffle, normalise, use_radiomics, transferring_model!=None, loss.__name__)\n",
    "\n",
    "    #bounding box values derived from NumpyAnalysis.ipynb\n",
    "    if use_radiomics:\n",
    "        segmentation_input_shape = (32, 32, 32, 16)\n",
    "    else:\n",
    "        segmentation_input_shape = (32, 32, 32, 2)\n",
    "    segmentation_output_shape = (32, 32, 32)\n",
    "    \n",
    "    train_generator = mask_bounding_box_predictor_generator(\n",
    "        'train', dataset, shuffle=shuffle, shift_centres=shift_centres, shift_sizes=shift_sizes, \n",
    "        normalise=normalise, use_radiomics=use_radiomics, batch_size=batch_size)\n",
    "    validation_generator = mask_bounding_box_predictor_generator(\n",
    "        'test', dataset, shuffle=False, shift_centres=False, shift_sizes=False, \n",
    "        normalise=normalise, use_radiomics=use_radiomics, batch_size=batch_size)\n",
    "\n",
    "\n",
    "    #copy all layer data over until the avg-pool/dense layer in the detection model; indices are found with model.summary()\n",
    "    #this is terrible :|\n",
    "    #for detect_layer, segment_layer in zip(detection_model.layers[1:-4], segmentation_model.layers[1:]):\n",
    "    #    segment_layer.set_weights(detect_layer.get_weights())\n",
    "\n",
    "    if transferring_model is not None:\n",
    "        segmentation_model = transferring_model\n",
    "        #for source_layer, target_layer in zip(transferring_model.layers, segmentation_model.layers):\n",
    "        #    target_layer.set_weights(source_layer.get_weights())\n",
    "    else:\n",
    "        segmentation_model = Resnet3DBuilder.build_resnet_18(segmentation_input_shape, segmentation_output_shape, mode='segmentation')\n",
    "        segmentation_model.compile(optimizer='adam',\n",
    "                      loss=[loss],\n",
    "                        metrics=[dice_coef, dice_metric])\n",
    "        \n",
    "    es = EarlyStopping(monitor='val_loss', patience=es_patience, restore_best_weights=True)\n",
    "    history = segmentation_model.fit_generator(train_generator, validation_data=validation_generator, verbose=1,\n",
    "                                             validation_steps=math.ceil(len(ct_data['test'][dataset])/batch_size), \n",
    "                                             epochs=total_epochs, steps_per_epoch=math.ceil(len(ct_data['train'][dataset])/batch_size), \n",
    "                                            callbacks=[es])\n",
    "    \n",
    "    loss_hist = history.history['loss']\n",
    "    val_hist = history.history['val_loss']\n",
    "    dice_index,_,_,_ = evaluate_segmentation_model(dataset, segmentation_model)\n",
    "    \n",
    "    with np.printoptions(precision=2):\n",
    "        print('loss:', np.array(loss_hist))\n",
    "        print('val_loss:', np.array(val_hist))\n",
    "        print('dice_index:', dice_index)\n",
    "        \n",
    "    with open('Models/keras models/training_log.txt', 'a+') as log:\n",
    "        log.write('{} \\ndice index: {} \\ntraining: {}\\nvalidation: {}\\n\\n'.format(get_model_name(), dice_index, loss_hist, val_hist))\n",
    "    \n",
    "    if save_model:\n",
    "        model_name = get_model_name()\n",
    "        segmentation_model.save('Models/keras models/' + model_name)\n",
    "        \n",
    "        \n",
    "    return segmentation_model\n",
    "\n",
    "'''segmentation_model = train_segmentation_model('HeadNeckCancer', batch_size=4, total_epochs=100,\n",
    "                                              shift_centres=False, shift_sizes=False, normalise=True, shuffle=True, use_radiomics=False)'''\n",
    "\n",
    "segmentation_model = train_segmentation_model('BreastCancer', batch_size=4, total_epochs=100,\n",
    "                                              shift_centres=False, shift_sizes=False, normalise=True, shuffle=True, use_radiomics=False)\n",
    "\n",
    "'''segmentation_model = train_segmentation_model('BreastCancer', batch_size=4, total_epochs=100,\n",
    "                                              shift_centres=False, shift_sizes=False, normalise=True, shuffle=True, use_radiomics=False,\n",
    "    transferring_model=load_model('/home/jzhe0882/Models/keras models/segmentation HeadNeckCancer batch=4 cycles=100 shift_centre=False shift_size=False shuffle=True normalise=True use_radiomics=False transferred=False jaccard_distance_loss.h5', \n",
    "                                custom_objects={'jaccard_distance_loss' : jaccard_distance_loss,\n",
    "                                    'dice_coef_loss': dice_coef_loss,\n",
    "                                          'dice_metric': dice_metric,\n",
    "                                          'dice_coef': dice_coef,\n",
    "                                       'conv4d' : conv4d.conv4d,\n",
    "                                    'l2': l2,\n",
    "                                    'Resnet3DBuilder': Resnet3DBuilder,\n",
    "                                    'tf' : tf,\n",
    "                                    'upsample3d_helper' : resnet3d3.upsample3d_helper},\n",
    "                                 compile=True))'''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING:tensorflow:From /opt/anaconda/lib/python3.7/site-packages/tensorflow/python/ops/math_ops.py:3066: to_int32 (from tensorflow.python.ops.math_ops) is deprecated and will be removed in a future version.\n",
      "Instructions for updating:\n",
      "Use tf.cast instead.\n"
     ]
    }
   ],
   "source": [
    "segmentation_model = load_model('/home/jzhe0882/Models/keras models/segmentation BreastCancer batch=8 cycles=100 shift_centre=False shift_size=False shuffle=True normalise=True use_radiomics=False transferred=True jaccard_distance_loss.h5',\n",
    "                               custom_objects={'jaccard_distance_loss' : jaccard_distance_loss,\n",
    "                                    'dice_coef_loss': dice_coef_loss,\n",
    "                                          'dice_metric': dice_metric,\n",
    "                                          'dice_coef': dice_coef})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING:tensorflow:From /opt/anaconda/lib/python3.7/site-packages/tensorflow/python/ops/math_ops.py:3066: to_int32 (from tensorflow.python.ops.math_ops) is deprecated and will be removed in a future version.\n",
      "Instructions for updating:\n",
      "Use tf.cast instead.\n"
     ]
    }
   ],
   "source": [
    "from keras.models import load_model\n",
    "\n",
    "segmentation_model=load_model('Models/keras models/segmentation BreastCancer False 50 8 True True.h5',\n",
    "                           custom_objects={'dice_coef_loss': dice_coef_loss,\n",
    "                                          'dice_metric': dice_metric,\n",
    "                                          'dice_coef': dice_coef})\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "F"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Playground"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def blyat(arr):\n",
    "    return np.roll(arr, [2,1], [0,1])\n",
    "    \n",
    "test = np.arange(20).reshape(4,5)\n",
    "print(test)\n",
    "test=blyat(test)\n",
    "print(test)\n",
    "test[-1:, :] = 0\n",
    "print(test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test = np.arange(1,31).reshape(2,3,5)\n",
    "print(test)\n",
    "centre = [1,1,1]\n",
    "print(test[centre[0], centre[1], centre[2]])\n",
    "test = get_bounding_box(test, np.array(centre), np.array([2,2,2]))\n",
    "print(test)\n",
    "test = align_bounding_box(test, np.array(centre), np.array([3,3,3]))\n",
    "print(test)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(detection_model.get_layer(index=-4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2, 3, 4)\n",
      "(4, 6, 8)\n",
      "[[[2. 1. 1. 2. 2. 4. 4. 0.]\n",
      "  [2. 1. 1. 2. 2. 4. 4. 0.]\n",
      "  [2. 1. 1. 2. 2. 4. 4. 0.]\n",
      "  [2. 1. 1. 2. 2. 4. 4. 0.]\n",
      "  [2. 1. 1. 2. 2. 4. 4. 0.]\n",
      "  [0. 0. 0. 0. 0. 0. 0. 0.]]\n",
      "\n",
      " [[2. 1. 1. 2. 2. 4. 4. 0.]\n",
      "  [2. 1. 1. 3. 3. 4. 4. 0.]\n",
      "  [2. 1. 1. 3. 3. 4. 4. 0.]\n",
      "  [2. 1. 1. 2. 2. 4. 4. 0.]\n",
      "  [2. 1. 1. 2. 2. 4. 4. 0.]\n",
      "  [0. 0. 0. 0. 0. 0. 0. 0.]]\n",
      "\n",
      " [[2. 1. 1. 2. 2. 4. 4. 0.]\n",
      "  [2. 1. 1. 3. 3. 4. 4. 0.]\n",
      "  [2. 1. 1. 3. 3. 4. 4. 0.]\n",
      "  [2. 1. 1. 2. 2. 4. 4. 0.]\n",
      "  [2. 1. 1. 2. 2. 4. 4. 0.]\n",
      "  [0. 0. 0. 0. 0. 0. 0. 0.]]\n",
      "\n",
      " [[0. 0. 0. 0. 0. 0. 0. 0.]\n",
      "  [0. 0. 0. 0. 0. 0. 0. 0.]\n",
      "  [0. 0. 0. 0. 0. 0. 0. 0.]\n",
      "  [0. 0. 0. 0. 0. 0. 0. 0.]\n",
      "  [0. 0. 0. 0. 0. 0. 0. 0.]\n",
      "  [0. 0. 0. 0. 0. 0. 0. 0.]]]\n"
     ]
    }
   ],
   "source": [
    "#resample bounding_box to target_size while keeping aspect ratios\n",
    "def resize_bounding_box2(bounding_box, target_size, mode):\n",
    "    size_ratio = np.amin(np.divide(target_size, bounding_box.shape))\n",
    "    size_ratio_index = np.argmin(np.divide(target_size, bounding_box.shape))\n",
    "    \n",
    "    #the resampling here is slightly inaccurate since new_size is rounded to the nearest int\n",
    "    bounding_box = sitk.GetImageFromArray(np.transpose(bounding_box))\n",
    "    new_size = np.rint(size_ratio * np.array(bounding_box.GetSize()))\n",
    "    new_spacing = tuple(np.divide(bounding_box.GetSpacing(), size_ratio))\n",
    "            \n",
    "    assert new_size[size_ratio_index] == target_size[size_ratio_index]\n",
    "    \n",
    "    resampler = sitk.ResampleImageFilter()\n",
    "    resampler.SetReferenceImage(bounding_box)\n",
    "    resampler.SetSize((int(new_size[0]), int(new_size[1]), int(new_size[2])))\n",
    "    resampler.SetOutputSpacing(new_spacing)\n",
    "\n",
    "    if mode == 'mask':\n",
    "        resampler.SetInterpolator(sitk.sitkNearestNeighbor)\n",
    "    else:\n",
    "        resampler.SetInterpolator(sitk.sitkLinear)\n",
    "        \n",
    "    new_box = np.transpose(sitk.GetArrayFromImage(resampler.Execute(bounding_box)))\n",
    "    \n",
    "    #pads the surroundings with 0's\n",
    "    return get_bounding_box(new_box, 0.5 * new_size, 0.5 * target_size)\n",
    "\n",
    "#resample bounding_box to target_size while keeping aspect ratios\n",
    "#alternative to resize_bounding_box that doesnt centre the volume\n",
    "def resize_bounding_box3(bounding_box, target_size, mode):\n",
    "    size_ratio = np.amin(np.divide(target_size, bounding_box.shape))\n",
    "    size_ratio_index = np.argmin(np.divide(target_size, bounding_box.shape))\n",
    "    \n",
    "    #the resampling here is slightly inaccurate since new_size is rounded to the nearest int\n",
    "    bounding_box = sitk.GetImageFromArray(np.transpose(bounding_box))\n",
    "    new_size = target_size\n",
    "    new_spacing = tuple(np.divide(bounding_box.GetSpacing(), size_ratio))\n",
    "        \n",
    "    assert new_size[size_ratio_index] == target_size[size_ratio_index]\n",
    "    \n",
    "    resampler = sitk.ResampleImageFilter()\n",
    "    resampler.SetReferenceImage(bounding_box)\n",
    "    resampler.SetSize((int(new_size[0]), int(new_size[1]), int(new_size[2])))\n",
    "    resampler.SetOutputSpacing(new_spacing)\n",
    "\n",
    "    if mode == 'mask':\n",
    "        resampler.SetInterpolator(sitk.sitkNearestNeighbor)\n",
    "    else:\n",
    "        resampler.SetInterpolator(sitk.sitkLinear)\n",
    "        \n",
    "    new_box = resampler.Execute(bounding_box)\n",
    "    new_box = np.transpose(sitk.GetArrayFromImage(resampler.Execute(new_box)))\n",
    "        \n",
    "    \n",
    "    return new_box\n",
    "\n",
    "\n",
    "test = np.array([[[2,1,2,4],\n",
    "                  [2,1,2,4],\n",
    "                  [2,1,2,4]], \n",
    "                 \n",
    "                 [[2,1,2,4],\n",
    "                  [2,1,3,4],\n",
    "                  [2,1,2,4]]], dtype=float)\n",
    "\n",
    "#test = np.arange(1,1001).reshape((10,10,10)).astype(float)\n",
    "print(test.shape)\n",
    "test = resize_bounding_box3(test, np.array([4,6,8]), 'mask')\n",
    "#test = resize_bounding_box3(test, np.array([2,3,4]), 'mask')\n",
    "print(test.shape)\n",
    "print(test)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "#old validation generator for resizing stuff\n",
    "maximal_bounding_volume = np.array([32,32,32])\n",
    "\n",
    "#generator for (ct bounding box, pet bounding_box)-> mask bounding box prediction models\n",
    "def mask_bounding_box_predictor_generator(mode='train', dataset='HeadNeckCancer', shift_centres=False, shift_sizes=False, shift_weight=0.7, \n",
    "                                          shuffle=False, normalise=True, use_radiomics=True, batch_size=4, region_radius=[16,16,16]):\n",
    "\n",
    "    def shift_centre(centre):\n",
    "        centre = np.random.randn(3) * shift_weight * np.array([20.4,4.6,21.2]) + centre #std taken from NumpyAnalysis\n",
    "        centre = np.maximum(centre, [0,0,0]) #keep centre within image bounds\n",
    "        centre = np.minimum(centre, volume_size)\n",
    "        return centre.astype(int)\n",
    "    \n",
    "    def shift_size(size):\n",
    "        size = np.random.randn(3) * shift_weight *  np.array([3.6,3,5.7]) + size #std taken from NumpyAnalysis\n",
    "        size = np.maximum(size, [1,1,1]) #keep size to > 1\n",
    "        size = np.minimum(size, maximal_bounding_volume)\n",
    "        return size\n",
    "        \n",
    "    def get_resized_bounding_box(source_volume, centre, size, target_size, mode):\n",
    "        bounding_box = get_bounding_box(source_volume, centre, 0.5 * size)\n",
    "        return resize_bounding_box(bounding_box, target_size, mode)\n",
    "            \n",
    "    def helper(pet_files, ct_files, pet_radiomics_files, ct_radiomics_files, mask_files, centres, sizes, \n",
    "               shift_centres, shift_sizes, shuffle, normalise, use_radiomics, batch_size):\n",
    "        while True:\n",
    "\n",
    "            if shuffle:\n",
    "                z = list(zip(pet_files, ct_files, pet_radiomics_files, ct_radiomics_files, centres, sizes, mask_files))\n",
    "                random.shuffle(z)\n",
    "                pet_files, ct_files, pet_radiomics_files, ct_radiomics_files, centres, sizes, mask_files = zip(*z)\n",
    "                                \n",
    "            ct_pet_batch = []\n",
    "            mask_batch = []\n",
    "            \n",
    "            for i in range(len(pet_files)):\n",
    "                \n",
    "                #randomly sample centres/sizes around the true centre/sizes to generate bounding boxes\n",
    "                if shift_centres:\n",
    "                    centre = shift_centre(centres[i])\n",
    "                else:\n",
    "                    centre = centres[i]\n",
    "\n",
    "                if shift_sizes:\n",
    "                    size = shift_size(sizes[i])\n",
    "                else:\n",
    "                    size = sizes[i]\n",
    "  \n",
    "\n",
    "                ct = get_resized_bounding_box(np.load(ct_files[i]), centre, size, maximal_bounding_volume, 'ct')\n",
    "                pet = get_resized_bounding_box(np.load(pet_files[i]), centre, size, maximal_bounding_volume, 'pet')\n",
    "                mask = get_resized_bounding_box(np.load(mask_files[i]), centre, size, maximal_bounding_volume, 'mask')\n",
    "                \n",
    "                ct_pet = zip_np_volumes(ct, pet)\n",
    "                \n",
    "                if use_radiomics:\n",
    "                    ct_radiomics = np.load(ct_radiomics_files[i])\n",
    "                    pet_radiomics = np.load(pet_radiomics_files[i])\n",
    "                    \n",
    "                    for axis in range(ct_radiomics.shape[-1]):\n",
    "                        ct_radiomic_feature = get_resized_bounding_box(ct_radiomics[:,:,:,axis], centre, size, maximal_bounding_volume, 'ct')\n",
    "                        pet_radiomic_feature = get_resized_bounding_box(pet_radiomics[:,:,:,axis], centre, size, maximal_bounding_volume, 'pet')\n",
    "                        ct_radiomic_feature = np.expand_dims(ct_radiomic_feature, axis=-1)\n",
    "                        pet_radiomic_feature = np.expand_dims(pet_radiomic_feature, axis=-1)\n",
    "                        \n",
    "                        ct_pet = np.concatenate((ct_pet, ct_radiomic_feature, pet_radiomic_feature), axis=-1)\n",
    "                \n",
    "                if normalise:\n",
    "                    for axis in range(ct_pet.shape[-1]):\n",
    "                        ct_pet[:,:,:,axis] = normalise_volume(ct_pet[:,:,:,axis])\n",
    "                \n",
    "                ct_pet_batch.append(ct_pet)\n",
    "                mask_batch.append(mask)\n",
    "\n",
    "                if len(ct_pet_batch) == batch_size:\n",
    "                    yield np.array(ct_pet_batch), np.array(mask_batch)\n",
    "                    ct_pet_batch.clear()\n",
    "                    mask_batch.clear()\n",
    "\n",
    "            if len(ct_pet_batch) > 0:\n",
    "                yield np.array(ct_pet_batch), np.array(mask_batch)\n",
    "                \n",
    "    return helper(pet_data[mode][dataset], ct_data[mode][dataset], pet_radiomics_data[mode][dataset], ct_radiomics_data[mode][dataset],\n",
    "                  mask_data[mode][dataset], centre_data[mode][dataset], size_data[mode][dataset],\n",
    "                 shift_centres, shift_sizes, shuffle, normalise, use_radiomics, batch_size)"
   ]
  },
  {
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}