{"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":"markdown","source":"# Preprocessing of training studies\n\nThis notebook preprocesses all training studies by cropping the brain from MRI images and storing the remaining slices in tensors.","metadata":{}},{"cell_type":"code","source":"import os\nimport sys \nimport json\nimport glob\nimport random\nimport re\nimport collections\nimport time\n\nimport numpy as np\nimport pandas as pd\nimport pydicom\nimport cv2\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\nfrom skimage.color import rgb2gray\nfrom skimage import data\nfrom skimage.filters import gaussian\nfrom skimage.segmentation import active_contour\nfrom skimage import measure\n\nfrom sklearn.model_selection import StratifiedKFold\nfrom sklearn.metrics import roc_auc_score\n\nfrom tqdm import tqdm","metadata":{"papermill":{"duration":1.048295,"end_time":"2021-07-14T20:26:46.309722","exception":false,"start_time":"2021-07-14T20:26:45.261427","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2021-09-12T18:18:30.123364Z","iopub.execute_input":"2021-09-12T18:18:30.123855Z","iopub.status.idle":"2021-09-12T18:18:33.884915Z","shell.execute_reply.started":"2021-09-12T18:18:30.123741Z","shell.execute_reply":"2021-09-12T18:18:33.883814Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pip install --upgrade ../input/mclahenumpy/mclahe-numpy","metadata":{"execution":{"iopub.status.busy":"2021-09-12T18:18:33.886527Z","iopub.execute_input":"2021-09-12T18:18:33.886837Z","iopub.status.idle":"2021-09-12T18:19:04.331473Z","shell.execute_reply.started":"2021-09-12T18:18:33.886805Z","shell.execute_reply":"2021-09-12T18:19:04.330307Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from mclahe import mclahe","metadata":{"execution":{"iopub.status.busy":"2021-09-12T18:19:04.333792Z","iopub.execute_input":"2021-09-12T18:19:04.334098Z","iopub.status.idle":"2021-09-12T18:19:04.342449Z","shell.execute_reply.started":"2021-09-12T18:19:04.334063Z","shell.execute_reply":"2021-09-12T18:19:04.3414Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"All MRI studies are stored in tensors of size (SIZE, SIZE, min(NUM_IMAGES, len(slices))\n\nParameter FIRST_HALF indicates whether first part of the dataset should be processed. In case of False, second half is prepared. This is to overcome the problem of a limited output data size (20gb).","metadata":{}},{"cell_type":"code","source":"mri_types = ['FLAIR', 'T1w', 'T1wCE', 'T2w']\nSIZE = 160\nNUM_IMAGES = 160\nFIRST_HALF = True","metadata":{"lines_to_end_of_cell_marker":2,"lines_to_next_cell":2,"papermill":{"duration":0.05565,"end_time":"2021-07-14T20:26:46.486521","exception":false,"start_time":"2021-07-14T20:26:46.430871","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2021-09-12T18:19:04.343935Z","iopub.execute_input":"2021-09-12T18:19:04.344285Z","iopub.status.idle":"2021-09-12T18:19:04.358487Z","shell.execute_reply.started":"2021-09-12T18:19:04.344253Z","shell.execute_reply":"2021-09-12T18:19:04.357292Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"output_dir = '/kaggle/tmp/dataset'\ndata_directory = '/kaggle/input/rsna-miccai-brain-tumor-radiogenomic-classification'","metadata":{"execution":{"iopub.status.busy":"2021-09-12T18:19:04.359631Z","iopub.execute_input":"2021-09-12T18:19:04.360028Z","iopub.status.idle":"2021-09-12T18:19:04.369528Z","shell.execute_reply.started":"2021-09-12T18:19:04.359995Z","shell.execute_reply":"2021-09-12T18:19:04.368687Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!mkdir -p {output_dir}/train","metadata":{"execution":{"iopub.status.busy":"2021-09-12T18:19:04.370527Z","iopub.execute_input":"2021-09-12T18:19:04.370793Z","iopub.status.idle":"2021-09-12T18:19:05.120845Z","shell.execute_reply.started":"2021-09-12T18:19:04.370767Z","shell.execute_reply":"2021-09-12T18:19:05.119403Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Functions to load images","metadata":{}},{"cell_type":"code","source":"def load_file(path):\n    return pydicom.read_file(path).pixel_array / 255.\n\n\ndef crop_dicom_image(data, crop_coordinates, img_size=SIZE):\n    dim0_min = 0 if crop_coordinates[0] < 0 or crop_coordinates[0] >= data.shape[0]-1 else crop_coordinates[0]\n    dim0_max = data.shape[0] if crop_coordinates[1] - dim0_min <= 0 or crop_coordinates[1] - dim0_min > data.shape[0] else crop_coordinates[1]\n    dim1_min = 0 if crop_coordinates[2] < 0 or crop_coordinates[2] >= data.shape[1]-1 else crop_coordinates[2]\n    dim1_max = data.shape[1] if crop_coordinates[3] - dim1_min <= 0 or crop_coordinates[3] - dim0_min > data.shape[0] else crop_coordinates[3]\n    data = data[dim0_min:dim0_max, dim1_min:dim1_max]\n    data = cv2.resize(data, (img_size, img_size))\n\n    return data\n\n\ndef natural_sort(l): \n    convert = lambda text: int(text) if text.isdigit() else text.lower()\n    alphanum_key = lambda key: [convert(c) for c in re.split('([0-9]+)', key)]\n    return sorted(l, key=alphanum_key)\n\n\ndef normalise(f):\n    f = mclahe(f)\n    return f\n\n\ndef find_biggest_snake(data_list, equal_to_prev, patience):\n    current_index = len(data_list) // 2\n    data = data_list[current_index]\n    s = np.linspace(0, 2 * np.pi, max(data.shape[0], data.shape[1]))\n    r = data.shape[0] / 2 + data.shape[0] / 2 * np.sin(s)\n    c = data.shape[1] / 2 + data.shape[1] / 2 * np.cos(s)\n    initial_snake = np.array([r, c]).T\n    \n    left_coordinates = find_biggest_snake_recursively(initial_snake, data_list, equal_to_prev, current_index, None, -1, patience=patience)\n    return find_biggest_snake_recursively(initial_snake, data_list, equal_to_prev, current_index + 1, left_coordinates, 1, patience=patience)\n    \n    \ndef find_biggest_snake_recursively(initial_snake, data_list, equal_to_prev, current_index, current_coordinates, index_diff, patience):\n    #print(\"find_biggest_snake_recursively\", current_index, current_coordinates, index_diff)\n    \n    if current_index < 0 or current_index >= len(data_list):\n        return current_coordinates\n    \n    if patience == 0:\n        #print(\"end of patience\")\n        return current_coordinates\n    \n    new_coordinates = [] \n    found_new_coordinates = False\n    \n    if current_coordinates is not None:\n        if index_diff == 1 and equal_to_prev[current_index]:\n            #print(\"found equal for diff 1\")\n            found_new_coordinates = True\n            new_coordinates = [coord for coord in current_coordinates]\n        elif index_diff == -1 and current_index < len(data_list) - 1 and equal_to_prev[current_index+1]:\n            found_new_coordinates = True\n            new_coordinates = [coord for coord in current_coordinates]\n            #print(\"found equal for diff -1\")\n\n    if not found_new_coordinates:\n        # Find contours at a constant value of 0.5\n        contours = measure.find_contours(data_list[current_index], 0.5)\n        \n        snake_min = [np.inf, np.inf]\n        snake_max = [0, 0]\n        for c in contours:\n            c_min = np.min(c, axis=0)\n            for i in [0, 1]:\n                if int(c_min[i]) < snake_min[i]:\n                    snake_min[i] = int(c_min[i])\n            c_max = np.max(c, axis=0)\n            for i in [0, 1]:\n                if int(c_max[i]) > snake_max[i]:\n                    snake_max[i] = int(c_max[i])\n        \n        coordinates = (snake_min[0], snake_max[0], snake_min[1], snake_max[1])\n\n        for i in range(4):\n            if (current_coordinates is None) or (i % 2 == 0 and coordinates[i] < current_coordinates[i]) or (i % 2 == 1 and coordinates[i] > current_coordinates[i]):\n                new_coordinates.append(coordinates[i])\n                found_new_coordinates = True\n            else:\n                new_coordinates.append(current_coordinates[i])\n    \n    if not found_new_coordinates:\n        patience -= 1\n        \n    return find_biggest_snake_recursively(initial_snake, data_list, equal_to_prev, current_index + index_diff, new_coordinates, index_diff, patience)\n\n\ndef load_dicom_images_3d_prime(scan_id, num_imgs=NUM_IMAGES, img_size=SIZE, mri_type=\"FLAIR\", split=\"train\"):\n    files = natural_sort(glob.glob(f\"{data_directory}/{split}/{scan_id}/{mri_type}/*.dcm\"))\n    \n    # excluding blank slices\n    for index in [0, -1]:\n        blank_slice_found = True\n        removed_blank = 0\n        while blank_slice_found:\n            data = load_file(files[index])\n            if np.count_nonzero(data - np.min(data)) / (data.shape[0] * data.shape[1]) < 0.003:\n                blank_slice_found = True\n                removed_blank += 1\n                del files[index]\n            else:\n                blank_slice_found = False\n    \n    if len(files) > num_imgs:\n        every_nth = len(files) / num_imgs\n        indexes = [min(int(round(i*every_nth)), len(files)-1) for i in range(0,num_imgs)]\n    else:\n        every_nth = 1\n        indexes = [i for i in range(0, len(files))]\n    \n    unique_indices = np.unique(np.array(indexes))\n    \n    files_to_load = [files[i] for i in indexes]\n    \n    images = []\n    equal_to_prev = []\n    prev_file = ''\n    for f in files_to_load:\n        if f == prev_file:\n            images.append(images[-1])\n            equal_to_prev.append(True)\n        else:\n            image = load_file(f)\n            images.append(image)\n            equal_to_prev.append(False)\n        \n        prev_file = f\n        \n    patience = int(every_nth * 2) + 1\n        \n    crop_coordinates = find_biggest_snake(images, equal_to_prev, patience)\n    \n    for i in range(len(images)):\n        if equal_to_prev[i]:\n            images[i] = images[i - 1]\n        else:\n            images[i] = crop_dicom_image(images[i], crop_coordinates)\n    \n    img3d = np.stack(images).T\n    \n    img3d = normalise(img3d)\n    \n    return np.expand_dims(img3d,0)","metadata":{"papermill":{"duration":0.035761,"end_time":"2021-07-14T20:26:46.726756","exception":false,"start_time":"2021-07-14T20:26:46.690995","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2021-09-12T18:19:05.122694Z","iopub.execute_input":"2021-09-12T18:19:05.123053Z","iopub.status.idle":"2021-09-12T18:19:05.165899Z","shell.execute_reply.started":"2021-09-12T18:19:05.123013Z","shell.execute_reply":"2021-09-12T18:19:05.164419Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#load_dicom_images_3d_prime(\"00113\", mri_type=\"FLAIR\")","metadata":{"execution":{"iopub.status.busy":"2021-09-12T18:19:05.172497Z","iopub.execute_input":"2021-09-12T18:19:05.172941Z","iopub.status.idle":"2021-09-12T18:19:05.177787Z","shell.execute_reply.started":"2021-09-12T18:19:05.172893Z","shell.execute_reply":"2021-09-12T18:19:05.176692Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"samples_to_exclude = [109, 123, 709]\n\ndf = pd.read_csv(f\"{data_directory}/train_labels.csv\")\nprint(\"original shape\", df.shape)\ndf = df[~df.BraTS21ID.isin(samples_to_exclude)]\nprint(\"new shape\", df.shape)\ndisplay(df)","metadata":{"execution":{"iopub.status.busy":"2021-09-12T18:19:05.179461Z","iopub.execute_input":"2021-09-12T18:19:05.179828Z","iopub.status.idle":"2021-09-12T18:19:05.232208Z","shell.execute_reply.started":"2021-09-12T18:19:05.179796Z","shell.execute_reply":"2021-09-12T18:19:05.231188Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"split = 'train'\n\n\ndf = df.iloc[0:df.shape[0] // 2] if FIRST_HALF else df.iloc[df.shape[0] // 2:]\n\nfor index, row in tqdm(df.iterrows(), total=df.shape[0]):\n    scan_id = str(row['BraTS21ID']).zfill(5)\n    out_dir = f'{output_dir}/{split}/{scan_id}'\n    os.makedirs(out_dir)\n    for mri_type in mri_types:\n        np.save(f'{out_dir}/{mri_type}.npy', load_dicom_images_3d_prime(scan_id, mri_type=mri_type))","metadata":{"execution":{"iopub.status.busy":"2021-09-12T18:19:05.23361Z","iopub.execute_input":"2021-09-12T18:19:05.233933Z","iopub.status.idle":"2021-09-12T18:19:31.251352Z","shell.execute_reply.started":"2021-09-12T18:19:05.233901Z","shell.execute_reply":"2021-09-12T18:19:31.249516Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!zip -r /kaggle/working/rsna-miccai-brain-tumor.zip {output_dir}","metadata":{"execution":{"iopub.status.busy":"2021-09-12T18:19:31.252349Z","iopub.status.idle":"2021-09-12T18:19:31.252805Z"},"trusted":true},"execution_count":null,"outputs":[]}]}