{"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":"In this notebook we offer solution for the Vesuvis Ink detection problem. \n\nUpon analysing [ink detection tutorial](https://www.kaggle.com/code/jpposma/vesuvius-challenge-ink-detection-tutorial) we noticed several segments where improvements could be made. \n\n- Random sampling of points for training set is not good.\n- We should use a CNN architecture that is more appropriate for 3D convolutions","metadata":{}},{"cell_type":"markdown","source":"### Setup","metadata":{}},{"cell_type":"code","source":"import torch\nimport torch.nn as nn\nimport torch.optim as optim\nimport numpy as np\nimport glob\nimport random\nimport PIL.Image as Image\nimport torch.utils.data as data\nfrom scipy.ndimage import convolve\nimport matplotlib.pyplot as plt\nimport matplotlib.patches as patches\nfrom tqdm import tqdm\nfrom ipywidgets import interact, fixed\nimport sys\nimport cv2\n\nPREFIX = '/kaggle/input/vesuvius-challenge-ink-detection/train/'\nBUFFER = 50  # Buffer size in x and y direction\nZ_START = 27 # First slice in the z direction to use\nZ_DIM = 10   # Number of slices in the z direction\nTRAINING_STEPS = 42000\nLEARNING_RATE = 0.05\nBATCH_SIZE = 32\nDEVICE = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\n\nplt.imshow(Image.open(PREFIX+\"1/ir.png\"), cmap=\"gray\")","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-06-25T19:36:23.551947Z","iopub.execute_input":"2023-06-25T19:36:23.552386Z","iopub.status.idle":"2023-06-25T19:36:28.135773Z","shell.execute_reply.started":"2023-06-25T19:36:23.552352Z","shell.execute_reply":"2023-06-25T19:36:28.131739Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Loading of binary images:\n\n* **mask.png**: a mask of which pixels contain data, and which pixels we should ignore.\n* **inklabels.png**: our label data: whether a pixel contains ink or no ink (which has been hand-labeled based on the infrared photo).","metadata":{}},{"cell_type":"code","source":"mask_1 = np.array(Image.open(PREFIX+\"1/mask.png\").convert('1'))\nlabel_1 = torch.from_numpy(np.array(Image.open(PREFIX+\"1/inklabels.png\"))).gt(0).float().to(DEVICE)\n\nmask_2 = np.array(Image.open(PREFIX+\"2/mask.png\").convert('1'))\nlabel_2 = torch.from_numpy(np.array(Image.open(PREFIX+\"2/inklabels.png\"))).gt(0).float().to(DEVICE)\n\nmask_3 = np.array(Image.open(PREFIX+\"3/mask.png\").convert('1'))\nlabel_3 = torch.from_numpy(np.array(Image.open(PREFIX+\"3/inklabels.png\"))).gt(0).float().to(DEVICE)","metadata":{"execution":{"iopub.status.busy":"2023-06-25T19:36:28.139414Z","iopub.execute_input":"2023-06-25T19:36:28.139763Z","iopub.status.idle":"2023-06-25T19:36:31.533012Z","shell.execute_reply.started":"2023-06-25T19:36:28.139733Z","shell.execute_reply":"2023-06-25T19:36:31.531979Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### What percentage of papirus is covered with ink?\n\nWe would also like to know what percentage of papirus is covered with ink. This information is usefull, since it will tell us how likely we will be to sample ink pixel in training dataset.","metadata":{}},{"cell_type":"code","source":"def ink_coverage_percentage(mask, label_cpu):\n    # We want to calculate amount of papirus that is covered with ink. This will be useful for the purpose of sampling data for training and testing...\n    return label_cpu.numpy().sum() / mask.sum()","metadata":{"execution":{"iopub.status.busy":"2023-06-25T18:07:08.294920Z","iopub.execute_input":"2023-06-25T18:07:08.297485Z","iopub.status.idle":"2023-06-25T18:07:08.304646Z","shell.execute_reply.started":"2023-06-25T18:07:08.297450Z","shell.execute_reply":"2023-06-25T18:07:08.303652Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ink_coverage_percentage(mask_1, label_1.cpu())","metadata":{"execution":{"iopub.status.busy":"2023-06-25T18:07:10.035311Z","iopub.execute_input":"2023-06-25T18:07:10.035738Z","iopub.status.idle":"2023-06-25T18:07:10.332793Z","shell.execute_reply.started":"2023-06-25T18:07:10.035701Z","shell.execute_reply":"2023-06-25T18:07:10.331793Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Surprisingly, less than 20% of fragments are covered with ink. This poses a problem. ","metadata":{}},{"cell_type":"markdown","source":"### The problem with random sampling of points for training dataset\n\n<a href=\"https://user-images.githubusercontent.com/22727759/224853655-3fad9edb-c798-452e-94d0-f74efe71c08e.mp4\"><img src=\"https://user-images.githubusercontent.com/22727759/224853385-ed190d89-f466-469c-82a9-499881759d57.gif\"/></a>\n\nSince ink pixels amount to less than 20% of the surface of scroll fragments, by randomly sampling points for our training dataset we are less likely to select points that contain ink that those non-ink ones. By doing so, we would end up with a model that is better in detecting non-ink pixels, which is not our initial idea. \n\n### Alternative:\n\nOur training set should have equal amount of ink and non-ink pixels. Therefore, our training dataset will have 50% ink and 50% non-ink datapoints. Additionally, we will make a categorization of ink pixels into border pixels and middle (middle of letter) pixels. We believe that those pixels that are in the middle of the letter contain more information that those on the edge. We will sample those gropus equally as well.\n\nSo, our training dataset will have 50% non-ink (blank) datapoints, 25% \"edge of the letter\" pixels and 25% middle of the letter pixels.","metadata":{}},{"cell_type":"markdown","source":"### Edge detection\n\nWhat we notice from the ink-label images is that the line stroke of the ink is actually very wide. In order to find out what is the line width on this image, we need to perform edge detection. \n","metadata":{}},{"cell_type":"markdown","source":"### Finding \"inner-letter\" pixels\n\nUpon expecting line stroke, and estimating that average ink line is 174 pixels wide, we make a decision that the **inner-letter** pixels are those pixels that contain ink and that have 35 or more pixels in all directions (left, right, above, below) that also contain ink. \n\nHere is the code for finding such pixels\n","metadata":{}},{"cell_type":"code","source":"def detect_inner_ink_pixels(image, pixels_num=35):\n    # Assuming 'image' is your input image with ink and non-ink pixels\n\n    # Binarize the image to obtain a binary mask\n    _, binary = cv2.threshold(image, 0, 1, cv2.THRESH_BINARY)\n\n    # Apply morphological erosion operation in order to extract inner pixels\n    kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (3, 3))\n    eroded = cv2.erode(image, kernel, iterations=pixels_num)\n\n    # Obtain the middle ink pixels\n    inner_ink_pixels = np.logical_and(binary, eroded)\n\n    return inner_ink_pixels.astype(np.uint8)","metadata":{"execution":{"iopub.status.busy":"2023-06-25T18:07:15.691319Z","iopub.execute_input":"2023-06-25T18:07:15.691685Z","iopub.status.idle":"2023-06-25T18:07:15.697862Z","shell.execute_reply.started":"2023-06-25T18:07:15.691656Z","shell.execute_reply":"2023-06-25T18:07:15.696656Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# TODO show image of middle of the letter pixels\ninner_pixels_1 = detect_inner_ink_pixels(label_1.cpu().numpy(), pixels_num=35)\ninner_pixels_2 = detect_inner_ink_pixels(label_2.cpu().numpy(), pixels_num=35)\ninner_pixels_3 = detect_inner_ink_pixels(label_3.cpu().numpy(), pixels_num=35)\n\n\n# inner_pixels_masked = np.ma.masked_where(inner_pixels_2 == 0, inner_pixels_2)\n# fig, ax = plt.subplots()\n# ax.imshow(label_2.cpu(), cmap='gray')\n# ax.imshow(inner_pixels_masked, cmap='autumn')\n# ax.set_title('Inner pixel')\n# plt.show()\n# del inner_pixels_masked","metadata":{"execution":{"iopub.status.busy":"2023-06-25T18:07:18.446174Z","iopub.execute_input":"2023-06-25T18:07:18.447125Z","iopub.status.idle":"2023-06-25T18:07:23.120868Z","shell.execute_reply.started":"2023-06-25T18:07:18.447081Z","shell.execute_reply":"2023-06-25T18:07:23.119844Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Inner pixels are marked with red color in the above image","metadata":{}},{"cell_type":"markdown","source":"Next, we'll load the 3d x-ray of the fragment. This is represented as a .tif image stack. The image stack is an array of 16-bit grayscale images. Each image represents a \"slice\" in the z-direction, going from below the papyrus, to above the papyrus. We'll convert it to a 4D tensor of 32-bit floats. We'll also convert the pixel values to the range [0, 1].\n\nTo save memory, we'll only load the innermost slices (`Z_DIM` of them). Let's look at them when we're done.","metadata":{}},{"cell_type":"code","source":"# Load the 3d x-ray scan, one slice at a time\nimages_1 = [np.array(Image.open(filename), dtype=np.float32)/65535.0 for filename in tqdm(sorted(glob.glob(PREFIX+\"1/surface_volume/*.tif\"))[Z_START:Z_START+Z_DIM])]\nimage_stack_1 = torch.stack([torch.from_numpy(image) for image in images_1], dim=0).to(DEVICE)\ndel images_1\n\nimages_2 = [np.array(Image.open(filename), dtype=np.float32)/65535.0 for filename in tqdm(sorted(glob.glob(PREFIX+\"2/surface_volume/*.tif\"))[Z_START:Z_START+Z_DIM])]\nimage_stack_2 = torch.stack([torch.from_numpy(image) for image in images_2], dim=0).to(DEVICE)\ndel images_2\n\nimages_3 = [np.array(Image.open(filename), dtype=np.float32)/65535.0 for filename in tqdm(sorted(glob.glob(PREFIX+\"3/surface_volume/*.tif\"))[Z_START:Z_START+Z_DIM])]\nimage_stack_3 = torch.stack([torch.from_numpy(image) for image in images_3], dim=0).to(DEVICE)\ndel images_3\n","metadata":{"execution":{"iopub.status.busy":"2023-06-25T18:07:23.128682Z","iopub.execute_input":"2023-06-25T18:07:23.131168Z","iopub.status.idle":"2023-06-25T18:09:00.139866Z","shell.execute_reply.started":"2023-06-25T18:07:23.131131Z","shell.execute_reply":"2023-06-25T18:09:00.138848Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Can you see the ink in these slices of the 3d x-ray scan..? Neither can we.\n\nNow we'll create a dataset of subvolumes. We use a small rectangle around the letter \"P\" for our evaluation, and we'll exclude those pixels from the training set. (It's actually a Greek letter \"rho\", which looks similar to our \"P\".)","metadata":{}},{"cell_type":"code","source":"rect_1 = (1100, 3500, 700, 950) # these are the coordinates for the fragment 1\nrect_2 = (2500, 6000, 1500, 1000) # these are the coordinates for the fragment 2\nrect_3 = (2100, 4400, 800, 950) # these are the coordinates for the fragment 3\n\nrect = rect_1\nfig, ax = plt.subplots()\nax.imshow(label_1.cpu())\npatch = patches.Rectangle((rect[0], rect[1]), rect[2], rect[3], linewidth=1, edgecolor='r', facecolor='none')\nax.add_patch(patch)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-25T18:09:38.468167Z","iopub.execute_input":"2023-06-25T18:09:38.468679Z","iopub.status.idle":"2023-06-25T18:09:40.414573Z","shell.execute_reply.started":"2023-06-25T18:09:38.468640Z","shell.execute_reply":"2023-06-25T18:09:40.413518Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class SubvolumeDataset(data.Dataset):\n    def __init__(self, image_stack, label, pixels):\n        self.image_stack = image_stack\n        self.label = label\n        self.pixels = pixels\n    def __len__(self):\n        return len(self.pixels)\n    def __getitem__(self, index):\n        y, x = self.pixels[index]\n        subvolume = self.image_stack[:, y-BUFFER:y+BUFFER+1, x-BUFFER:x+BUFFER+1].view(1, Z_DIM, BUFFER*2+1, BUFFER*2+1)\n        inklabel = self.label[y, x].view(1)\n        return subvolume, inklabel","metadata":{"execution":{"iopub.status.busy":"2023-06-25T18:09:44.205668Z","iopub.execute_input":"2023-06-25T18:09:44.206025Z","iopub.status.idle":"2023-06-25T18:09:44.214264Z","shell.execute_reply.started":"2023-06-25T18:09:44.205996Z","shell.execute_reply":"2023-06-25T18:09:44.212987Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class MyInkDetector(torch.nn.Module):\n    def __init__(self):\n        super().__init__()\n\n        self.conv1 = nn.Sequential(\n            nn.Conv3d(1, 16, kernel_size=5, stride=(1, 2, 2), padding=1),\n            nn.BatchNorm3d(16),\n            nn.LeakyReLU(0.01, inplace=True)\n        )\n\n        self.conv2 = nn.Sequential(\n            nn.Conv3d(16, 32, kernel_size=5, stride=2, padding=1),\n            nn.BatchNorm3d(32),\n            nn.LeakyReLU(0.01, inplace=True)\n        )\n\n        self.conv3 = nn.Sequential(\n            nn.Conv3d(32, 64, kernel_size=5, stride=2, padding=1),\n            nn.BatchNorm3d(64),\n            nn.LeakyReLU(0.01, inplace=True)\n        )\n\n        self.middle_layer = nn.Sequential(\n            nn.AdaptiveAvgPool3d(1),\n            nn.Flatten()\n        )\n\n        self.fc_layer = nn.Sequential(\n            nn.Linear(64, 128),\n            nn.ReLU(inplace=True),\n            nn.Linear(128, 128),\n            nn.ReLU(inplace=True),\n            nn.Linear(128, 1)\n        )\n\n    def forward(self, x):\n        x = self.conv1(x)\n        x = self.conv2(x)\n        x = self.conv3(x)\n        x = self.middle_layer(x)\n        x = self.fc_layer(x)\n        return x","metadata":{"execution":{"iopub.status.busy":"2023-06-25T18:09:48.489928Z","iopub.execute_input":"2023-06-25T18:09:48.490270Z","iopub.status.idle":"2023-06-25T18:09:48.501997Z","shell.execute_reply.started":"2023-06-25T18:09:48.490242Z","shell.execute_reply":"2023-06-25T18:09:48.500939Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model = MyInkDetector().to(DEVICE)","metadata":{"execution":{"iopub.status.busy":"2023-06-25T18:09:51.465374Z","iopub.execute_input":"2023-06-25T18:09:51.465751Z","iopub.status.idle":"2023-06-25T18:09:51.518474Z","shell.execute_reply.started":"2023-06-25T18:09:51.465721Z","shell.execute_reply":"2023-06-25T18:09:51.517579Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"num_px_samples_03 = int(28444028 * 0.3)\nnum_px_samples_04 = int(28444028 * 0.4)","metadata":{"execution":{"iopub.status.busy":"2023-06-25T18:09:53.359281Z","iopub.execute_input":"2023-06-25T18:09:53.359950Z","iopub.status.idle":"2023-06-25T18:09:53.364791Z","shell.execute_reply.started":"2023-06-25T18:09:53.359917Z","shell.execute_reply":"2023-06-25T18:09:53.363715Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Generating pixel lists from fragments","metadata":{}},{"cell_type":"code","source":"# proportions = [0.2, 0.2, 0.6]\ndef generate_pixel_lists(mask, rect, label, inner_pixels, num_samples):\n    \n    print(\"Generating pixel lists...\")\n    # Split our dataset into train and val. The pixels inside the rect are the\n    # val set, and the pixels outside the rect are the train set.\n    # Adapted from https://www.kaggle.com/code/jamesdavey/100x-faster-pixel-coordinate-generator-1s-runtime\n    # Create a Boolean array of the same shape as the bitmask, initially all True\n    not_border = np.zeros(mask.shape, dtype=bool)\n    not_border[BUFFER:mask.shape[0]-BUFFER, BUFFER:mask.shape[1]-BUFFER] = True\n    arr_mask = np.array(mask) * not_border\n    inside_rect = np.zeros(mask.shape, dtype=bool) * arr_mask\n    # Sets all indexes with inside_rect array to True\n    inside_rect[rect[1]:rect[1]+rect[3]+1, rect[0]:rect[0]+rect[2]+1] = True\n    # Set the pixels within the inside_rect to False\n    outside_rect = np.ones(mask.shape, dtype=bool) * arr_mask\n    outside_rect[rect[1]:rect[1]+rect[3]+1, rect[0]:rect[0]+rect[2]+1] = False\n    pixels_inside_rect = np.argwhere(inside_rect)\n    pixels_outside_rect = np.argwhere(outside_rect)\n\n    del not_border\n    del arr_mask\n\n    del pixels_outside_rect\n    combined_pixels_list = []\n\n\n    inner_letter_pixels_list = np.argwhere(inner_pixels * outside_rect)\n    num_samples_inner_letter_px = int(num_samples * 0.25)\n    samples_list_inner_letter = random.choices(inner_letter_pixels_list,\n                                               k=num_samples_inner_letter_px)\n    del inner_letter_pixels_list\n    combined_pixels_list = combined_pixels_list + samples_list_inner_letter\n    del samples_list_inner_letter\n    print('Finished sampling list #1')\n\n    border_letter_pixels_list = np.argwhere(label.cpu().numpy() - inner_pixels)\n    num_samples_border_letter_px = int(num_samples * 0.25)\n    samples_list_border_letter = random.choices(border_letter_pixels_list,\n                                                k=num_samples_border_letter_px)\n    del border_letter_pixels_list\n    combined_pixels_list = combined_pixels_list + samples_list_border_letter\n    del samples_list_border_letter\n    print('Finished sampling list #2')\n\n\n    blank_paper_pixels_list = np.argwhere((mask - label.cpu().numpy()) * outside_rect)\n    num_samples_blank_paper_px = int(num_samples * 0.5)\n    samples_list_blank_paper = random.choices(blank_paper_pixels_list,\n                                              k=num_samples_blank_paper_px)\n    del blank_paper_pixels_list\n    combined_pixels_list = combined_pixels_list + samples_list_blank_paper\n    \n    del samples_list_blank_paper\n    print('Finished sampling list #3')\n    print('Pixel list successfully generated!')\n    del inside_rect\n    del outside_rect\n    \n    return np.array(combined_pixels_list)\n    ","metadata":{"execution":{"iopub.status.busy":"2023-06-25T18:09:58.625974Z","iopub.execute_input":"2023-06-25T18:09:58.626444Z","iopub.status.idle":"2023-06-25T18:09:58.643639Z","shell.execute_reply.started":"2023-06-25T18:09:58.626409Z","shell.execute_reply":"2023-06-25T18:09:58.642613Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pixel_list_fragment_1 = generate_pixel_lists(mask_1, rect_1, label_1,\n                                             inner_pixels_1, num_samples=num_px_samples_03)\npixel_list_fragment_2 = generate_pixel_lists(mask_2, rect_2, label_2,\n                                             inner_pixels_2, num_samples=num_px_samples_04)\npixel_list_fragment_3 = generate_pixel_lists(mask_3, rect_3, label_3,\n                                             inner_pixels_3, num_samples=num_px_samples_03)","metadata":{"execution":{"iopub.status.busy":"2023-06-25T18:10:22.235523Z","iopub.execute_input":"2023-06-25T18:10:22.235895Z","iopub.status.idle":"2023-06-25T18:11:11.603721Z","shell.execute_reply.started":"2023-06-25T18:10:22.235867Z","shell.execute_reply":"2023-06-25T18:11:11.602680Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# TRAINING_STEPS = 3000\nprint(\"Training...\")\ndataset_fragment_1 = SubvolumeDataset(image_stack_1, label_1, pixel_list_fragment_1)\ndataset_fragment_2 = SubvolumeDataset(image_stack_2, label_2, pixel_list_fragment_2)\ndataset_fragment_3 = SubvolumeDataset(image_stack_3, label_3, pixel_list_fragment_3)\n\ntrain_dataset = data.ConcatDataset([dataset_fragment_1, dataset_fragment_2, dataset_fragment_3])\ndel dataset_fragment_1\ndel dataset_fragment_2\ndel dataset_fragment_3\ntrain_loader = data.DataLoader(train_dataset, batch_size=BATCH_SIZE, shuffle=True)\n\ncriterion = nn.BCEWithLogitsLoss()\noptimizer = optim.SGD(model.parameters(), lr=LEARNING_RATE)\nscheduler = torch.optim.lr_scheduler.OneCycleLR(optimizer, max_lr=LEARNING_RATE, total_steps=TRAINING_STEPS)\nmodel.train()\nrunning_loss = 0.0\nfor i, (subvolumes, inklabels) in tqdm(enumerate(train_loader), total=TRAINING_STEPS):\n    if i >= TRAINING_STEPS:\n        break\n    optimizer.zero_grad()\n    outputs = model(subvolumes.to(DEVICE))\n    loss = criterion(outputs, inklabels.to(DEVICE))\n    loss.backward()\n    optimizer.step()\n    scheduler.step()\n    running_loss += loss.item()\n    if i % 3000 == 3000-1:\n        print(\"Loss:\", running_loss / 3000)\n        running_loss = 0.0","metadata":{"execution":{"iopub.status.busy":"2023-06-25T18:11:11.605460Z","iopub.execute_input":"2023-06-25T18:11:11.605900Z","iopub.status.idle":"2023-06-25T18:34:27.585761Z","shell.execute_reply.started":"2023-06-25T18:11:11.605865Z","shell.execute_reply":"2023-06-25T18:34:27.584797Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"del train_loader","metadata":{"execution":{"iopub.status.busy":"2023-06-25T18:40:42.040821Z","iopub.execute_input":"2023-06-25T18:40:42.041197Z","iopub.status.idle":"2023-06-25T18:40:42.045647Z","shell.execute_reply.started":"2023-06-25T18:40:42.041166Z","shell.execute_reply":"2023-06-25T18:40:42.044614Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Finally, we'll generate a prediction image. We'll use the model to predict the presence of ink for each pixel in our rectangle (the val set). Conceptually it looks like this:\n\n<a href=\"https://user-images.githubusercontent.com/22727759/224853653-7cffd0a4-c6fa-49a2-93c1-e3c820863a51.mp4\"><img src=\"https://user-images.githubusercontent.com/22727759/224853379-09ae991e-02be-4ecc-a652-313165b3005c.gif\"/></a>\n\n\nThis should take about a minute.\n\nRemember that the model has never seen the label data within the rectangle before!\n\nWe'll plot it side-by-side with the label image. Are you able to recognize the letter \"P\" in it?","metadata":{}},{"cell_type":"markdown","source":"### Evaluating the model on a cut out rect","metadata":{}},{"cell_type":"code","source":"rect = rect_1\nnot_border = np.zeros(mask_1.shape, dtype=bool)\nnot_border[BUFFER:mask_1.shape[0]-BUFFER, BUFFER:mask_1.shape[1]-BUFFER] = True\narr_mask = np.array(mask_1) * not_border\ninside_rect = np.zeros(mask_1.shape, dtype=bool) * arr_mask\ninside_rect[rect[1]:rect[1]+rect[3]+1, rect[0]:rect[0]+rect[2]+1] = True\n\n\ndel not_border\ndel arr_mask\n\npx_inside_rect_1 = np.argwhere(inside_rect)\ndel inside_rect\n\n\n\neval_dataset = SubvolumeDataset(image_stack_1, label_1, px_inside_rect_1)\neval_loader = data.DataLoader(eval_dataset, batch_size=BATCH_SIZE, shuffle=False)\noutput = torch.zeros_like(label_1).float()\nmodel.eval()\nwith torch.no_grad():\n    for i, (subvolumes, _) in enumerate(tqdm(eval_loader)):\n        for j, value in enumerate(model(subvolumes.to(DEVICE))):\n            output[tuple(px_inside_rect_1[i*BATCH_SIZE+j])] = value\n\nfig, (ax1, ax2) = plt.subplots(1, 2)\nax1.imshow(output.cpu(), cmap='gray')\nax2.imshow(label_1.cpu(), cmap='gray')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-25T19:53:42.484130Z","iopub.execute_input":"2023-06-25T19:53:42.485220Z","iopub.status.idle":"2023-06-25T19:56:22.665321Z","shell.execute_reply.started":"2023-06-25T19:53:42.485185Z","shell.execute_reply":"2023-06-25T19:56:22.664143Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Since our output has to be binary, we have to choose a threshold, say 40% confidence.","metadata":{}},{"cell_type":"code","source":"THRESHOLD = 0.4\nfig, (ax1, ax2) = plt.subplots(1, 2)\nax1.imshow(output.gt(THRESHOLD).cpu(), cmap='gray')\nax2.imshow(label_1.cpu(), cmap='gray')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-25T19:58:30.714757Z","iopub.execute_input":"2023-06-25T19:58:30.715142Z","iopub.status.idle":"2023-06-25T19:58:33.357926Z","shell.execute_reply.started":"2023-06-25T19:58:30.715111Z","shell.execute_reply":"2023-06-25T19:58:33.356855Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Measures of predictive performance\n\n- #### Classification accuracy\n- #### Sensitivity\n- #### Specificity\n- #### Log-loss","metadata":{}},{"cell_type":"code","source":"def ca(a: torch.Tensor, b: torch.Tensor):\n\n    \"\"\"\n    ca is short for Classification accuracy, it is the percentage of correctly classified examples\n    The assumption is that we have binary classification problem, and two classes 1 and 0.\n    :param a: Tensor\n    :param b: Tensor\n    :return:\n    \"\"\"\n\n    score = torch.sum(a.eq(b)) / (a.size()[0] * a.size()[1])\n    return float(score.numpy())\n\ndef calc_metrics(a: torch.Tensor, b: torch.Tensor):\n    \"\"\"\n    Sensitivity is a measure of correctly classified examples of the positive class\n\n    :param a: Tensor\n    :param b: Tensor\n    :return:\n    \"\"\"\n    \n    # Convert tensors to Boolean tensors\n    a = a.bool()\n    b = b.bool()\n    \n    # True Positives (TP): elements where both a and b are 1\n    TP = torch.logical_and(a, b).sum().item()\n\n    # True Negatives (TN): elements where both a and b are 0\n    TN = torch.logical_not(torch.logical_or(a, b)).sum().item()\n\n    # False Positives (FP): elements where a is 1 and b is 0\n    FP = torch.logical_and(torch.logical_not(a), b).sum().item()\n\n    # False Negatives (FN): elements where a is 0 and b is 1\n    FN = torch.logical_and(a, torch.logical_not(b)).sum().item()\n\n    # Calculate sensitivity and specificity\n    sensitivity = TP / (TP + FN)\n    specificity = TN / (TN + FP)\n    precision = TP / (TP + FP)\n    recall = sensitivity\n\n    return sensitivity, specificity, precision, recall\n\ndef f05_score(precision: float, recall: float, beta: float):\n    return (1 + beta**2) * precision * recall / (beta**2 * precision + recall)\n    ","metadata":{"execution":{"iopub.status.busy":"2023-06-25T19:58:40.387863Z","iopub.execute_input":"2023-06-25T19:58:40.388218Z","iopub.status.idle":"2023-06-25T19:58:40.399266Z","shell.execute_reply.started":"2023-06-25T19:58:40.388189Z","shell.execute_reply":"2023-06-25T19:58:40.398188Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Results of the model","metadata":{}},{"cell_type":"code","source":"predictions = output.gt(THRESHOLD).float()\n# rect_1 = (1100, 3500, 700, 950)\n# rect_2 = (2500, 6000, 1500, 1000)\n# rect_3 = (2100, 4400, 800, 950)\npredictions_rect = predictions[3500:4450, 1100:1800]\nval_rect = label_1[3500:4450, 1100:1800]\n# predictions_rect = predictions[4400:5350, 2100:2900]\n# val_rect = label_3[4400:5350, 2100:2900]\nacc = ca(predictions_rect.cpu(), val_rect.cpu())\nsens, spec, precision, recall = calc_metrics(predictions_rect.cpu(), val_rect.cpu())\nf_score = f05_score(precision, recall, beta=0.5)\n\nlog_loss_criterion = torch.nn.BCELoss()\nlog_loss = log_loss_criterion(predictions_rect.cpu(), val_rect.cpu()).numpy()\n\nprint('F0.5 Score: ' + str(f_score))\nprint('Accuracy: ' + str(acc))\nprint('Sensitivity: ' + str(sens))\nprint('Specificity: ' + str(spec))\nprint('Log loss: ' + str(log_loss))","metadata":{"execution":{"iopub.status.busy":"2023-06-25T19:58:58.791027Z","iopub.execute_input":"2023-06-25T19:58:58.791402Z","iopub.status.idle":"2023-06-25T19:58:58.829618Z","shell.execute_reply.started":"2023-06-25T19:58:58.791370Z","shell.execute_reply":"2023-06-25T19:58:58.828632Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"And now the submission code","metadata":{}},{"cell_type":"code","source":"# Adapted from https://www.kaggle.com/code/stainsby/fast-tested-rle/notebook\n# and https://www.kaggle.com/code/kotaiizuka/faster-rle/notebook\ndef rle(output):\n    pixels = np.where(output.flatten().cpu() > THRESHOLD, 1, 0).astype(np.uint8)\n    pixels[0] = 0\n    pixels[-1] = 0\n    runs = np.where(pixels[1:] != pixels[:-1])[0] + 2\n    runs[1::2] = runs[1::2] - runs[:-1:2]\n    return ' '.join(str(x) for x in runs)\n","metadata":{"execution":{"iopub.status.busy":"2023-06-25T14:29:14.040288Z","iopub.execute_input":"2023-06-25T14:29:14.040763Z","iopub.status.idle":"2023-06-25T14:29:14.047789Z","shell.execute_reply.started":"2023-06-25T14:29:14.040638Z","shell.execute_reply":"2023-06-25T14:29:14.046654Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Test set evaluation","metadata":{}},{"cell_type":"code","source":"import pandas as pd\nfrom pathlib import Path\nfrom collections import defaultdict\ntest_path = Path('/kaggle/input/vesuvius-challenge-ink-detection/test/')\n\n# mask_1 = np.array(Image.open(PREFIX+\"1/mask.png\").convert('1'))\n# label_1 = torch.from_numpy(np.array(Image.open(PREFIX+\"1/inklabels.png\"))).gt(0).float().to(DEVICE)\n# train_dataset = SubvolumeDataset(image_stack_1, label_1, pixels_outside_rect)\n\n\ntest_fragments = [test_path / fragment_name for fragment_name in test_path.iterdir()]\ntest_fragments = sorted(map(lambda path: path.resolve(), test_fragments))\nprint(\"All fragments:\", test_fragments)\n\nsubmission = defaultdict(list)\n\nfor test_fragment in test_fragments:\n    mask = np.array(Image.open(str(test_fragment  / \"mask.png\")).convert('1'))\n    surface_volume_paths = sorted((test_fragment / \"surface_volume\").rglob(\"*.tif\"))\n    label = torch.zeros(mask.shape)\n    \n    not_border = np.zeros(mask.shape, dtype=bool)\n    not_border[BUFFER:mask.shape[0]-BUFFER, BUFFER:mask.shape[1]-BUFFER] = True\n    arr_mask = np.array(mask) * not_border\n\n    # we don't convert to torch since it doesn't support uint16\n    # Image loadings\n    print('Loading test image stack')\n    images = [np.array(Image.open(fn), dtype=np.float32)/65535.0 for fn in tqdm(surface_volume_paths[Z_START:Z_START+Z_DIM])]\n    image_stack = torch.stack([torch.from_numpy(image) for image in images], dim=0).to(DEVICE)\n    print('Finished loading the images')\n                             \n    # decide which pixels you are going to use\n    pixels_list = np.argwhere(arr_mask)    \n    \n    eval_dataset = SubvolumeDataset(image_stack, label, pixels_list)\n    eval_loader = data.DataLoader(eval_dataset, batch_size=BATCH_SIZE, shuffle=False)\n    output = torch.zeros_like(label).float()\n    model.eval()\n    with torch.no_grad():\n        for i, (subvolumes, _) in enumerate(tqdm(eval_loader)):\n            for j, value in enumerate(model(subvolumes.to(DEVICE))):\n                output[tuple(pixels_list[i*BATCH_SIZE+j])] = value\n                             \n    \n    prediction_image = output.gt(THRESHOLD).cpu()\n    submission[\"Id\"].append(test_fragment.name)\n    submission[\"Predicted\"].append(rle(prediction_image))\n                             \npd.DataFrame.from_dict(submission).to_csv(\"/kaggle/working/submission.csv\", index=False)","metadata":{"execution":{"iopub.status.busy":"2023-06-25T15:08:52.726996Z","iopub.execute_input":"2023-06-25T15:08:52.727429Z","iopub.status.idle":"2023-06-25T15:08:56.729921Z","shell.execute_reply.started":"2023-06-25T15:08:52.727394Z","shell.execute_reply":"2023-06-25T15:08:56.728462Z"},"trusted":true},"execution_count":null,"outputs":[]}]}