{"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":"## Biomed DataChallenge\n\n#### Este ejercicio consiste en clasificar tumores cerebrales de dos tipos diferentes a partir de características de radiómica que se extraen de imágenes de resonancia magnética. Para ello se utilizarán las imágenes del data challenge RSNA MICCAI Brain Tumor Radiogenomic classification. \n","metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","execution":{"iopub.execute_input":"2022-11-02T13:26:07.811286Z","iopub.status.busy":"2022-11-02T13:26:07.810440Z","iopub.status.idle":"2022-11-02T13:26:07.833856Z","shell.execute_reply":"2022-11-02T13:26:07.833179Z","shell.execute_reply.started":"2022-11-02T13:26:07.811191Z"}}},{"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)","metadata":{"execution":{"iopub.status.busy":"2023-05-27T23:08:07.791025Z","iopub.execute_input":"2023-05-27T23:08:07.791769Z","iopub.status.idle":"2023-05-27T23:08:07.795986Z","shell.execute_reply.started":"2023-05-27T23:08:07.791733Z","shell.execute_reply":"2023-05-27T23:08:07.795239Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install pyradiomics","metadata":{"scrolled":true,"execution":{"iopub.status.busy":"2023-05-27T23:07:10.295056Z","iopub.execute_input":"2023-05-27T23:07:10.295537Z","iopub.status.idle":"2023-05-27T23:08:02.815727Z","shell.execute_reply.started":"2023-05-27T23:07:10.295496Z","shell.execute_reply":"2023-05-27T23:08:02.814257Z"},"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 radiomics import glcm, imageoperations, setVerbosity\nfrom skimage.measure import label, regionprops\nimport scipy\n#from skimage import morphology\nfrom scipy import ndimage\nsetVerbosity(60)\n\ndataset_path = '../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/'","metadata":{"execution":{"iopub.status.busy":"2023-05-27T23:08:02.823279Z","iopub.execute_input":"2023-05-27T23:08:02.823595Z","iopub.status.idle":"2023-05-27T23:08:07.788782Z","shell.execute_reply.started":"2023-05-27T23:08:02.823560Z","shell.execute_reply":"2023-05-27T23:08:07.787723Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# carga lista ordenada de los archivos\ndataset_dirs = sorted(os.listdir(dataset_path))","metadata":{"execution":{"iopub.status.busy":"2023-05-27T23:08:07.797507Z","iopub.execute_input":"2023-05-27T23:08:07.798119Z","iopub.status.idle":"2023-05-27T23:08:07.993783Z","shell.execute_reply.started":"2023-05-27T23:08:07.798070Z","shell.execute_reply":"2023-05-27T23:08:07.992611Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"reader = sitk.ImageSeriesReader()\nreader.LoadPrivateTagsOn()","metadata":{"execution":{"iopub.status.busy":"2023-05-15T08:59:22.668753Z","iopub.execute_input":"2023-05-15T08:59:22.669084Z","iopub.status.idle":"2023-05-15T08:59:22.675957Z","shell.execute_reply.started":"2023-05-15T08:59:22.669056Z","shell.execute_reply":"2023-05-15T08:59:22.674966Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Registro y normalización","metadata":{}},{"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    resampler.SetTransform(sitk.AffineTransform(image.GetDimension()))\n    resampler.SetOutputSpacing(ref_image.GetSpacing())\n    resampler.SetSize(ref_image.GetSize())\n    resampler.SetOutputDirection(ref_image.GetDirection())\n    resampler.SetOutputOrigin(ref_image.GetOrigin())\n    resampler.SetDefaultPixelValue(image.GetPixelIDValue())\n    resamped_image = resampler.Execute(image)\n    \n    return resamped_image","metadata":{"execution":{"iopub.status.busy":"2023-05-15T08:59:22.677175Z","iopub.execute_input":"2023-05-15T08:59:22.677610Z","iopub.status.idle":"2023-05-15T08:59:22.688987Z","shell.execute_reply.started":"2023-05-15T08:59:22.677580Z","shell.execute_reply":"2023-05-15T08:59:22.688152Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def normalizeImage(data):\n    return (data - np.min(data)) / (np.max(data) - np.min(data))","metadata":{"execution":{"iopub.status.busy":"2023-05-15T08:59:22.690113Z","iopub.execute_input":"2023-05-15T08:59:22.690590Z","iopub.status.idle":"2023-05-15T08:59:22.700990Z","shell.execute_reply.started":"2023-05-15T08:59:22.690561Z","shell.execute_reply":"2023-05-15T08:59:22.699982Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ndef get_img(index):\n    \n    ind = dataset_dirs[index]\n    \n    if ind not in rem_cases:\n    \n        filenamesDICOM = reader.GetGDCMSeriesFileNames(f'{dataset_path}/{dataset_dirs[index]}/T1w')\n        reader.SetFileNames(filenamesDICOM)\n        t1_sitk = reader.Execute()\n\n        filenamesDICOM = reader.GetGDCMSeriesFileNames(f'{dataset_path}/{dataset_dirs[index]}/FLAIR')\n        reader.SetFileNames(filenamesDICOM)\n        flair_sitk = reader.Execute()\n\n        filenamesDICOM = reader.GetGDCMSeriesFileNames(f'{dataset_path}/{dataset_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 = normalizeImage(sitk.GetArrayFromImage(t1_sitk))\n        flair_resampled_array = normalizeImage(sitk.GetArrayFromImage(flair_resampled))\n        t1wce_resampled_array = normalizeImage(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, ind\n    \n    else: \n        return -1","metadata":{"execution":{"iopub.status.busy":"2023-05-15T08:59:22.702160Z","iopub.execute_input":"2023-05-15T08:59:22.703024Z","iopub.status.idle":"2023-05-15T08:59:22.717372Z","shell.execute_reply.started":"2023-05-15T08:59:22.702982Z","shell.execute_reply":"2023-05-15T08:59:22.716332Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Segmentación","metadata":{}},{"cell_type":"markdown","source":"Para segmentar los tumores, vamos a utilizar una red preentrenada. Podéis encontrar la documentación de esta red en https://github.com/mateuszbuda/brain-segmentation-pytorch: ","metadata":{}},{"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":"2023-05-15T08:59:22.719022Z","iopub.execute_input":"2023-05-15T08:59:22.719899Z","iopub.status.idle":"2023-05-15T08:59:23.315832Z","shell.execute_reply.started":"2023-05-15T08:59:22.719866Z","shell.execute_reply":"2023-05-15T08:59:23.314511Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Para segmentar una imagen, podemos ejecutar el código siguiente. Fijaros en la necesidad de cambiar los ejes. Una vez la red nos devuelve la ROI, debemos binarizarla. En este ejemplo se ha elegido un umbral de 0.5. Para cambiar las imágenes de formato sitk imagen a array y viceversa utilizamos `sitk.GetArrayFromImage` o `sitk.GetImageFromArray`. ","metadata":{}},{"cell_type":"markdown","source":"Prueba con casos diferentes, por ejemplo 1,2,5 y 48 y visualiza los resultados de la segmentación.","metadata":{}},{"cell_type":"code","source":"rem_cases = ['00003','00009','00011','00014','00044','00054','00061','00078','00085','00098','00105','00107',\n             '00109','00110','00121','00140','00148','00149','00170','00186','00194','00204','00206','00212',\n             '00218','00222','00236','00240','00243','00251','00283','00288','00293','00296','00305','00311',\n             '00313','00316','00317','00321','00329','00332','00340','00344','00350','00351','00353','00364',\n             '00379','00390','00413','00433','00440','00441','00442','00445','00451','00464','00495','00496',\n             '00500','00502','00506','00514','00518','00538','00540','00547','00550','00555','00568','00572',\n             '00578','00582','00587','00596','00602','00604','00613','00615','00641','00645','00650','00651',\n             '00654','00659','00674','00679','00708','00715','00728','00732','00736','00737','00750','00758',\n             '00767','00773','00782','00784','00811']","metadata":{"execution":{"iopub.status.busy":"2023-05-15T09:00:34.469569Z","iopub.execute_input":"2023-05-15T09:00:34.469960Z","iopub.status.idle":"2023-05-15T09:00:34.477861Z","shell.execute_reply.started":"2023-05-15T09:00:34.469931Z","shell.execute_reply":"2023-05-15T09:00:34.476730Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"case=1\n\nim, ind = get_img(case)\ntest_img = np.array([np.moveaxis(np.array(im.resize((256, 256))), -1, 0)])\ntest_res = model(torch.Tensor(test_img))\n\nimage = sitk.GetImageFromArray(test_img[0,1,:,:].reshape(1, 256, 256))\nmask = sitk.GetImageFromArray(np.array([\n        test_res[0][0].detach().cpu().numpy() > 0.5\n    ]).astype(np.uint8))\n\nim = sitk.GetArrayFromImage(image)\nmsk = sitk.GetArrayFromImage(mask)\n\nfig, (ax1, ax2) = plt.subplots(1,2)\nax1.imshow(im[0, :, :])\nax1.set_title('Flair')\nax1.axis('off')\nax2.imshow(msk[0, :, :])\nax2.set_title('Tumor segmentation')\nax2.axis('off')","metadata":{"execution":{"iopub.status.busy":"2023-05-15T09:01:14.877791Z","iopub.execute_input":"2023-05-15T09:01:14.879001Z","iopub.status.idle":"2023-05-15T09:01:20.640249Z","shell.execute_reply.started":"2023-05-15T09:01:14.878963Z","shell.execute_reply":"2023-05-15T09:01:20.639225Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Realizamos una limpieza de la máscara quitando los \"blobs\" con áreas muy pequeñas sueltas. En el caso de segmentaciones con varios \"blobs\" grandes, nos quedamos con el de mayor área. Además rellenamos agujeros de dentro de la segmentación con la función `binary_fill_holes`. Finalmente realizamos un \"crop\" a la imagen, quedándonos únicamente con la zona de interés (tumor). Los casos con mala segmentació (i.e. caso 2) los quitaremos más adelante.","metadata":{}},{"cell_type":"code","source":"# Ha dicho que esta no la utilicemos porque es la que da error\ndef remove_small_blobs(binary_mask: np.ndarray, min_size: int = 0):\n    \"\"\"\n    Removes from the input mask all the blobs having less than N adjacent pixels.\n    We set the small objects to the background label 0.\n    \"\"\"\n    if min_size > 0:\n        dtype = binary_mask.dtype\n        binary_mask = morphology.remove_small_objects(binary_mask.astype(bool), min_size=min_size)\n        binary_mask = binary_mask.astype(dtype)\n    return binary_mask","metadata":{"execution":{"iopub.status.busy":"2023-05-15T09:01:20.641870Z","iopub.execute_input":"2023-05-15T09:01:20.642177Z","iopub.status.idle":"2023-05-15T09:01:20.649109Z","shell.execute_reply.started":"2023-05-15T09:01:20.642150Z","shell.execute_reply":"2023-05-15T09:01:20.647983Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def keep_bigger_blob(binary_mask: np.ndarray):\n    \n    labels_mask = label(binary_mask)     \n    regions = regionprops(labels_mask)\n    regions.sort(key=lambda x: x.area, reverse=True)\n    if len(regions) > 1:\n        for rg in regions[1:]:\n            labels_mask[rg.coords[:,0],rg.coords[:,1],rg.coords[:,2]] = 0\n    labels_mask[labels_mask!=0] = 1\n    mask = labels_mask\n    \n    return mask","metadata":{"execution":{"iopub.status.busy":"2023-05-15T09:01:20.650185Z","iopub.execute_input":"2023-05-15T09:01:20.650497Z","iopub.status.idle":"2023-05-15T09:01:20.660315Z","shell.execute_reply.started":"2023-05-15T09:01:20.650462Z","shell.execute_reply":"2023-05-15T09:01:20.659498Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"image0 = sitk.GetImageFromArray(test_img[0,1,:,:].reshape(1, 256, 256))\n\nmask = np.array([test_res[0][0].detach().cpu().numpy() > 0.5]).astype(np.uint8)\n#mask_clean = remove_small_blobs(mask, 100)\nmask_clean2 = keep_bigger_blob(mask)\nmask_cl = ndimage.binary_fill_holes(mask_clean2[0,:,:]).astype(int)\nmask_clean2[0,:,:] = mask_cl\nmask_clean2 = sitk.GetImageFromArray(mask_clean2)\n# Crop the image\n# bb is the bounding box, upon which the image and mask are cropped\nbb, correctedMask = imageoperations.checkMask(image0, mask_clean2, label=1)\nif correctedMask is not None:\n    mask_clean2 = correctedMask\nimage, mask = imageoperations.cropToTumorMask(image0, mask_clean2, bb)","metadata":{"execution":{"iopub.status.busy":"2023-05-15T09:01:20.662484Z","iopub.execute_input":"2023-05-15T09:01:20.662807Z","iopub.status.idle":"2023-05-15T09:01:20.730926Z","shell.execute_reply.started":"2023-05-15T09:01:20.662779Z","shell.execute_reply":"2023-05-15T09:01:20.729775Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"msk0 = mask\nmsk = sitk.GetArrayFromImage(mask)\n\nT1 = test_img[0,0,:,:].reshape(1, 256, 256)\nFlair = test_img[0,1,:,:].reshape(1, 256, 256)\nT1wce = test_img[0,2,:,:].reshape(1, 256, 256)\n\nf, axarr = plt.subplots(1,4, figsize=(20, 20))\naxarr[0].imshow(T1[0,:,:])\naxarr[0].set_title('T1 %s' %(ind))\naxarr[0].axis('off')\naxarr[1].imshow(Flair[0,:,:])\naxarr[1].set_title('Flair %s' %(ind))\naxarr[1].axis('off')\naxarr[2].imshow(T1wce[0,:,:])\naxarr[2].set_title('T1wce %s' %(ind))\naxarr[2].axis('off')\n#axarr[3].imshow(mask_clean[0,:,:])\n#axarr[3].set_title('Mask1 %s' %(ind))\n#axarr[3].axis('off')\naxarr[3].imshow(msk[0,:,:])\naxarr[3].set_title('Mask2 %s' %(ind))\naxarr[3].axis('off')","metadata":{"execution":{"iopub.status.busy":"2023-05-15T09:01:20.732062Z","iopub.execute_input":"2023-05-15T09:01:20.732372Z","iopub.status.idle":"2023-05-15T09:01:21.300314Z","shell.execute_reply.started":"2023-05-15T09:01:20.732345Z","shell.execute_reply":"2023-05-15T09:01:21.299438Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Extraemos las características de radiómica.","metadata":{}},{"cell_type":"markdown","source":"Para extraer las características, podemos elegir cualquiera de las 3 secuencias disponibles y la máscara de segmentación obtenida. En este ejemplo vamos a extraer las características de primer orden de una imagen utilizando la secuencia *Flair*. Fíjate que la librería pyRadiomics necesita imágenes sitk como entrada. Revisa la librería pyradiomics para ver todas las características que puedes extraer (https://pyradiomics.readthedocs.io/en/latest/features.html).","metadata":{}},{"cell_type":"code","source":"firstOrderFeatures = radiomics.firstorder.RadiomicsFirstOrder(image, mask)\n# Extraemos todas las características de la categoría First Order:\nfirstOrderFeatures.enableAllFeatures()\n# Calcula las características y hace el print\nprint('Calculating first order features...',)\nresult = firstOrderFeatures.execute()\nprint('done')\nprint(result)","metadata":{"execution":{"iopub.status.busy":"2023-05-15T09:01:21.301768Z","iopub.execute_input":"2023-05-15T09:01:21.302349Z","iopub.status.idle":"2023-05-15T09:01:21.316006Z","shell.execute_reply.started":"2023-05-15T09:01:21.302318Z","shell.execute_reply":"2023-05-15T09:01:21.314915Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Creamos un dataframe dónde guardaremos todas las características.","metadata":{}},{"cell_type":"code","source":"data = pd.DataFrame(list(result.values()))\ndata = data.T  \ndata.columns = list(result.keys())\ndata","metadata":{"execution":{"iopub.status.busy":"2023-05-15T09:01:21.317412Z","iopub.execute_input":"2023-05-15T09:01:21.318038Z","iopub.status.idle":"2023-05-15T09:01:21.361064Z","shell.execute_reply.started":"2023-05-15T09:01:21.317999Z","shell.execute_reply":"2023-05-15T09:01:21.360083Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# TAREA","metadata":{}},{"cell_type":"markdown","source":"1.- Extraer todas las características posibles con la librería pyradiomics.\n\n2.- Implementar estrategias de selección de características para reducir el número de características.\n\n3.- Aplicar un clasificador que permita distinguir entre los dos tipos de tumores.\n\n4.- Guardar las predicciones del dataset de test (los 100 últimos casos) en un csv y subirlo a la plataforma cuyo enlace se encuentra en el AV.\n\n5.- Escribir un documento en formato artículo científico (template en el AV) con los siguientes apartados: 1.Introducción, 2.Material y métodos, 3.Resultados, 4.Discusión y 5.Conclusiones.\n\n6.- Subir el artículo y el notebook al AV.","metadata":{}},{"cell_type":"markdown","source":"Lista casos que se deben excluir debido a una mala segmentación.","metadata":{}},{"cell_type":"code","source":"rem_cases = ['00003','00009','00011','00014','00044','00054','00061','00078','00085','00098','00105','00107',\n             '00109','00110','00121','00140','00148','00149','00170','00186','00194','00204','00206','00212',\n             '00218','00222','00236','00240','00243','00251','00283','00288','00293','00296','00305','00311',\n             '00313','00316','00317','00321','00329','00332','00340','00344','00350','00351','00353','00364',\n             '00379','00390','00413','00433','00440','00441','00442','00445','00451','00464','00495','00496',\n             '00500','00502','00506','00514','00518','00538','00540','00547','00550','00555','00568','00572',\n             '00578','00582','00587','00596','00602','00604','00613','00615','00641','00645','00650','00651',\n             '00654','00659','00674','00679','00708','00715','00728','00732','00736','00737','00750','00758',\n             '00767','00773','00782','00784','00811']","metadata":{"execution":{"iopub.status.busy":"2023-05-15T09:01:21.362345Z","iopub.execute_input":"2023-05-15T09:01:21.362705Z","iopub.status.idle":"2023-05-15T09:01:21.372178Z","shell.execute_reply.started":"2023-05-15T09:01:21.362676Z","shell.execute_reply":"2023-05-15T09:01:21.371109Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\n# carga lista ordenada de los archivos\ndataset_dirs = sorted(os.listdir(dataset_path))\ndataset_dirs2 = [i for i in dataset_dirs if i not in rem_cases]","metadata":{"execution":{"iopub.status.busy":"2023-05-15T09:01:21.373485Z","iopub.execute_input":"2023-05-15T09:01:21.373778Z","iopub.status.idle":"2023-05-15T09:01:21.389166Z","shell.execute_reply.started":"2023-05-15T09:01:21.373752Z","shell.execute_reply":"2023-05-15T09:01:21.387985Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_img(index):\n    \n    ind = dataset_dirs2[index]\n    \n    \n    filenamesDICOM = reader.GetGDCMSeriesFileNames(f'{dataset_path}/{dataset_dirs2[index]}/T1w')\n    reader.SetFileNames(filenamesDICOM)\n    t1_sitk = reader.Execute()\n\n    filenamesDICOM = reader.GetGDCMSeriesFileNames(f'{dataset_path}/{dataset_dirs2[index]}/FLAIR')\n    reader.SetFileNames(filenamesDICOM)\n    flair_sitk = reader.Execute()\n\n    filenamesDICOM = reader.GetGDCMSeriesFileNames(f'{dataset_path}/{dataset_dirs2[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 = normalizeImage(sitk.GetArrayFromImage(t1_sitk))\n    flair_resampled_array = normalizeImage(sitk.GetArrayFromImage(flair_resampled))\n    t1wce_resampled_array = normalizeImage(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\n","metadata":{"execution":{"iopub.status.busy":"2023-05-15T09:01:21.392147Z","iopub.execute_input":"2023-05-15T09:01:21.392525Z","iopub.status.idle":"2023-05-15T09:01:21.402928Z","shell.execute_reply.started":"2023-05-15T09:01:21.392492Z","shell.execute_reply":"2023-05-15T09:01:21.402111Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"imagenes = []\ndfs =  []\n\nfor i in range(len(dataset_dirs2)):\n    imagen = get_img(i)\n    test_img = np.array([np.moveaxis(np.array(imagen.resize((256, 256))), -1, 0)])\n    test_res = model(torch.Tensor(test_img))\n\n    image = sitk.GetImageFromArray(test_img[0,1,:,:].reshape(1, 256, 256))\n    mask = sitk.GetImageFromArray(np.array([\n            test_res[0][0].detach().cpu().numpy() > 0.5\n        ]).astype(np.uint8))\n\n    im = sitk.GetArrayFromImage(image)\n    msk = sitk.GetArrayFromImage(mask)\n\n    image0 = sitk.GetImageFromArray(test_img[0,1,:,:].reshape(1, 256, 256))\n\n    mask = np.array([test_res[0][0].detach().cpu().numpy() > 0.5]).astype(np.uint8)\n    #mask_clean = remove_small_blobs(mask, 100)\n    mask_clean2 = keep_bigger_blob(mask)\n    mask_cl = ndimage.binary_fill_holes(mask_clean2[0,:,:]).astype(int)\n    mask_clean2[0,:,:] = mask_cl\n    mask_clean2 = sitk.GetImageFromArray(mask_clean2)\n    # Crop the image\n    # bb is the bounding box, upon which the image and mask are cropped\n    bb, correctedMask = imageoperations.checkMask(image0, mask_clean2, label=1)\n    if correctedMask is not None:\n        mask_clean2 = correctedMask\n    image, mask = imageoperations.cropToTumorMask(image0, mask_clean2, bb)\n\n    msk0 = mask\n    msk = sitk.GetArrayFromImage(mask)\n\n    T1 = test_img[0,0,:,:].reshape(1, 256, 256)\n    Flair = test_img[0,1,:,:].reshape(1, 256, 256)\n    T1wce = test_img[0,2,:,:].reshape(1, 256, 256)\n\n    # hacer todo excpto2D\n    firstOrderFeatures = radiomics.firstorder.RadiomicsFirstOrder(image, mask)\n    # Extraemos todas las características de la categoría First Order:\n    firstOrderFeatures.enableAllFeatures()\n    # Calcula las características y hace el print\n    result = firstOrderFeatures.execute()\n    \n    data = pd.DataFrame(list(result.values()))\n    data = data.T  \n    data.columns = list(result.keys())\n    \n    dfs.append(data)\n    print(i)\ndfs = pd.concat(dfs) # de lista a pd/df\n    \n    ","metadata":{"execution":{"iopub.status.busy":"2023-05-15T09:01:21.404100Z","iopub.execute_input":"2023-05-15T09:01:21.405155Z","iopub.status.idle":"2023-05-15T09:55:19.262900Z","shell.execute_reply.started":"2023-05-15T09:01:21.405121Z","shell.execute_reply":"2023-05-15T09:55:19.261099Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dfs","metadata":{"execution":{"iopub.status.busy":"2023-05-15T09:55:19.265590Z","iopub.execute_input":"2023-05-15T09:55:19.265989Z","iopub.status.idle":"2023-05-15T09:55:19.312194Z","shell.execute_reply.started":"2023-05-15T09:55:19.265948Z","shell.execute_reply":"2023-05-15T09:55:19.311067Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"np.array([np.moveaxis(np.array(get_img(0).resize((256, 256))), -1, 0)]).shape","metadata":{"execution":{"iopub.status.busy":"2023-05-15T11:14:14.319353Z","iopub.execute_input":"2023-05-15T11:14:14.320137Z","iopub.status.idle":"2023-05-15T11:14:17.526657Z","shell.execute_reply.started":"2023-05-15T11:14:14.320092Z","shell.execute_reply":"2023-05-15T11:14:17.525446Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"imagenes = []\ndfs1st =  []\ndfs3d = []\ndfsGLCM = []\ndfsGLSZM = []\ndfsGLRLM = []\ndfsNGTDM = []\ndfsGLDM = []\n\ndf_final = []\n\n\nfor i in range(len(dataset_dirs2)):\n# for i in range(2):\n    imagen = get_img(i)\n    test_img = np.array([np.moveaxis(np.array(imagen.resize((256, 256))), -1, 0)])\n#     test_res = model(torch.Tensor(test_img))\n\n#     image = sitk.GetImageFromArray(test_img[0,1,:,:].reshape(1, 256, 256))\n#     mask = sitk.GetImageFromArray(np.array([\n#             test_res[0][0].detach().cpu().numpy() > 0.5\n#         ]).astype(np.uint8))\n\n#     im = sitk.GetArrayFromImage(image)\n#     msk = sitk.GetArrayFromImage(mask)\n\n#     image0 = sitk.GetImageFromArray(test_img[0,1,:,:].reshape(1, 256, 256))\n\n#     mask = np.array([test_res[0][0].detach().cpu().numpy() > 0.5]).astype(np.uint8)\n#     #mask_clean = remove_small_blobs(mask, 100)\n#     mask_clean2 = keep_bigger_blob(mask)\n#     mask_cl = ndimage.binary_fill_holes(mask_clean2[0,:,:]).astype(int)\n#     mask_clean2[0,:,:] = mask_cl\n#     mask_clean2 = sitk.GetImageFromArray(mask_clean2)\n#     # Crop the image\n#     # bb is the bounding box, upon which the image and mask are cropped\n#     bb, correctedMask = imageoperations.checkMask(image0, mask_clean2, label=1)\n#     if correctedMask is not None:\n#         mask_clean2 = correctedMask\n#     image, mask = imageoperations.cropToTumorMask(image0, mask_clean2, bb)\n\n#     msk0 = mask\n#     msk = sitk.GetArrayFromImage(mask)\n\n    T1 = test_img[0,0,:,:].reshape(1, 256, 256)\n    Flair = test_img[0,1,:,:].reshape(1, 256, 256)\n    T1wce = test_img[0,2,:,:].reshape(1, 256, 256)\n    \n    seqs = [T1, Flair, T1wce]\n    \n    dfs1st = []\n    dfs3d = []\n    dfsGLCM = []\n    dfsGLSZM = []\n    dfsGLRLM = []\n    dfsNGTDM = []\n    dfsGLDM = []\n\n    for j in range(3):\n    \n        test_res = model(torch.Tensor(test_img))\n\n        image = sitk.GetImageFromArray(seqs[j])\n        mask = sitk.GetImageFromArray(np.array([\n                test_res[0][0].detach().cpu().numpy() > 0.5\n            ]).astype(np.uint8))\n\n        im = sitk.GetArrayFromImage(image)\n        msk = sitk.GetArrayFromImage(mask)\n\n        image0 = sitk.GetImageFromArray(seqs[j])\n\n        mask = np.array([test_res[0][0].detach().cpu().numpy() > 0.5]).astype(np.uint8)\n        #mask_clean = remove_small_blobs(mask, 100)\n        mask_clean2 = keep_bigger_blob(mask)\n        mask_cl = ndimage.binary_fill_holes(mask_clean2[0,:,:]).astype(int)\n        mask_clean2[0,:,:] = mask_cl\n        mask_clean2 = sitk.GetImageFromArray(mask_clean2)\n        # Crop the image\n        # bb is the bounding box, upon which the image and mask are cropped\n        bb, correctedMask = imageoperations.checkMask(image0, mask_clean2, label=1)\n        if correctedMask is not None:\n            mask_clean2 = correctedMask\n        image, mask = imageoperations.cropToTumorMask(image0, mask_clean2, bb)\n\n        msk0 = mask\n        msk = sitk.GetArrayFromImage(mask)\n\n\n        ##### ------  CARACTERISTICAS ------------------\n\n        # -------------- firstOrderFeatures\n\n        firstOrderFeatures = radiomics.firstorder.RadiomicsFirstOrder(image, mask)\n        # Extraemos todas las características de la categoría First Order:\n        firstOrderFeatures.enableAllFeatures()\n        # Calcula las características y hace el print\n        result1st = firstOrderFeatures.execute()\n\n        data1st = pd.DataFrame(list(result1st.values()))\n        data1st = data1st.T  \n        data1st.columns = list(result1st.keys())\n\n        dfs1st.append(data1st)\n\n        # -------------- shape 3d\n\n        shape3d = radiomics.shape.RadiomicsShape(image, mask)\n\n        shape3d.enableAllFeatures()\n        result3d = shape3d.execute()\n\n        data3d = pd.DataFrame(list(result3d.values()))\n        data3d = data3d.T  \n        data3d.columns = list(result3d.keys())\n\n        dfs3d.append(data3d)\n\n        # -------------- GLCM\n\n        GLCM = radiomics.glcm.RadiomicsGLCM(image, mask)\n\n        GLCM.enableAllFeatures()\n        resultGLCM = GLCM.execute()\n\n        dataGLCM = pd.DataFrame(list(resultGLCM.values()))\n        dataGLCM = dataGLCM.T  \n        dataGLCM.columns = list(resultGLCM.keys())\n\n        dfsGLCM.append(dataGLCM)\n\n        # -------------- GLSZM\n        GLSZM = radiomics.glszm.RadiomicsGLSZM(image, mask)\n\n        GLSZM.enableAllFeatures()\n        resultGLSZM = GLSZM.execute()\n\n        dataGLSZM = pd.DataFrame(list(resultGLSZM.values()))\n        dataGLSZM = dataGLSZM.T  \n        dataGLSZM.columns = list(resultGLSZM.keys())\n\n        dfsGLSZM.append(dataGLSZM)\n\n        # -------------- GLRLM\n\n        GLRLM = radiomics.glrlm.RadiomicsGLRLM(image, mask)\n\n        GLRLM.enableAllFeatures()\n        resultGLRLM = GLRLM.execute()\n\n        dataGLRLM = pd.DataFrame(list(resultGLRLM.values()))\n        dataGLRLM = dataGLRLM.T  \n        dataGLRLM.columns = list(resultGLRLM.keys())\n\n        dfsGLRLM.append(dataGLRLM)\n\n        # -------------- NGTDM\n\n        NGTDM = radiomics.ngtdm.RadiomicsNGTDM(image, mask)\n\n        NGTDM.enableAllFeatures()\n        resultNGTDM = NGTDM.execute()\n\n        dataNGTDM = pd.DataFrame(list(resultNGTDM.values()))\n        dataNGTDM = dataNGTDM.T  \n        dataNGTDM.columns = list(resultNGTDM.keys())\n\n        dfsNGTDM.append(dataNGTDM)\n\n        # -------------- GLDM\n\n\n        GLDM = radiomics.gldm.RadiomicsGLDM(image, mask)\n\n        GLDM.enableAllFeatures()\n        resultGLDM = GLDM.execute()\n\n        dataGLDM = pd.DataFrame(list(resultGLDM.values()))\n        dataGLDM = dataGLDM.T  \n        dataGLDM.columns = list(resultGLDM.keys())\n\n        dfsGLDM.append(dataGLDM)\n\n\n        print('seq',j)\n    \n    dfs1st = pd.concat(dfs1st)\n    dfs3d = pd.concat(dfs3d)\n    dfsGLCM = pd.concat(dfsGLCM)\n    dfsGLSZM = pd.concat(dfsGLSZM)\n    dfsGLRLM = pd.concat(dfsGLRLM)\n    dfsNGTDM = pd.concat(dfsNGTDM)\n    dfsGLDM = pd.concat(dfsGLDM)\n    DFS = [dfs1st, dfs3d, dfsGLCM, dfsGLSZM, dfsGLRLM, dfsNGTDM, dfsGLDM]\n    DFS = pd.concat(DFS, axis=1)\n    df_final.append(DFS)\n\n    print('img', i)","metadata":{"execution":{"iopub.status.busy":"2023-05-15T12:20:16.662854Z","iopub.execute_input":"2023-05-15T12:20:16.663297Z","iopub.status.idle":"2023-05-15T13:34:32.723455Z","shell.execute_reply.started":"2023-05-15T12:20:16.663263Z","shell.execute_reply":"2023-05-15T13:34:32.721471Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"len(df_final)","metadata":{"execution":{"iopub.status.busy":"2023-05-15T13:34:32.730144Z","iopub.execute_input":"2023-05-15T13:34:32.730597Z","iopub.status.idle":"2023-05-15T13:34:32.746485Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import pandas as pd\n\nsufijos = ['_T1', '_Flair', '_t1wce']\n\ndataframes_f = []\nfor i in range(len(df_final)):\n    # Concatenar todas las filas del dataframe en una sola fila\n    fila_concatenada = pd.concat([df_final[i].iloc[j] for j in range(len(df_final[i]))], axis=0)\n\n    # Obtener los nombres de las variables originales\n    nombres_originales = df_final[0].columns\n\n    # Crear los nuevos nombres de las columnas\n    nuevos_nombres = []\n    for nombre in nombres_originales:\n        nuevos_nombres.extend([nombre + sufijo for sufijo in sufijos])\n\n    # Asignar los nuevos nombres de las columnas\n    fila_concatenada.index = nuevos_nombres\n\n    # Crear un nuevo dataframe con una sola fila\n    nuevo_dataframe = pd.DataFrame(fila_concatenada)\n    \n    dataframes_f.append(nuevo_dataframe.T)\n","metadata":{"execution":{"iopub.status.busy":"2023-05-15T13:34:32.769079Z","iopub.execute_input":"2023-05-15T13:34:32.769664Z","iopub.status.idle":"2023-05-15T13:34:34.194685Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_final = pd.concat(dataframes_f, axis=0)\ndf_final\n","metadata":{"execution":{"iopub.status.busy":"2023-05-15T13:34:34.195942Z","iopub.execute_input":"2023-05-15T13:34:34.196270Z","iopub.status.idle":"2023-05-15T13:34:34.323685Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_final.to_csv('/kaggle/working/df_radiomica.csv', index=False)\n","metadata":{"execution":{"iopub.status.busy":"2023-05-15T13:40:54.121491Z","iopub.execute_input":"2023-05-15T13:40:54.122722Z","iopub.status.idle":"2023-05-15T13:40:54.397664Z","shell.execute_reply.started":"2023-05-15T13:40:54.122673Z","shell.execute_reply":"2023-05-15T13:40:54.396217Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"First Order Statistics (19 features) <br>\nShape-based (3D) (16 features)<br>\nGray Level Cooccurence Matrix (24 features)<br>\nGray Level Run Length Matrix (16 features)<br>\nGray Level Size Zone Matrix (16 features)<br>\nNeighbouring Gray Tone Difference Matrix (5 features)<br>\nGray Level Dependence Matrix (14 features)","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from __future__ import print_function\n\nimport logging\n\nimport numpy as np\nimport pywt\nimport SimpleITK as sitk\nimport six\nfrom six.moves import range\n","metadata":{"execution":{"iopub.status.busy":"2023-05-15T15:43:55.410071Z","iopub.execute_input":"2023-05-15T15:43:55.411070Z","iopub.status.idle":"2023-05-15T15:43:55.417989Z","shell.execute_reply.started":"2023-05-15T15:43:55.411017Z","shell.execute_reply":"2023-05-15T15:43:55.416598Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"logger = logging.getLogger(__name__)","metadata":{"execution":{"iopub.status.busy":"2023-05-15T15:43:57.302008Z","iopub.execute_input":"2023-05-15T15:43:57.303169Z","iopub.status.idle":"2023-05-15T15:43:57.308586Z","shell.execute_reply.started":"2023-05-15T15:43:57.303123Z","shell.execute_reply":"2023-05-15T15:43:57.307298Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def getWaveletImage(inputImage, inputMask, **kwargs):\n  \"\"\"\n  Applies wavelet filter to the input image and yields the decompositions and the approximation.\n\n  Following settings are possible:\n\n  - start_level [0]: integer, 0 based level of wavelet which should be used as first set of decompositions\n    from which a signature is calculated\n  - level [1]: integer, number of levels of wavelet decompositions from which a signature is calculated.\n  - wavelet [\"coif1\"]: string, type of wavelet decomposition. Enumerated value, validated against possible values\n    present in the ``pyWavelet.wavelist()``. Current possible values (pywavelet version 0.4.0) (where an\n    aditional number is needed, range of values is indicated in []):\n\n    - haar\n    - dmey\n    - sym[2-20]\n    - db[1-20]\n    - coif[1-5]\n    - bior[1.1, 1.3, 1.5, 2.2, 2.4, 2.6, 2.8, 3.1, 3.3, 3.5, 3.7, 3.9, 4.4, 5.5, 6.8]\n    - rbio[1.1, 1.3, 1.5, 2.2, 2.4, 2.6, 2.8, 3.1, 3.3, 3.5, 3.7, 3.9, 4.4, 5.5, 6.8]\n\n  Returned filter name reflects wavelet type:\n  wavelet[level]-<decompositionName>\n\n  N.B. only levels greater than the first level are entered into the name.\n\n  :return: Yields each wavelet decomposition and final approximation, corresponding imaget type name and ``kwargs``\n    (customized settings).\n  \"\"\"\n  global logger\n\n  logger.debug('Generating Wavelet images')\n\n  Nd = inputImage.GetDimension()\n  axes = list(range(Nd - 1, -1, -1))\n  if kwargs.get('force2D', False):\n    axes.remove(kwargs.get('force2Ddimension', 0))\n\n  approx, ret = _swt3(inputImage, tuple(axes), **kwargs)\n\n  for idx, wl in enumerate(ret, start=1):\n    for decompositionName, decompositionImage in wl.items():\n      logger.info('Computing Wavelet %s', decompositionName)\n\n      if idx == 1:\n        inputImageName = 'wavelet-%s' % (decompositionName)\n      else:\n        inputImageName = 'wavelet%s-%s' % (idx, decompositionName)\n      logger.debug('Yielding %s image', inputImageName)\n      return decompositionImage, inputImageName, kwargs\n\n  if len(ret) == 1:\n    inputImageName = 'wavelet-%s' % ('L' * len(axes))\n  else:\n    inputImageName = 'wavelet%s-%s' % (len(ret), ('L' * len(axes)))\n  logger.debug('Yielding approximation (%s) image', inputImageName)\n  return approx, inputImageName, kwargs","metadata":{"execution":{"iopub.status.busy":"2023-05-15T15:43:58.661118Z","iopub.execute_input":"2023-05-15T15:43:58.661529Z","iopub.status.idle":"2023-05-15T15:43:58.672799Z","shell.execute_reply.started":"2023-05-15T15:43:58.661496Z","shell.execute_reply":"2023-05-15T15:43:58.670912Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def _swt3(inputImage, axes, **kwargs):  # Stationary Wavelet Transform 3D\n  wavelet = kwargs.get('wavelet', 'coif1')\n  level = kwargs.get('level', 1)\n  start_level = kwargs.get('start_level', 0)\n\n  matrix = sitk.GetArrayFromImage(inputImage)  # This function gets a numpy array from the SimpleITK Image \"inputImage\"\n  matrix = np.asarray(matrix) # The function np.asarray converts \"matrix\" (which could be also a tuple) into an array.\n\n  original_shape = matrix.shape\n  # original_shape becomes a tuple (?,?,?) containing the number of rows, columns, and slices of the image\n  # this is of course dependent on the number of dimensions, but the same principle holds\n  padding = tuple([(0, 1 if dim % 2 != 0 else 0) for dim in original_shape])\n  # padding is necessary because of pywt.swtn (see function Notes)\n  data = matrix.copy()  # creates a modifiable copy of \"matrix\" and we call it \"data\"\n  data = np.pad(data, padding, 'wrap')  # padding the tuple \"padding\" previously computed\n\n  if not isinstance(wavelet, pywt.Wavelet):\n    wavelet = pywt.Wavelet(wavelet)\n\n  for i in range(0, start_level):  # if start_level = 0 (default) this for loop never gets executed\n    # compute all decompositions and saves them in \"dec\" dict\n    dec = pywt.swtn(data, wavelet, level=1, start_level=0, axes=axes)[0]\n    # copies in \"data\" just the \"aaa\" decomposition (i.e. approximation; No of consecutive 'a's = len(axes))\n    data = dec['a' * len(axes)].copy()\n\n  ret = []  # initialize empty list\n  for i in range(start_level, start_level + level):\n    # compute the n-dimensional stationary wavelet transform\n    dec = pywt.swtn(data, wavelet, level=1, start_level=0, axes=axes)[0]\n    # Copy the approximation into data (approximation in output / input for next levels)\n    data = dec['a' * len(axes)].copy()\n\n    dec_im = {}  # initialize empty dict\n    for decName, decImage in six.iteritems(dec):\n      # Returning the approximiation is done only for the last loop,\n      # and is handled separately below (by building it from `data`)\n      # There for, skip it here\n      if decName == 'a' * len(axes):\n        continue\n      decTemp = decImage.copy()\n      decTemp = decTemp[tuple(slice(None, -1 if dim % 2 != 0 else None) for dim in original_shape)]\n      sitkImage = sitk.GetImageFromArray(decTemp)\n      sitkImage.CopyInformation(inputImage)\n      dec_im[str(decName).replace('a', 'L').replace('d', 'H')] = sitkImage\n      # modifies 'a' with 'L' (Low-pass filter) and 'd' with 'H' (High-pass filter)\n\n    ret.append(dec_im)  # appending all the filtered sitk images (stored in \"dec_im\") to the \"ret\" list\n\n  data = data[tuple(slice(None, -1 if dim % 2 != 0 else None) for dim in original_shape)]\n  approximation = sitk.GetImageFromArray(data)\n  approximation.CopyInformation(inputImage)\n\n  return approximation, ret  # returns the approximation and the detail (ret) coefficients of the stationary wavelet decomposition\n","metadata":{"execution":{"iopub.status.busy":"2023-05-15T15:44:00.484593Z","iopub.execute_input":"2023-05-15T15:44:00.485015Z","iopub.status.idle":"2023-05-15T15:44:00.498611Z","shell.execute_reply.started":"2023-05-15T15:44:00.484973Z","shell.execute_reply":"2023-05-15T15:44:00.497277Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"imagenes = []\ndfs1st =  []\ndfs3d = []\ndfsGLCM = []\ndfsGLSZM = []\ndfsGLRLM = []\ndfsNGTDM = []\ndfsGLDM = []\n\n\nIMG_wavelets = []\ndf_final_wavelets = []\n\n\nfor i in range(len(dataset_dirs2)):\n# for i in range(2):\n    imagen = get_img(i)\n    test_img = np.array([np.moveaxis(np.array(imagen.resize((256, 256))), -1, 0)])\n#     test_res = model(torch.Tensor(test_img))\n\n#     image = sitk.GetImageFromArray(test_img[0,1,:,:].reshape(1, 256, 256))\n#     mask = sitk.GetImageFromArray(np.array([\n#             test_res[0][0].detach().cpu().numpy() > 0.5\n#         ]).astype(np.uint8))\n\n#     im = sitk.GetArrayFromImage(image)\n#     msk = sitk.GetArrayFromImage(mask)\n\n#     image0 = sitk.GetImageFromArray(test_img[0,1,:,:].reshape(1, 256, 256))\n\n#     mask = np.array([test_res[0][0].detach().cpu().numpy() > 0.5]).astype(np.uint8)\n#     #mask_clean = remove_small_blobs(mask, 100)\n#     mask_clean2 = keep_bigger_blob(mask)\n#     mask_cl = ndimage.binary_fill_holes(mask_clean2[0,:,:]).astype(int)\n#     mask_clean2[0,:,:] = mask_cl\n#     mask_clean2 = sitk.GetImageFromArray(mask_clean2)\n#     # Crop the image\n#     # bb is the bounding box, upon which the image and mask are cropped\n#     bb, correctedMask = imageoperations.checkMask(image0, mask_clean2, label=1)\n#     if correctedMask is not None:\n#         mask_clean2 = correctedMask\n#     image, mask = imageoperations.cropToTumorMask(image0, mask_clean2, bb)\n\n#     msk0 = mask\n#     msk = sitk.GetArrayFromImage(mask)\n\n    T1 = test_img[0,0,:,:].reshape(1, 256, 256)\n    Flair = test_img[0,1,:,:].reshape(1, 256, 256)\n    T1wce = test_img[0,2,:,:].reshape(1, 256, 256)\n    \n    seqs = [T1, Flair, T1wce]\n    \n    dfs1st = []\n    dfs3d = []\n    dfsGLCM = []\n    dfsGLSZM = []\n    dfsGLRLM = []\n    dfsNGTDM = []\n    dfsGLDM = []\n    wavelets = []\n\n    for j in range(3):\n    \n        test_res = model(torch.Tensor(test_img))\n\n        image = sitk.GetImageFromArray(seqs[j])\n        mask = sitk.GetImageFromArray(np.array([\n                test_res[0][0].detach().cpu().numpy() > 0.5\n            ]).astype(np.uint8))\n        \n        \n\n        im = sitk.GetArrayFromImage(image)\n        msk = sitk.GetArrayFromImage(mask)\n\n        image0 = sitk.GetImageFromArray(seqs[j])\n\n        mask = np.array([test_res[0][0].detach().cpu().numpy() > 0.5]).astype(np.uint8)\n        #mask_clean = remove_small_blobs(mask, 100)\n        mask_clean2 = keep_bigger_blob(mask)\n        mask_cl = ndimage.binary_fill_holes(mask_clean2[0,:,:]).astype(int)\n        mask_clean2[0,:,:] = mask_cl\n        mask_clean2 = sitk.GetImageFromArray(mask_clean2)\n        # Crop the image\n        # bb is the bounding box, upon which the image and mask are cropped\n        bb, correctedMask = imageoperations.checkMask(image0, mask_clean2, label=1)\n        if correctedMask is not None:\n            mask_clean2 = correctedMask\n        image, mask = imageoperations.cropToTumorMask(image0, mask_clean2, bb)\n\n        msk0 = mask\n        msk = sitk.GetArrayFromImage(mask)\n\n        wavelet = getWaveletImage(image, mask)\n        \n        ##### ------  CARACTERISTICAS ------------------\n    \n        \n        \n        # -------------- firstOrderFeatures\n\n        firstOrderFeatures = radiomics.firstorder.RadiomicsFirstOrder(wavelet[0], mask)\n        # Extraemos todas las características de la categoría First Order:\n        firstOrderFeatures.enableAllFeatures()\n        # Calcula las características y hace el print\n        result1st = firstOrderFeatures.execute()\n\n        data1st = pd.DataFrame(list(result1st.values()))\n        data1st = data1st.T  \n        data1st.columns = list(result1st.keys())\n\n        dfs1st.append(data1st)\n\n        # -------------- shape 3d\n\n        shape3d = radiomics.shape.RadiomicsShape(wavelet[0], mask)\n\n        shape3d.enableAllFeatures()\n        result3d = shape3d.execute()\n\n        data3d = pd.DataFrame(list(result3d.values()))\n        data3d = data3d.T  \n        data3d.columns = list(result3d.keys())\n\n        dfs3d.append(data3d)\n\n        # -------------- GLCM\n\n        GLCM = radiomics.glcm.RadiomicsGLCM(wavelet[0], mask)\n\n        GLCM.enableAllFeatures()\n        resultGLCM = GLCM.execute()\n\n        dataGLCM = pd.DataFrame(list(resultGLCM.values()))\n        dataGLCM = dataGLCM.T  \n        dataGLCM.columns = list(resultGLCM.keys())\n\n        dfsGLCM.append(dataGLCM)\n\n        # -------------- GLSZM\n        GLSZM = radiomics.glszm.RadiomicsGLSZM(wavelet[0], mask)\n\n        GLSZM.enableAllFeatures()\n        resultGLSZM = GLSZM.execute()\n\n        dataGLSZM = pd.DataFrame(list(resultGLSZM.values()))\n        dataGLSZM = dataGLSZM.T  \n        dataGLSZM.columns = list(resultGLSZM.keys())\n\n        dfsGLSZM.append(dataGLSZM)\n\n        # -------------- GLRLM\n\n        GLRLM = radiomics.glrlm.RadiomicsGLRLM(wavelet[0], mask)\n\n        GLRLM.enableAllFeatures()\n        resultGLRLM = GLRLM.execute()\n\n        dataGLRLM = pd.DataFrame(list(resultGLRLM.values()))\n        dataGLRLM = dataGLRLM.T  \n        dataGLRLM.columns = list(resultGLRLM.keys())\n\n        dfsGLRLM.append(dataGLRLM)\n\n        # -------------- NGTDM\n\n        NGTDM = radiomics.ngtdm.RadiomicsNGTDM(wavelet[0], mask)\n\n        NGTDM.enableAllFeatures()\n        resultNGTDM = NGTDM.execute()\n\n        dataNGTDM = pd.DataFrame(list(resultNGTDM.values()))\n        dataNGTDM = dataNGTDM.T  \n        dataNGTDM.columns = list(resultNGTDM.keys())\n\n        dfsNGTDM.append(dataNGTDM)\n\n        # -------------- GLDM\n\n\n        GLDM = radiomics.gldm.RadiomicsGLDM(wavelet[0], mask)\n\n        GLDM.enableAllFeatures()\n        resultGLDM = GLDM.execute()\n\n        dataGLDM = pd.DataFrame(list(resultGLDM.values()))\n        dataGLDM = dataGLDM.T  \n        dataGLDM.columns = list(resultGLDM.keys())\n\n        dfsGLDM.append(dataGLDM)\n\n\n        print('seq',j)\n    \n    dfs1st = pd.concat(dfs1st)\n    dfs3d = pd.concat(dfs3d)\n    dfsGLCM = pd.concat(dfsGLCM)\n    dfsGLSZM = pd.concat(dfsGLSZM)\n    dfsGLRLM = pd.concat(dfsGLRLM)\n    dfsNGTDM = pd.concat(dfsNGTDM)\n    dfsGLDM = pd.concat(dfsGLDM)\n    DFS = [dfs1st, dfs3d, dfsGLCM, dfsGLSZM, dfsGLRLM, dfsNGTDM, dfsGLDM]\n    DFS = pd.concat(DFS, axis=1)\n    df_final_wavelets.append(DFS)\n    IMG_wavelets.append(wavelets)\n\n    print('img', i)","metadata":{"_kg_hide-output":false,"execution":{"iopub.status.busy":"2023-05-15T16:10:29.652022Z","iopub.execute_input":"2023-05-15T16:10:29.652595Z","iopub.status.idle":"2023-05-15T17:14:14.857691Z","shell.execute_reply.started":"2023-05-15T16:10:29.652553Z","shell.execute_reply":"2023-05-15T17:14:14.849596Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import pandas as pd\n\nsufijos = ['_T1_wavelet', '_Flair_wavelet', '_t1wce_wavelet']\n\ndataframes_f = []\nfor i in range(len(df_final_wavelets)):\n    # Concatenar todas las filas del dataframe en una sola fila\n    fila_concatenada = pd.concat([df_final_wavelets[i].iloc[j] for j in range(len(df_final_wavelets[i]))], axis=0)\n\n    # Obtener los nombres de las variables originales\n    nombres_originales = df_final_wavelets[0].columns\n\n    # Crear los nuevos nombres de las columnas\n    nuevos_nombres = []\n    for nombre in nombres_originales:\n        nuevos_nombres.extend([nombre + sufijo for sufijo in sufijos])\n\n    # Asignar los nuevos nombres de las columnas\n    fila_concatenada.index = nuevos_nombres\n\n    # Crear un nuevo dataframe con una sola fila\n    nuevo_dataframe = pd.DataFrame(fila_concatenada)\n    \n    dataframes_f.append(nuevo_dataframe.T)","metadata":{"execution":{"iopub.status.busy":"2023-05-15T17:14:14.869205Z","iopub.execute_input":"2023-05-15T17:14:14.869784Z","iopub.status.idle":"2023-05-15T17:14:16.266121Z","shell.execute_reply.started":"2023-05-15T17:14:14.869682Z","shell.execute_reply":"2023-05-15T17:14:16.264934Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_final_wavelets = pd.concat(dataframes_f, axis=0)\ndf_final_wavelets","metadata":{"execution":{"iopub.status.busy":"2023-05-15T17:14:16.268815Z","iopub.execute_input":"2023-05-15T17:14:16.269598Z","iopub.status.idle":"2023-05-15T17:14:16.409170Z","shell.execute_reply.started":"2023-05-15T17:14:16.269548Z","shell.execute_reply":"2023-05-15T17:14:16.408220Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_final_wavelets.to_csv('/kaggle/working/df_wavelets_radiomica.csv', index=False)\n","metadata":{"execution":{"iopub.status.busy":"2023-05-15T17:14:16.411917Z","iopub.execute_input":"2023-05-15T17:14:16.412298Z","iopub.status.idle":"2023-05-15T17:14:16.717167Z","shell.execute_reply.started":"2023-05-15T17:14:16.412265Z","shell.execute_reply":"2023-05-15T17:14:16.715523Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df1 = pd.read_csv('/kaggle/working/df_radiomica.csv')\ndf2 = pd.read_csv('/kaggle/working/df_wavelets_radiomica.csv')\n\nDF_FINAL_RADIOMICA = pd.concat([df1, df2], axis=1)\n\nDF_FINAL_RADIOMICA.to_csv('/kaggle/working/DF_FINAL_RADIOMICA.csv', index=False)\nDF_FINAL_RADIOMICA","metadata":{"execution":{"iopub.status.busy":"2023-05-15T17:23:01.977144Z","iopub.execute_input":"2023-05-15T17:23:01.978877Z","iopub.status.idle":"2023-05-15T17:23:02.600709Z","shell.execute_reply.started":"2023-05-15T17:23:01.978820Z","shell.execute_reply":"2023-05-15T17:23:02.599374Z"},"trusted":true},"execution_count":null,"outputs":[]}]}