{"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":"code","source":"import numpy as np\nimport glob\nimport PIL.Image as Image\nimport matplotlib.patches as patches\nfrom tqdm import tqdm\nfrom ipywidgets import interact, fixed\nimport torch\nimport torch.nn as nn\nimport torch.optim as optim\nimport torch.utils.data as data\nimport matplotlib.pyplot as plt\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 = 30000\nLEARNING_RATE = 0.03\nBATCH_SIZE = 32\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-05-09T08:29:35.590834Z","iopub.execute_input":"2023-05-09T08:29:35.591689Z","iopub.status.idle":"2023-05-09T08:29:42.248281Z","shell.execute_reply.started":"2023-05-09T08:29:35.591652Z","shell.execute_reply":"2023-05-09T08:29:42.247194Z"},"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(\"m.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-05-09T08:29:42.250810Z","iopub.execute_input":"2023-05-09T08:29:42.255067Z","iopub.status.idle":"2023-05-09T08:29:48.671063Z","shell.execute_reply.started":"2023-05-09T08:29:42.255022Z","shell.execute_reply":"2023-05-09T08:29:48.669934Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"images = [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-05-09T08:29:48.672895Z","iopub.execute_input":"2023-05-09T08:29:48.673561Z","iopub.status.idle":"2023-05-09T08:30:08.996083Z","shell.execute_reply.started":"2023-05-09T08:29:48.673520Z","shell.execute_reply":"2023-05-09T08:30:08.995068Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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-05-09T08:30:08.998695Z","iopub.execute_input":"2023-05-09T08:30:08.999694Z","iopub.status.idle":"2023-05-09T08:30:10.877087Z","shell.execute_reply.started":"2023-05-09T08:30:08.999657Z","shell.execute_reply":"2023-05-09T08:30:10.876084Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now we'll define a PyTorch dataset and (super simple) 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\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":{"_uuid":"9cf6abaa-b06b-4d05-a6d8-b7f012cf2e2c","_cell_guid":"25930dc1-1883-4fcd-9d40-f4aaa2ca152b","collapsed":false,"jupyter":{"outputs_hidden":false},"_kg_hide-output":true,"execution":{"iopub.status.busy":"2023-05-09T08:30:10.878841Z","iopub.execute_input":"2023-05-09T08:30:10.879562Z","iopub.status.idle":"2023-05-09T08:30:10.904479Z","shell.execute_reply.started":"2023-05-09T08:30:10.879522Z","shell.execute_reply":"2023-05-09T08:30:10.903201Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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()\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()\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":{"_uuid":"29fb8028-1ec5-42fb-b61a-658e1181c951","_cell_guid":"790108d2-5260-4276-9105-9da5a27c28f5","collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2023-05-09T08:30:10.906279Z","iopub.execute_input":"2023-05-09T08:30:10.906635Z","iopub.status.idle":"2023-05-09T08:38:38.623329Z","shell.execute_reply.started":"2023-05-09T08:30:10.906600Z","shell.execute_reply":"2023-05-09T08:38:38.622198Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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-05-09T08:38:38.625067Z","iopub.execute_input":"2023-05-09T08:38:38.625876Z","iopub.status.idle":"2023-05-09T08:40:05.638992Z","shell.execute_reply.started":"2023-05-09T08:38:38.625832Z","shell.execute_reply":"2023-05-09T08:40:05.637720Z"},"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-05-09T08:40:05.640561Z","iopub.execute_input":"2023-05-09T08:40:05.640985Z","iopub.status.idle":"2023-05-09T08:40:08.327691Z","shell.execute_reply.started":"2023-05-09T08:40:05.640949Z","shell.execute_reply":"2023-05-09T08:40:08.326445Z"},"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":"\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)\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-05-09T08:40:08.329530Z","iopub.execute_input":"2023-05-09T08:40:08.329903Z","iopub.status.idle":"2023-05-09T08:40:08.796741Z","shell.execute_reply.started":"2023-05-09T08:40:08.329867Z","shell.execute_reply":"2023-05-09T08:40:08.795707Z"},"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}}]}