{"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":"#https://www.kaggle.com/boojum/connecting-voxel-spaces/","metadata":{"execution":{"iopub.status.busy":"2021-08-26T19:16:52.908262Z","iopub.execute_input":"2021-08-26T19:16:52.909013Z","iopub.status.idle":"2021-08-26T19:16:52.914227Z","shell.execute_reply.started":"2021-08-26T19:16:52.908911Z","shell.execute_reply":"2021-08-26T19:16:52.912801Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"I thought that all pipelines presented on public notebooks are giving random output, so I decided to think a little bit differently. The idea is to extract features using MRI images and tumor mask. Here I'm giving an example for one image, maybe I'm doing something wrong, but I'll be glad to get feedback from the community. Also I'm worried that it will be out of time on the whole test set\n\n### Registration \nhttps://www.kaggle.com/boojum/connecting-voxel-spaces/\n![](https://sun9-45.userapi.com/impg/Flbnug2OUli1ecXsoIKeUasIGXGj_5hqjX4cRg/z2nfz8-b3a0.jpg?size=2560x1153&quality=96&sign=0536543610f1655d967af88dbc775e98&type=album)\n\n### Segmentation \nUnet\nhttps://pytorch.org/hub/mateuszbuda_brain-segmentation-pytorch_unet/\n\n### Features\n\nhttps://pyradiomics.readthedocs.io/en/latest/features.html#module-radiomics.shape2D","metadata":{}},{"cell_type":"code","source":"!pip install pyradiomics","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2021-08-26T19:16:52.916183Z","iopub.execute_input":"2021-08-26T19:16:52.916913Z","iopub.status.idle":"2021-08-26T19:17:06.819519Z","shell.execute_reply.started":"2021-08-26T19:16:52.916806Z","shell.execute_reply":"2021-08-26T19:17:06.818234Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os\nimport sys \nfrom tqdm import tqdm\nimport numpy as np\nimport pandas as pd\nfrom PIL import Image\nimport pydicom\nimport torch\nimport nibabel as nib\nimport matplotlib.pyplot as plt\nimport SimpleITK as sitk\nimport radiomics\n\ntrain_path = '../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/'","metadata":{"execution":{"iopub.status.busy":"2021-08-26T19:17:06.821181Z","iopub.execute_input":"2021-08-26T19:17:06.821480Z","iopub.status.idle":"2021-08-26T19:17:09.158076Z","shell.execute_reply.started":"2021-08-26T19:17:06.821450Z","shell.execute_reply":"2021-08-26T19:17:09.156917Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_dirs = sorted(os.listdir(train_path))","metadata":{"execution":{"iopub.status.busy":"2021-08-26T19:17:09.159982Z","iopub.execute_input":"2021-08-26T19:17:09.160359Z","iopub.status.idle":"2021-08-26T19:17:09.215641Z","shell.execute_reply.started":"2021-08-26T19:17:09.160326Z","shell.execute_reply":"2021-08-26T19:17:09.214561Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"reader = sitk.ImageSeriesReader()\nreader.LoadPrivateTagsOn()","metadata":{"execution":{"iopub.status.busy":"2021-08-26T19:17:09.228062Z","iopub.execute_input":"2021-08-26T19:17:09.228396Z","iopub.status.idle":"2021-08-26T19:17:09.247219Z","shell.execute_reply.started":"2021-08-26T19:17:09.228368Z","shell.execute_reply":"2021-08-26T19:17:09.246053Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def resample(image, ref_image):\n\n    resampler = sitk.ResampleImageFilter()\n    resampler.SetReferenceImage(ref_image)\n    resampler.SetInterpolator(sitk.sitkLinear)\n    \n    resampler.SetTransform(sitk.AffineTransform(image.GetDimension()))\n\n    resampler.SetOutputSpacing(ref_image.GetSpacing())\n\n    resampler.SetSize(ref_image.GetSize())\n\n    resampler.SetOutputDirection(ref_image.GetDirection())\n\n    resampler.SetOutputOrigin(ref_image.GetOrigin())\n\n    resampler.SetDefaultPixelValue(image.GetPixelIDValue())\n\n    resamped_image = resampler.Execute(image)\n    \n    return resamped_image","metadata":{"execution":{"iopub.status.busy":"2021-08-26T19:17:09.250112Z","iopub.execute_input":"2021-08-26T19:17:09.250635Z","iopub.status.idle":"2021-08-26T19:17:09.257876Z","shell.execute_reply.started":"2021-08-26T19:17:09.250595Z","shell.execute_reply":"2021-08-26T19:17:09.256987Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def normalize(data):\n    return (data - np.min(data)) / (np.max(data) - np.min(data))","metadata":{"execution":{"iopub.status.busy":"2021-08-26T19:17:09.259537Z","iopub.execute_input":"2021-08-26T19:17:09.259919Z","iopub.status.idle":"2021-08-26T19:17:09.271332Z","shell.execute_reply.started":"2021-08-26T19:17:09.259886Z","shell.execute_reply":"2021-08-26T19:17:09.270257Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ndef get_img(index):\n    filenamesDICOM = reader.GetGDCMSeriesFileNames(f'{train_path}/{train_dirs[index]}/T1w')\n    reader.SetFileNames(filenamesDICOM)\n    t1_sitk = reader.Execute()\n\n    filenamesDICOM = reader.GetGDCMSeriesFileNames(f'{train_path}/{train_dirs[index]}/FLAIR')\n    reader.SetFileNames(filenamesDICOM)\n    flair_sitk = reader.Execute()\n\n    filenamesDICOM = reader.GetGDCMSeriesFileNames(f'{train_path}/{train_dirs[index]}/T1wCE')\n    reader.SetFileNames(filenamesDICOM)\n    t1wce_sitk = reader.Execute()\n\n    flair_resampled = resample(flair_sitk, t1_sitk)\n    t1wce_resampled = resample(t1wce_sitk, t1_sitk)\n\n    t1_sitk_array = normalize(sitk.GetArrayFromImage(t1_sitk))\n    flair_resampled_array = normalize(sitk.GetArrayFromImage(flair_resampled))\n    t1wce_resampled_array = normalize(sitk.GetArrayFromImage(t1wce_resampled))\n\n    stacked = np.stack([t1_sitk_array, flair_resampled_array, t1wce_resampled_array])\n\n    to_rgb = stacked[:,t1_sitk_array.shape[0]//2,:,:].transpose(1,2,0)\n    im = Image.fromarray((to_rgb * 255).astype(np.uint8))\n    return im","metadata":{"execution":{"iopub.status.busy":"2021-08-26T19:17:09.272841Z","iopub.execute_input":"2021-08-26T19:17:09.273200Z","iopub.status.idle":"2021-08-26T19:17:09.290410Z","shell.execute_reply.started":"2021-08-26T19:17:09.273146Z","shell.execute_reply":"2021-08-26T19:17:09.289089Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model = torch.hub.load('mateuszbuda/brain-segmentation-pytorch', 'unet',\n    in_channels=3, out_channels=1, init_features=32, pretrained=True)","metadata":{"execution":{"iopub.status.busy":"2021-08-26T19:17:09.292071Z","iopub.execute_input":"2021-08-26T19:17:09.292378Z","iopub.status.idle":"2021-08-26T19:17:15.479822Z","shell.execute_reply.started":"2021-08-26T19:17:09.292350Z","shell.execute_reply":"2021-08-26T19:17:15.478940Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"im = get_img(3)","metadata":{"execution":{"iopub.status.busy":"2021-08-26T19:17:15.481446Z","iopub.execute_input":"2021-08-26T19:17:15.482196Z","iopub.status.idle":"2021-08-26T19:17:24.240066Z","shell.execute_reply.started":"2021-08-26T19:17:15.482148Z","shell.execute_reply":"2021-08-26T19:17:24.239211Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_img = np.array([np.moveaxis(np.array(im.resize((256, 256))), -1, 0)])\ntest_res = model(torch.Tensor(test_img))","metadata":{"execution":{"iopub.status.busy":"2021-08-26T19:17:24.241630Z","iopub.execute_input":"2021-08-26T19:17:24.242385Z","iopub.status.idle":"2021-08-26T19:17:24.676560Z","shell.execute_reply.started":"2021-08-26T19:17:24.242336Z","shell.execute_reply":"2021-08-26T19:17:24.675617Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.imshow(test_img[0, 1, :, :])","metadata":{"execution":{"iopub.status.busy":"2021-08-26T19:17:24.678307Z","iopub.execute_input":"2021-08-26T19:17:24.679143Z","iopub.status.idle":"2021-08-26T19:17:24.908265Z","shell.execute_reply.started":"2021-08-26T19:17:24.679089Z","shell.execute_reply":"2021-08-26T19:17:24.907264Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.imshow(test_res[0][0].detach().cpu().numpy() > 0.5)","metadata":{"execution":{"iopub.status.busy":"2021-08-26T19:17:24.910056Z","iopub.execute_input":"2021-08-26T19:17:24.910464Z","iopub.status.idle":"2021-08-26T19:17:25.105241Z","shell.execute_reply.started":"2021-08-26T19:17:24.910417Z","shell.execute_reply":"2021-08-26T19:17:25.103839Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"shape = radiomics.shape2D.RadiomicsShape2D(\n    sitk.GetImageFromArray(test_img), \n    sitk.GetImageFromArray(np.array([\n        test_res[0][0].detach().cpu().numpy() > 0.5\n    ]).astype(np.uint8)),\n    force2D=True\n)","metadata":{"execution":{"iopub.status.busy":"2021-08-26T19:17:25.106835Z","iopub.execute_input":"2021-08-26T19:17:25.107253Z","iopub.status.idle":"2021-08-26T19:17:25.129999Z","shell.execute_reply.started":"2021-08-26T19:17:25.107216Z","shell.execute_reply":"2021-08-26T19:17:25.128907Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"shape.getMeshSurfaceFeatureValue(), \\\nshape.getPixelSurfaceFeatureValue(),\\\nshape.getPerimeterFeatureValue(), \\\nshape.getPerimeterSurfaceRatioFeatureValue(), \\\nshape.getSphericityFeatureValue(), \\\nshape.getSphericalDisproportionFeatureValue(), \\\nshape.getMaximumDiameterFeatureValue(), \\\nshape.getMajorAxisLengthFeatureValue(), \\\nshape.getMinorAxisLengthFeatureValue(), \\\nshape.getElongationFeatureValue(), \\","metadata":{"execution":{"iopub.status.busy":"2021-08-26T19:17:25.131235Z","iopub.execute_input":"2021-08-26T19:17:25.131664Z","iopub.status.idle":"2021-08-26T19:17:25.140442Z","shell.execute_reply.started":"2021-08-26T19:17:25.131615Z","shell.execute_reply":"2021-08-26T19:17:25.139361Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## P.S. Warning\nThe segmentation doesn't not working properly on all images due to the different tumor modalities. The model was trained on low-grade tumors, so be aware. Let's see the examples","metadata":{}},{"cell_type":"markdown","source":"#### For this type of broken segmentation like in '00009' we could perform a little trick to fix. \nWe need to multiply mask to the empty space of the original image","metadata":{}},{"cell_type":"code","source":"im = get_img(6)\ntest_img = np.array([np.moveaxis(np.array(im.resize((256, 256))), -1, 0)])\ntest_res = model(torch.Tensor(test_img))\n\nf, axarr = plt.subplots(1,5, figsize=(20, 20))\naxarr[0].imshow(test_img[0, 0])\naxarr[1].imshow(test_img[0, 1])\naxarr[2].imshow(test_img[0, 2])\naxarr[3].imshow(test_res[0][0].detach().cpu().numpy() > 0.5)\naxarr[4].imshow((test_res[0][0].detach().cpu().numpy() > 0.5) * (test_img[0, 1] != 0))","metadata":{"execution":{"iopub.status.busy":"2021-08-26T19:36:52.201266Z","iopub.execute_input":"2021-08-26T19:36:52.201673Z","iopub.status.idle":"2021-08-26T19:36:57.843503Z","shell.execute_reply.started":"2021-08-26T19:36:52.201641Z","shell.execute_reply":"2021-08-26T19:36:57.842247Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### For the \"00003\" image the broken segmentation could be due to the bad registration and this trick won't work","metadata":{}},{"cell_type":"code","source":"im = get_img(2)\ntest_img = np.array([np.moveaxis(np.array(im.resize((256, 256))), -1, 0)])\ntest_res = model(torch.Tensor(test_img))\n\nf, axarr = plt.subplots(1,5, figsize=(20, 20))\naxarr[0].imshow(test_img[0, 0])\naxarr[1].imshow(test_img[0, 1])\naxarr[2].imshow(test_img[0, 2])\naxarr[3].imshow(test_res[0][0].detach().cpu().numpy() > 0.5)\naxarr[4].imshow((test_res[0][0].detach().cpu().numpy() > 0.5) * (test_img[0, 1] != 0))","metadata":{"execution":{"iopub.status.busy":"2021-08-26T19:36:57.846273Z","iopub.execute_input":"2021-08-26T19:36:57.846767Z","iopub.status.idle":"2021-08-26T19:37:01.231544Z","shell.execute_reply.started":"2021-08-26T19:36:57.846718Z","shell.execute_reply":"2021-08-26T19:37:01.230464Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Let's see more examples of broken segmentation","metadata":{}},{"cell_type":"code","source":"im = get_img(9)\ntest_img = np.array([np.moveaxis(np.array(im.resize((256, 256))), -1, 0)])\ntest_res = model(torch.Tensor(test_img))\n\nf, axarr = plt.subplots(1,5, figsize=(20, 20))\naxarr[0].imshow(test_img[0, 0])\naxarr[1].imshow(test_img[0, 1])\naxarr[2].imshow(test_img[0, 2])\naxarr[3].imshow(test_res[0][0].detach().cpu().numpy() > 0.5)\naxarr[4].imshow((test_res[0][0].detach().cpu().numpy() > 0.5) * (test_img[0, 1] != 0))","metadata":{"execution":{"iopub.status.busy":"2021-08-26T19:34:57.110825Z","iopub.execute_input":"2021-08-26T19:34:57.111255Z","iopub.status.idle":"2021-08-26T19:35:17.608005Z","shell.execute_reply.started":"2021-08-26T19:34:57.111218Z","shell.execute_reply":"2021-08-26T19:35:17.606358Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"im = get_img(20)\ntest_img = np.array([np.moveaxis(np.array(im.resize((256, 256))), -1, 0)])\ntest_res = model(torch.Tensor(test_img))\n\nf, axarr = plt.subplots(1,5, figsize=(20, 20))\naxarr[0].imshow(test_img[0, 0])\naxarr[1].imshow(test_img[0, 1])\naxarr[2].imshow(test_img[0, 2])\naxarr[3].imshow(test_res[0][0].detach().cpu().numpy() > 0.5)\naxarr[4].imshow((test_res[0][0].detach().cpu().numpy() > 0.5) * (test_img[0, 1] != 0))","metadata":{"execution":{"iopub.status.busy":"2021-08-26T19:38:13.785860Z","iopub.execute_input":"2021-08-26T19:38:13.786297Z","iopub.status.idle":"2021-08-26T19:38:19.344937Z","shell.execute_reply.started":"2021-08-26T19:38:13.786260Z","shell.execute_reply":"2021-08-26T19:38:19.343920Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"im = get_img(24)\ntest_img = np.array([np.moveaxis(np.array(im.resize((256, 256))), -1, 0)])\ntest_res = model(torch.Tensor(test_img))\n\nf, axarr = plt.subplots(1,5, figsize=(20, 20))\naxarr[0].imshow(test_img[0, 0])\naxarr[1].imshow(test_img[0, 1])\naxarr[2].imshow(test_img[0, 2])\naxarr[3].imshow(test_res[0][0].detach().cpu().numpy() > 0.5)\naxarr[4].imshow((test_res[0][0].detach().cpu().numpy() > 0.5) * (test_img[0, 1] != 0))","metadata":{"execution":{"iopub.status.busy":"2021-08-26T19:39:36.249586Z","iopub.execute_input":"2021-08-26T19:39:36.250011Z","iopub.status.idle":"2021-08-26T19:39:43.748464Z","shell.execute_reply.started":"2021-08-26T19:39:36.249980Z","shell.execute_reply":"2021-08-26T19:39:43.747664Z"},"trusted":true},"execution_count":null,"outputs":[]}]}