{"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\nimport pydicom\nimport numpy as np \nimport cv2\nimport os \nimport shutil\nimport pandas as pd\nfrom tqdm import tqdm\nimport seaborn as sns\nimport matplotlib.pyplot as plt\nimport torch\nimport torch.nn.functional as F\n\nCONFIG = {\n    'REMOVE_ZERO_SLICES' : True,\n    'INPUT_DIM' : (128,128,128)\n}\n\nTRAIN_DIR = '../input/rsna-miccai-brain-tumor-radiogenomic-classification/train'\nTEST_DIR = '../input/rsna-miccai-brain-tumor-radiogenomic-classification/test'","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2021-08-20T01:52:18.414476Z","iopub.execute_input":"2021-08-20T01:52:18.415073Z","iopub.status.idle":"2021-08-20T01:52:21.176835Z","shell.execute_reply.started":"2021-08-20T01:52:18.414973Z","shell.execute_reply":"2021-08-20T01:52:21.175965Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def trim3d(arr_3d, plot=False):\n    slices = arr_3d.shape[2]\n    largest_area = 0\n    largest_idx = 0\n    x,y,w,h = 0,0,arr_3d.shape[0],arr_3d.shape[1]\n    for i in range(slices):\n        s = arr_3d[:,:,i]\n        if np.sum(s) > 0.0:\n        \n            _,thresh = cv2.threshold(s,0,255,cv2.THRESH_BINARY)\n            contours,hierarchy = cv2.findContours(thresh,cv2.RETR_EXTERNAL,cv2.CHAIN_APPROX_SIMPLE)\n            areas = [cv2.contourArea(c) for c in contours]\n            max_idx = np.argmax(areas)\n\n            if areas[max_idx] >= largest_area:\n                largest_area = areas[max_idx]\n                largest_idx = i\n                x,y,w,h = cv2.boundingRect(contours[max_idx])\n    \n        \n    result = arr_3d[y:y+h, x:x+w, :]\n    \n    if plot:\n        fig, axes = plt.subplots(nrows=1, ncols=2)\n        axes[0].imshow(arr_3d[:,:,largest_idx])\n        axes[1].imshow(result[:,:,largest_idx])\n        plt.show()\n    \n    return result\n\ndef dicom_to_img(ds):\n    data = pydicom.pixel_data_handlers.apply_voi_lut(ds.pixel_array, ds)\n    if ds.PhotometricInterpretation == \"MONOCHROME1\":\n        data = np.amax(data) - data\n    data = data - np.min(data)\n    data = data / np.max(data)\n    data = (data * 255).astype(np.uint8)\n    return data\n\ndef same_size(img):\n    img = img.astype(np.float32) \n    img = torch.from_numpy(img) \n    img = img.unsqueeze(0).unsqueeze(0) \n    img = F.interpolate(img, size=CONFIG['INPUT_DIM'], mode='trilinear', align_corners=True) \n    img = img.type(torch.uint8) \n    img = img.squeeze(0).squeeze(0)\n    img = img.numpy()\n    return img","metadata":{"execution":{"iopub.status.busy":"2021-08-20T01:52:21.178431Z","iopub.execute_input":"2021-08-20T01:52:21.178830Z","iopub.status.idle":"2021-08-20T01:52:21.194250Z","shell.execute_reply.started":"2021-08-20T01:52:21.178790Z","shell.execute_reply":"2021-08-20T01:52:21.192870Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Look at Data Dimensions","metadata":{}},{"cell_type":"code","source":"flair_df, t1w_df, t2w_df, t1wCE_df = [],[],[],[]\n\npatients = os.listdir(TRAIN_DIR)\nfor patient in tqdm(patients):\n    patient_dir = os.path.join(TRAIN_DIR, patient)\n    img_types = os.listdir(patient_dir)\n    for im_type in img_types:\n        files = os.listdir(os.path.join(patient_dir, im_type)) #input\n        \n        ds = pydicom.dcmread(os.path.join(patient_dir, im_type, files[0]))\n        \n        entry = [patient, ds.pixel_array.shape[0], ds.pixel_array.shape[1], len(files)]\n        if im_type == 'FLAIR':\n            flair_df.append(entry)\n        elif im_type == 'T1wCE':\n            t1wCE_df.append(entry)\n        elif im_type == 'T1w':\n            t1w_df.append(entry)\n        else:\n            t2w_df.append(entry)\n            \nflair = pd.DataFrame(flair_df, columns=[\"patient\",\"height\", \"width\", \"depth\"])\nt1w = pd.DataFrame(t1w_df, columns=[\"patient\",\"height\", \"width\", \"depth\"])\nt1wCE = pd.DataFrame(t1wCE_df, columns=[\"patient\",\"height\", \"width\", \"depth\"])\nt2w = pd.DataFrame(t2w_df, columns=[\"patient\",\"height\", \"width\", \"depth\"])","metadata":{"execution":{"iopub.status.busy":"2021-08-20T01:52:21.196843Z","iopub.execute_input":"2021-08-20T01:52:21.197399Z","iopub.status.idle":"2021-08-20T01:52:21.527293Z","shell.execute_reply.started":"2021-08-20T01:52:21.197220Z","shell.execute_reply":"2021-08-20T01:52:21.525822Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"datalist = [flair, t1w, t1wCE, t2w]\nnames = [\"FLAIR\", \"T1w\", \"T1wCE\", \"T2w\"]\nfig, axes = plt.subplots(nrows=4, ncols=3, figsize=(15,10))\nfor i, x in enumerate(datalist):\n    sns.histplot(data=x['height'],binwidth=10, ax=axes[i, 0], multiple='stack')\n    axes[i,0].set_title(f'{names[i]} Height')\n    sns.histplot(data=x['width'],binwidth=10, ax=axes[i, 1], multiple='stack')\n    axes[i,1].set_title(f'{names[i]} Width')\n    sns.histplot(data=x['depth'],binwidth=10, ax=axes[i, 2], multiple='stack')\n    axes[i,2].set_title(f'{names[i]} Depth')\nplt.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-08-20T01:52:21.528872Z","iopub.execute_input":"2021-08-20T01:52:21.529162Z","iopub.status.idle":"2021-08-20T01:52:23.153435Z","shell.execute_reply.started":"2021-08-20T01:52:21.529136Z","shell.execute_reply":"2021-08-20T01:52:23.152431Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Crop Images to Object of Interest","metadata":{}},{"cell_type":"code","source":"save_dir = os.getcwd()\nsave_train_dir = os.path.join(save_dir, 'train')\nsave_test_dir = os.path.join(save_dir, 'test')\n\nif not os.path.exists(save_train_dir):\n    os.makedirs(save_train_dir)\n\nif not os.path.exists(save_test_dir):\n    os.makedirs(save_test_dir)","metadata":{"execution":{"iopub.status.busy":"2021-08-20T01:52:23.154669Z","iopub.execute_input":"2021-08-20T01:52:23.154971Z","iopub.status.idle":"2021-08-20T01:52:23.162007Z","shell.execute_reply.started":"2021-08-20T01:52:23.154944Z","shell.execute_reply":"2021-08-20T01:52:23.160072Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def combine_scans(input_dir, output_dir):\n    flair_df, t1w_df, t2w_df, t1wCE_df = [],[],[],[]\n    patients = os.listdir(input_dir)\n    for patient in tqdm(patients):\n        img_types = os.listdir(os.path.join(input_dir, patient))\n        img_types.sort() #FFLAIR -> T1w -> T1wCE -> T2w\n        patient_dir = os.path.join(output_dir, patient)\n        if not os.path.exists(patient_dir):\n            os.makedirs(patient_dir) #output\n            \n        # For each image type\n\n        for i,img_type in enumerate(img_types):\n            img_dir = os.path.join(input_dir, patient, img_type) #input\n            img_files = os.listdir(img_dir)\n            img_files = [os.path.join(img_dir, x) for x in img_files]\n            \n            ds = [pydicom.dcmread(x) for x in img_files]\n            ds.sort(key = lambda x: int(x.ImagePositionPatient[2]))\n            imgs = [dicom_to_img(x) for x in ds]\n            if CONFIG['REMOVE_ZERO_SLICES']:\n                imgs = [x for x in imgs if np.sum(x) > 0]\n                \n            if len(imgs) > 0: # if there are only slices of blanks. \n                img3d = np.stack(imgs, axis=-1)\n                if np.sum(img3d) > 0:\n                    img3d = trim3d(img3d)\n                    # img3d = same_size(img3d)\n                    entry = [patient, img3d.shape[0], img3d.shape[1], img3d.shape[2]]\n                    \n                    if img_type == 'FLAIR':\n                        flair_df.append(entry)\n                    elif img_type == 'T1wCE':\n                        t1wCE_df.append(entry)\n                    elif img_type == 'T1w':\n                        t1w_df.append(entry)\n                    else:\n                        t2w_df.append(entry)\n                        \n                    save_img_dir = os.path.join(patient_dir, img_type)\n                    if not os.path.exists(save_img_dir):\n                        os.makedirs(save_img_dir)\n                    file_path = os.path.join(save_img_dir, f'{patient}_{img_type}_Image.npz')\n                    np.savez_compressed(file_path, img3d)\n                else:\n                    print(f'All Blanks for {patient} - {img_type}')\n            else:\n                print(f\"No Images for patient {patient} - {img_type}\")\n    \n    flair = pd.DataFrame(flair_df, columns=[\"patient\",\"height\", \"width\", \"depth\"])\n    t1w = pd.DataFrame(t1w_df, columns=[\"patient\",\"height\", \"width\", \"depth\"])\n    t1wCE = pd.DataFrame(t1wCE_df, columns=[\"patient\",\"height\", \"width\", \"depth\"])\n    t2w = pd.DataFrame(t2w_df, columns=[\"patient\",\"height\", \"width\", \"depth\"])\n    \n    return [flair, t1w, t1wCE, t2w]\n        \ndatalist_train = combine_scans(TRAIN_DIR, save_train_dir)\nprint(\"Done with Train Dataset\")\nprint(\"Start Test Dataset\")\ndatalist_test = combine_scans(TEST_DIR, save_test_dir)","metadata":{"execution":{"iopub.status.busy":"2021-08-20T01:52:23.163299Z","iopub.execute_input":"2021-08-20T01:52:23.163624Z","iopub.status.idle":"2021-08-20T01:52:47.706153Z","shell.execute_reply.started":"2021-08-20T01:52:23.163596Z","shell.execute_reply":"2021-08-20T01:52:47.705171Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"names = [\"FLAIR\", \"T1w\", \"T1wCE\", \"T2w\"]\nfig, axes = plt.subplots(nrows=4, ncols=3, figsize=(15,10))\nfor i, x in enumerate(datalist_train):\n    sns.histplot(data=x['height'],binwidth=10, ax=axes[i, 0], multiple='stack')\n    axes[i,0].set_title(f'{names[i]} Height')\n    sns.histplot(data=x['width'],binwidth=10, ax=axes[i, 1], multiple='stack')\n    axes[i,1].set_title(f'{names[i]} Width')\n    sns.histplot(data=x['depth'],binwidth=10, ax=axes[i, 2], multiple='stack')\n    axes[i,2].set_title(f'{names[i]} Depth')\nplt.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-08-20T01:52:47.709333Z","iopub.execute_input":"2021-08-20T01:52:47.709631Z","iopub.status.idle":"2021-08-20T01:52:49.265261Z","shell.execute_reply.started":"2021-08-20T01:52:47.709602Z","shell.execute_reply":"2021-08-20T01:52:49.264282Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# shutil.rmtree('./')\n# os.listdir('./')","metadata":{"execution":{"iopub.status.busy":"2021-08-20T01:52:49.267066Z","iopub.execute_input":"2021-08-20T01:52:49.267373Z","iopub.status.idle":"2021-08-20T01:52:49.271020Z","shell.execute_reply.started":"2021-08-20T01:52:49.267345Z","shell.execute_reply":"2021-08-20T01:52:49.270107Z"},"trusted":true},"execution_count":null,"outputs":[]}]}