{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\ninput_shape  = (224,224,3)\n\nimport glob\nimport pickle\nimport random\nimport pydicom\n\nimport pandas as pd\nimport numpy as np\nfrom tqdm import tqdm\nfrom numba import jit, prange\n\nimport matplotlib.pyplot as plt\nimport matplotlib\nimport seaborn as sns\n\nimport tensorflow as tf\nprint(tf.__version__)\n\nimport tensorflow.keras as keras\nfrom tensorflow.keras.models import Sequential, Model\nfrom tensorflow.keras.layers import Dense, concatenate, Input\nfrom scipy.spatial import distance\nfrom tensorflow.keras.preprocessing import image\nfrom tensorflow.keras.applications.vgg19 import VGG19, preprocess_input\n#from tensorflow.keras.applications.inception_v3 import InceptionV3, preprocess_input\nfrom tensorflow.keras.layers import Conv2D, GlobalMaxPooling2D, Conv3D, MaxPool3D, BatchNormalization, Dropout, GlobalAveragePooling3D, Flatten, LeakyReLU\nfrom tensorflow.keras.initializers import HeUniform\nfrom tensorflow.keras.models import load_model\nfrom tensorflow.keras.optimizers import Adam\n\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.manifold import TSNE\n\nimport os\nimport cv2\nfrom pydicom.pixel_data_handlers.util import apply_voi_lut\n\n\nimport time\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\nimport re\n# from gluoncv.utils import viz\n# import mxnet as mx\nimport multiprocessing as mp\n\n# Input data files are available in the read-only \"../input/\" directory\n\n\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2021-09-06T17:23:50.900183Z","iopub.execute_input":"2021-09-06T17:23:50.900518Z","iopub.status.idle":"2021-09-06T17:23:58.443835Z","shell.execute_reply.started":"2021-09-06T17:23:50.900449Z","shell.execute_reply":"2021-09-06T17:23:58.442994Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"! dir '../input/rsna-miccai-brain-tumor-radiogenomic-classification/'\n\nroot_dir = '../input/rsna-miccai-brain-tumor-radiogenomic-classification/'\ndf_train = pd.read_csv(root_dir+'train_labels.csv')\nsns.countplot(data=df_train, x='MGMT_value')","metadata":{"execution":{"iopub.status.busy":"2021-09-06T17:23:58.445424Z","iopub.execute_input":"2021-09-06T17:23:58.445818Z","iopub.status.idle":"2021-09-06T17:23:59.383174Z","shell.execute_reply.started":"2021-09-06T17:23:58.445779Z","shell.execute_reply":"2021-09-06T17:23:59.382227Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Add the full paths for each id for different types of sequences to the csv \ndef full_ids(data):\n    zeros = 5 - len(str(data))\n    if zeros > 0:\n        prefix = ''.join(['0' for i in range(zeros)])\n    \n    return prefix+str(data)\n        \n\ndf_train['PatientID'] = df_train['BraTS21ID'].apply(full_ids)\n\n# Add all the paths to the df for easy access\ndf_train['flair'] = df_train['PatientID'].apply(lambda file_id : root_dir+'train/'+file_id+'/FLAIR/')\ndf_train['t1w']   = df_train['PatientID'].apply(lambda file_id : root_dir+'train/'+file_id+'/T1w/')\ndf_train['t1wce'] = df_train['PatientID'].apply(lambda file_id : root_dir+'train/'+file_id+'/T1wCE/')\ndf_train['t2w']   = df_train['PatientID'].apply(lambda file_id : root_dir+'train/'+file_id+'/T2w/')\n\ndf_train = df_train.drop(['BraTS21ID'], axis=1)\n# Delete column\ndf_train","metadata":{"execution":{"iopub.status.busy":"2021-09-06T17:23:59.384864Z","iopub.execute_input":"2021-09-06T17:23:59.385189Z","iopub.status.idle":"2021-09-06T17:23:59.424585Z","shell.execute_reply.started":"2021-09-06T17:23:59.385145Z","shell.execute_reply":"2021-09-06T17:23:59.423506Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def _dicom2array(path, voi_lut=True, fix_monochrome=True, resize=True, img_size=256):\n    dicom = pydicom.read_file(path)\n    data  = apply_voi_lut(dicom.pixel_array, dicom)\n    if dicom.PhotometricInterpretation == \"MONOCHROME1\":\n        data = np.amax(data) - data\n    #Normalize the data: subtract off the minimum, divide by the maximum, convert to 256 uint8\n    data = data - np.min(data)\n    data = data/np.max(data)\n    data = (data * 255).astype(np.uint8)\n    \n    #Resize images to target value\n    data = cv2.resize(data, (img_size, img_size))\n    data = cv2.equalizeHist(data)\n    return data\n\n\ndef _circumscriber(img: np.array) -> np.array:\n    #First is vertical, second is horizontal, third is slices\n    vmin = 0\n    vlimit = img.shape[0]\n    hmin = 0\n    hlimit = img.shape[1]\n    \n    for i in range(vlimit):\n        if np.max(img[i, :]) == 0:\n            vmin += 1\n        else:\n            break\n    vmax = vmin + 1\n    for i in range(vmin+1, vlimit):\n        if np.max(img[i, :]) > 0:\n            vmax += 1\n        else:\n            break\n    \n    for j in range(hlimit):\n        if np.max(img[:, j]) == 0:\n            hmin += 1\n        else:\n            break\n    hmax = hmin + 1\n    for j in range(hmin+1, hlimit):\n        if np.max(img[:, j]) > 0:\n            hmax += 1\n        else:\n            break\n    return img[vmin:vmax, hmin:hmax]\n\ndef crop_image(img,tol=100):\n    # img is 2D image data\n    # tol  is tolerance\n    mask = img>tol\n    return img[np.ix_(mask.any(1),mask.any(0))]\n\nimport operator\ndef cropND(img, bounding):\n    start = tuple(map(lambda a, da: a//2-da//2, img.shape, bounding))\n    end   = tuple(map(operator.add, start, bounding))\n    slices = tuple(map(slice, start, end))\n    return img[slices]\n\n\ndcm_path = '../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00014/T1w/Image-190.dcm'\n#print(data)\nimg1  = _dicom2array(dcm_path)\nprint(img1.shape)\n\n#print(img.shape)\n\nfig = plt.figure(figsize=(128, 128))\n\nfig.add_subplot(15, 15, 1)\nplt.imshow(img1, cmap=\"bone\")\n\n","metadata":{"execution":{"iopub.status.busy":"2021-09-06T19:05:08.901572Z","iopub.execute_input":"2021-09-06T19:05:08.902205Z","iopub.status.idle":"2021-09-06T19:05:09.539228Z","shell.execute_reply.started":"2021-09-06T19:05:08.902156Z","shell.execute_reply":"2021-09-06T19:05:09.538124Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import scipy\n\ntrain_path = '/kaggle/input/rsna-miccai-brain-tumor-radiogenomic-classification'\n\n\ndef save_slices(root_dir, patient, modality, mgmt, slices):\n    filename = '{}\\\\{}_{}_{}.pickle'.format(root_dir, patient, modality, mgmt)\n    with open(filename, 'wb') as handle:\n        pickle.dump(slices, handle, protocol=pickle.HIGHEST_PROTOCOL)\n    \n    \n    \ndef load_slices(patient, modality, mgmt):\n    filename = '{}\\\\{}_{}_{}.pickle'.format(root_dir, patient, modality, mgmt)\n    with open(filename, 'rb') as handle:\n        slices = pickle.load(handle)\n    return slices  \n\n\n\ndef get_ordered_slices(path:str):\n    \"\"\"\n    get a list of ordered slices\n    \"\"\"\n    images = [os.path.basename(f) for f in glob.glob(f'{path}/*.dcm')]\n    images.sort(key=lambda f: int(re.sub('\\D', '', f)))\n    return images\n      \n\ndef get_ordered_slices_restrict(path:str):\n    \"\"\"\n    get a list of ordered slices\n    \"\"\"\n    images   = [os.path.basename(f) for f in glob.glob(f'{path}/*.dcm')]\n    images.sort(key=lambda f: int(re.sub('\\D', '', f)))\n    len_imgs = len(images)\n    fourth   = len_imgs//4\n    images   = images[fourth:-fourth]\n   \n    return images    \n        \n\n\n\ndef load_FULL_brain_restricted(dcm_path):\n    \"\"\"\n    send all of the images in the chosen modality, in order, as a single 3D np array\n    \"\"\"\n    \n    images = get_ordered_slices(dcm_path)\n    real_images = [_dicom2array(dcm_path + f) for f in images]\n    # np.count_nonzero(anchor) < 10000)\n    good_images = np.array([im for im in real_images if (np.count_nonzero(im) > 7000)]).T\n    \n    #final_image = _circumscriber(good_images)\n    return good_images\n\n\ndef load_FULL_brain(dcm_path):\n    \"\"\"\n    send all of the images in the chosen modality, in order, as a single 3D np array\n    \"\"\"\n    \n    images = get_ordered_slices(dcm_path)\n    real_images = [_dicom2array(dcm_path + f) for f in images]\n    # np.count_nonzero(anchor) < 10000)\n    good_images = np.array([im for im in real_images if (np.max(im) > 0)]).T\n    #final_image = _circumscriber(good_images)\n    return good_images\n\n\ndef get_nr_slices(path: str):\n    \"\"\"\n    return nr of slices\n    \"\"\"\n    return len(glob.glob(f'{path}/*.dcm'))\n    \n# Old implementaiton   \n# def get_ordered_slice(path:str, index:int):\n#     \"\"\"\n#     get specific slice\n#     \"\"\"\n#     images = [os.path.basename(f) for f in glob.glob(f'{path}/*.dcm')]\n#     images.sort(key=lambda f: int(re.sub('\\D', '', f)))\n#     dicom_path = \"{}{}\".format(path, images[index])\n#     return _dicom2array(dicom_path)\n\n\n\ndef get_ordered_slice(slice_arr:str, slice_path:int, slice_idx: int):\n    \"\"\"\n    get specific slice\n    \"\"\"\n    #print(slice_arr[slice_idx])\n    #print(slice_path)\n   \n    #print(\"-------------\")\n    slice_file      = slice_arr[slice_idx]\n    slice_fullpath  = \"{}{}\".format(slice_path, slice_file)\n    return np.array(_dicom2array(slice_fullpath))\n\n\n\"\"\"\n'Accent', 'Accent_r', 'Blues', 'Blues_r', 'BrBG', 'BrBG_r', 'BuGn', 'BuGn_r', 'BuPu', 'BuPu_r', 'CMRmap', 'CMRmap_r', 'Dark2', 'Dark2_r', 'GnBu', 'GnBu_r', 'Greens', 'Greens_r', 'Greys', 'Greys_r', 'OrRd', 'OrRd_r', 'Oranges', 'Oranges_r', 'PRGn', 'PRGn_r', 'Paired', 'Paired_r', 'Pastel1', 'Pastel1_r', 'Pastel2', 'Pastel2_r', 'PiYG', 'PiYG_r', 'PuBu', 'PuBuGn', 'PuBuGn_r', 'PuBu_r', 'PuOr', 'PuOr_r', 'PuRd', 'PuRd_r', 'Purples', 'Purples_r', 'RdBu', 'RdBu_r', 'RdGy', 'RdGy_r', 'RdPu', 'RdPu_r', 'RdYlBu', 'RdYlBu_r', 'RdYlGn', 'RdYlGn_r', 'Reds', 'Reds_r', 'Set1', 'Set1_r', 'Set2', 'Set2_r', 'Set3', 'Set3_r', 'Spectral', 'Spectral_r', 'Wistia', 'Wistia_r', 'YlGn', 'YlGnBu', 'YlGnBu_r', 'YlGn_r', 'YlOrBr', 'YlOrBr_r', 'YlOrRd', 'YlOrRd_r', 'afmhot', 'afmhot_r', 'autumn', 'autumn_r', 'binary', 'binary_r', 'bone', 'bone_r', 'brg', 'brg_r', 'bwr', 'bwr_r', 'cividis', 'cividis_r', 'cool', 'cool_r', 'coolwarm', 'coolwarm_r', 'copper', 'copper_r', 'crest', 'crest_r', 'cubehelix', 'cubehelix_r', 'flag', 'flag_r', 'flare', 'flare_r', 'gist_earth', 'gist_earth_r', 'gist_gray', 'gist_gray_r', 'gist_heat', 'gist_heat_r', 'gist_ncar', 'gist_ncar_r', 'gist_rainbow', 'gist_rainbow_r', 'gist_stern', 'gist_stern_r', 'gist_yarg', 'gist_yarg_r', 'gnuplot', 'gnuplot2', 'gnuplot2_r', 'gnuplot_r', 'gray', 'gray_r', 'hot', 'hot_r', 'hsv', 'hsv_r', 'icefire', 'icefire_r', 'inferno', 'inferno_r', 'jet', 'jet_r', 'magma', 'magma_r', 'mako', 'mako_r', 'nipy_spectral', 'nipy_spectral_r', 'ocean', 'ocean_r', 'pink', 'pink_r', 'plasma', 'plasma_r', 'prism', 'prism_r', 'rainbow', 'rainbow_r', 'rocket', 'rocket_r', 'seismic', 'seismic_r', 'spring', 'spring_r', 'summer', 'summer_r', 'tab10', 'tab10_r', 'tab20', 'tab20_r', 'tab20b', 'tab20b_r', 'tab20c', 'tab20c_r', 'terrain', 'terrain_r', 'turbo', 'turbo_r', 'twilight', 'twilight_r', 'twilight_shifted', 'twilight_shifted_r', 'viridis', 'viridis_r', 'vlag', 'vlag_r', 'winter', 'winter_r'\n\"\"\"\n\ndef plot_imgs(imgs, cols=4, size=7, is_rgb=True, title=\"\", cmap='nipy_spectral', img_size=(224,224)):\n    rows = imgs.shape[2]//cols + 2\n    fig = plt.figure(figsize=(cols*size, rows*size))\n    for i in range(imgs.shape[2]-1):\n        img = imgs[:,:,i]\n        fig.add_subplot(rows, cols, i+1)\n        plt.imshow(img)\n    plt.suptitle(title)\n    plt.show()\n\n\n    \ndef plot_img_list(imgs, cols=4, size=7, is_rgb=True, title=\"\", cmap='bone', img_size=(224,224)):\n    rows = len(imgs)\n    fig = plt.figure(figsize=(cols*size, rows*size))\n    for i in range(rows):\n        fig.add_subplot(rows, cols, i+1)\n        plt.imshow(imgs[i])\n    plt.suptitle(title)\n    plt.show()\n    \n\n    \ndef resample(image, scan, new_spacing=[1,1,1]):\n    # Determine current pixel spacing\n    spacing       = np.array([scan[0].SliceThickness] + list(scan[0].PixelSpacing), dtype=np.float32)\n    resize_factor = spacing / new_spacing\n    new_shape     = np.round(image.shape * resize_factor)\n    \n    # recompute the resize factor and spacing such that we match the rounded new shape above\n    rounded_resize_factor = new_shape / image.shape\n    rounded_new_spacing   = spacing / rounded_resize_factor\n    print(rounded_resize_factor)\n    # zoom with resize factor\n    image = scipy.ndimage.interpolation.zoom(image, (1,1,rounded_resize_factor[2]), mode='nearest')\n    \n    return image, rounded_new_spacing\n\n\n\ndef load_scans(dcm_path):\n    \n    slices = [pydicom.dcmread(file) for file in glob.glob(dcm_path+\"/*\")]\n    slices.sort(key = lambda x: float(x.ImagePositionPatient[2]))\n#     else: #basepath == \"../input/osic-pulmonary-fibrosis-progression/\":    \n#         # This competition shows missing values in ImagePosition,\n#         # this is why we are sorting by filename number\n#         files = listdir(dcm_path)\n#         file_nums = [np.int(file.split(\".\")[0]) for file in files]\n#         sorted_file_nums = np.sort(file_nums)[::-1]\n#        slices = [pydicom.dcmread(dcm_path + \"/\" + str(file_num) + \".dcm\" ) for file_num in sorted_file_nums]\n    return slices\n\n\n\n\npath2 = \"../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00018/FLAIR/\"\n#path2      = \"../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00122/T1wCE/\"\n\n\nboxtest1   = load_FULL_brain(path2)\nscans      = load_scans(path2)\n\nprint(boxtest1.shape)\n\n#get a volume for 3d cnn\nvol = boxtest1[:,:,60:60+32]\nprint(vol.shape)\nplot_imgs(vol)\nt0 = time.time()\nimg_resampled, spacing = resample(boxtest1, scans, [1,1,1])\nt1 = time.time()\nprint(t1-t0)\nprint(img_resampled.shape)\nplot_imgs(img_resampled)\n\n#boxtest2 = load_FULL_brain_restricted(path2)\n#print(boxtest2.shape)\n\n#plot_imgs(boxtest2[:,:,-30:-10])","metadata":{"execution":{"iopub.status.busy":"2021-09-06T20:05:24.902932Z","iopub.execute_input":"2021-09-06T20:05:24.903470Z","iopub.status.idle":"2021-09-06T20:05:34.482479Z","shell.execute_reply.started":"2021-09-06T20:05:24.903431Z","shell.execute_reply":"2021-09-06T20:05:34.480422Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\ndef get_image_plane(data):\n    '''\n    Returns the MRI's plane from the dicom data.\n    \n    '''\n    x1,y1,_,x2,y2,_ = [round(j) for j in data.ImageOrientationPatient]\n    cords = [x1,y1,x2,y2]\n\n    if cords == [1,0,0,0]:\n        return 'coronal'\n    if cords == [1,0,0,1]:\n        return 'axial'\n    if cords == [0,1,0,0]:\n        return 'sagittal'\n    \ndef get_image(data):\n    '''\n    Returns the image data as a numpy array.\n    '''  \n    if np.max(data.pixel_array)==0:\n        img = data.pixel_array\n    else:\n        img = data.pixel_array/np.max(data.pixel_array)\n        img = (img * 255).astype(np.uint8)\n        \n    return img\n\ndef sorted_image_dirs(path: str):\n    '''\n    Sorts the list of image directories by image number in a path\n    '''\n    dirs = glob.glob(path+'*')\n    dirs.sort(key=lambda x: int(x.split('/')[-1].split('-')[-1].split('.')[0]))\n    \n    return dirs\n\n\ndef get_all_images(path: str):\n    '''\n    Returns a list of (non blank) images from a given path (of shape [non_blank_image_count, 512, 512])\n    '''\n    image_dirs = sorted_image_dirs(path)\n    images = []\n    \n    for directory in image_dirs:\n        data = pydicom.dcmread(directory)\n        \n        img = get_image(data)\n        \n        # Exclude the blank images\n        if np.max(img)!=0:\n            images.append(img)\n        else:\n            pass\n    \n    return images\n\ndef get_plane(path: str, plane=False):\n    '''\n    Returns the middle image from the path\n    \n    if plane=True returns the plane\n    '''\n    image_dirs = sorted_image_dirs(path)\n\n    dicom = pydicom.dcmread(image_dirs[len(image_dirs)//2])\n    img   = get_image(dicom)\n    plane = get_image_plane(dicom)\n  \n    return plane\n\n","metadata":{"execution":{"iopub.status.busy":"2021-09-06T17:24:00.189762Z","iopub.execute_input":"2021-09-06T17:24:00.190033Z","iopub.status.idle":"2021-09-06T17:24:00.204745Z","shell.execute_reply.started":"2021-09-06T17:24:00.189996Z","shell.execute_reply":"2021-09-06T17:24:00.204142Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"seq_axes = {'flair':('flair_axis', []), 't1w':('t1w_axis', []), 't1wce':('t1wce_axis',[]), 't2w':('t2w_axis',[])}\n\nfor key, value in seq_axes.items():\n    df_train[value[0]] = df_train[key].apply(lambda dicom_dir : get_plane(dicom_dir, plane=True))\n","metadata":{"execution":{"iopub.status.busy":"2021-09-06T17:24:00.205823Z","iopub.execute_input":"2021-09-06T17:24:00.206270Z","iopub.status.idle":"2021-09-06T17:25:53.575151Z","shell.execute_reply.started":"2021-09-06T17:24:00.206243Z","shell.execute_reply":"2021-09-06T17:25:53.574321Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Exclude bad samples\nto_exclude = [\"00109\", \"00123\", \"00709\"]\nprint(df_train.shape[0])\ndf_train = df_train[~df_train['PatientID'].isin(to_exclude)]\nprint(df_train.shape[0])\n\ndf_train.head()","metadata":{"execution":{"iopub.status.busy":"2021-09-06T17:25:53.578007Z","iopub.execute_input":"2021-09-06T17:25:53.578397Z","iopub.status.idle":"2021-09-06T17:25:53.597187Z","shell.execute_reply.started":"2021-09-06T17:25:53.578359Z","shell.execute_reply":"2021-09-06T17:25:53.596389Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Create TEST dataset","metadata":{}},{"cell_type":"code","source":"df_test = df_train.sample(frac=0.17, random_state=2022)","metadata":{"execution":{"iopub.status.busy":"2021-09-06T17:25:53.598991Z","iopub.execute_input":"2021-09-06T17:25:53.599457Z","iopub.status.idle":"2021-09-06T17:25:53.610955Z","shell.execute_reply.started":"2021-09-06T17:25:53.599420Z","shell.execute_reply":"2021-09-06T17:25:53.610202Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(len(df_train.index))\n\nprint(len(df_test.index))\n\nsns.countplot(data=df_test, x='MGMT_value')","metadata":{"execution":{"iopub.status.busy":"2021-09-06T17:25:53.612278Z","iopub.execute_input":"2021-09-06T17:25:53.612728Z","iopub.status.idle":"2021-09-06T17:25:53.725920Z","shell.execute_reply.started":"2021-09-06T17:25:53.612684Z","shell.execute_reply":"2021-09-06T17:25:53.725129Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"patient = \"00811\"\nslices_path = df_train['flair'][df_train['PatientID'] == patient ].iloc[0]\nprint(slices_path)\nflair_images = load_FULL_brain_restricted(slices_path)\n\nprint('No of images:', flair_images.shape[2])\nprint('MGMT: ', df_train['MGMT_value'][df_train['PatientID'] == patient ].iloc[0])\nprint(flair_images.shape)\n# np_images =  np.array([im for im in list(flair_images.values())]).T\n\n# Expand the dimensions of the image, this is so that it fits the expected model input format\n# The shape of the input should be (batch_size, width, height, nr_channels) = (1, 224, 224, 3)\nimgs_expanded = [np.expand_dims(flair_images[:,:,idx], axis = 2) for idx  in range(flair_images.shape[2])]\nimgs_tensor   = [tf.convert_to_tensor(img) for img in imgs_expanded]\nrgb_converted = [tf.image.grayscale_to_rgb(img) for img in imgs_tensor]\n\n# rgb_converted = tf.image.grayscale_to_rgb(img_tensor)\n\nimg_expanded = np.expand_dims(rgb_converted[12], axis = 0)\nprint(\"Input shape:\",img_expanded.shape)\nfig, ax = plt.subplots(1,2, figsize=(30,10))\n\nindex = -2\nax[0].imshow(rgb_converted[index], cmap='icefire')\nax[1].imshow(flair_images[:,:,index], cmap='icefire')\n\n# Pre-process the img in the same way original images were\nimg_ready = preprocess_input(img_expanded)\nprint(\"Input shape <after preproc>:\", img_ready.shape)\n# result = embedding_model.predict(img_ready)\n# print(\"Output shape:\", result.shape)","metadata":{"execution":{"iopub.status.busy":"2021-09-06T17:25:53.726934Z","iopub.execute_input":"2021-09-06T17:25:53.727185Z","iopub.status.idle":"2021-09-06T17:25:54.661753Z","shell.execute_reply.started":"2021-09-06T17:25:53.727162Z","shell.execute_reply":"2021-09-06T17:25:54.661116Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Extract FLAIR \ndf_flair_mgmt_on  = df_train[df_train[\"MGMT_value\"]== 1][['PatientID', 'flair', 'flair_axis']]\ndf_flair_mgmt_off = df_train[df_train[\"MGMT_value\"]== 0][['PatientID', 'flair', 'flair_axis']]\n\n# Extract t1w \ndf_t1w_mgmt_on  = df_train[df_train[\"MGMT_value\"]== 1][['PatientID', 't1w', 't1w_axis']]\ndf_t1w_mgmt_off = df_train[df_train[\"MGMT_value\"]== 0][['PatientID', 't1w', 't1w_axis']]\n\n# Extract t1wce\ndf_t1wce_mgmt_on  = df_train[df_train[\"MGMT_value\"]== 1][['PatientID', 't1wce', 't1wce_axis']]\ndf_t1wce_mgmt_off = df_train[df_train[\"MGMT_value\"]== 0][['PatientID', 't1wce', 't1wce_axis']]\n\n# Extract t2w\ndf_t2w_mgmt_on  = df_train[df_train[\"MGMT_value\"]== 1][['PatientID', 't2w', 't2w_axis']]\ndf_t2w_mgmt_off = df_train[df_train[\"MGMT_value\"]== 0][['PatientID', 't2w', 't2w_axis']]\n\n# Test slicing\ndf_flair_mgmt_on.head()","metadata":{"execution":{"iopub.status.busy":"2021-09-06T17:25:54.662636Z","iopub.execute_input":"2021-09-06T17:25:54.662962Z","iopub.status.idle":"2021-09-06T17:25:54.699959Z","shell.execute_reply.started":"2021-09-06T17:25:54.662937Z","shell.execute_reply":"2021-09-06T17:25:54.699074Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Cleaning up\ndf_flair_mgmt_on  = df_flair_mgmt_on.rename(columns={\"flair_axis\": \"axis\", 'flair': 'slices'})\ndf_flair_mgmt_off = df_flair_mgmt_off.rename(columns={\"flair_axis\": \"axis\", 'flair': 'slices'})\n\ndf_t1w_mgmt_on  = df_t1w_mgmt_on.rename(columns={\"t1w_axis\": \"axis\", 't1w': 'slices'})\ndf_t1w_mgmt_off = df_t1w_mgmt_off.rename(columns={\"t1w_axis\": \"axis\",'t1w': 'slices'})\n\ndf_t2w_mgmt_on  = df_t2w_mgmt_on.rename(columns={\"t2w_axis\": \"axis\", 't2w': 'slices'})\ndf_t2w_mgmt_off = df_t2w_mgmt_off.rename(columns={\"t2w_axis\": \"axis\", 't2w': 'slices'})\n\ndf_t1wce_mgmt_on  = df_t1wce_mgmt_on.rename(columns={\"t1wce_axis\": \"axis\", 't1wce': 'slices'})\ndf_t1wce_mgmt_off = df_t1wce_mgmt_off.rename(columns={\"t1wce_axis\": \"axis\",  't1wce': 'slices'})","metadata":{"execution":{"iopub.status.busy":"2021-09-06T17:25:54.701074Z","iopub.execute_input":"2021-09-06T17:25:54.701306Z","iopub.status.idle":"2021-09-06T17:25:54.711139Z","shell.execute_reply.started":"2021-09-06T17:25:54.701285Z","shell.execute_reply":"2021-09-06T17:25:54.710437Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, ax = plt.subplots(1,4, figsize=(30,10))\n\nsns.countplot(x='axis', data=df_flair_mgmt_on, ax=ax[0])\nsns.countplot(x='axis', data=df_t1w_mgmt_on, ax=ax[1])\nsns.countplot(x='axis', data=df_t2w_mgmt_on, ax=ax[2])\nsns.countplot(x='axis', data=df_t1wce_mgmt_on, ax=ax[3])\nplt.show()\n\nfig, ax = plt.subplots(1,4, figsize=(30,10))\nsns.countplot(x='axis', data=df_flair_mgmt_off, ax=ax[0])\nsns.countplot(x='axis', data=df_t1w_mgmt_off, ax=ax[1])\nsns.countplot(x='axis', data=df_t2w_mgmt_off, ax=ax[2])\nsns.countplot(x='axis', data=df_t1wce_mgmt_off, ax=ax[3])\n\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-09-06T17:25:54.712341Z","iopub.execute_input":"2021-09-06T17:25:54.712723Z","iopub.status.idle":"2021-09-06T17:25:55.522683Z","shell.execute_reply.started":"2021-09-06T17:25:54.712688Z","shell.execute_reply":"2021-09-06T17:25:55.522115Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_train[\"flair_arr\"] = df_train['flair'].apply(lambda slices_path : get_ordered_slices(slices_path))\ndf_train[\"t1w_arr\"]   = df_train['t1w'].apply(lambda slices_path : get_ordered_slices(slices_path))\ndf_train[\"t1wce_arr\"] = df_train['t1wce'].apply(lambda slices_path : get_ordered_slices(slices_path))\ndf_train[\"t2w_arr\"]   = df_train['t2w'].apply(lambda slices_path : get_ordered_slices(slices_path))","metadata":{"execution":{"iopub.status.busy":"2021-09-06T17:38:39.911956Z","iopub.execute_input":"2021-09-06T17:38:39.912382Z","iopub.status.idle":"2021-09-06T17:38:46.598356Z","shell.execute_reply.started":"2021-09-06T17:38:39.912348Z","shell.execute_reply":"2021-09-06T17:38:46.597186Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"entry = 11\nprint(df_train[\"flair\"].iloc[entry])\n\nprint(len(df_train[\"flair_arr\"].iloc[entry]))\nprint(len(df_train[\"t1w_arr\"].iloc[entry]))\nprint(len(df_train[\"t1wce_arr\"].iloc[entry]))\nprint(len(df_train[\"t2w_arr\"].iloc[entry]))\n\n\ndef get_len(data):\n    lens = []\n    for xi in data:\n        lens.append(len(xi))\n    return np.array(lens)\n    \n    \ndef plot_hist(np_arr,axd):\n    sns.histplot(np_arr, bins = [0,50,100,150,200,250, 300, 350, 400, 450, 500, 550, 600], ax=axd)\n    \n\ndf_train.head()\n# 240\n# 32\n# 360\n# 288\n\n\nfig, ax   = plt.subplots(1,4, figsize=(30,10))\nflair_arr = get_len(np.array(df_train[\"flair_arr\"]))\nt1w_arr   = get_len(np.array(df_train[\"t1w_arr\"]))\nt2w_arr   = get_len(np.array(df_train[\"t2w_arr\"]))\nt1wce_arr = get_len(np.array(df_train[\"t1wce_arr\"]))\n\nplot_hist(flair_arr, ax[0])\nplot_hist(t1w_arr, ax[1])\nplot_hist(t1wce_arr, ax[2])\nplot_hist(t2w_arr, ax[3])\n\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-09-06T20:03:46.063699Z","iopub.execute_input":"2021-09-06T20:03:46.064255Z","iopub.status.idle":"2021-09-06T20:03:46.853652Z","shell.execute_reply.started":"2021-09-06T20:03:46.064211Z","shell.execute_reply":"2021-09-06T20:03:46.852001Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_train.head(15)","metadata":{"execution":{"iopub.status.busy":"2021-09-06T19:59:28.361676Z","iopub.execute_input":"2021-09-06T19:59:28.362124Z","iopub.status.idle":"2021-09-06T19:59:28.486171Z","shell.execute_reply.started":"2021-09-06T19:59:28.362089Z","shell.execute_reply":"2021-09-06T19:59:28.484250Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# generate training images\n# an image contain all 4 modalities disposed in a square\n# df_mgmt_on  = [df_flair_mgmt_on,\n#             df_t1w_mgmt_on,\n#             df_t1wce_mgmt_on,\n#             df_t2w_mgmt_on,\n#            ]\n\n# df_mgmt_off = [df_flair_mgmt_off, \n#             df_t1w_mgmt_off, \n#             df_t1wce_mgmt_off, \n#             df_t2w_mgmt_off \n#            ]\n","metadata":{"execution":{"iopub.status.busy":"2021-09-06T17:25:55.557891Z","iopub.execute_input":"2021-09-06T17:25:55.558190Z","iopub.status.idle":"2021-09-06T17:25:55.562651Z","shell.execute_reply.started":"2021-09-06T17:25:55.558165Z","shell.execute_reply":"2021-09-06T17:25:55.561723Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def _circumscriber(img: np.array) -> np.array:\n    #First is vertical, second is horizontal, third is slices\n    vmin = 0\n    vlimit = img.shape[0]\n    hmin = 0\n    hlimit = img.shape[1]\n    \n    for i in range(vlimit):\n        if np.max(img[i, :]) == 0:\n            vmin += 1\n        else:\n            break\n    vmax = vmin + 1\n    for i in range(vmin+1, vlimit):\n        if np.max(img[i, :]) > 0:\n            vmax += 1\n        else:\n            break\n    \n    for j in range(hlimit):\n        if np.max(img[:, j]) == 0:\n            hmin += 1\n        else:\n            break\n    hmax = hmin + 1\n    for j in range(hmin+1, hlimit):\n        if np.max(img[:, j]) > 0:\n            hmax += 1\n        else:\n            break\n    return img[vmin:vmax, hmin:hmax]\n\n\n\ndef load_slices_meta(path: str):\n    # load the DICOM files\n    images = [os.path.basename(f) for f in glob.glob(f'{path}/*.dcm')]\n    # images.sort(key=lambda f: int(re.sub('\\D', '', f)))\n    images = [pydicom.read_file(\"{}/{}\".format(path,f)) for f in images]\n    slices = sorted(images, key=lambda s: s.SliceLocation)\n    try:\n        slice_thickness = np.abs(slices[0].ImagePositionPatient[2] - slices[1].ImagePositionPatient[2])\n    except:\n        slice_thickness = np.abs(slices[0].SliceLocation - slices[1].SliceLocation)\n        \n    for s in slices:\n        s.SliceThickness = slice_thickness\n  \n    return slices\n\n\n    \ndef _dicom_meta2array(dicom, voi_lut=True, fix_monochrome=True, resize=True, img_size=512):\n    data  = apply_voi_lut(dicom.pixel_array, dicom)\n    #data = dicom.pixel_array\n    if dicom.PhotometricInterpretation == \"MONOCHROME1\":\n        data = np.amax(data) - data\n    #Normalize the data: subtract off the minimum, divide by the maximum, convert to 256 uint8\n    data = data - np.min(data)\n    data = data/np.max(data)\n    data = (data * 255).astype(np.uint8)\n    # Resize images to target value\n    data = cv2.resize(data, (img_size, img_size))\n    data = cv2.equalizeHist(data)\n    return data\n\n\n\ndef reslice(slices):\n    # pixel aspects, assuming all slices are the same\n    # ps = slices[0].PixelSpacing\n    # ss = slices[0].SliceThickness\n    # ax_aspect  = ps[1]/ps[0]\n    # sag_aspect = ps[1]/ss\n    # cor_aspect = ss/ps[0]\n\n    # create 3D array\n    img_shape = list(_dicom_meta2array(slices[0]).shape)\n    img_shape.append(len(slices))\n    img3d = np.zeros(img_shape)\n\n    # fill 3D array with the images from the files\n    for i, s in enumerate(slices):\n        img2d = _dicom_meta2array(s)\n        img3d[:,:,i] = img2d\n    \n    print(img3d.shape)\n    # img3d = _circumscriber(img3d)\n    # print(img3d.shape)\n    \n    # plot 3 orthogonal slices\n    fig, ax = plt.subplots(1,3, figsize=(30,10))\n    \n    p_slice = img3d[:,:,img_shape[2]//3]\n    #p_slice = _circumscriber(p_slice)\n    print(p_slice.shape)\n    ax[0].imshow(p_slice)\n    #ax[0].set_aspect(ax_aspect)\n    \n    img_size = 512\n    #sag_aspect = sag_aspect * (img_shape[0]/img_shape[2] + 2)\n    p_slice = np.rot90(img3d[:,img_shape[1]//3,:])\n    p_slice = cv2.resize(p_slice, (img_size, img_size))\n    #p_slice = _circumscriber(p_slice)\n    # shape of 3d: (512, 512, 352)\n    ax[1].imshow(p_slice)\n    # ax[1].set_aspect(sag_aspect)\n    # if(len(slices) < 50):\n    #    ax[1].set_aspect(sag_aspect*20)\n    \n\n    # cor_aspect = cor_aspect * (img_shape[2]/img_shape[0])\n    p_slice = np.rot90(img3d[img_shape[0]//3,:,:])\n    p_slice = cv2.resize(p_slice, (img_size, img_size))\n    #p_slice = _circumscriber(p_slice)\n    print(p_slice.shape)\n    ax[2].imshow(p_slice)\n    #ax[2].set_aspect(cor_aspect)\n    #if(len(slices) < 50): \n    #    ax[2].set_aspect(cor_aspect*20)\n\n\nslices_path_t1wce = \"../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00014/T1wCE/\"\nslices_path_t1w   = \"../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00014/T1w/\"\nslices_path_flair = \"../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00014/FLAIR/\"\nslices_path_t2w   = \"../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00014/T2w/\"\n\n\nslices_t1wce = load_slices_meta(slices_path_t1wce)\nslices_flair = load_slices_meta(slices_path_flair)\nslices_t1w   = load_slices_meta(slices_path_t1w)\nslices_t2w   = load_slices_meta(slices_path_t2w)\n\n\nprint(len(slices_t1wce))\nprint(len(slices_flair))\nprint(len(slices_t1w))\nprint(len(slices_t2w))\n\n\n# reslice(slices_t1wce)\n# reslice(slices_t1w)\n# reslice(slices_t2w)\n# reslice(slices_flair)","metadata":{"execution":{"iopub.status.busy":"2021-09-06T17:25:55.564239Z","iopub.execute_input":"2021-09-06T17:25:55.564628Z","iopub.status.idle":"2021-09-06T17:26:08.862093Z","shell.execute_reply.started":"2021-09-06T17:25:55.564590Z","shell.execute_reply":"2021-09-06T17:26:08.861204Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import scipy\n\n\ndef resample(slices, new_spacing=[1,1,1]):\n    # Determine current pixel spacing\n    new_spacing = map(float, (new_spacing))\n    image = scipy.ndimage.interpolation.zoom(slices, new_spacing)\n    return image\n    \n    \ndef reslice_hu(slices):\n    # create 3D array\n    img_shape = list(slices.T.shape)\n    img_shape.append(len(slices))\n    img3d=np.zeros(slices.T.shape)\n    img_size = 512\n    # fill 3D array with the images from the files\n    for i, d_slice in enumerate(slices):\n        #print(d_slice.shape)\n        img3d[:,:,i] = d_slice\n\n    # plot 3 orthogonal slices\n    fig, ax = plt.subplots(1,3, figsize=(30,10))\n    ax[0].imshow(img3d[:,:,img_shape[2]//3])\n    \n    p_slice = np.rot90(img3d[:,img_shape[1]//3,:])\n    p_slice = cv2.resize(p_slice, (img_size, img_size))\n    ax[1].imshow(p_slice)\n   \n    p_slice = np.rot90(img3d[img_shape[0]//3,:,:])\n    p_slice = cv2.resize(p_slice, (img_size, img_size))\n    ax[2].imshow(p_slice)","metadata":{"execution":{"iopub.status.busy":"2021-09-06T17:33:19.867768Z","iopub.execute_input":"2021-09-06T17:33:19.868273Z","iopub.status.idle":"2021-09-06T17:33:19.882324Z","shell.execute_reply.started":"2021-09-06T17:33:19.868236Z","shell.execute_reply":"2021-09-06T17:33:19.881234Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"patient = \"00005\"\n\nslices_path_t1wce = df_train['t1wce'][df_train['PatientID'] == patient ].iloc[0]\nprint(slices_path_t1wce)\n#flair_images = load_FULL_brain_restricted(slices_path)\n#print(flair_images.shape)\n#plt.imshow(flair_images[:,:,7])\n\nslices_twice  = load_slices_meta(slices_path_t1wce)\nslices_pixels = get_pixels_hu(slices_twice)\n\nprint(slices_pixels.shape)\n#pix_resampled, spacing = resample(slices_pixels, slices, [4,1,1])\nprint(slices_pixels.shape)\n\n\n# reslice(slices_twice)\n# reslice_hu(slices_pixels)\n\npath_00014 = \"../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00014/FLAIR/\"\npath2 = \"../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00122/T1wCE/\"\n#boxtest1   = load_FULL_brain(path_00014)\n#print(boxtest1.shape)\n\n#boxtest2 = load_FULL_brain_restricted(path2)\n#print(boxtest2.shape)\n\n#plot_imgs(boxtest2[:,:,-30:-10])\ndf_train.head()\n","metadata":{"execution":{"iopub.status.busy":"2021-09-06T17:33:39.072094Z","iopub.execute_input":"2021-09-06T17:33:39.072444Z","iopub.status.idle":"2021-09-06T17:33:39.908491Z","shell.execute_reply.started":"2021-09-06T17:33:39.072415Z","shell.execute_reply":"2021-09-06T17:33:39.907485Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Create model","metadata":{}},{"cell_type":"code","source":"def create_3d_conv_layer(input_mod):\n    x = Conv3D(filters=32, kernel_size=3, activation=\"relu\",padding='same')(input_mod)\n    x = Conv3D(filters=32, kernel_size=3, activation=\"relu\",padding='same')(x)\n    x = MaxPool3D(pool_size=2)(x)\n    x = BatchNormalization()(x)\n\n    x = Conv3D(filters=64, kernel_size=3, activation=\"relu\")(x)\n    x = Conv3D(filters=64, kernel_size=3, activation=\"relu\")(x)\n    x = MaxPool3D(pool_size=2)(x)\n    x = BatchNormalization()(x)\n\n    x = Conv3D(filters=128, kernel_size=3, activation=\"relu\")(x)\n    x = Conv3D(filters=128, kernel_size=3, activation=\"relu\")(x)\n    x = MaxPool3D(pool_size=2)(x)\n    x = BatchNormalization()(x)\n\n    x = Conv3D(filters=256, kernel_size=3, activation=\"relu\", padding='same')(x)\n    x = MaxPool3D(pool_size=2, padding='same')(x)\n    x = BatchNormalization()(x)\n     \n    \n    x = GlobalAveragePooling3D()(x)\n    x = Dense(units=512, activation=\"relu\")(x)\n    x = Dropout(0.3)(x)\n    \n    return x\n\n    \n    \ndef get_model_3d(width=256, height=256, depth=32):\n    \"\"\"Build a 3D convolutional neural network model.\"\"\"\n\n    # inputs = keras.Input((width, height, depth, 1))\n    input_flair = Input((width, height, depth, 1))\n    input_t1w   = Input((width, height, depth, 1))\n    input_t1wce = Input((width, height, depth, 1))\n    input_t2w   = Input((width, height, depth, 1))\n    \n    # models\n    mod_flair = create_3d_conv_layer(input_flair)\n    mod_t1w   = create_3d_conv_layer(input_t1w)\n    mod_t1wce = create_3d_conv_layer(input_t1wce)    \n    mod_t2w   = create_3d_conv_layer(input_t2w)\n    \n    conv = concatenate([mod_flair, mod_t1w, mod_t1wce, mod_t2w])\n    \n    conv   = Flatten()(conv)\n    dense  = Dense(512)(conv)\n    dense  = LeakyReLU(alpha=0.1)(dense)\n    dense  = Dropout(0.5)(dense)\n    output = Dense(units=1, activation='sigmoid')(dense)\n\n    model  = Model(inputs=[input_flair, input_t1w, input_t1wce, input_t2w], outputs=output, name=\"3dcnn\")\n    return model\n\n\n\nfrom cloud_tpu_client import Client\n# c = Client()\n# c.configure_tpu_version(tf.__version__, restart_type=\"always\")\n\n#detect and init the TPU\ntpu = tf.distribute.cluster_resolver.TPUClusterResolver.connect()\nprint(\"Running on TPU \", tpu.cluster_spec().as_dict()[\"worker\"])\n\n# #instantiate a distribution strategy\ntpu_strategy = tf.distribute.experimental.TPUStrategy(tpu)\nwith tpu_strategy.scope():\n    # Build model.\n    model = get_model(width=256, height=256, depth=32)\n    # Compile model.\n    initial_learning_rate = 0.0001\n    lr_schedule = tf.keras.optimizers.schedules.ExponentialDecay(\n        initial_learning_rate, decay_steps=100000, decay_rate=0.96, staircase=True\n    )\n\n    model.compile(\n        loss=\"binary_crossentropy\",\n        optimizer=keras.optimizers.Adam(learning_rate=lr_schedule),\n        metrics=[\"acc\"],\n    )\n    \n    \nmodel.summary()\nfrom keras.utils.vis_utils import plot_model\nplot_model(model, to_file='model_plot.png', show_shapes=True, show_layer_names=True)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def create_3d_conv_layer(input_mod):\n    x = Conv3D(filters=32, kernel_size=3, activation=\"relu\",padding='same')(input_mod)\n    x = Conv3D(filters=32, kernel_size=3, activation=\"relu\",padding='same')(x)\n    x = MaxPool3D(pool_size=2)(x)\n    x = BatchNormalization()(x)\n\n    x = Conv3D(filters=64, kernel_size=3, activation=\"relu\")(x)\n    x = Conv3D(filters=64, kernel_size=3, activation=\"relu\")(x)\n    x = MaxPool3D(pool_size=2)(x)\n    x = BatchNormalization()(x)\n\n    x = Conv3D(filters=128, kernel_size=3, activation=\"relu\")(x)\n    x = Conv3D(filters=128, kernel_size=3, activation=\"relu\")(x)\n    x = MaxPool3D(pool_size=2)(x)\n    x = BatchNormalization()(x)\n\n    x = Conv3D(filters=256, kernel_size=3, activation=\"relu\", padding='same')(x)\n    x = MaxPool3D(pool_size=2, padding='same')(x)\n    x = BatchNormalization()(x)\n     \n    \n    x = GlobalAveragePooling3D()(x)\n    x = Dense(units=512, activation=\"relu\")(x)\n    x = Dropout(0.3)(x)\n    \n    return x\n\n    \n    \ndef get_model_3d(width=256, height=256, depth=32):\n    \"\"\"Build a 3D convolutional neural network model.\"\"\"\n\n    # inputs = keras.Input((width, height, depth, 1))\n    input_flair = Input((width, height, depth, 1))\n    input_t1w   = Input((width, height, depth, 1))\n    input_t1wce = Input((width, height, depth, 1))\n    input_t2w   = Input((width, height, depth, 1))\n    \n    # models\n    mod_flair = create_3d_conv_layer(input_flair)\n    mod_t1w   = create_3d_conv_layer(input_t1w)\n    mod_t1wce = create_3d_conv_layer(input_t1wce)    \n    mod_t2w   = create_3d_conv_layer(input_t2w)\n    \n    conv = concatenate([mod_flair, mod_t1w, mod_t1wce, mod_t2w])\n    \n    conv   = Flatten()(conv)\n    dense  = Dense(512)(conv)\n    dense  = LeakyReLU(alpha=0.1)(dense)\n    dense  = Dropout(0.5)(dense)\n    output = Dense(units=1, activation='sigmoid')(dense)\n\n    model  = Model(inputs=[input_flair, input_t1w, input_t1wce, input_t2w], outputs=output, name=\"3dcnn\")\n    return model\n\n\n\nfrom cloud_tpu_client import Client\n# c = Client()\n# c.configure_tpu_version(tf.__version__, restart_type=\"always\")\n\n#detect and init the TPU\ntpu = tf.distribute.cluster_resolver.TPUClusterResolver.connect()\nprint(\"Running on TPU \", tpu.cluster_spec().as_dict()[\"worker\"])\n\n# #instantiate a distribution strategy\ntpu_strategy = tf.distribute.experimental.TPUStrategy(tpu)\nwith tpu_strategy.scope():\n    # Build model.\n    model = get_model(width=256, height=256, depth=32)\n    # Compile model.\n    initial_learning_rate = 0.0001\n    lr_schedule = tf.keras.optimizers.schedules.ExponentialDecay(\n        initial_learning_rate, decay_steps=100000, decay_rate=0.96, staircase=True\n    )\n\n    model.compile(\n        loss=\"binary_crossentropy\",\n        optimizer=keras.optimizers.Adam(learning_rate=lr_schedule),\n        metrics=[\"acc\"],\n    )\n    \n    \nmodel.summary()\nfrom keras.utils.vis_utils import plot_model\nplot_model(model, to_file='model_plot.png', show_shapes=True, show_layer_names=True)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#print(df_train.head())\nprint(len(df_train.index))\ndf_train[(df_train['flair_axis'] == \"axial\") |  (df_train['t1w_axis'] == \"axial\") |  (df_train['t2w_axis'] == \"axial\") |  (df_train['t1wce_axis'] == \"axial\")]\n","metadata":{"execution":{"iopub.status.busy":"2021-09-06T04:42:40.043651Z","iopub.execute_input":"2021-09-06T04:42:40.044053Z","iopub.status.idle":"2021-09-06T04:42:40.070576Z","shell.execute_reply.started":"2021-09-06T04:42:40.044019Z","shell.execute_reply":"2021-09-06T04:42:40.069868Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"path_00014 = \"../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00014/FLAIR/\"\n#boxtest1   = load_FULL_brain(path_00014)\n#print(boxtest1.shape)\n\nt0= time.time()\nboxtest2 = load_FULL_brain_restricted(path_00014)\nt1 = time.time()\nprint(t1-t0)\n# returns  (512, 512, 82)","metadata":{"execution":{"iopub.status.busy":"2021-09-06T04:42:40.071794Z","iopub.execute_input":"2021-09-06T04:42:40.072081Z","iopub.status.idle":"2021-09-06T04:42:41.130574Z","shell.execute_reply.started":"2021-09-06T04:42:40.072053Z","shell.execute_reply":"2021-09-06T04:42:41.129504Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"entry = df_train.sample(1).iloc[0]\n\nentry[\"flair_axis\"]","metadata":{"execution":{"iopub.status.busy":"2021-09-06T04:42:41.132007Z","iopub.execute_input":"2021-09-06T04:42:41.132596Z","iopub.status.idle":"2021-09-06T04:42:41.140633Z","shell.execute_reply.started":"2021-09-06T04:42:41.132540Z","shell.execute_reply":"2021-09-06T04:42:41.139669Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"modalities = { 0 : \"flair\",\n               1 : \"t1w\",\n               2 : \"t1wce\",\n               3 : \"t2w\"\n             }\n\ndef train_preprocessing(volume, label):\n    \"\"\"Process training data by rotating and adding a channel.\"\"\"\n    # Rotate volume\n    #volume = rotate(volume)\n    volume = tf.expand_dims(volume, axis=3)\n    return volume, label\n\n\ndef validation_preprocessing(volume, label):\n    \"\"\"Process validation data by only adding a channel.\"\"\"\n    volume = tf.expand_dims(volume, axis=3)\n    return volume, label\n\n\ndef get_ordered_slices_subset(path:str, proportion:float, slices_count:int):\n    \"\"\"\n    get a list of ordered slices\n    \"\"\"\n    images = [os.path.basename(f) for f in glob.glob(f'{path}/*.dcm')]\n    images.sort(key=lambda f: int(re.sub('\\D', '', f)))\n    \n    len_imgs = len(images)\n    fifth    = len_imgs//5\n    # Most of images have a lot of empty space, get rid of 2/5 of them\n    images   = images[fifth:-fifth]\n    len_imgs = len(images)\n    \n    if len_imgs < slices_count:\n        return None\n    elif len_imgs == slices_count:\n        return images\n    \n    index_start = int(proportion * len_imgs)                   \n    if index_start + slices_count > len_imgs:\n        subset = images[-slices_count:]\n    else:\n        subset = images[index_start:index_start+slices_count]\n        \n    return subset\n\n\ndef is_bad_sample(selection):\n    bad_slices = [\"00584/T1w\", \"00569/T2w\", '00998/T1wCE', '00584/T1w', '00137/T1w', '00998/T2w', '00589/T1w' ]\n    for badness in bad_slices:\n            if (badness in selection):\n                return 1\n    return 0\n\n\ndef load_3d_volume(df, slices_count=32):\n    \"\"\"\n    send all of the images in the chosen modality, in order, as a single 3D np array\n    \"\"\"\n    coin  = np.random.choice([0,1])\n    proportion = random.random()\n    #modality   = np.random.choice(list(modalities.values()))\n    \n    \n    patient  = patients.sample(1).iloc[0]\n    mgmt_val = patient[\"MGMT_value\"].iloc[0]\n    \n    dcm_path_flair = patient[modalities[0]]\n    dcm_path_t1w   = patient[modalities[1]]\n    dcm_path_t1wce = patient[modalities[2]]\n    dcm_path_t2w   = patient[modalities[3]]\n    \n    images_flair = get_ordered_slices_subset(dcm_path_flair, proportion, slices_count)\n    images_t1w   = get_ordered_slices_subset(dcm_path_t1w, proportion, slices_count)\n    images_t1wce = get_ordered_slices_subset(dcm_path_t1wce, proportion, slices_count)\n    images_t2w   = get_ordered_slices_subset(dcm_path_t2w, proportion, slices_count)\n    \n    while images_flair == None or images_t1w == None or images_t1wce == None or images_t2w == None:\n        patient  = patients.sample(1).iloc[0]\n        mgmt_val = patient[\"MGMT_value\"].iloc[0]\n        \n        dcm_path_flair = patient[modalities[0]]\n        dcm_path_t1w   = patient[modalities[1]]\n        dcm_path_t1wce = patient[modalities[2]]\n        dcm_path_t2w   = patient[modalities[3]]\n\n        images_flair = get_ordered_slices_subset(dcm_path_flair, proportion, slices_count)\n        images_t1w   = get_ordered_slices_subset(dcm_path_t1w, proportion, slices_count)\n        images_t1wce = get_ordered_slices_subset(dcm_path_t1wce, proportion, slices_count)\n        images_t2w   = get_ordered_slices_subset(dcm_path_t2w, proportion, slices_count)\n        \n    \n    # Load meta from files\n    pixels_flair = [_dicom2array(dcm_path_flair + f) for f in images_flair]\n    pixels_t1w   = [_dicom2array(dcm_path_t1w + f) for f in images_t1w]\n    pixels_t1wce = [_dicom2array(dcm_path_t1wce + f) for f in images_t1wce]\n    pixels_t2w   = [_dicom2array(dcm_path_t2w + f) for f in images_t2w]\n    \n    # np.count_nonzero(anchor) < 10000)\n    vol_flair = np.array([im for im in pixels_flair]).T\n    vol_flair = tf.convert_to_tensor(tf.expand_dims(vol_flair, axis=3))\n    \n    vol_t1w   = np.array([im for im in pixels_t1w]).T\n    vol_t1w   = tf.convert_to_tensor(tf.expand_dims(vol_t1w, axis=3))\n    \n    vol_t1wce = np.array([im for im in pixels_t1wce]).T\n    vol_t1wce = tf.convert_to_tensor(tf.expand_dims(vol_t1wce, axis=3))\n    \n    vol_t2w   = np.array([im for im in pixels_t2w]).T\n    vol_t2w   = tf.convert_to_tensor(tf.expand_dims(vol_t2w, axis=3))\n    \n    \n    #final_image = _circumscriber(good_images)\n    volumes = [vol_flair,vol_t1w, vol_t1wce, vol_t2w]\n    return volumes,  mgmt_val\n\n\ndef generate_batches(df, total_items):\n    #patients_axial = df_train[(df_train['flair_axis'] == \"axial\") &  (df_train['t1w_axis'] == \"axial\") &  (df_train['t2w_axis'] == \"axial\") &  (df_train['t1wce_axis'] == \"axial\")]\n    slices_count = 32\n    i = 0\n    while True:\n        volumes, label = load_3d_volume(df, slices_count)\n        yield volumes, tf.convert_to_tensor(label)\n        #i += 1\n    \n\n    \ndef generate_batches_2(df, total_items):\n    #patients_axial = df_train[(df_train['flair_axis'] == \"axial\") &  (df_train['t1w_axis'] == \"axial\") &  (df_train['t2w_axis'] == \"axial\") &  (df_train['t1wce_axis'] == \"axial\")]\n    slices_count = 32\n    i            = 0\n    flair  = np.zeros((total_items, 256, 256, slices_count, 1))\n    t1w    = np.zeros((total_items, 256, 256, slices_count, 1))\n    t1wce  = np.zeros((total_items, 256, 256, slices_count, 1))\n    t2w    = np.zeros((total_items, 256, 256, slices_count, 1))\n    labels = np.zeros((total_items, ))\n    \n    for i in range(total_items):\n        volumes, label = load_3d_volume(df, slices_count)\n        flair[i]  = volumes[0]\n        t1w[i]    = volumes[1]\n        t1wce[i]  = volumes[2]\n        t2w[i]    = volumes[3]\n        labels[i] = label\n    return [flair, t1w, t1wce, t2w], labels\n\n        \n# TEST        \npatient_entry = df_train.iloc[1]\naxis     = \"flair_axis\"\nmodality = \"flair\"\n\nt0 = time.time()\nsample_generator = generate_batches_2(df_test, 4)\nsample_test = sample_generator\nt1 = time.time()\nprint(t1-t0)\nprint(\"--------\")\n\nprint(len(sample_test[0]))\nprint(len(sample_test))\n\nprint(sample_test[0][0].shape)\nprint(sample_test[0][1].shape)\nprint(sample_test[0][2].shape)\nprint(sample_test[0][3].shape)\n\n\nprint(sample_test[1].shape)\n#print(sample_test[0])\nprint(sample_test[1])\n\n","metadata":{"execution":{"iopub.status.busy":"2021-09-06T04:42:41.141767Z","iopub.execute_input":"2021-09-06T04:42:41.142127Z","iopub.status.idle":"2021-09-06T04:42:47.036600Z","shell.execute_reply.started":"2021-09-06T04:42:41.142093Z","shell.execute_reply":"2021-09-06T04:42:47.035728Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# TEST  \nx = model.predict(sample_test[0])\n","metadata":{"execution":{"iopub.status.busy":"2021-09-06T04:42:47.037960Z","iopub.execute_input":"2021-09-06T04:42:47.038365Z","iopub.status.idle":"2021-09-06T04:43:03.520464Z","shell.execute_reply.started":"2021-09-06T04:42:47.038322Z","shell.execute_reply":"2021-09-06T04:43:03.519510Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#train_generator = tf.data.Dataset.from_generator(generate_batches_train, output_types={ \"slices\":tf.uint8, \"cls_lbl\":tf.uint8})\n# validation_generator = tf.data.Dataset.from_generator(generate_batches_test, output_types={ \"slices\":tf.uint8, \"cls_lbl\":tf.uint8})\n\n# train_generator      = next(generate_batches_train())\n# validation_generator = next(generate_batches_test()) \n\n\n#@tf.function\ndef input_fn(df, batch_size, epochs):\n    features_shape = [256, 256, 32, 1]\n    labels_shape   = [] \n    dataset = tf.data.Dataset.from_generator(lambda: generate_batches(df, batch_size),\n                                             output_types=((tf.uint8,\n                                                            tf.uint8,\n                                                            tf.uint8,\n                                                            tf.uint8), tf.int64),\n                                             output_shapes=((tf.TensorShape(features_shape),\n                                                             tf.TensorShape(features_shape),\n                                                             tf.TensorShape(features_shape),\n                                                             tf.TensorShape(features_shape)), tf.TensorShape(labels_shape))\n                                            )\n                                            \n\n    #dataset = dataset.repeat(epochs)\n    dataset = dataset.batch(batch_size)\n    dataset = dataset.prefetch(2)\n    return dataset\n\n\n\ndef input_fn_2(df, batch_size):\n    sample_generator = generate_batches_2(df, batch_size)\n                                    \n    dataset = tf.data.Dataset.from_tensor_slices(sample_generator)\n    #dataset = dataset.batch(batch_size)\n    dataset = dataset.prefetch(2)\n    return dataset\n\n\n# TEST\n# batch_size = 8\n# epochs    = 10\n# train_generator = input_fn_2(df_train, batch_size)\n# data = train_generator.take(1)\n# images, labels = list(data)[0]\n# print(labels)\n# print(images.shape)\n\n# for next_element in train_generator:\n#     tf.print(next_element[1].shape)\n#     tf.print(next_element[0].shape)    ","metadata":{"execution":{"iopub.status.busy":"2021-09-06T04:43:03.521565Z","iopub.execute_input":"2021-09-06T04:43:03.521828Z","iopub.status.idle":"2021-09-06T04:43:03.528880Z","shell.execute_reply.started":"2021-09-06T04:43:03.521803Z","shell.execute_reply":"2021-09-06T04:43:03.528010Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\n# Define callbacks.\nload_locally = tf.saved_model.LoadOptions(experimental_io_device='/job:localhost')\ncheckpoint_cb = tf.keras.callbacks.ModelCheckpoint(\n    \"3d_image_classification.h5\", save_best_only=True\n)\nearly_stopping_cb = tf.keras.callbacks.EarlyStopping(monitor=\"val_acc\", patience=15)\n\n\n# Train the model, doing validation at the end of each epoch\n\nepochs      = 10\nbatch_size  = 16\nnum_batches = len(df_train.index)//batch_size * 2\n\n# epochs = epochs,\n# steps_per_epoch = num_batches,\n# train_generator      = input_fn(df_train, batch_size, epochs)\n# validation_generator = input_fn(df_test, batch_size, epochs)\n\n\n#_train_generator      = generate_batches_2(df_train, batch_size)\n#_validation_generator = generate_batches_2(df_test, batch_size)\n#print(_train_generator)\n\n# model.fit(_train_generator,\n#     validation_data=_validation_generator,\n#     epochs=epochs,\n#     steps_per_epoch = num_batches,\n#     validation_steps = num_batches,\n#     shuffle=True,\n#     verbose=1,\n#     callbacks=[checkpoint_cb,early_stopping_cb])\nlosses     = []\nval_losses = []\nstop = 5\nfor i in tqdm(range(10000)):\n    X_train, y_train = generate_batches_2(df_train, batch_size)\n    X_test, y_test   = generate_batches_2(df_test, batch_size)\n    \n    h_callback = model.fit(X_train, y_train,\n        validation_data=(X_test, y_test),\n        epochs=epochs,\n        shuffle=True,\n        verbose=1,\n        callbacks=[checkpoint_cb,early_stopping_cb])\n    losses     = losses     + h_callback.history['loss']\n    val_losses = val_losses + h_callback.history['val_loss']\n    if i == stop:\n        break\n    \n    \n    \n# 1\n# 144/144 [==============================] -                loss: 0.6855 - acc: 0.5620 - val_loss: 0.6783 - val_acc: 0.5894\n# 144/144 [==============================] - 728s 5s/step - loss: 0.6856 - acc: 0.5675 - val_loss: 0.9366 - val_acc: 0.5043\n# 144/144 [==============================] - 729s 5s/step - loss: 0.6770 - acc: 0.5692 - val_loss: 0.7183 - val_acc: 0.5694\n# 144/144 [==============================] - 720s 5s/step - loss: 0.6839 - acc: 0.5593 - val_loss: 0.6938 - val_acc: 0.5495\n# 144/144 [==============================] - 715s 5s/step - loss: 0.6919 - acc: 0.5324 - val_loss: 0.6973 - val_acc: 0.5000\n# 144/144 [==============================] - 730s 5s/step - loss: 0.6812 - acc: 0.5609 - val_loss: 0.7013 - val_acc: 0.5582\n\n\n#2","metadata":{"scrolled":true,"execution":{"iopub.status.busy":"2021-09-06T10:42:44.080704Z","iopub.execute_input":"2021-09-06T10:42:44.081121Z","iopub.status.idle":"2021-09-06T10:51:07.820180Z","shell.execute_reply.started":"2021-09-06T10:42:44.081078Z","shell.execute_reply":"2021-09-06T10:51:07.819121Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print( h_callback.history['loss'])\n\nplt.plot(losses)\nplt.xlabel('Epochs')\nplt.ylabel('Loss')\n","metadata":{"execution":{"iopub.status.busy":"2021-09-06T10:39:57.964863Z","iopub.execute_input":"2021-09-06T10:39:57.965332Z","iopub.status.idle":"2021-09-06T10:39:58.116573Z","shell.execute_reply.started":"2021-09-06T10:39:57.965291Z","shell.execute_reply":"2021-09-06T10:39:58.115494Z"},"trusted":true},"execution_count":null,"outputs":[]}]}