{"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":"import numpy as np\nimport pandas as pd\nimport tensorflow as tf\nimport pydicom\nimport cv2\nimport matplotlib.pyplot as plt\nfrom scipy.interpolate import interp1d\nimport ipywidgets as widgets\n\nimport os\nfrom pathlib import Path\n\nplt.style.use('ggplot')","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2021-08-14T05:59:33.506929Z","iopub.execute_input":"2021-08-14T05:59:33.507428Z","iopub.status.idle":"2021-08-14T05:59:33.513572Z","shell.execute_reply.started":"2021-08-14T05:59:33.507396Z","shell.execute_reply":"2021-08-14T05:59:33.512865Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Resources Path\nBefore, we have put all dataset information into well structure by using panda dataframe.<br>\nJust for making life easier!","metadata":{}},{"cell_type":"code","source":"base_dir = Path(\"/kaggle/input\")\nds_dir = base_dir.joinpath(\"rsna-miccai-brain-tumor-radiogenomic-classification\")\nlabels = ds_dir.joinpath(\"train_labels.csv\")\ntrain_dir = ds_dir.joinpath(\"train\")\ntest_dir = ds_dir.joinpath(\"test\")\nclean_tain_dir = Path(\"/kaggle/working/brain\")\nIMCAT = [\"FLAIR\",\"T1w\",\"T1wCE\",\"T2w\"]\ntest_pre = os.listdir(test_dir)\ntest_pre = [int(pre)for pre in test_pre]\ntest_pre.sort()\ntrain_labels_df = pd.read_csv(labels)\ntrain_df = pd.read_pickle(str(base_dir.joinpath(\"bt-in-nutshell/aggregate.pkl\")))\ntrain_df.head()","metadata":{"execution":{"iopub.status.busy":"2021-08-14T05:59:33.719025Z","iopub.execute_input":"2021-08-14T05:59:33.719639Z","iopub.status.idle":"2021-08-14T05:59:43.384450Z","shell.execute_reply.started":"2021-08-14T05:59:33.719589Z","shell.execute_reply":"2021-08-14T05:59:43.383491Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Preprocessing","metadata":{}},{"cell_type":"markdown","source":"### Normalization method\n#### Z-score","metadata":{}},{"cell_type":"markdown","source":"$$\n\\frac{im_i - \\mu_i}{\\sigma_i}\n$$","metadata":{}},{"cell_type":"code","source":"def z_score(im):\n    \"\"\"\n    z-score nomalization\n    \"\"\"\n    mask_im = im>im.mean()\n    logical_mask = mask_im>0.\n    mean = im[logical_mask].mean()\n    std = im[logical_mask].std()\n    return (im-mean)/std","metadata":{"execution":{"iopub.status.busy":"2021-08-14T05:59:43.386323Z","iopub.execute_input":"2021-08-14T05:59:43.386737Z","iopub.status.idle":"2021-08-14T05:59:43.393351Z","shell.execute_reply.started":"2021-08-14T05:59:43.386694Z","shell.execute_reply":"2021-08-14T05:59:43.392096Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Nyul","metadata":{}},{"cell_type":"code","source":"def train_nyul(images,i_min=1,i_max=99,i_s_min = 1,i_s_max=100,l_perc=10,u_perc = 90,n_land = 10):\n    percs = np.concatenate(([i_min],np.arange(l_perc,u_perc+1,n_land),[i_max]))\n    standard_scale = np.zeros(len(percs))\n    \n    for _im in images:\n        mask_data = _im>np.mean(_im)\n        masked = _im[mask_data>0]\n        landmarks = np.percentile(masked, percs)\n        min_p = np.percentile(masked, i_min)\n        max_p = np.percentile(masked, i_max)\n        f = interp1d([min_p, max_p], [i_s_min, i_s_max])\n        landmarks = np.array(f(landmarks))\n        standard_scale += landmarks\n    standard_scale = standard_scale / len(images)\n    return standard_scale, percs\n\ndef go_on_hist(image,standard_scale,landmark_percs):\n    mask_data = image>image.mean()\n    masked = image[mask_data>0]\n    landmarks = np.percentile(masked, landmark_percs)\n    f = interp1d(landmarks, standard_scale, fill_value='extrapolate')\n    normed = f(image)\n    return normed","metadata":{"execution":{"iopub.status.busy":"2021-08-14T05:59:43.395153Z","iopub.execute_input":"2021-08-14T05:59:43.395461Z","iopub.status.idle":"2021-08-14T05:59:43.407482Z","shell.execute_reply.started":"2021-08-14T05:59:43.395432Z","shell.execute_reply":"2021-08-14T05:59:43.406775Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Tool function","metadata":{}},{"cell_type":"code","source":"def image_state(im):\n    \"\"\"\n    get image mean and std\n    \"\"\"\n    non_zero_pixels = im[np.nonzero(im)]\n    mean = np.mean(non_zero_pixels)\n    std = np.std(non_zero_pixels)\n    return mean,std\n\ndef level(mu,s,scale_factor=1.7):\n    \"\"\"\n    determine level for creating mask\n    \"\"\"\n    return mu+scale_factor*s\n\ndef calc_idx(image):\n    mean,std = image_state(image)\n    non_zero_pixels = np.count_nonzero(image>level(mean,std))\n    return non_zero_pixels\n\ndef top_valuable(images):\n    indices = [calc_idx(image)for image in images]\n    top_image = np.argsort(indices)[::-1][0]\n    return top_image\n\ndef top_valuable_line(image,axis):\n    mean,std = image_state(image)\n    non_cero_pixels = np.count_nonzero(image > level(mean,std),  axis=axis)\n    top_line = np.argsort(non_cero_pixels)[::-1][0]\n    return top_line\n\ndef calc_center(center_pos,r,c,t):\n    if t == \"Axial\":\n        center=[\n              center_pos[0] + center_pos[0] * c,\n              center_pos[1],\n              center_pos[2] - center_pos[1] * r\n        ]\n    if t == \"Saggital\":\n        center = [\n              center_pos[0],\n              center_pos[1] + center_pos[0] * c,\n              center_pos[2] - center_pos[1] * r\n        ]\n    if t==\"Coronal\":\n        center = [\n              center_pos[0] + center_pos[0] * c,\n              center_pos[1] + center_pos[0] * r,\n              center_pos[2]\n        ]\n    return center\n\ndef find_nearest(_centeres,_center,_ori):\n    axis_move = {'Sagittal': 0, 'Coronal': 1, 'Axial': 2}\n    scan = np.argsort(np.abs(_centeres - _center),axis=0)[0][axis_move[_ori]]\n    return scan\n\ndef read_dicom_images(dframe,case,study,orientation):\n    s_dframe = dframe.loc[dframe[\"patient_id\"]==case]\n    if study is not None:\n        s_dframe = s_dframe.loc[dframe[\"study\"]==study]\n    if orientation is not None:\n        s_dframe = s_dframe.loc[dframe[\"orientation\"]==orientation]\n    ims_paths = s_dframe[\"path\"].to_list()\n    sorted_paths = sorted(ims_paths,key=lambda x:int(x.split(\"Image-\")[1].split(\".\")[0]))\n    _ans = list()\n    _ans_path = list()\n    for p in sorted_paths:\n        p_im = pydicom.read_file(p).pixel_array\n        if np.all((p_im==0)):\n            continue\n        _ans.append(p_im)\n        _ans_path.append(p)\n    return _ans,_ans_path\n\ndef clean_zero_images(row):\n    im_obj = Image(row[\"path\"])\n    if np.all((im_obj.image==0)):\n        return False\n    if np.count_nonzero(im_obj.image)/(im_obj.image.shape[0]*im_obj.image.shape[1])<0.1:\n        return False\n    return True\n\ndef get_bounding_box(image):\n    mins = np.min(np.nonzero(image),axis=1)\n    maxs = np.max(np.nonzero(image),axis=1)\n    return mins[0],mins[1],maxs[0],maxs[1]\n\ndef extract_bounding_boxes(images):\n    bb = []\n    for im in images:\n        bb.append([*get_bounding_box(im)])\n    \n    return np.array(bb)\n\ndef extract_stuff(images,bb):\n    cropped = list()\n    \n    for idx,im in enumerate(images):\n        x_min,y_min,x_max,y_max = bb[idx]\n        cropped.append(im[x_min:x_max,y_min:y_max])\n    return cropped","metadata":{"execution":{"iopub.status.busy":"2021-08-14T06:10:39.459116Z","iopub.execute_input":"2021-08-14T06:10:39.459683Z","iopub.status.idle":"2021-08-14T06:10:39.484982Z","shell.execute_reply.started":"2021-08-14T06:10:39.459635Z","shell.execute_reply":"2021-08-14T06:10:39.484143Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def plot_observation(im):\n    \"\"\"\n    take an image and plot and make hist of it\n    \"\"\"\n    pixels = im.ravel()\n    non_z_pix = pixels[np.nonzero(pixels)]\n    mean,std = image_state(im)\n    threshold = np.count_nonzero(non_z_pix>level(mean,std))\n    fig, (axi, axh,axk) = plt.subplots(1, 3, figsize = (20,4), gridspec_kw={'width_ratios': [1, 4, 4]})\n    normal_non_z = z_score(non_z_pix)\n    normal_non_mean,normal_non_std = image_state(normal_non_z)\n    normal_threshold = np.count_nonzero(normal_non_z>level(normal_non_mean,normal_non_std))\n    \n    \n    fig.suptitle(f\"# over threshold normal-{threshold} & zero_score-{normal_threshold}\")\n    \n    axk.hist(normal_non_z, 200)\n    axk.set_title(\"Zero Score\")\n    ax_limits = axk.get_ylim()\n    axk.vlines(normal_non_mean, ymin=ax_limits[0], ymax=ax_limits[1], colors='b', label = \"mean\")\n    axk.vlines(normal_non_mean+normal_non_std, ymin=ax_limits[0], ymax=ax_limits[1], colors='b', linestyles='dotted',label=\"mean+std\")\n    axk.vlines(level(normal_non_mean,normal_non_std), ymin=ax_limits[0], ymax=ax_limits[1], colors='g', linestyles='dashed',label=\"threshold\")\n    axk.set_xlim(-6,6)\n    axk.legend(loc=\"upper left\")\n    axk.grid(False)\n    \n    axh.hist(non_z_pix, 200)\n    axh.set_title(\"Original\")\n    ax_limits = axh.get_ylim()\n    axh.vlines(mean, ymin=ax_limits[0], ymax=ax_limits[1], colors='b', label = \"mean\")\n    axh.vlines(mean+std, ymin=ax_limits[0], ymax=ax_limits[1], colors='b', linestyles='dotted',label=\"mean+std\")\n    axh.vlines(level(mean,std), ymin=ax_limits[0], ymax=ax_limits[1], colors='g', linestyles='dashed',label=\"threshold\")\n    axh.grid(False)\n    \n    axi.imshow(im, cmap = plt.cm.gray)\n    axi.grid(False)\n    axi.axis('off')\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2021-08-14T05:59:43.430630Z","iopub.execute_input":"2021-08-14T05:59:43.430940Z","iopub.status.idle":"2021-08-14T05:59:43.446665Z","shell.execute_reply.started":"2021-08-14T05:59:43.430908Z","shell.execute_reply":"2021-08-14T05:59:43.445507Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Some MIR sample\nwe have **three** study for each of observation:\n* **FLAIR**\n* **T1w**\n* **T1wCE**\n* **T2w**\n\n<p>As we can see, for each study, we have series of images, after that when we have ploted them, We saw some cases to mention:</p>\n\n* **Contrast**: <p> each study has specific contrast, and he can use techniques to enhance them</p>\n* **Some empty images**: <p> there were some empty images in folder, and for having a cleaner dataset we can explicitly remove and put them away.</p>\n* **Series**: <p> we have series of images for each study, but a subset of these images are valuable. As [mentioned](https://www.kaggle.com/josecarmona/btrc-eda-final), we can score images based on non zero voxels, and histogram.</p>\n* **Various shape**: <p>some images in **256x256** and **512x512**</p>\n* **Plates**:<p>sampled MRI have different posisioned and plates</p>\n    1. **Saggital**\n    2. **Coronal**\n    3. **Axial**\n\n\n","metadata":{}},{"cell_type":"markdown","source":"#### Some plot","metadata":{}},{"cell_type":"code","source":"cols_name = [\"FLAIR\",\"T1w\",\"T1wCE\",\"T2w\"]\nrows_name = [\"Saggital\",\"Coronal\",\"Axial\"]\n\nsample1 = \"00386\"\nsample_flair_images,_ = read_dicom_images(train_df,sample1,\"FLAIR\",None)\nsample_tw1_images,_=read_dicom_images(train_df,sample1,\"T1w\",None)\nsample_t1wce_image,_=read_dicom_images(train_df,sample1,\"T1wCE\",None)\nsample_t2w_images,_=read_dicom_images(train_df,sample1,\"T2w\",None)\n\nfig, axes = plt.subplots(nrows=1, ncols=4, figsize=(12, 8))\n\nfor ax, col in zip(axes, cols_name):\n    ax.set_title(col)\n\naxes[0].imshow(sample_flair_images[20],cmap = plt.cm.gray)\naxes[0].grid(False)\naxes[0].axis('off')\n\naxes[1].imshow(sample_tw1_images[20],cmap = plt.cm.gray)\naxes[1].grid(False)\naxes[1].axis('off')\n\naxes[2].imshow(sample_t2w_images[20],cmap = plt.cm.gray)\naxes[2].grid(False)\naxes[2].axis('off')\n\naxes[3].imshow(sample_t1wce_image[20],cmap = plt.cm.gray)\naxes[3].grid(False)\naxes[3].axis('off')\n\n    \nfig.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-08-14T05:59:43.447985Z","iopub.execute_input":"2021-08-14T05:59:43.448566Z","iopub.status.idle":"2021-08-14T05:59:46.221425Z","shell.execute_reply.started":"2021-08-14T05:59:43.448527Z","shell.execute_reply":"2021-08-14T05:59:46.220227Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Compare\n#### Zero score","metadata":{}},{"cell_type":"code","source":"patient_id = \"00386\"\nflair_images,_ = read_dicom_images(train_df,patient_id,\"FLAIR\",\"Axial\")\nfor im in flair_images:\n    plot_observation(im)","metadata":{"execution":{"iopub.status.busy":"2021-08-14T05:59:46.222896Z","iopub.execute_input":"2021-08-14T05:59:46.223339Z","iopub.status.idle":"2021-08-14T06:00:43.387857Z","shell.execute_reply.started":"2021-08-14T05:59:46.223293Z","shell.execute_reply":"2021-08-14T06:00:43.386896Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Nyul \n> we use a more cleaned data to create a trained landmark on nyul algorithm","metadata":{}},{"cell_type":"code","source":"c_train_df = pd.read_pickle(str(base_dir.joinpath(\"bt-in-nutshell/cleaned_data.pkl\")))\nc_train_df.head()","metadata":{"execution":{"iopub.status.busy":"2021-08-14T06:00:43.390213Z","iopub.execute_input":"2021-08-14T06:00:43.390541Z","iopub.status.idle":"2021-08-14T06:00:43.434017Z","shell.execute_reply.started":"2021-08-14T06:00:43.390513Z","shell.execute_reply":"2021-08-14T06:00:43.432828Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"flaired_images = c_train_df.loc[c_train_df[\"study\"]==\"FLAIR\"]\nflaired_axial_images = flaired_images.loc[flaired_images[\"orientation\"]==\"Axial\"][\"path\"]\nflaired_axial_images = [pydicom.read_file(p).pixel_array for p in flaired_axial_images]","metadata":{"execution":{"iopub.status.busy":"2021-08-14T06:00:43.435799Z","iopub.execute_input":"2021-08-14T06:00:43.436104Z","iopub.status.idle":"2021-08-14T06:00:44.612070Z","shell.execute_reply.started":"2021-08-14T06:00:43.436076Z","shell.execute_reply":"2021-08-14T06:00:44.611015Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### nyul hyperparameter","metadata":{}},{"cell_type":"code","source":"i_min = widgets.IntSlider(\n    value=1,\n    min=1,\n    max=99,\n    step=5,\n    description='Minimum Percentil:',\n    disabled=False,\n    continuous_update=False,\n    orientation='horizontal',\n    readout=True,\n    readout_format='d'\n)\n\ni_max = widgets.IntSlider(\n    value=99,\n    min=1,\n    max=99,\n    step=4,\n    description='Maximum Percentil:',\n    disabled=False,\n    continuous_update=False,\n    orientation='horizontal',\n    readout=True,\n    readout_format='d'\n)\n\ni_s_min = widgets.IntSlider(\n    value=1,\n    min=1,\n    max=100,\n    step=5,\n    description='Minimum standard Percentil:',\n    disabled=False,\n    continuous_update=False,\n    orientation='horizontal',\n    readout=True,\n    readout_format='d'\n)\n\ni_s_max = widgets.IntSlider(\n    value=100,\n    min=1,\n    max=100,\n    step=5,\n    description='Maximum Standard Percentil:',\n    disabled=False,\n    continuous_update=False,\n    orientation='horizontal',\n    readout=True,\n    readout_format='d'\n)\n\nl_perc = i_s_max = widgets.IntSlider(\n    value=10,\n    min=1,\n    max=100,\n    step=5,\n    description='Low Middle Percentil:',\n    disabled=False,\n    continuous_update=False,\n    orientation='horizontal',\n    readout=True,\n    readout_format='d'\n)\n\nu_perc = i_s_max = widgets.IntSlider(\n    value=90,\n    min=1,\n    max=100,\n    step=5,\n    description='Upper Middle Percentil:',\n    disabled=False,\n    continuous_update=False,\n    orientation='horizontal',\n    readout=True,\n    readout_format='d'\n)\n\nn_land = widgets.IntText(\n    value=10,\n    description='Number of landmarks:',\n    disabled=False\n)\n\ncontainer = widgets.VBox([i_min,i_max,i_s_min,i_s_max,l_perc,u_perc,n_land])\ncontainer","metadata":{"execution":{"iopub.status.busy":"2021-08-14T06:00:44.613225Z","iopub.execute_input":"2021-08-14T06:00:44.613565Z","iopub.status.idle":"2021-08-14T06:00:44.702099Z","shell.execute_reply.started":"2021-08-14T06:00:44.613535Z","shell.execute_reply":"2021-08-14T06:00:44.701381Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"standard_scale,landmarks = train_nyul(flaired_axial_images,i_min.value,i_max.value,i_s_min.value,i_s_max.value,l_perc.value,u_perc.value,n_land.value)\nfig,ax = plt.subplots(1,1,figsize = (15,4))\nax.plot(standard_scale,landmarks)\nax_limits = ax.get_ylim()\nax.grid(False)\nax.set_xlabel(\"landmark\")\nfor land in landmarks:\n    ax.vlines(land, ymin=ax_limits[0], ymax=ax_limits[1], colors='b', linestyles='dotted')\nax.set_ylabel(\"standard\")\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2021-08-14T06:00:44.703127Z","iopub.execute_input":"2021-08-14T06:00:44.703422Z","iopub.status.idle":"2021-08-14T06:00:45.627700Z","shell.execute_reply.started":"2021-08-14T06:00:44.703395Z","shell.execute_reply":"2021-08-14T06:00:45.626657Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"normal_flaired_axial_images = [go_on_hist(im,standard_scale,landmarks) for im in flaired_axial_images]\nbboxes = extract_bounding_boxes(flaired_axial_images)\ncropped_flaired_axial_images = extract_stuff(flaired_axial_images,bboxes)\ncropped_normal_flaired_axial_images = extract_stuff(normal_flaired_axial_images,bboxes)","metadata":{"execution":{"iopub.status.busy":"2021-08-14T06:14:20.000325Z","iopub.execute_input":"2021-08-14T06:14:20.000881Z","iopub.status.idle":"2021-08-14T06:14:22.954976Z","shell.execute_reply.started":"2021-08-14T06:14:20.000834Z","shell.execute_reply":"2021-08-14T06:14:22.953976Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, axes = plt.subplots(10, 2, figsize = (20,100), gridspec_kw={'width_ratios': [50, 50]})\naxes[0][0].set_title(\"Original\")\naxes[0][1].set_title(\"Nyul Normalized\")\n\nfor i,(im,norm_im) in enumerate(zip(cropped_flaired_axial_images[:10],cropped_normal_flaired_axial_images[:10])):\n    axes[i][0].imshow(im,cmap=\"gray\")\n    axes[i][0].grid(False)\n    axes[i][0].axis(\"off\")\n    axes[i][1].imshow(norm_im,cmap=\"gray\")\n    axes[i][1].grid(False)\n    axes[i][1].axis(\"off\")\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2021-08-14T06:15:26.268008Z","iopub.execute_input":"2021-08-14T06:15:26.268549Z","iopub.status.idle":"2021-08-14T06:15:29.115331Z","shell.execute_reply.started":"2021-08-14T06:15:26.268517Z","shell.execute_reply":"2021-08-14T06:15:29.114019Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}