{"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":"### Notebook\n- This is a notebook explaining the [Ink Detection progress prize on Kaggle](https://www.kaggle.com/competitions/vesuvius-challenge), which is part of the larger [Vesuvius Challenge](https://scrollprize.org).\n\n- For more background on the process of ink detection, be sure to check out [Tutorial 4: Ink Detection](https://scrollprize.org/tutorial4) on the Vesuvius Challenge website.\n\n### Goal of the Competition\n- Join the $1,000,000+ Vesuvius Challenge to resurrect an ancient library from the ashes of a volcano. In this competition you are tasked with detecting ink from 3D X-ray scans and reading the contents. Thousands of scrolls were part of a library located in a Roman villa in Herculaneum, a town next to Pompeii. This villa was buried by the Vesuvius eruption nearly 2000 years ago. Due to the heat of the volcano, the scrolls were carbonized, and are now impossible to open without breaking them. These scrolls were discovered a few hundred years ago and have been waiting to be read using modern techniques.\n- Detect ink from 3d x-ray scans of fragments of papyrus which became detached from some of the excavated scrolls.\n- The ink used in the Herculaneum scrolls does not show up readily in X-ray scans but  machine learning models can detect it. \n\n<img src=\"https://user-images.githubusercontent.com/177461/224853397-3cf86dc2-45b4-4e7c-9ec2-28a733791a75.jpg\" width=\"200\"/>\n\n- It's an infrared photo, since the ink is better visible in infrared light.\n\n\n### Ink-Detecttion Tutorial\n- Ink detection is the task of taking data from a 3D X-ray scan around a papyrus surface, and identifying the locations of the inked parts of the papyrus.\n- Ink seems to be radiolucent, making it hard to detect on 3D X-ray scans.\n- Campfire & En-Gedi scrolls: Ink shows up as brighter voxels in 3D X-ray scans, so ink detection can be done by taking the brightest pixel in some voxel region.\nHerculaneum scrolls & fragments: Ink does not seem directly visible in 3D X-ray scans, but there does seem to be data there, since machine learning models can detect it.","metadata":{"_uuid":"748d3554-7759-4443-b1b8-916e19ca50ed","_cell_guid":"a18c7d0f-a17a-4603-956d-2b531016c536","trusted":true}},{"cell_type":"markdown","source":"### Source \n - https://pytorch.org/docs/stable/generated/torch.optim.SGD.html\n - https://pytorch.org/docs/stable/optim.html\n - https://pytorch.org/docs/stable/generated/torch.optim.lr_scheduler.PolynomialLR.html#torch.optim.lr_scheduler.PolynomialLR\n - https://pytorch.org/docs/stable/generated/torch.nn.Conv3d.html\n - https://pytorch.org/docs/stable/nn.functional.html\n - https://stackoverflow.com/questions/68606661/what-is-difference-between-nn-module-and-nn-sequential\n - https://pytorch.org/docs/stable/generated/torch.optim.lr_scheduler.ConstantLR.html#torch.optim.lr_scheduler.ConstantLR\n - https://pytorch.org/docs/stable/nn.html\n - https://pytorch.org/docs/stable/generated/torch.nn.L1Loss.html#torch.nn.L1Loss\n - https://pytorch.org/docs/stable/generated/torch.nn.KLDivLoss.html#torch.nn.KLDivLoss\n - https://discuss.pytorch.org/t/runtimeerror-non-empty-3d-or-4d-batch-mode-tensor-expected-for-input/95422/2","metadata":{}},{"cell_type":"markdown","source":"### TORCH.OPTIM\n- torch.optim is a package implementing various optimization algorithms\n\n### torch.optim.lr_scheduler\n- lr_scheduler provides several methods to adjust the learning rate based on the number of epochs. torch.optim.lr_scheduler.ReduceLROnPlateau allows dynamic learning rate reducing based on some validation measurements.\n- lr_scheduler.OneCycleLR . Sets the learning rate of each parameter group according to the 1cycle learning rate policy.\n-  - lr_scheduler.CosineAnnealingLR .Set the learning rate of each parameter group using a cosine annealing schedule, where \n\n### CONSTANTLR\n- Decays the learning rate of each parameter group by a small constant factor until the number of epoch reaches a pre-defined milestone: total_iters.\n\n### L1LOSS\n - Creates a criterion that measures the mean absolute error (MAE) between each element in the input.\n \n ### KLDIVLOSS\n - The Kullback-Leibler divergence loss.\n ### RMSprop\n - Implements RMSprop algorithm.\n \n ### CHAINEDSCHEDULER\n  - Chains list of learning rate schedulers.\n  \n  ### - Rprop\n  - Implements the resilient backpropagation algorithm.","metadata":{}},{"cell_type":"code","source":"import torch\nimport torch.nn as nn\nimport torch.optim as optim\nimport numpy as np\nimport glob\nimport PIL.Image as Image\nimport torch.utils.data as data\nimport matplotlib.pyplot as plt\nimport matplotlib.patches as patches\nfrom tqdm import tqdm\nfrom ipywidgets import interact, fixed\n\nPREFIX = '/kaggle/input/vesuvius-challenge/train/1/'\nBUFFER = 30  # 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 = 20000\nLEARNING_RATE = 0.05\nBATCH_SIZE = 24\nDEVICE = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\n\nplt.imshow(Image.open(PREFIX+\"ir.png\"), cmap=\"gray\")","metadata":{"_uuid":"91db1348-a896-4607-8686-f6c6df6419ed","_cell_guid":"ff3c9fb9-0c86-4acf-9162-c741c46e53a4","collapsed":false,"jupyter":{"outputs_hidden":false},"_kg_hide-output":false,"execution":{"iopub.status.busy":"2023-06-01T04:51:11.466993Z","iopub.execute_input":"2023-06-01T04:51:11.467988Z","iopub.status.idle":"2023-06-01T04:51:17.174548Z","shell.execute_reply.started":"2023-06-01T04:51:11.467933Z","shell.execute_reply":"2023-06-01T04:51:17.173365Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's load these binary images:\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":{"_uuid":"ee5f3f31-6cc0-4dbf-b4e1-e9afb18bc0ea","_cell_guid":"7607926b-5a4d-4107-8935-955288d53ebd","trusted":true}},{"cell_type":"code","source":"mask = np.array(Image.open(PREFIX+\"mask.png\").convert('1'))\nlabel = torch.from_numpy(np.array(Image.open(PREFIX+\"inklabels.png\"))).gt(0).float().to(DEVICE)\nfig, (ax1, ax2) = plt.subplots(1, 2)\nax1.set_title(\"mask.png\")\nax1.imshow(mask, cmap='gray')\nax2.set_title(\"inklabels.png\")\nax2.imshow(label.cpu(), cmap='gray')\nplt.show()","metadata":{"_uuid":"351fefa9-30dd-4e8d-bf3b-1aaa0fb33905","_cell_guid":"43a9acbe-f2c5-4976-b9ee-027e62c27a83","collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2023-06-01T04:51:26.026658Z","iopub.execute_input":"2023-06-01T04:51:26.027523Z","iopub.status.idle":"2023-06-01T04:51:32.593917Z","shell.execute_reply.started":"2023-06-01T04:51:26.027483Z","shell.execute_reply":"2023-06-01T04:51:32.592700Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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\n- To save memory, we'll only load the innermost slices (`Z_DIM` of them). Let's look at them when we're done.","metadata":{"_uuid":"74a6ae91-2aa2-4eb3-9f3c-956dc54cf017","_cell_guid":"09e95f98-439b-49c3-aae2-550fadd533df","trusted":true}},{"cell_type":"code","source":"# Load the 3d x-ray scan, one slice at a time\nimages = [np.array(Image.open(filename), dtype=np.float32)/65535.0 for filename in tqdm(sorted(glob.glob(PREFIX+\"surface_volume/*.tif\"))[Z_START:Z_START+Z_DIM])]\nimage_stack = torch.stack([torch.from_numpy(image) for image in images], dim=0).to(DEVICE)\n\nfig, axes = plt.subplots(1, len(images), figsize=(15, 3))\nfor image, ax in zip(images, axes):\n  ax.imshow(np.array(Image.fromarray(image).resize((image.shape[1]//20, image.shape[0]//20)), dtype=np.float32), cmap='gray')\n  ax.set_xticks([]); ax.set_yticks([])\nfig.tight_layout()\nplt.show()","metadata":{"_uuid":"929b9e7a-d30c-4462-a7b9-0765a76cdb4e","_cell_guid":"f75858f9-06ad-43cc-be5f-ab09738a58c1","collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2023-06-01T04:51:37.594955Z","iopub.execute_input":"2023-06-01T04:51:37.595672Z","iopub.status.idle":"2023-06-01T04:51:58.532025Z","shell.execute_reply.started":"2023-06-01T04:51:37.595631Z","shell.execute_reply":"2023-06-01T04:51:58.530911Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Load the 3d x-ray scan, one slice at a time\nimages = [np.array(Image.open(filename), dtype=np.float32)/65535.0 for filename in tqdm(sorted(glob.glob(PREFIX+\"surface_volume/*.tif\"))[Z_START:Z_START+Z_DIM])]\nimage_stack = torch.stack([torch.from_numpy(image) for image in images], dim=0).to(DEVICE)\nimage_stack1 = torch.stack([torch.from_numpy(image) for image2 in images], dim=1).to(DEVICE)\n\nfig, axes = plt.subplots(1, len(images), figsize=(15, 3))\nfor image, ax in zip(images, axes):\n   ax.imshow(np.array(Image.fromarray(image).resize((image.shape[1]//20, image.shape[0]//20)), dtype=np.float32), cmap='prism')\n   ax.set_xticks([]); ax.set_yticks([])\nfig.tight_layout()\nplt.show()\nfor image1, ax in zip(images, axes):\n       ax.imshow(np.array(Image.fromarray(image1).resize((image1.shape[1]//20, image1.shape[1]//20)), dtype=np.float32), cmap='gray')\n\n       ax.set_xticks([]); ax.set_yticks([])\nfig.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-01T04:52:43.745433Z","iopub.execute_input":"2023-06-01T04:52:43.745923Z","iopub.status.idle":"2023-06-01T04:52:59.967327Z","shell.execute_reply.started":"2023-06-01T04:52:43.745877Z","shell.execute_reply":"2023-06-01T04:52:59.965848Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- 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":{"_uuid":"18a40c3b-4241-4f8d-b808-e1cab35b8f6e","_cell_guid":"173e0034-8ce4-4a12-9a0b-43cff0022249","trusted":true}},{"cell_type":"code","source":"rect = (1100, 3500, 700, 950)\nfig, ax = plt.subplots()\nax.imshow(label.cpu())\npatch = patches.Rectangle((rect[0], rect[1]), rect[2], rect[3], linewidth=2, edgecolor='r', facecolor='none')\nax.add_patch(patch)\nplt.show()","metadata":{"_uuid":"39a855d9-38f5-4f2d-ac7f-1f65529d1a3e","_cell_guid":"62d60bbe-85f1-4df4-b1e6-f4337268b11c","collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2023-06-01T04:54:29.605959Z","iopub.execute_input":"2023-06-01T04:54:29.606342Z","iopub.status.idle":"2023-06-01T04:54:31.453058Z","shell.execute_reply.started":"2023-06-01T04:54:29.606307Z","shell.execute_reply":"2023-06-01T04:54:31.451988Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- Define a PyTorch dataset and a model.","metadata":{"_uuid":"4061cd69-3dd7-495c-838f-17280c1041b2","_cell_guid":"456fbd67-b8ed-4622-925a-3efc9fb7eb4a","trusted":true}},{"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\n\n# model = nn.Sequential(\n#     nn.Conv3d(1, 16, 3, 1, 1), nn.AvgPool3d(kernel_size=3, stride=2, padding=1),\n#     nn.Conv3d(16, 32, 3, 1, 1), nn.AvgPool3d(kernel_size=3, stride=2, padding=1),\n#     nn.Conv3d(32, 64, 3, 1, 1), nn.AvgPool3d(kernel_size=3, stride=2, padding=1),\n#     nn.Flatten(start_dim=1),\n#     nn.LazyLinear(512), nn.ReLU(),\n#     nn.LazyLinear(1), nn.LogSoftmax()\n# ).to(DEVICE)","metadata":{"execution":{"iopub.status.busy":"2023-05-31T13:02:22.025497Z","iopub.execute_input":"2023-05-31T13:02:22.025964Z","iopub.status.idle":"2023-05-31T13:02:22.032102Z","shell.execute_reply.started":"2023-05-31T13:02:22.025925Z","shell.execute_reply":"2023-05-31T13:02:22.030876Z"},"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\n\n# model = nn.Sequential(\n#         nn.Conv3d(1, 16, 3, 1, 1),\n#         nn.MaxPool3d(2, 2),\n#         nn.Conv3d(16, 32, 3, 1, 1),\n#         nn.MaxPool3d(2, 2),\n#         nn.Conv3d(32, 64, 3, 1, 1),\n#         nn.MaxPool3d(2, 2),\n#         nn.Flatten(start_dim=1),\n#         nn.LazyLinear(128), nn.ReLU(),\n#         nn.Linear(64*2*2*2, 120),\n#         nn.Linear(120, 10),\n#         nn.LazyLinear(1), nn.Sigmoid()\n# ).to(DEVICE)","metadata":{"execution":{"iopub.status.busy":"2023-05-31T13:02:22.033600Z","iopub.execute_input":"2023-05-31T13:02:22.034536Z","iopub.status.idle":"2023-05-31T13:02:22.047111Z","shell.execute_reply.started":"2023-05-31T13:02:22.034495Z","shell.execute_reply":"2023-05-31T13:02:22.045989Z"},"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\n\nmodel = nn.Sequential(\n    nn.Conv3d(1, 16, 3, 1, 1), nn.MaxPool3d(2, 2),\n    nn.Conv3d(16, 32, 3, 1, 1), nn.MaxPool3d(2, 2),\n    nn.Conv3d(32, 64, 3, 1, 1), nn.MaxPool3d(2, 2),\n    nn.Flatten(start_dim=1),\n    nn.LazyLinear(128), nn.ReLU(),\n    nn.LazyLinear(1), nn.Sigmoid()\n).to(DEVICE)","metadata":{"execution":{"iopub.status.busy":"2023-06-01T04:54:34.226708Z","iopub.execute_input":"2023-06-01T04:54:34.227535Z","iopub.status.idle":"2023-06-01T04:54:34.254824Z","shell.execute_reply.started":"2023-06-01T04:54:34.227495Z","shell.execute_reply":"2023-06-01T04:54:34.253673Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- Model training Conceptually:\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>","metadata":{"_uuid":"cfbc9413-b6ff-4a69-ae60-22dddef18f11","_cell_guid":"e7b84e1e-7903-4290-9e09-3bf2cb62b222","trusted":true}},{"cell_type":"code","source":"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\nnot_border = np.zeros(mask.shape, dtype=bool)\nnot_border[BUFFER:mask.shape[0]-BUFFER, BUFFER:mask.shape[1]-BUFFER] = True\narr_mask = np.array(mask) * not_border\ninside_rect = np.zeros(mask.shape, dtype=bool) * arr_mask\n# Sets all indexes with inside_rect array to True\ninside_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\noutside_rect = np.ones(mask.shape, dtype=bool) * arr_mask\noutside_rect[rect[1]:rect[1]+rect[3]+1, rect[0]:rect[0]+rect[2]+1] = False\npixels_inside_rect = np.argwhere(inside_rect)\npixels_outside_rect = np.argwhere(outside_rect)\n\nprint(\"Training...\")\ntrain_dataset = SubvolumeDataset(image_stack, label, pixels_outside_rect)\ntrain_loader = data.DataLoader(train_dataset, batch_size=BATCH_SIZE, shuffle=True)\ncriterion = nn.BCELoss() \n# criterion =nn.TripletMarginWithDistanceLoss(distance_function=nn.PairwiseDistance())\n# criterion = nn.TripletMarginLoss(margin=1.0, p=2)\n# criterion = nn.CosineEmbeddingLoss()\n# criterion = nn.L1Loss()\n# criterion = nn.KLDivLoss(reduction=\"batchmean\")\n# optimizer = optim.Adadelta(model.parameters(), lr=LEARNING_RATE)\n# optimizer = optim.RMSprop(model.parameters(), lr=LEARNING_RATE)\noptimizer = optim.ASGD(model.parameters(), lr=LEARNING_RATE) # Ink detected\n# optimizer = optim.Rprop(model.parameters(), lr=LEARNING_RATE)\n# optimizer = optim.Adagrad(model.parameters(), lr=LEARNING_RATE)\n# scheduler = torch.optim.lr_scheduler.CyclicLR(optimizer, base_lr=0.01, max_lr=0.1)\n# scheduler = torch.optim.lr_scheduler.ConstantLR(optimizer, factor=0.5, total_iters=4)\nscheduler1 = torch.optim.lr_scheduler.ConstantLR(optimizer, factor=0.1, total_iters=2)\nscheduler2 = torch.optim.lr_scheduler.ExponentialLR(optimizer, gamma=0.9)\nscheduler3 = torch.optim.lr_scheduler.StepLR(optimizer, step_size=30, gamma=0.1)\nlmbda = lambda epoch: 0.95\nscheduler4 = torch.optim.lr_scheduler.MultiplicativeLR(optimizer, lr_lambda=lmbda)\nscheduler5 = torch.optim.lr_scheduler.LinearLR(optimizer, start_factor=0.5, total_iters=4)\nscheduler6 = torch.optim.lr_scheduler.CosineAnnealingWarmRestarts(optimizer, T_0=30, T_mult=1)\nINITIAL_LEARNING_RATE = 0.01\nmin_lr = 0.0001\nlambda1 = lambda epoch: max(0.99 ** epoch, min_lr / INITIAL_LEARNING_RATE)\nscheduler7 = torch.optim.lr_scheduler.LambdaLR(optimizer, lr_lambda=[lambda1])\nscheduler8 = torch.optim.lr_scheduler.PolynomialLR(optimizer, total_iters=4, power=1.0)\nscheduler9 = optim.lr_scheduler.CyclicLR(optimizer,base_lr=0.001,max_lr=0.1,mode='triangular2',cycle_momentum=False)\n# scheduler10 = torch.optim.lr_scheduler.OneCycleLR(optimizer, max_lr=0.01, steps_per_epoch=128, epochs=10,cycle_momentum=False)\nscheduler = torch.optim.lr_scheduler.ChainedScheduler([scheduler1, scheduler2,scheduler3,scheduler4,scheduler5,scheduler6,scheduler7,scheduler8,scheduler9])\nmodel.train()\n# running_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-01T05:39:40.588428Z","iopub.execute_input":"2023-06-01T05:39:40.588832Z","iopub.status.idle":"2023-06-01T05:44:07.747332Z","shell.execute_reply.started":"2023-06-01T05:39:40.588794Z","shell.execute_reply":"2023-06-01T05:44:07.746249Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- 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\n\n- Model has never seen the label data within the rectangle before!\n- Plot it side-by-side with the label image. Are you able to recognize the letter \"P\" in it?","metadata":{"_uuid":"8edfa121-b2ef-419e-acbe-6d430fb50133","_cell_guid":"3c6ad763-f47c-4aaa-8c12-cb95c2d28d74","trusted":true}},{"cell_type":"code","source":"eval_dataset = SubvolumeDataset(image_stack, label, pixels_inside_rect)\neval_loader = data.DataLoader(eval_dataset, batch_size=BATCH_SIZE, shuffle=False)\noutput = torch.zeros_like(label).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(pixels_inside_rect[i*BATCH_SIZE+j])] = value\n\nfig, (ax1, ax2) = plt.subplots(1, 2)\nax1.imshow(output.cpu(), cmap='gray')\nax2.imshow(label.cpu(), cmap='gray')\nplt.show()","metadata":{"_uuid":"68013338-2e52-4f0b-b15b-b18e071aa5da","_cell_guid":"56812c41-904d-4645-bddf-49b19fe2685d","collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2023-06-01T05:44:57.706413Z","iopub.execute_input":"2023-06-01T05:44:57.706815Z","iopub.status.idle":"2023-06-01T05:46:25.701940Z","shell.execute_reply.started":"2023-06-01T05:44:57.706772Z","shell.execute_reply":"2023-06-01T05:46:25.700890Z"},"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":{"_uuid":"15c0a510-b4d1-4e14-974e-cbe5a7ac6b8e","_cell_guid":"49b12c62-7143-4d27-be79-e20e8cd9f5fe","trusted":true}},{"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.cpu(), cmap='gray')\nplt.show()","metadata":{"_uuid":"3eb989fa-2823-49f7-a55f-1e178884e344","_cell_guid":"4f1a3ed4-43c8-4a9c-9f49-dab4d3047048","collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2023-06-01T05:48:00.406362Z","iopub.execute_input":"2023-06-01T05:48:00.407191Z","iopub.status.idle":"2023-06-01T05:48:03.078677Z","shell.execute_reply.started":"2023-06-01T05:48:00.407146Z","shell.execute_reply":"2023-06-01T05:48:03.077486Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Finally, Kaggle expects a runlength-encoded submission.csv file, so let's output that.","metadata":{"_uuid":"04dc6e9a-5178-4ffa-95b3-a64783d4cf1a","_cell_guid":"df73edb8-b09d-4e84-b33a-b384dbe486fe","trusted":true}},{"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)\nrle_output = rle(output)\n# This doesn't make too much sense, but let's just output in the required format\n# so notebook works as a submission. :-)\nprint(\"Id,Predicted\\na,\" + rle_output + \"\\nb,\" + rle_output, file=open('submission.csv', 'w'))","metadata":{"_uuid":"512cd6ab-2794-4ad3-87cc-fc240561f286","_cell_guid":"79b4b495-2e93-49ba-b9d8-434dccf49907","collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2023-06-01T05:48:35.027338Z","iopub.execute_input":"2023-06-01T05:48:35.027726Z","iopub.status.idle":"2023-06-01T05:48:35.493459Z","shell.execute_reply.started":"2023-06-01T05:48:35.027691Z","shell.execute_reply":"2023-06-01T05:48:35.492353Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Hurray! We've detected ink! Now, can you do better? :-) For example, you could start with this [example submission](https://www.kaggle.com/code/danielhavir/vesuvius-challenge-example-submission).","metadata":{"_uuid":"e84d0aa9-a297-4a90-b4f2-a8afb7c389c5","_cell_guid":"799ddbd9-6862-4a43-9ae4-3c2d2f01da85","trusted":true}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}