{"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":"markdown","source":"# EJERCICIO RADIÓMICA\n## Ángel Guevara y Arturo Sirvent","metadata":{}},{"cell_type":"code","source":"!pip install pyradiomics","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-05-02T20:29:14.359000Z","iopub.execute_input":"2022-05-02T20:29:14.359556Z","iopub.status.idle":"2022-05-02T20:29:27.261117Z","shell.execute_reply.started":"2022-05-02T20:29:14.359519Z","shell.execute_reply":"2022-05-02T20:29:27.260168Z"},"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\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.preprocessing import MinMaxScaler\nfrom sklearn.linear_model import LogisticRegression \nfrom sklearn.metrics import plot_roc_curve, classification_report\n\n\ntrain_path = '../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/'\ntrain_label_path = '../input/rsna-miccai-brain-tumor-radiogenomic-classification/train_labels.csv'","metadata":{"execution":{"iopub.status.busy":"2022-05-02T20:29:28.640591Z","iopub.execute_input":"2022-05-02T20:29:28.640961Z","iopub.status.idle":"2022-05-02T20:29:29.838807Z","shell.execute_reply.started":"2022-05-02T20:29:28.640925Z","shell.execute_reply":"2022-05-02T20:29:29.837098Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_dirs = sorted(os.listdir(train_path))","metadata":{"execution":{"iopub.status.busy":"2022-05-02T20:29:32.286316Z","iopub.execute_input":"2022-05-02T20:29:32.286741Z","iopub.status.idle":"2022-05-02T20:29:32.357983Z","shell.execute_reply.started":"2022-05-02T20:29:32.286707Z","shell.execute_reply":"2022-05-02T20:29:32.357110Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"reader = sitk.ImageSeriesReader()\nreader.LoadPrivateTagsOn()","metadata":{"execution":{"iopub.status.busy":"2022-05-02T20:29:34.716239Z","iopub.execute_input":"2022-05-02T20:29:34.716647Z","iopub.status.idle":"2022-05-02T20:29:34.722619Z","shell.execute_reply.started":"2022-05-02T20:29:34.716610Z","shell.execute_reply":"2022-05-02T20:29:34.721747Z"},"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":"2022-05-02T20:29:36.338649Z","iopub.execute_input":"2022-05-02T20:29:36.339342Z","iopub.status.idle":"2022-05-02T20:29:36.346777Z","shell.execute_reply.started":"2022-05-02T20:29:36.339293Z","shell.execute_reply":"2022-05-02T20:29:36.346015Z"},"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":"2022-05-02T20:29:40.970138Z","iopub.execute_input":"2022-05-02T20:29:40.970713Z","iopub.status.idle":"2022-05-02T20:29:40.976073Z","shell.execute_reply.started":"2022-05-02T20:29:40.970655Z","shell.execute_reply":"2022-05-02T20:29:40.974994Z"},"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":"2022-05-02T20:29:42.584883Z","iopub.execute_input":"2022-05-02T20:29:42.585265Z","iopub.status.idle":"2022-05-02T20:29:42.593856Z","shell.execute_reply.started":"2022-05-02T20:29:42.585233Z","shell.execute_reply":"2022-05-02T20:29:42.592582Z"},"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":"2022-05-02T20:29:47.499084Z","iopub.execute_input":"2022-05-02T20:29:47.500018Z","iopub.status.idle":"2022-05-02T20:29:52.736508Z","shell.execute_reply.started":"2022-05-02T20:29:47.499952Z","shell.execute_reply":"2022-05-02T20:29:52.735354Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Carga y visualización de los datos","metadata":{}},{"cell_type":"code","source":"# Nos quedamos con las primeras 80 imágenes\nim = []\nfor i in range(80):\n    im.append(get_img(i))\n","metadata":{"execution":{"iopub.status.busy":"2022-05-02T20:30:13.553333Z","iopub.execute_input":"2022-05-02T20:30:13.553749Z","iopub.status.idle":"2022-05-02T20:37:45.278997Z","shell.execute_reply.started":"2022-05-02T20:30:13.553715Z","shell.execute_reply":"2022-05-02T20:37:45.274452Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Cargamos las etiquetas\nlabels_train = pd.read_csv(train_label_path)[0:80]\nlabels_train=labels_train.set_index('BraTS21ID')\nlabels_train","metadata":{"execution":{"iopub.status.busy":"2022-05-02T20:40:33.881028Z","iopub.execute_input":"2022-05-02T20:40:33.881707Z","iopub.status.idle":"2022-05-02T20:40:33.968094Z","shell.execute_reply.started":"2022-05-02T20:40:33.881634Z","shell.execute_reply":"2022-05-02T20:40:33.966835Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Aplicamos el modelo de segmentación a las imágenes que nos hemos guardado\ntest_img=[]\ntest_res=[]\nfor i in range(80):\n    test_img.append(np.array([np.moveaxis(np.array(im[i].resize((256, 256))), -1, 0)]))\n    test_res.append(model(torch.Tensor(test_img[i])))","metadata":{"execution":{"iopub.status.busy":"2022-05-02T20:40:37.667410Z","iopub.execute_input":"2022-05-02T20:40:37.667808Z","iopub.status.idle":"2022-05-02T20:41:19.552184Z","shell.execute_reply.started":"2022-05-02T20:40:37.667777Z","shell.execute_reply":"2022-05-02T20:41:19.550910Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Mostramos las primeras 20 imágenes. Las 3 primeras columnas son las imagénes de las tres pruebas que se realizan mientras \n# que la última es la segmentación a partir de esas imágenes.\n# En cada fila tenemos una muestra.\nplt.figure(figsize=(20,70))\nfor i in range(20):\n    plt.subplot(20,4,1+i*4)\n    plt.imshow(test_img[i][0, 0])\n    plt.title('Prueba1')\n    plt.xticks([])\n    plt.yticks([])\n    plt.subplot(20,4,2+i*4)\n    plt.imshow(test_img[i][0, 1])\n    plt.title('Prueba2')\n    plt.xticks([])\n    plt.yticks([])\n    plt.subplot(20,4,3+i*4)\n    plt.imshow(test_img[i][0, 2])\n    plt.title('Prueba3')\n    plt.xticks([])\n    plt.yticks([])\n    plt.subplot(20,4,4+i*4)\n    plt.imshow(test_res[i].detach().numpy()[0,0])\n    plt.title('Segmentación')\n    plt.xticks([])\n    #plt.title()\n    plt.yticks([])\n\nplt.tight_layout()","metadata":{"execution":{"iopub.status.busy":"2022-05-02T20:41:43.458526Z","iopub.execute_input":"2022-05-02T20:41:43.458930Z","iopub.status.idle":"2022-05-02T20:41:50.947830Z","shell.execute_reply.started":"2022-05-02T20:41:43.458897Z","shell.execute_reply":"2022-05-02T20:41:50.946945Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Extracción de características","metadata":{}},{"cell_type":"markdown","source":"Ahora, extraemos diferentes características de las imágenes para poder construir un clasificador. ","metadata":{}},{"cell_type":"code","source":"# Características de Forma 2D\nshape=[]\nfor i in range(80):\n    aux=radiomics.shape2D.RadiomicsShape2D(\n        sitk.GetImageFromArray(test_img[i]), \n        sitk.GetImageFromArray(np.array([\n            test_res[i][0][0].detach().cpu().numpy() > 0.5\n        ]).astype(np.uint8)),\n        force2D=True\n    )\n    shape.append(aux)","metadata":{"execution":{"iopub.status.busy":"2022-05-02T20:42:00.602675Z","iopub.execute_input":"2022-05-02T20:42:00.603329Z","iopub.status.idle":"2022-05-02T20:42:00.794580Z","shell.execute_reply.started":"2022-05-02T20:42:00.603291Z","shell.execute_reply":"2022-05-02T20:42:00.793425Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"MeshSurface=[]\nPixelSurface=[]\nPerimeter=[]\nPerimeterSurfaceRatio=[]\nSpericity=[]\nSphericalDisproportion=[]\nMaximumDiameter=[]\nMajorAxisLength=[]\nMinorAxisLenth=[]\nElongation=[]\n\nfor i in range(80):\n    MeshSurface.append(shape[i].getMeshSurfaceFeatureValue())\n    PixelSurface.append(shape[i].getPixelSurfaceFeatureValue())\n    Perimeter.append(shape[i].getPerimeterFeatureValue())\n    PerimeterSurfaceRatio.append(shape[i].getPerimeterSurfaceRatioFeatureValue())\n    Spericity.append(shape[i].getSphericityFeatureValue())\n    SphericalDisproportion.append(shape[i].getSphericalDisproportionFeatureValue())\n    MaximumDiameter.append(shape[i].getMaximumDiameterFeatureValue())\n    MajorAxisLength.append(shape[i].getMajorAxisLengthFeatureValue())\n    MinorAxisLenth.append(shape[i].getMinorAxisLengthFeatureValue())\n    Elongation.append(shape[i].getElongationFeatureValue())\n","metadata":{"execution":{"iopub.status.busy":"2022-05-02T20:42:04.743591Z","iopub.execute_input":"2022-05-02T20:42:04.744194Z","iopub.status.idle":"2022-05-02T20:42:04.759366Z","shell.execute_reply.started":"2022-05-02T20:42:04.744156Z","shell.execute_reply":"2022-05-02T20:42:04.757041Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Construimos un dataframe con las características de forma\ndf = pd.DataFrame()\ndf['ID'] = [int(i) for i in train_dirs][0:80]\ndf = df.set_index('ID')\ndf['MeshSurface'] = MeshSurface\ndf['PixelSurface'] = PixelSurface\ndf['Perimeter'] = Perimeter\ndf['PerimeterSurfaceRatio'] = PerimeterSurfaceRatio\ndf['Spericity'] = Spericity\ndf['SphericalDisproportion'] = SphericalDisproportion\ndf['MaximumDiameter'] = MaximumDiameter\ndf['MajorAxisLength'] = MajorAxisLength\ndf['MinorAxisLenth'] = MinorAxisLenth\ndf['Elongation'] = Elongation\n#df = df.join(labels_train)\ndf","metadata":{"execution":{"iopub.status.busy":"2022-05-02T20:42:07.947621Z","iopub.execute_input":"2022-05-02T20:42:07.948105Z","iopub.status.idle":"2022-05-02T20:42:08.001630Z","shell.execute_reply.started":"2022-05-02T20:42:07.948063Z","shell.execute_reply":"2022-05-02T20:42:07.999247Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Características de Textura (GLCM)\ntexturas=[]\nfor i in range(80):\n    aux=radiomics.glcm.RadiomicsGLCM(\n        sitk.GetImageFromArray(test_img[i][0,0,:,:].reshape(1, 256, 256)), \n        sitk.GetImageFromArray(np.array([\n            test_res[i][0][0].detach().cpu().numpy() > 0.5\n        ]).astype(np.uint8)),\n        force2D=True\n    )\n    texturas.append(aux)","metadata":{"execution":{"iopub.status.busy":"2022-05-02T20:42:11.305008Z","iopub.execute_input":"2022-05-02T20:42:11.305675Z","iopub.status.idle":"2022-05-02T20:42:11.453270Z","shell.execute_reply.started":"2022-05-02T20:42:11.305638Z","shell.execute_reply":"2022-05-02T20:42:11.451901Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Extraemos las características con la función `.execute()` y las vamos almacenando en la lista de resultados.\nresults=[]\nfor i in range(80):\n    texturas[i].enableAllFeatures()\n    results.append(texturas[i].execute())","metadata":{"execution":{"iopub.status.busy":"2022-05-02T20:42:15.554980Z","iopub.execute_input":"2022-05-02T20:42:15.555731Z","iopub.status.idle":"2022-05-02T20:42:16.247876Z","shell.execute_reply.started":"2022-05-02T20:42:15.555689Z","shell.execute_reply":"2022-05-02T20:42:16.246744Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Construimos un dataframe con las características de texturas\ndf2 = pd.DataFrame(results)\ndf2['ID'] = [int(i) for i in train_dirs][0:80]\ndf2 = df2.set_index('ID')\ndf2","metadata":{"execution":{"iopub.status.busy":"2022-05-02T20:42:20.051323Z","iopub.execute_input":"2022-05-02T20:42:20.051781Z","iopub.status.idle":"2022-05-02T20:42:20.103234Z","shell.execute_reply.started":"2022-05-02T20:42:20.051742Z","shell.execute_reply":"2022-05-02T20:42:20.102307Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Características de primer orden\norden1=[]\nfor i in range(80):\n    aux=radiomics.firstorder.RadiomicsFirstOrder(\n        sitk.GetImageFromArray(test_img[i][0,0,:,:].reshape(1, 256, 256)), \n        sitk.GetImageFromArray(np.array([\n            test_res[i][0][0].detach().cpu().numpy() > 0.5\n        ]).astype(np.uint8)),\n        force2D=True\n    )\n    orden1.append(aux)","metadata":{"execution":{"iopub.status.busy":"2022-05-02T20:42:23.113746Z","iopub.execute_input":"2022-05-02T20:42:23.114262Z","iopub.status.idle":"2022-05-02T20:42:23.272966Z","shell.execute_reply.started":"2022-05-02T20:42:23.114220Z","shell.execute_reply":"2022-05-02T20:42:23.272014Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Extraemos las características con la función `.execute()` y las vamos almacenando en la lista de resultados2.\nresults2=[]\nfor i in range(80):\n    orden1[i].enableAllFeatures()\n    results2.append(orden1[i].execute())","metadata":{"execution":{"iopub.status.busy":"2022-05-02T20:42:28.349728Z","iopub.execute_input":"2022-05-02T20:42:28.350343Z","iopub.status.idle":"2022-05-02T20:42:28.679063Z","shell.execute_reply.started":"2022-05-02T20:42:28.350288Z","shell.execute_reply":"2022-05-02T20:42:28.677919Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Construimos un dataframe con las características de primer orden\ndf3 = pd.DataFrame(results2)\ndf3['ID'] = [int(i) for i in train_dirs][0:80]\ndf3 = df3.set_index('ID')\ndf3","metadata":{"execution":{"iopub.status.busy":"2022-05-02T20:42:32.015315Z","iopub.execute_input":"2022-05-02T20:42:32.015783Z","iopub.status.idle":"2022-05-02T20:42:32.061552Z","shell.execute_reply.started":"2022-05-02T20:42:32.015746Z","shell.execute_reply":"2022-05-02T20:42:32.060241Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Una vez hemos extraído todas las características que queríamos, las juntamos todas en un solo dataframe con el que podamos trabajar. También normalizaremos los datos antes de aplicar ningún modelo.","metadata":{}},{"cell_type":"code","source":"# Unimos los 3 dataframes que contienen las diferentes características de las imágenes\ndf = df.join(df2)\ndf = df.join(df3)\ndf","metadata":{"execution":{"iopub.status.busy":"2022-05-02T20:42:35.216148Z","iopub.execute_input":"2022-05-02T20:42:35.216542Z","iopub.status.idle":"2022-05-02T20:42:35.271149Z","shell.execute_reply.started":"2022-05-02T20:42:35.216500Z","shell.execute_reply":"2022-05-02T20:42:35.269952Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Normalizamos los datos \nscaler = MinMaxScaler()\ndf_norm  = scaler.fit_transform(df)\n","metadata":{"execution":{"iopub.status.busy":"2022-05-02T20:42:38.449330Z","iopub.execute_input":"2022-05-02T20:42:38.449956Z","iopub.status.idle":"2022-05-02T20:42:38.463979Z","shell.execute_reply.started":"2022-05-02T20:42:38.449901Z","shell.execute_reply":"2022-05-02T20:42:38.462709Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Modelos de clasificación","metadata":{}},{"cell_type":"markdown","source":"A continuación, probamos varios clasificadores y los comparamos. ","metadata":{}},{"cell_type":"code","source":"# Dividimos en train y test\n\nx_tr, x_ts, y_tr, y_ts = train_test_split(df_norm, labels_train.to_numpy().squeeze(), test_size=0.2, random_state=0 ) ","metadata":{"execution":{"iopub.status.busy":"2022-05-02T20:42:41.775221Z","iopub.execute_input":"2022-05-02T20:42:41.775816Z","iopub.status.idle":"2022-05-02T20:42:41.788844Z","shell.execute_reply.started":"2022-05-02T20:42:41.775764Z","shell.execute_reply":"2022-05-02T20:42:41.787421Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Regresión Logística\nclf = LogisticRegression(random_state=0)\nclf.fit(x_tr, y_tr)\npreds = clf.predict(x_ts)\n\nprint(classification_report(y_ts, preds))","metadata":{"execution":{"iopub.status.busy":"2022-05-02T20:44:03.497826Z","iopub.execute_input":"2022-05-02T20:44:03.498243Z","iopub.status.idle":"2022-05-02T20:44:03.539965Z","shell.execute_reply.started":"2022-05-02T20:44:03.498208Z","shell.execute_reply":"2022-05-02T20:44:03.538609Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Dibujamos la curva ROC\nplot_roc_curve(clf,x_ts,y_ts)\nplt.plot([0,1],[0,1])","metadata":{"execution":{"iopub.status.busy":"2022-05-02T20:44:23.616240Z","iopub.execute_input":"2022-05-02T20:44:23.616747Z","iopub.status.idle":"2022-05-02T20:44:23.811901Z","shell.execute_reply.started":"2022-05-02T20:44:23.616597Z","shell.execute_reply":"2022-05-02T20:44:23.810365Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Random Forest. \n\nfrom sklearn.ensemble import RandomForestClassifier \n\nRF = RandomForestClassifier(n_estimators=100, max_depth=2, random_state=0)\nRF.fit(x_tr, y_tr)\npreds_RF = RF.predict(x_ts)\n\nprint(classification_report(y_ts, preds_RF))","metadata":{"execution":{"iopub.status.busy":"2022-05-02T20:43:47.718040Z","iopub.execute_input":"2022-05-02T20:43:47.718428Z","iopub.status.idle":"2022-05-02T20:43:47.964984Z","shell.execute_reply.started":"2022-05-02T20:43:47.718396Z","shell.execute_reply":"2022-05-02T20:43:47.963855Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Dibujamos la curva ROC del RandomForest\nplot_roc_curve(RF,x_ts,y_ts)\nplt.plot([0,1],[0,1])","metadata":{"execution":{"iopub.status.busy":"2022-05-02T20:44:35.347243Z","iopub.execute_input":"2022-05-02T20:44:35.347816Z","iopub.status.idle":"2022-05-02T20:44:35.565336Z","shell.execute_reply.started":"2022-05-02T20:44:35.347764Z","shell.execute_reply":"2022-05-02T20:44:35.563864Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# SVC\n\nfrom sklearn.svm import SVC \n\nsvc = SVC(kernel='rbf', gamma=1, random_state=0)\nsvc.fit(x_tr, y_tr)\npreds_svc = svc.predict(x_ts)\n\nprint(classification_report(y_ts, preds_svc))","metadata":{"execution":{"iopub.status.busy":"2022-05-02T20:44:47.730342Z","iopub.execute_input":"2022-05-02T20:44:47.730767Z","iopub.status.idle":"2022-05-02T20:44:47.746190Z","shell.execute_reply.started":"2022-05-02T20:44:47.730735Z","shell.execute_reply":"2022-05-02T20:44:47.745369Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Dibujamos la curva ROC del SVC\nplot_roc_curve(svc,x_ts,y_ts)\nplt.plot([0,1],[0,1])","metadata":{"execution":{"iopub.status.busy":"2022-05-02T20:44:49.906781Z","iopub.execute_input":"2022-05-02T20:44:49.907436Z","iopub.status.idle":"2022-05-02T20:44:50.102670Z","shell.execute_reply.started":"2022-05-02T20:44:49.907385Z","shell.execute_reply":"2022-05-02T20:44:50.101723Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Vemos cuáles son las características más importantes según el modelo de Regresión Logística, es decir, cuáles son los biomarcadores más relevantes para la clasificación de las imágenes. Lo hacemos solo sobre este clasificador porque parecen el más relevante atendiendo a los resultados obtenidos de precisión y recall.","metadata":{}},{"cell_type":"code","source":"# Regresión Logística : Importancia de características\n# Usamos la función `permutation_importance` para calcular la importancia de cada característica. Básicamente lo que hace esta \n# función es barajar una de las características y ver como afecta eso al rendimiento del clasificador. Esto lo hace para cada \n# variable.\nfrom sklearn.inspection import permutation_importance\n\nimportancias = permutation_importance(clf, x_tr, y_tr, random_state=0)\n\n# Nos quedamos solo con las características más relevantes.\nindx = np.where(importancias['importances_mean'] >= abs(0.006))\nplt.bar(x= df.columns[indx], height=importancias['importances_mean'][indx])\nplt.xticks(rotation=90)\n","metadata":{"execution":{"iopub.status.busy":"2022-05-02T20:45:02.602991Z","iopub.execute_input":"2022-05-02T20:45:02.603401Z","iopub.status.idle":"2022-05-02T20:45:02.964145Z","shell.execute_reply.started":"2022-05-02T20:45:02.603367Z","shell.execute_reply":"2022-05-02T20:45:02.963083Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Parece que las características más importantes a la hora de la clasificación son el Perímetro, el Diámetro Máximo y el Coeficiente de Correlación Máximo (MCC). Después de estas, tenemos otras variables como el Contraste, la Entropía y la Diferencia de Varianzas que están un nivel por debajo en cuánto a la importancia.","metadata":{}},{"cell_type":"markdown","source":"## Conclusiones","metadata":{}},{"cell_type":"markdown","source":"A nosotros nos interesa que el modelo sea capaz de acertar el máximo de casos positivos, es decir, que acierte el máximo de casos en los que hay un tumor en el órgano sobre el total de casos con tumor. Esto se corresponde con que el modelo tenga un recall alto. Además también querremos que de todos los casos en los que decimos que hay un tumor acertemos el máximo, es decir que tenga una precisión alta también. \n\nEn principio, nos interesa más lo dicho anteriormente sobre la clase 1, pero tampoco podemos dejar de lado la clase 0. Por tanto, lo interesante será un recall y una precisión alta en ambas clases. \n\nObservando las métricas de los tres modelos propuestos juntos con sus curvas ROC, el modelo de Regresión Logística nos parece el mejor puesto que tiene los mejores valores de recall y precisión para la clase 1 (que es la más importante de detectar a priori) y, además, también tiene los mejores resultados para la clase 0. \n\nComparando con el resto de modelos, vemos que el SVC es peor en general tanto para la clase 0 como para la clase 1, y que el Random Forest tiene los peores resultados para la clase 1. Por esto y por lo anterior, la regresión logística funcionaría mejor en este caso de clasificación.\n\n","metadata":{}},{"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.78586Z","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.78626Z","shell.execute_reply":"2021-08-26T19:38:19.34392Z"},"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.24998Z","shell.execute_reply":"2021-08-26T19:39:43.747664Z"},"trusted":true},"execution_count":null,"outputs":[]}]}