{"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":"# Install MONAI package \nThe MONAI framework is the open-source foundation being created by Project MONAI. MONAI is a freely available, community-supported, PyTorch-based framework for deep learning in healthcare imaging. \nYou can learn more in [MONAI's website](https://monai.io/index.html).","metadata":{}},{"cell_type":"code","source":"!pip install monai","metadata":{"execution":{"iopub.status.busy":"2022-04-24T00:38:58.557247Z","iopub.execute_input":"2022-04-24T00:38:58.557612Z","iopub.status.idle":"2022-04-24T00:39:12.591297Z","shell.execute_reply.started":"2022-04-24T00:38:58.557522Z","shell.execute_reply":"2022-04-24T00:39:12.590069Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Import needed packages","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport torch\nimport matplotlib.pyplot as plt\nimport os\nfrom sklearn.model_selection import StratifiedKFold\nfrom glob import glob\nimport pydicom as dcm\n\nfrom monai.transforms import(\n    \n    Compose,\n    AddChanneld,\n    Resized,\n    ToTensord,\n    Spacingd,\n    Orientationd,\n    ScaleIntensityRanged,\n    CropForegroundd,\n)\nfrom monai.data import Dataset, DataLoader\nfrom monai.utils import set_determinism\n\nimport pandas as pd\n","metadata":{"execution":{"iopub.status.busy":"2022-04-24T00:39:12.594648Z","iopub.execute_input":"2022-04-24T00:39:12.594922Z","iopub.status.idle":"2022-04-24T00:39:22.158679Z","shell.execute_reply.started":"2022-04-24T00:39:12.594889Z","shell.execute_reply":"2022-04-24T00:39:22.157420Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Let's take a look at a .dcm file\nDICOM (Digital Imaging and Communications in Medicine) is a standard protocol for the management and transmission of medical images and related data and is used in many healthcare facilities.\nFiles in DICOM format have the .dcm suffix. \nFor CT Scans, each .dcm file is a single slice from the exam.\nIn the next cell, we will show the structure of a .dcm file in this dataset.","metadata":{}},{"cell_type":"code","source":"\"\"\"\nIn this cell, we define the input_dir, where all files are in, and the imaging_dir,\nwhere we have all our dicom files.\n\"\"\"\n\ninput_dir = '/kaggle/input/unifesp-fatty-liver'\nimaging_dir = os.path.join(input_dir, 'fatty-liver-dataset', 'd_2')","metadata":{"execution":{"iopub.status.busy":"2022-04-24T00:39:22.160057Z","iopub.execute_input":"2022-04-24T00:39:22.160672Z","iopub.status.idle":"2022-04-24T00:39:22.168438Z","shell.execute_reply.started":"2022-04-24T00:39:22.160635Z","shell.execute_reply":"2022-04-24T00:39:22.165514Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\"\"\"\nHere, we open a dicom file by calling the dcmread function inside the pydicom library.\n\"\"\"\n\ni = 2022  # change this to open a different dicom file\nfile = sorted(glob(imaging_dir + '/*' + '/*' + '/*' + '/*'))[i]\ndicom = dcm.dcmread(file)\nprint(dicom)","metadata":{"execution":{"iopub.status.busy":"2022-04-24T00:39:22.171382Z","iopub.execute_input":"2022-04-24T00:39:22.171799Z","iopub.status.idle":"2022-04-24T00:40:00.048718Z","shell.execute_reply.started":"2022-04-24T00:39:22.171751Z","shell.execute_reply":"2022-04-24T00:40:00.047440Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\"\"\"\nA dicom file has many tags. We can acess a specific tag, for example, the tag (0010,0020), which \ncontains the patient ID, can be acessed like this:\ndicom[0x0010, 0x0020]\nWe can also acess the tag value and its name.\n\n\"\"\"\nprint(dicom[0x0010,0x0020])\nprint(dicom[0x0010,0x0020].value)\nprint(dicom[0x0010,0x0020].name)","metadata":{"execution":{"iopub.status.busy":"2022-04-24T00:40:00.050274Z","iopub.execute_input":"2022-04-24T00:40:00.051023Z","iopub.status.idle":"2022-04-24T00:40:00.059205Z","shell.execute_reply.started":"2022-04-24T00:40:00.050969Z","shell.execute_reply":"2022-04-24T00:40:00.057858Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\"\"\"\nLast, but not least, the image pixels can also be acessed in our dicom file.\n\"\"\"\nprint(dicom.pixel_array)\nprint(dicom.pixel_array.shape)","metadata":{"execution":{"iopub.status.busy":"2022-04-24T00:40:00.061404Z","iopub.execute_input":"2022-04-24T00:40:00.062824Z","iopub.status.idle":"2022-04-24T00:40:00.075211Z","shell.execute_reply.started":"2022-04-24T00:40:00.062763Z","shell.execute_reply":"2022-04-24T00:40:00.074522Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Now, we are going to split the patients into training, cross-validation and testing sets\nIt is very important that the split is done at the patient level (the first subdirectory inside 'd_2'), because we don't want the same patient to have studies in different sets.","metadata":{}},{"cell_type":"code","source":"\"\"\"\nFirst, we split all patients between Chest and Abdomen. \nWe get this information by iterating through the exam directories, opening a dicom file \nand reading the Study Description Tag. \n\"\"\"\n\nlist_chest = []\nlist_abdomen = []\nfor file in (sorted(glob(imaging_dir + '/*' + '/*' + '/*'))):  \n    dicom = dcm.dcmread(os.path.join(file, os.listdir(file)[0]), stop_before_pixels=True)\n    \n    # Here we iterate through each series directory (the last subdirectory) and read the first dicom file\n    # (doesn't matter which one you read). We stop before the pixel array, so the process takes less time.\n    if dicom[0x0008,0x1030].value == 'Chest':\n        list_chest.append(file.split('/')[6])\n    elif dicom[0x0008,0x1030].value == 'Abdomen':\n        list_abdomen.append(file.split('/')[6])\n\n","metadata":{"execution":{"iopub.status.busy":"2022-04-24T00:40:00.076949Z","iopub.execute_input":"2022-04-24T00:40:00.077627Z","iopub.status.idle":"2022-04-24T00:40:05.351289Z","shell.execute_reply.started":"2022-04-24T00:40:00.077593Z","shell.execute_reply":"2022-04-24T00:40:05.350298Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\"\"\"\nGet the unique values from Chest and Abdomen lists and convert to integers\n\"\"\"\nlist_chest = np.unique(list_chest).astype(int)\nlist_abdomen = np.unique(list_abdomen).astype(int)","metadata":{"execution":{"iopub.status.busy":"2022-04-24T00:40:05.352752Z","iopub.execute_input":"2022-04-24T00:40:05.353641Z","iopub.status.idle":"2022-04-24T00:40:05.359947Z","shell.execute_reply.started":"2022-04-24T00:40:05.353594Z","shell.execute_reply":"2022-04-24T00:40:05.358986Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\"\"\"\nFind the index of the patients who have both Chest and Abdomen CTs, and remove from Chest List.\nThis is important because we are gonna use the Chest List to create the validation set, which should\nonly contain Chest Scans, because our final goal is to create a model to perform well on these studies.\n\"\"\"\n\n_, idx, __ = np.intersect1d(list_chest, list_abdomen, assume_unique=True, return_indices=True)\nlist_chest = np.delete(list_chest, idx)\n","metadata":{"execution":{"iopub.status.busy":"2022-04-24T00:40:05.361568Z","iopub.execute_input":"2022-04-24T00:40:05.362682Z","iopub.status.idle":"2022-04-24T00:40:05.382040Z","shell.execute_reply.started":"2022-04-24T00:40:05.362640Z","shell.execute_reply":"2022-04-24T00:40:05.380790Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now that we have our chest and abdomen lists, we remove from the chest group the patients from the testing set, which are already established by the competition.","metadata":{}},{"cell_type":"code","source":"# First we list the patients from the testing set by reading the submission file.\ntest_split = pd.read_csv(os.path.join(input_dir,'sample_submission.csv'))['Id'].to_numpy()\n# Then, we certify that all these patients are chest patients.\ncommon, idx, _ = np.intersect1d(list_chest, test_split, assume_unique=True, return_indices=True)\n\nif all(common == test_split):\n    print('All patients from the test set are chest patients.')\nelse:\n    print('There are patients from the test set which are not chest patients.')\n# Finally, we remove the testing patients.\nlist_chest = np.delete(list_chest, idx)","metadata":{"execution":{"iopub.status.busy":"2022-04-24T00:40:05.386382Z","iopub.execute_input":"2022-04-24T00:40:05.386749Z","iopub.status.idle":"2022-04-24T00:40:05.415142Z","shell.execute_reply.started":"2022-04-24T00:40:05.386711Z","shell.execute_reply":"2022-04-24T00:40:05.413804Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now, we are going to divide the chest group into n cross validation groups, so we can train multiple models and ensemble them. ","metadata":{}},{"cell_type":"code","source":"n_splits = 3  # change here the number of cross validation groups\nn_abd = len(list_abdomen)\nn_chest = len(list_chest)\nn_total = n_abd + n_chest\nn_val = int(n_chest/n_splits)\npercent_test = 100*n_val/n_total\n\nprint('CT Abdomen:' , n_abd)\nprint('CT Chest:' , n_chest)\nprint('Total:', n_total)\nprint('Validation Split:', n_val, 'Percentage of total:', f'{percent_test:.1f}')","metadata":{"execution":{"iopub.status.busy":"2022-04-24T00:40:05.417126Z","iopub.execute_input":"2022-04-24T00:40:05.417495Z","iopub.status.idle":"2022-04-24T00:40:05.428680Z","shell.execute_reply.started":"2022-04-24T00:40:05.417428Z","shell.execute_reply":"2022-04-24T00:40:05.427392Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Based on the above split sizes, we are gonna split the chest group into three cross validation sets, in a stratified manner, so that we don't get an unbalenced group. Here, we choose this number of splits so that we get a validation set size similar to the testing set size. Feel free to play around with this number.","metadata":{}},{"cell_type":"code","source":"labels = pd.read_csv(os.path.join(input_dir, 'train.csv'))  # load the traning csv\nlabels","metadata":{"execution":{"iopub.status.busy":"2022-04-24T00:40:05.430618Z","iopub.execute_input":"2022-04-24T00:40:05.431426Z","iopub.status.idle":"2022-04-24T00:40:05.471534Z","shell.execute_reply.started":"2022-04-24T00:40:05.431370Z","shell.execute_reply":"2022-04-24T00:40:05.470335Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"labels.index = labels['Id']  # change the index so we can use the loc function to retrieve the Ids\nchest_labels = np.array(labels.loc[list_chest,'ground_truth']) # create an array with only our chest group\nchest_labels.shape","metadata":{"execution":{"iopub.status.busy":"2022-04-24T00:40:05.474050Z","iopub.execute_input":"2022-04-24T00:40:05.474418Z","iopub.status.idle":"2022-04-24T00:40:05.488701Z","shell.execute_reply.started":"2022-04-24T00:40:05.474373Z","shell.execute_reply":"2022-04-24T00:40:05.487147Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now, let's create the cross validation and training sets and save them in txt files.","metadata":{}},{"cell_type":"code","source":"splits_dir = '/kaggle/working/splits'\nkf = StratifiedKFold(n_splits=3, shuffle=True, random_state=2022)\ni=0\ntry:\n    os.mkdir(splits_dir)\nexcept FileExistsError:\n    print('\"splits\" directory already exists.')\nfor train_idx, val_idx in kf.split(list_chest, chest_labels):\n    \n        \n    np.savetxt(os.path.join(splits_dir, f'train_split_{i}.txt'), \n               np.concatenate((list_chest[train_idx], list_abdomen)), delimiter=',', fmt='%d')\n    np.savetxt(os.path.join(splits_dir,f'val_split_{i}.txt'), list_chest[val_idx], delimiter=',',\n               fmt='%d')\n    i+=1\n\n\nnp.savetxt(os.path.join(splits_dir,f'test_split.txt'), test_split, delimiter=',', fmt='%d')","metadata":{"execution":{"iopub.status.busy":"2022-04-24T00:40:05.490412Z","iopub.execute_input":"2022-04-24T00:40:05.490798Z","iopub.status.idle":"2022-04-24T00:40:05.512637Z","shell.execute_reply.started":"2022-04-24T00:40:05.490751Z","shell.execute_reply":"2022-04-24T00:40:05.511261Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Finally, let's create a function to load train and validation groups from .txt created above, and a function to load the test group.","metadata":{}},{"cell_type":"code","source":"def load_train_val(split, split_dir):\n    \"\"\"\n    This function loads a train and validation split as NumPy arrays, from the txt files we created before.\n    \n    'split': this parameter is the split number from the k-fold cross validation. It should be an integer\n    between 0 (inclusive) and the k value (not inclusive) you chose.\n    'split_dir': the directory where you saved your splits.\n    \"\"\"\n    n_splits = 3  # change this to your k-fold number\n    if split in range(n_splits):\n        train = np.loadtxt(os.path.join(split_dir, f'train_split_{split}.txt'), dtype=int, delimiter=',')\n        val = np.loadtxt(os.path.join(split_dir, f'val_split_{split}.txt'), dtype=int, delimiter=',')\n        return train, val\n    else:\n        raise ValueError('Please specify a valid split number.')\n\n        \ndef load_test(split_dir):\n    \"\"\"\n    This function loads the testing set into a NumPy array.\n    'split_dir': the directory where you saved your splits.\n    \"\"\"\n    return np.loadtxt(os.path.join(split_dir, 'test_split.txt'),dtype=int, delimiter=',')","metadata":{"execution":{"iopub.status.busy":"2022-04-24T00:40:05.514742Z","iopub.execute_input":"2022-04-24T00:40:05.515020Z","iopub.status.idle":"2022-04-24T00:40:05.523782Z","shell.execute_reply.started":"2022-04-24T00:40:05.514990Z","shell.execute_reply":"2022-04-24T00:40:05.522642Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# After we created our splits, let's load our images\n","metadata":{}},{"cell_type":"markdown","source":"We will create some functions to make our life easier. \nFirst, let's create a custom Dataset class to deal with Dicom images.","metadata":{}},{"cell_type":"code","source":"class DicomDataset(Dataset):\n    def __init__(self, data, transform):\n        self.data = data\n        self.transform = transform\n\n    def __len__(self):\n        return len(self.data)\n\n    def __getitem__(self, index):\n        item = {'img': None, 'label': None}\n        item['img'] = dcm.dcmread(self.data[index]['img'])\n        slope = float(item['img'][0x0028, 0x1053].value)\n        intercept = float(item['img'][0x0028, 0x1052].value)\n\n        item['img'] = np.array(item['img'].pixel_array, dtype=np.float32)\n        item['img'] *= slope\n        item['img'] += intercept\n        item['label'] = float(self.data[index]['label'])\n        if self.transform:\n            item = self.transform(item)\n\n        return item","metadata":{"execution":{"iopub.status.busy":"2022-04-24T00:40:05.525413Z","iopub.execute_input":"2022-04-24T00:40:05.525717Z","iopub.status.idle":"2022-04-24T00:40:05.542146Z","shell.execute_reply.started":"2022-04-24T00:40:05.525682Z","shell.execute_reply":"2022-04-24T00:40:05.541342Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now, let's create a function that takes our split and creates a traning and a validation loader.","metadata":{}},{"cell_type":"code","source":"def prepare_train(train, val, imaging_dir, labels_dir, batch_size, pixdim=(1.5, 1.5), a_min=-200, \n                  a_max=200, spatial_size=[128,128]):\n\n    \"\"\"\n    This function creates a training and a validation loader, with all the transforms you want for\n    preprocessing. Here, we include some basic transforms, but you can add more operations that you \n    find in the Monai documentation.\n    https://monai.io/docs.html\n    \"\"\"\n    \n                    \n    set_determinism(seed=0)\n    path_images = sorted(glob(imaging_dir + '/*' + '/*' + '/*' + '/*.dcm'))\n\n    \n    path_train_images = [image for image in path_images if int(image.split('/')[-4]) in train]\n    path_val_images = [image for image in path_images if int(image.split('/')[-4]) in val]    \n    \n    df = pd.read_csv(os.path.join(labels_dir, 'train.csv'), index_col=0)\n    \n    \n    \n    train_files = [{\"img\": image_name, \"label\": df.loc[label_name, 'ground_truth']} for image_name,\n                 label_name in zip(path_train_images, train)]\n    val_files = [{\"img\": image_name, \"label\": df.loc[label_name, 'ground_truth']} for image_name,\n                 label_name in zip(path_val_images, val)]\n    \n    train_transforms = Compose(\n        [\n            AddChanneld(keys=[\"img\"]),\n            Spacingd(keys=[\"img\"], pixdim=pixdim, mode=(\"bilinear\")),\n            ScaleIntensityRanged(keys=[\"img\"], a_min=a_min, a_max=a_max, b_min=0.0, b_max=1.0, clip=True), \n            Resized(keys=[\"img\"], spatial_size=spatial_size),\n            Orientationd(keys=['img'], axcodes='RA'),\n            ToTensord(keys=[\"img\", 'label']),\n\n        ]\n    )\n\n    val_transforms = Compose(\n        [\n            AddChanneld(keys=[\"img\"]),\n            Spacingd(keys=[\"img\"], pixdim=pixdim, mode=(\"bilinear\")),\n            ScaleIntensityRanged(keys=[\"img\"], a_min=a_min, a_max=a_max, b_min=0.0, b_max=1.0, clip=True), \n            Resized(keys=[\"img\"], spatial_size=spatial_size),\n            Orientationd(keys=['img'], axcodes='RA'),\n            ToTensord(keys=[\"img\", 'label']),\n\n            \n        ]\n    )\n\n    \n\n    train_ds = DicomDataset(data=train_files, transform=train_transforms)\n    train_loader = DataLoader(train_ds, batch_size=batch_size, shuffle=True)\n\n    val_ds = DicomDataset(data=val_files, transform=val_transforms)\n    val_loader = DataLoader(val_ds, batch_size=batch_size, shuffle=True)\n\n    return train_loader, val_loader\n","metadata":{"execution":{"iopub.status.busy":"2022-04-24T00:40:05.544161Z","iopub.execute_input":"2022-04-24T00:40:05.544405Z","iopub.status.idle":"2022-04-24T00:40:05.563986Z","shell.execute_reply.started":"2022-04-24T00:40:05.544376Z","shell.execute_reply":"2022-04-24T00:40:05.562548Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Our last step here, let's create a function to show one batch from a DataLoader.","metadata":{}},{"cell_type":"code","source":"def show_batch(dl):\n    \"\"\"\n    This function is to show one batch from your datasets, so that you can si if the it is okay or you need \n    to change/delete something.\n\n    `dl`: this parameter should take the data loader, which means you need to prepare first and apply \n    the transforms that you want. After that, pass it to this function so that you visualize the patients\n    with the transforms that you want.\n    \"\"\"\n\n    for batch_data in dl:\n        batch_size = len(batch_data['label'])\n        fig, ax = plt.subplots(batch_size//4,4,figsize=(16,batch_size))\n        \n        for i in range(len(batch_data['label'])):\n            label = batch_data['label'][i].item()\n            \n            ax[i//4,i%4].imshow(batch_data['img'][i,0], cmap='gray', interpolation='bilinear', aspect='auto')\n            ax[i//4,i%4].set_title(f'label: {label}')\n            ax[i//4,i%4].xaxis.set_ticklabels([])\n            ax[i//4,i%4].yaxis.set_ticklabels([])\n        plt.show()\n        break\n","metadata":{"execution":{"iopub.status.busy":"2022-04-24T00:40:05.565675Z","iopub.execute_input":"2022-04-24T00:40:05.565930Z","iopub.status.idle":"2022-04-24T00:40:05.579902Z","shell.execute_reply.started":"2022-04-24T00:40:05.565896Z","shell.execute_reply":"2022-04-24T00:40:05.578860Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Alright, now that everything is ready, let's print a batch of slices from the training set!","metadata":{}},{"cell_type":"code","source":"labels_dir = '/kaggle/input/unifesp-fatty-liver'\ntrain, val = load_train_val(0, '/kaggle/working/splits')\ntrain_loader, val_loader = prepare_train(train, val, imaging_dir, labels_dir, 16, a_min=-200, a_max=200)\nshow_batch(train_loader)\n","metadata":{"execution":{"iopub.status.busy":"2022-04-24T00:40:05.581300Z","iopub.execute_input":"2022-04-24T00:40:05.581855Z","iopub.status.idle":"2022-04-24T00:40:10.246499Z","shell.execute_reply.started":"2022-04-24T00:40:05.581819Z","shell.execute_reply":"2022-04-24T00:40:10.245383Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Now it's your turn! Try to create a Machine Learning algorithm that can predict liver steatosis!","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}