{"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":"# Things to be solved..\n- pt no. 109(no FLAIR) & 123(no T1w) & 709(no FLAIR) should be excluded!\n- choice :  \n1. Normalization(for 2D slice? or 3D bulk image?) : 2 options\n2. Again Normlization(for Patch? or for overall 2D(or 3D) image?) : * 2 options\n3. Normalization method (mean-norm, zero-mean, unit-varia~) : * n options\n#### 그런데! But all related to normalization.. (학습 속도를 빠르게 하느냐 차이니까...뭐가 되었던 normalization을 해주기만 하면 되는거 아닌지?)","metadata":{}},{"cell_type":"markdown","source":"## 0. Load required libraries & csv files","metadata":{}},{"cell_type":"code","source":"# Install required packages\n! pip install natsort","metadata":{"execution":{"iopub.status.busy":"2021-08-11T01:32:28.477827Z","iopub.execute_input":"2021-08-11T01:32:28.478531Z","iopub.status.idle":"2021-08-11T01:32:36.237553Z","shell.execute_reply.started":"2021-08-11T01:32:28.478486Z","shell.execute_reply":"2021-08-11T01:32:36.236572Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import pandas as pd \nimport matplotlib.pyplot as plt\n\nimport os \nimport cv2\nimport numpy as np\nimport glob\nimport pydicom\nfrom pydicom.pixel_data_handlers.util import apply_voi_lut\n\n#from torch.utils.data import Dataset, DataLoader\nimport SimpleITK as sitk\n\n# to sort file names by its order\nfrom natsort import natsorted\n%matplotlib inline\n\n# Global variables\nmodalities = ['FLAIR', 'T1w', 'T1wCE', 'T2w']","metadata":{"execution":{"iopub.status.busy":"2021-08-11T01:32:37.389064Z","iopub.execute_input":"2021-08-11T01:32:37.38943Z","iopub.status.idle":"2021-08-11T01:32:37.402596Z","shell.execute_reply.started":"2021-08-11T01:32:37.389398Z","shell.execute_reply":"2021-08-11T01:32:37.401474Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Load y : labels & preds\ntrain_labels = pd.read_csv('../input/rsna-miccai-brain-tumor-radiogenomic-classification/train_labels.csv')\npreds_labels = pd.read_csv('../input/rsna-miccai-brain-tumor-radiogenomic-classification/sample_submission.csv')\n#train_labels = pd.read_csv('../png_data/train_labels.csv')\n#preds_labels = pd.read_csv('../png_data/sample_submission.csv')\n\nprint(f'Number of patients = {len(train_labels)}'.format())\ntrain_labels.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 1. Show sample","metadata":{}},{"cell_type":"code","source":"# Reference: https://www.kaggle.com/ayuraj/brain-tumor-eda-and-interactive-viz-with-w-b\ndef ReadMRI(path, voi_lut = True, fix_monochrome = True):\n    # Original from: https://www.kaggle.com/raddar/convert-dicom-to-np-array-the-correct-way\n    dicom = pydicom.read_file(path)\n    \n    # VOI LUT (if available by DICOM device) is used to transform raw DICOM data to \n    # \"human-friendly\" view\n    if voi_lut:\n        data = apply_voi_lut(dicom.pixel_array, dicom)\n    else:\n        data = dicom.pixel_array\n               \n    # depending on this value, X-ray may look inverted - fix that:\n    if fix_monochrome and dicom.PhotometricInterpretation == \"MONOCHROME1\":\n        data = np.amax(data) - data\n        \n    # min-max normalization\n    data = data - np.min(data)\n    data = data / np.max(data)\n    data = (data * 255).astype(np.uint8)\n    \n    return data","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def load_imgs(idx, train=True): \n    images_idx = {}\n    for modal in modalities:\n        images_modal = []\n        if train:\n            file_path_list = glob.glob('../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/'+idx+'/'+modal+'/*')\n        else:\n            file_path_list = glob.glob('../input/rsna-miccai-brain-tumor-radiogenomic-classification/test/'+idx+'/'+modal+'/*')\n            \n        # Should be sorted again. \n        file_path_list = natsorted(file_path_list)\n        \n        for file_path in file_path_list:\n            image = ReadMRI(file_path)\n            # In case, direc is empty!\n            if len(image) == 0:\n                print('yes')\n                image = np.zeros((1,256,256)) # pt no. 109(no FLAIR) & 123(no T1w) & 709(no FLAIR) >> excluded! \n            images_modal.append(image)\n        images_idx[modal] = np.array(images_modal)\n        \n    return images_idx","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n# Just to check how long would \"load_imgs\" be. \n# %%time 위에 주석 달면 time 체크 불가....!?\nfor i in range(1):\n    idx = str(train_labels.BraTS21ID[i]).zfill(5)\n    print(idx)\n    imgs = load_imgs(idx)\n    img_ = imgs[modalities[0]]#.mean(axis=0)\n    print(img_.shape)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Again, just to check sample\nfig, ax = plt.subplots(2,2, figsize=(10,10))\nfor i in range(2):\n    for j in range(2):\n        m = ax[i,j].imshow(imgs[modalities[2*i+j]].mean(axis=0))\n        ax[i,j].set_title(modalities[2*i+j])\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 2. Preprocessing each image (per pt, modality)","metadata":{}},{"cell_type":"code","source":"# 1. Create a data loader with N4biasFieldCorrectionImageFilter\nclass PreprocessedImage2Dver():\n    def __init__(self, list_BraTS21ID, list_labels=None,\n                dim=(512, 512), n_modals=len(modalities), n_classes=2,\n                 num_slices_from_center:int=3,\n                 is_train=True, transform = None): # For single pt. \n        self.dim = dim\n        ### ????? 1. How to do batch norm for this task..?\n        #self.batch_size = batch_size\n        self.list_labels = list_labels\n        self.is_train = (list_labels is not None)\n        self.list_BraTS21ID = list_BraTS21ID\n        self.n_modals = n_modals # number of modalities\n        self.num_slices_from_center = num_slices_from_center\n        \n    \n    def __getitem__(self, index):\n        BraTS21ID_temp = self.list_BraTS21ID[index] #self.list_BraTS21ID[index*self.batch_size:(index+1)*self.batch_size] # index로 적어준 batch 만 진행!\n        \n        X = self.__data_generation(BraTS21ID_temp)\n        \n        if self.is_train:\n            y = self.list_labels[index] #self.list_labels[index*self.batch_size:(index+1)*self.batch_size]\n            return np.array(X), np.array(y)\n        else:\n            return np.array(X)\n    \n    def __data_generation(self, BraTS21ID_temp):\n        new_imgs = np.zeros((self.num_slices_from_center*2,  *self.dim, self.n_modals))\n        #print(new_imgs.shape)\n        #new_imgs = None\n        \n        idx = str(BraTS21ID_temp).zfill(5)\n        imgs = load_imgs(idx, train=self.is_train) # imgs = {'FLAIR': ~, 'T1w': ~, 'T1wCE': ~} lib다. \n        \n        index_modal = 0\n        for modal in modalities:\n            corrct_norm_imgs_modal = []\n\n            img_modal = imgs[modal] #.shape :  ex. (288, 256, 256) \n            \n            if img_modal.shape[0] < self.num_slices_from_center *2 :\n                print('The number of slice is smaller than total number of slices!!')\n                break\n            \n            central_slices = int(img_modal.shape[0]/2)\n            \n            list_slices = list(range(central_slices-self.num_slices_from_center, central_slices+self.num_slices_from_center))\n            for slice_i in list_slices:\n                img_slice_i = cv2.resize(img_modal[slice_i,:,:], dsize= self.dim, interpolation = cv2.INTER_LINEAR) #???\n                img_slice_i_ndarray = np.array(img_slice_i, dtype = 'float32')\n\n                # 1. Removing radiofrequency inhomogeneity using N4BiasFieldCorrection\n                # ref : https://www.kaggle.com/josepc/rsna-effnet\n                inputImage = sitk.GetImageFromArray(img_slice_i_ndarray)\n                maskImage = sitk.GetImageFromArray((img_slice_i_ndarray>0.1)*1) #sitk.OtsuThreshold(inputImage, 0,1,200) \n                inputImage = sitk.Cast(inputImage, sitk.sitkFloat32) #?? 왜 32로?\n                maskImage = sitk.Cast(maskImage, sitk.sitkUInt8) #?? 왜 8로?\n                corrector = sitk.N4BiasFieldCorrectionImageFilter()\n                numberFittingLevels = 4 # ??\n                maxIter = 100 # ??\n                if maxIter is not None:\n                    corrector.SetMaximumNumberOfIterations([maxIter]*numberFittingLevels) # ??\n                corrected_image = corrector.Execute(inputImage, maskImage)\n                \n                corrected_image = sitk.GetArrayFromImage(corrected_image)\n                \n                # 2. Normalization (근데 앞에 load_imgs에서 normalization term이 있긴 했음.)\n                # ????? 2. 학습 속도만 해결하면 되는거니, 어떻게든 normalization만 하면 되는걸지?\n                max_2D = np.amax(corrected_image) \n                min_2D = np.amin(corrected_image)\n                mean_2D = np.mean(corrected_image)\n                \n                normalized_corrected_image = (corrected_image-min_2D)/(max_2D-min_2D)\n                \n                \n                corrct_norm_imgs_modal.append(normalized_corrected_image)\n                \n            corrct_norm_imgs_modal = np.array(corrct_norm_imgs_modal)\n            \n            new_imgs[:, :, :, index_modal] = corrct_norm_imgs_modal\n            \n            index_modal += 1\n        \n        return new_imgs # shape = (num_slices_from_center*2, *dim, n_modals) !!! not 3D. 4D.","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 3. Check sample preprocessed image\n### - num_slices_from_center:int=3\n### - for \"pt no.00002\" \n### - for 2D slice image at center","metadata":{}},{"cell_type":"code","source":"# Do train_test_split. \n# - Train : Test = 8:2 / Stratified random sampling / random_state = 42\nfrom sklearn.model_selection import train_test_split\nX_train, X_val, y_train, y_val = train_test_split(train_labels.BraTS21ID, train_labels.MGMT_value,\n                                                 test_size = .2, random_state = 42, stratify = train_labels.MGMT_value)","metadata":{"execution":{"iopub.status.busy":"2021-08-11T14:12:57.579916Z","iopub.execute_input":"2021-08-11T14:12:57.580539Z","iopub.status.idle":"2021-08-11T14:12:58.576733Z","shell.execute_reply.started":"2021-08-11T14:12:57.580442Z","shell.execute_reply":"2021-08-11T14:12:58.575351Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ??????? 3. 원래의 W*H 사이즈를 보존해야?\ndim = (256, 256) \ntrain_set = PreprocessedImage2Dver(X_train, y_train, dim=dim)\nval_set = PreprocessedImage2Dver(X_val, y_val, dim=dim)\ntest_set = PreprocessedImage2Dver(preds_labels.BraTS21ID,  dim = dim)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sample_preprocessed_imgs = train_set[1][0] # pt no.00002\n\ncenter_slice_idx = int(sample_preprocessed_imgs.shape[0]/2)\n#print(center_slice_idx) # In this cas, returns 3. \nimg_at_center = sample_preprocessed_imgs[center_slice_idx] # slice at center","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Check sampled space of sample_img (FLAIR ver). \nprint(img_at_center[100:150, 100:130, 0])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Again, just to check sample preprocessed image.\nfig, ax = plt.subplots(2,2, figsize=(10,10))\nfor i in range(2):\n    for j in range(2):\n        m = ax[i,j].imshow(img_at_center[:, :, 2*i+j], vmin=0, vmax=1)\n        ax[i,j].set_title(modalities[2*i+j])\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Before preprocessing","metadata":{}},{"cell_type":"code","source":"sample_pt_ID = str(train_labels.BraTS21ID[1]).zfill(5)\n\nbefore_prepro_img = load_imgs(idx=sample_pt_ID)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('Before preprocessing - pt no.00002')\n# Again, just to check sample\nfig, ax = plt.subplots(2,2, figsize=(10,10))\nfor i in range(2):\n    for j in range(2):\n        center_slice_idx = int(before_prepro_img[modalities[2*i+j]].shape[0]/2)\n        #print(center_slice_idx)\n        m = ax[i,j].imshow(before_prepro_img[modalities[2*i+j]][center_slice_idx])\n        ax[i,j].set_title(modalities[2*i+j])\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}