{"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":"
\n

\n HuBMAP + HPA - Hacking the Human Body\n

\n \n
\n\n
\n\n
\n\n
","metadata":{}},{"cell_type":"markdown","source":"---\n\n

\n\n

\n\n

TABLE OF CONTENTS

\n\n---\n\n

0    IMPORTS & GENERAL SETTINGS

\n\n---\n\n

1    BACKGROUND INFORMATION

\n\n---\n\n\n

2    HELPER FUNCTIONS

\n\n---\n\n

3    DATASET CREATION AND EXPLORATION

\n\n---\n\n

4    MODELLING

\n\n---","metadata":{}},{"cell_type":"markdown","source":"
\n\n\n\n

0  IMPORTS & GENERAL SETTINGS     ⤒

","metadata":{}},{"cell_type":"code","source":"import matplotlib.pyplot as plt\nimport pandas as pd\nimport numpy as np\nfrom tqdm import tqdm\nimport os\n\nfrom PIL import Image\n\nimport torch\nfrom torch.utils.data import Dataset, DataLoader\nfrom torch.utils.tensorboard import SummaryWriter\n\nimport albumentations as A\nimport cv2\n\n!pip install monai\nimport monai\nfrom monai.inferers import sliding_window_inference\nfrom monai.metrics import DiceMetric\nfrom monai.visualize import plot_2d_or_3d_image\nfrom monai.data import decollate_batch\n\nfrom monai.transforms import (\n Activations,\n AddChannel,\n AsDiscrete,\n Compose,\n LoadImage,\n RandRotate90,\n RandSpatialCrop,\n ScaleIntensity,\n EnsureType,\n)","metadata":{"execution":{"iopub.status.busy":"2022-06-27T13:35:06.502518Z","iopub.execute_input":"2022-06-27T13:35:06.502872Z","iopub.status.idle":"2022-06-27T13:35:32.484425Z","shell.execute_reply.started":"2022-06-27T13:35:06.502842Z","shell.execute_reply":"2022-06-27T13:35:32.48337Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Global Settings\n\nINPUT_PATH = \"../input/hubmap-organ-segmentation\"\n\n# IMG Settings\n\nIMG_SHAPE = (256, 256)\nIMG_MEAN = [0.485, 0.456, 0.406]\nIMG_STD = [0.229, 0.224, 0.225]\n\n# Device\ndevice = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")","metadata":{"execution":{"iopub.status.busy":"2022-06-27T14:53:42.679611Z","iopub.execute_input":"2022-06-27T14:53:42.680128Z","iopub.status.idle":"2022-06-27T14:53:42.689236Z","shell.execute_reply.started":"2022-06-27T14:53:42.680083Z","shell.execute_reply":"2022-06-27T14:53:42.688057Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n\n\n\n

1  BACKGROUND INFORMATION    ⤒

\n\n---\n\n

1.1 BASIC COMPETITION INFORMATION

\n\n---\n\n
PRIMARY TASK DESCRIPTION\n\nThis is a placeholder copy/paste from the OVERVIEW>DESCRIPTION page of the Kaggle competition\n\nThe Human BioMolecular Atlas Program (HuBMAP) is working to create a Human Reference Atlas at the cellular level. Sponsored by the National Institutes of Health (NIH), HuBMAP and Indiana University’s Cyberinfrastructure for Network Science Center (CNS) have partnered with institutions across the globe for this endeavor. A major partner is the Human Protein Atlas (HPA), a Swedish research program aiming to map the protein expression in human cells, tissues, and organs, funded by the Knut and Alice Wallenberg Foundation.\n\nIn this competition, you’ll identify and segment functional tissue units (FTUs) across five human organs. You'll build your model using a dataset of tissue section images, with the best submissions segmenting FTUs as accurately as possible.\n\nIf successful, you'll help accelerate the world’s understanding of the relationships between cell and tissue organization. With a better idea of the relationship of cells, researchers will have more insight into the function of cells that impact human health. Further, the Human Reference Atlas constructed by HuBMAP will be freely available for use by researchers and pharmaceutical companies alike, potentially improving and prolonging human life.\n\n

1.2 COMPETITION EVALUATION

\n\n---\n\n
GENERAL EVALUATION INFORMATION\n\nThis competition is evaluated on the mean Dice coefficient. The Dice coefficient can be used to compare the pixel-wise agreement between a predicted segmentation and its corresponding ground truth. The formula is given by:\n\n
$Dice = \\large\\frac{2 * |X \\cap Y|}{|X| + |Y|}$

\n\nwhere X is the predicted set of pixels and Y is the ground truth. The Dice coefficient is defined to be 1 when both X and Y are empty. The leaderboard score is the mean of the Dice coefficients for each image in the test set.\n\n
SUBMISSION FILE INFORMATION\n\nIn order to reduce the submission file size, our metric uses run-length encoding on the pixel values. Instead of submitting an exhaustive list of indices for your segmentation, you will submit pairs of values that contain a start position and a run length. E.g. '1 3' implies starting at pixel 1 and running a total of 3 pixels (1,2,3).\n\nNote that, at the time of encoding, the mask should be binary, meaning the masks for all objects in an image are joined into a single large mask. A value of 0 should indicate pixels that are not masked, and a value of 1 will indicate pixels that are masked.\n\nThe competition format requires a space delimited list of pairs. For example, '1 3 10 5' implies pixels 1,2,3,10,11,12,13,14 are to be included in the mask. The metric checks that the pairs are sorted, positive, and the decoded pixel values are not duplicated. The pixels are numbered from top to bottom, then left to right: 1 is pixel (1,1), 2 is pixel (2,1), etc.\nThe file should contain a header and have the following format:\n\n```\nid,rle\n2,1 3\n5,3 4\n6,4 5 66 77\netc...\n```\n\n
\n\n

1.3 DATASET OVERVIEW

\n\n---\n\n
GENERAL INFORMATION\n\nIn this competition we are segmenting neuronal cells in images. The training annotations are provided as run length encoded masks, and the images are in PNG format. The number of images is small, but the number of annotated objects is quite high. The hidden test set is roughly 240 images.\n\n
FILES\n\n[train/test].csv Metadata for the train/test set. Only the first few rows of the test set are available for download.\n\n* `id` - The image ID.\n\n* `organ` - The organ that the biopsy sample was taken from.\n\n* `data_source` - Whether the image was provided by HuBMAP or HPA.\n\n* `img_height` - The height of the image in pixels.\n\n* `img_width` - The width of the image in pixels.\n \n* `pixel_size` - The height/width of a single pixel from this image in micrometers. All HPA images have a pixel size of 0.4 µm. For HuBMAP imagery the pixel size is 0.5 µm for kidney, 0.2290 µm for large intestine, 0.7562 µm for lung, 0.4945 µm for spleen, and 6.263 µm for prostate.\n \n* `tissue_thickness` - The thickness of the biopsy sample in micrometers. All HPA images have a thickness of 4 µm. The HuBMAP samples have tissue slice thicknesses 10 µm for kidney, 8 µm for large intestine, 4 µm for spleen, 5 µm for lung, and 5 µm for prostate.\n\n* `rle` - The target column. A run length encoded copy of the annotations. Provided for the training set only.\n\n* `age` - The patient's age in years. Provided for the training set only.\n\n* `sex` - The sex of the patient. Provided for the training set only.\nsample_submission.csv\n\n* `id` - The image ID.\n\n* `rle` - A run length encoded mask of the FTUs in the image.","metadata":{}},{"cell_type":"markdown","source":"
\n\n\n\n\n\n

\n 2  HELPER FUNCTION    ⤒\n

\n\n---","metadata":{}},{"cell_type":"code","source":"def rle_decode(mask_rle, shape, color=1):\n '''\n mask_rle: run-length as string formated (start length)\n shape: (height,width) of array to return \n Returns numpy array, 1 - mask, 0 - background\n '''\n s = mask_rle.split()\n starts, lengths = [np.asarray(x, dtype=int) for x in (s[0:][::2], s[1:][::2])]\n starts -= 1\n ends = starts + lengths\n img = np.zeros(shape[0] * shape[1], dtype=np.float32)\n for lo, hi in zip(starts, ends):\n img[lo : hi] = color\n return img.reshape(shape).T","metadata":{"execution":{"iopub.status.busy":"2022-06-27T13:35:32.562808Z","iopub.execute_input":"2022-06-27T13:35:32.563233Z","iopub.status.idle":"2022-06-27T13:35:32.572558Z","shell.execute_reply.started":"2022-06-27T13:35:32.563194Z","shell.execute_reply":"2022-06-27T13:35:32.571477Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n\n\n\n\n\n

\n 3  DATASET CREATION AND EXPLORATION    ⤒\n

\n\n---\n\n

3.0 Torch Dataset Class

\n\n---","metadata":{}},{"cell_type":"code","source":"class HubMAP_Dataset(Dataset):\n def __init__(self, train):\n self.dataframe = pd.read_csv(os.path.join(INPUT_PATH, \"train.csv\"))\n \n # Stratified dataframe division by organ trait\n self.input = self.dataframe.groupby(\"organ\", group_keys=False).apply(lambda x: x.sample(min(len(x), int(0.8*len(x)))))\n \n if not train:\n self.input = self.dataframe.drop(self.input.index)\n \n # Augmentations\n self.transform = A.Compose([\n A.Resize(height = IMG_SHAPE[0], width = IMG_SHAPE[1]),\n A.Normalize(mean=IMG_MEAN, std=IMG_STD)\n ])\n \n def __len__(self):\n return self.input.shape[0]\n \n def __getitem__(self, ndx):\n # Image properties\n img_path = str(self.input.iloc[ndx, 0]) + \".tiff\"\n shape = (int(self.input.iloc[ndx, 3]), int(self.input.iloc[ndx, 4]))\n \n # Image/Mask Loading\n img = Image.open(os.path.join(INPUT_PATH, \"train_images\", img_path))\n img = np.array(img)\n \n mask = rle_decode(self.input.iloc[ndx, -3], shape)\n \n aug_data = self.transform(image = img, mask = mask)\n \n return torch.tensor(np.moveaxis(aug_data[\"image\"], 2, 0), dtype=torch.float32), torch.tensor(np.expand_dims(aug_data[\"mask\"], axis=0), dtype=torch.int16)\n ","metadata":{"execution":{"iopub.status.busy":"2022-06-27T14:53:55.418978Z","iopub.execute_input":"2022-06-27T14:53:55.419408Z","iopub.status.idle":"2022-06-27T14:53:55.431653Z","shell.execute_reply.started":"2022-06-27T14:53:55.419373Z","shell.execute_reply":"2022-06-27T14:53:55.430539Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"

3.1 VISUALIZE THE TRAIN DATA

\n\n---","metadata":{}},{"cell_type":"code","source":"img = Image.open(\"../input/hubmap-organ-segmentation/train_images/10044.tiff\")\nimg","metadata":{"execution":{"iopub.status.busy":"2022-06-27T13:35:32.592439Z","iopub.execute_input":"2022-06-27T13:35:32.593175Z","iopub.status.idle":"2022-06-27T13:35:37.314017Z","shell.execute_reply.started":"2022-06-27T13:35:32.593136Z","shell.execute_reply":"2022-06-27T13:35:37.312913Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"img = np.array(img)\n\ndf = pd.read_csv(os.path.join(INPUT_PATH, \"train.csv\"))\nmask = rle_decode(df.iloc[0, -3], (3000, 3000))\n\nplt.figure(figsize=(15, 15))\nplt.imshow(img)\nplt.imshow(mask, alpha=0.5)","metadata":{"execution":{"iopub.status.busy":"2022-06-27T13:35:37.315259Z","iopub.execute_input":"2022-06-27T13:35:37.315778Z","iopub.status.idle":"2022-06-27T13:35:40.233083Z","shell.execute_reply.started":"2022-06-27T13:35:37.315716Z","shell.execute_reply":"2022-06-27T13:35:40.231968Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n\n\n\n\n\n

\n 4  MODELLING    ⤒\n

\n\n---\n\n

4.0 Train Loop

\n\n---","metadata":{}},{"cell_type":"code","source":"# Model Elements\n\n#model = monai.networks.nets.UNet(\n# spatial_dims=2,\n# in_channels=3,\n# out_channels=1,\n# channels=(16, 32, 64, 128, 256),\n# strides=(2, 2, 2, 2),\n# num_res_units=2,\n# ).to(device)\n\nmodel = monai.networks.nets.SegResNet(\n spatial_dims = 2,\n in_channels=3,\n out_channels=1,\n upsample_mode = \"deconv\").to(device)\n\nloss_function = monai.losses.DiceLoss(sigmoid=True)\noptimizer = torch.optim.Adam(model.parameters(), 1e-3)\n\npost_trans = Compose([EnsureType(), Activations(sigmoid=True), AsDiscrete(threshold=0.5)])\ndice_metric = DiceMetric(include_background=True, reduction=\"mean\", get_not_nans=False)\n\n# Dataset\ntrain_ds = HubMAP_Dataset(train=True)\ntrain_loader = DataLoader(train_ds, batch_size=16)\n\nval_ds = HubMAP_Dataset(train=False)\nval_loader = DataLoader(val_ds, batch_size=4)\n\n# Training Loop\n\nval_interval = 2\nbest_metric = -1\nbest_metric_epoch = -1\nnum_epochs = 10\nepoch_loss_values = list()\nmetric_values = list()\nwriter = SummaryWriter()\nfor epoch in range(num_epochs):\n print(\"-\" * 10)\n print(f\"epoch {epoch + 1}/{num_epochs}\")\n model.train()\n epoch_loss = 0\n step = 0\n for batch_data in tqdm(train_loader):\n step += 1\n inputs, labels = batch_data[0].to(device), batch_data[1].to(device)\n optimizer.zero_grad()\n outputs = model(inputs)\n loss = loss_function(outputs, labels)\n loss.backward()\n optimizer.step()\n epoch_loss += loss.item()\n epoch_len = len(train_ds) // train_loader.batch_size\n #print(f\"{step}/{epoch_len}, train_loss: {loss.item():.4f}\")\n writer.add_scalar(\"train_loss\", loss.item(), epoch_len * epoch + step)\n epoch_loss /= step\n epoch_loss_values.append(epoch_loss)\n print(f\"epoch {epoch + 1} average loss: {epoch_loss:.4f}\")\n\n if (epoch + 1) % val_interval == 0:\n model.eval()\n with torch.no_grad():\n val_images = None\n val_labels = None\n val_outputs = None\n for val_data in tqdm(val_loader):\n val_images, val_labels = val_data[0].to(device), val_data[1].to(device)\n roi_size = (96, 96)\n sw_batch_size = 4\n val_outputs = sliding_window_inference(val_images, roi_size, sw_batch_size, model)\n val_outputs = [post_trans(i) for i in decollate_batch(val_outputs)]\n # compute metric for current iteration\n dice_metric(y_pred=val_outputs, y=val_labels)\n # aggregate the final mean dice result\n metric = dice_metric.aggregate().item()\n # reset the status for next validation round\n dice_metric.reset()\n metric_values.append(metric)\n if metric > best_metric:\n best_metric = metric\n best_metric_epoch = epoch + 1\n torch.save(model.state_dict(), \"best_metric_model_segmentation2d_array.pth\")\n print(\"saved new best metric model\")\n print(\n \"current epoch: {} current mean dice: {:.4f} best mean dice: {:.4f} at epoch {}\".format(\n epoch + 1, metric, best_metric, best_metric_epoch\n )\n )\n writer.add_scalar(\"val_mean_dice\", metric, epoch + 1)\n # plot the last model output as GIF image in TensorBoard with the corresponding image and label\n plot_2d_or_3d_image(val_images, epoch + 1, writer, index=0, tag=\"image\")\n plot_2d_or_3d_image(val_labels, epoch + 1, writer, index=0, tag=\"label\")\n plot_2d_or_3d_image(val_outputs, epoch + 1, writer, index=0, tag=\"output\")","metadata":{"execution":{"iopub.status.busy":"2022-06-27T14:54:06.917606Z","iopub.execute_input":"2022-06-27T14:54:06.917987Z","iopub.status.idle":"2022-06-27T17:43:35.994869Z","shell.execute_reply.started":"2022-06-27T14:54:06.917955Z","shell.execute_reply":"2022-06-27T17:43:35.993816Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"

4.1 Trained Model Results

\n\n---","metadata":{}},{"cell_type":"code","source":"for idx,batch in enumerate(val_loader):\n if idx == 0:\n img, mask = batch\n else:\n break","metadata":{"execution":{"iopub.status.busy":"2022-06-27T17:43:35.997112Z","iopub.execute_input":"2022-06-27T17:43:35.998124Z","iopub.status.idle":"2022-06-27T17:43:37.218714Z","shell.execute_reply.started":"2022-06-27T17:43:35.99808Z","shell.execute_reply":"2022-06-27T17:43:37.217512Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#model = monai.networks.nets.UNet(\n# spatial_dims=2,\n# in_channels=3,\n# out_channels=1,\n# channels=(16, 32, 64, 128, 256),\n# strides=(2, 2, 2, 2),\n# num_res_units=2,\n# )\n\nmodel = monai.networks.nets.SegResNet(\n spatial_dims = 2,\n in_channels=3,\n out_channels=1,\n upsample_mode = \"deconv\")\n\nmodel.load_state_dict(torch.load(\"./best_metric_model_segmentation2d_array.pth\"))\n\nout = model(img)","metadata":{"execution":{"iopub.status.busy":"2022-06-27T17:43:37.283656Z","iopub.execute_input":"2022-06-27T17:43:37.284031Z","iopub.status.idle":"2022-06-27T17:43:37.777018Z","shell.execute_reply.started":"2022-06-27T17:43:37.283988Z","shell.execute_reply":"2022-06-27T17:43:37.775974Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, ax = plt.subplots(2, 2, figsize=(15, 15))\n\nfor i in range(2):\n for j in range(2):\n if i == 1:\n aux = 1\n else:\n aux = 0\n \n ax[i,j].imshow(img[0+i+j+aux, 0, ...], cmap=\"gray\")\n ax[i,j].imshow(out[0+i+j+aux, 0, ...].detach().numpy(), cmap=\"Greens\", alpha = 0.5)\n ax[i,j].imshow(mask[0+i+j+aux, 0, ...], cmap=\"Blues\", alpha = 0.3)","metadata":{"execution":{"iopub.status.busy":"2022-06-27T17:47:21.722109Z","iopub.execute_input":"2022-06-27T17:47:21.722513Z","iopub.status.idle":"2022-06-27T17:47:22.81488Z","shell.execute_reply.started":"2022-06-27T17:47:21.72247Z","shell.execute_reply":"2022-06-27T17:47:22.813919Z"},"trusted":true},"execution_count":null,"outputs":[]}]}