{"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":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport cv2\nimport matplotlib.pyplot as plt\n\nimport torch\nimport torch.nn as nn\nimport torch.optim as optim\nimport glob\nimport PIL.Image as Image\nimport torch.utils.data as data\nimport matplotlib.patches as patches\nfrom tqdm import tqdm\nfrom ipywidgets import interact, fixed\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 = 300\nLEARNING_RATE = 0.03\nBATCH_SIZE = 32\nDEVICE = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-output":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-18T19:30:14.369593Z","iopub.execute_input":"2023-03-18T19:30:14.370730Z","iopub.status.idle":"2023-03-18T19:30:15.666039Z","shell.execute_reply.started":"2023-03-18T19:30:14.370687Z","shell.execute_reply":"2023-03-18T19:30:15.664740Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## **Goal of the Competition**\n\nJoin the 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. There is a $150,000 grand prize available to the first team that can read these scrolls from a 3d x-ray scan.\n\n\n## **Dataset Description**\n\nYour challenge is to recover where ink is present from 3d x-ray scans of detached fragments of ancient papyrus scrolls. This is an important subproblem in the overall task of solving the Vesuvius Challenge.\n\nThis is a Code Competition. When your submitted notebook is scored, the actual test data will be made available to your notebook.\n\n**Files**\n\n**[train/test]/[fragment_id]/surface_volume/[image_id].tif** slices from the 3d x-ray surface volume. Each file contains a greyscale slice in the z-direction. Each fragment contains 65 slices. Combined this image stack gives us width * height * 65 number of voxels per fragment. You can expect two fragments in the hidden test set, which together are roughly the same size as a single training fragment. The sample slices available to download in the test folders are simply copied from training fragment one.\n\n## **Analyze Data**\n\nWe will first define a variable and display an image using Matplotlib\n\nThe **PREFIX** variable as the file path prefix for the data files. It assumes that the data is located in the Kaggle input directory for the \"Vesuvius Challenge Ink Detection\" competition, specifically in the \"train/1/\" subdirectory.\n\nThis line displays an image using Matplotlib. It opens the file \"ir.png\" in the \"PREFIX\" directory using the Image.open() function from the PIL library and displays it using plt.imshow(). The cmap=\"gray\" argument specifies that the image should be displayed in grayscale. The semicolon at the end of the line suppresses the output of the imshow() function, so it doesn't display any extra information in the notebook.","metadata":{}},{"cell_type":"code","source":"PREFIX = '/kaggle/input/vesuvius-challenge-ink-detection/train/1/'\nplt.imshow(Image.open(PREFIX+\"ir.png\"), cmap=\"gray\");","metadata":{"execution":{"iopub.status.busy":"2023-03-18T19:30:15.671153Z","iopub.execute_input":"2023-03-18T19:30:15.671479Z","iopub.status.idle":"2023-03-18T19:30:19.356454Z","shell.execute_reply.started":"2023-03-18T19:30:15.671444Z","shell.execute_reply":"2023-03-18T19:30:19.355519Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's load these binary images:\n\n**mask.png:** a mask of which pixels contain data, and which pixels we should ignore.\n\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 = 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":{"execution":{"iopub.status.busy":"2023-03-18T19:30:19.357560Z","iopub.execute_input":"2023-03-18T19:30:19.358565Z","iopub.status.idle":"2023-03-18T19:30:22.747553Z","shell.execute_reply.started":"2023-03-18T19:30:19.358528Z","shell.execute_reply":"2023-03-18T19:30:22.745611Z"},"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\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 = [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":{"execution":{"iopub.status.busy":"2023-03-18T19:30:22.751597Z","iopub.execute_input":"2023-03-18T19:30:22.752412Z","iopub.status.idle":"2023-03-18T19:30:33.775903Z","shell.execute_reply.started":"2023-03-18T19:30:22.752349Z","shell.execute_reply":"2023-03-18T19:30:33.774961Z"},"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":{"execution":{"iopub.status.busy":"2023-03-18T19:30:33.777282Z","iopub.execute_input":"2023-03-18T19:30:33.777993Z","iopub.status.idle":"2023-03-18T19:30:35.655908Z","shell.execute_reply.started":"2023-03-18T19:30:33.777954Z","shell.execute_reply":"2023-03-18T19:30:35.654763Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## **Pytorch Model**\n\nWe define a PyTorch dataset class called \"SubvolumeDataset\" that takes in an image stack, a label map, and a list of pixel coordinates as input. It is designed to be used with 3D image data, where the \"image_stack\" input is a 4D tensor of size (Z_DIM, Y_DIM, X_DIM, 1), representing a stack of 2D images along the Z-axis. The \"label\" input is a 2D tensor of size (Y_DIM, X_DIM) that specifies the ground-truth label for each pixel in the image stack. The \"pixels\" input is a list of 2-tuples, where each tuple contains the y and x coordinates of a pixel in the image stack.\n\nThe \"len\" method of the dataset class returns the number of pixels in the input \"pixels\" list, and the \"getitem\" method returns a tuple containing a subvolume tensor and its corresponding label tensor. The subvolume tensor is a 4D tensor of size (1, Z_DIM, BUFFER2+1, BUFFER2+1), where BUFFER is a constant value defined elsewhere in the code. The subvolume tensor is extracted from the image stack by taking a slice along the Z-axis centered around the pixel coordinates specified in the \"pixels\" list. The label tensor is a scalar tensor that specifies the ground-truth label for the pixel at the center of the subvolume.\n\nWe also build a PyTorch sequential neural network model that consists of several convolutional layers, dropout layers, max pooling layers, and linear layers. The model is designed to take in a 4D input tensor of size (1, Z_DIM, BUFFER2+1, BUFFER2+1) and output a scalar value between 0 and 1 using a sigmoid activation function. The model is initialized on a specified device (e.g. GPU) using the \".to()\" method.","metadata":{}},{"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        \n    def __len__(self):\n        return len(self.pixels)\n    \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.Dropout3d(0.2), nn.MaxPool3d(2, 2),\n    nn.Conv3d(16, 32, 3, 1, 1), nn.Dropout3d(0.2), nn.MaxPool3d(2, 2),\n    nn.Conv3d(32, 64, 3, 1, 1), nn.Dropout3d(0.2), nn.MaxPool3d(2, 2),\n    nn.Flatten(start_dim=1),\n    nn.LazyLinear(128), nn.ReLU(),\n    nn.Dropout(0.2),\n    nn.LazyLinear(1), nn.Sigmoid()\n).to(DEVICE)","metadata":{"execution":{"iopub.status.busy":"2023-03-18T19:36:12.063369Z","iopub.execute_input":"2023-03-18T19:36:12.063886Z","iopub.status.idle":"2023-03-18T19:36:12.078894Z","shell.execute_reply.started":"2023-03-18T19:36:12.063842Z","shell.execute_reply":"2023-03-18T19:36:12.077686Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We will now train the model.","metadata":{}},{"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.\npixels_inside_rect = []\npixels_outside_rect = []\nfor pixel in zip(*np.where(mask == 1)):\n    if pixel[1] < BUFFER or pixel[1] >= mask.shape[1]-BUFFER or pixel[0] < BUFFER or pixel[0] >= mask.shape[0]-BUFFER:\n        continue # Too close to the edge\n    if pixel[1] >= rect[0] and pixel[1] <= rect[0]+rect[2] and pixel[0] >= rect[1] and pixel[0] <= rect[1]+rect[3]:\n        pixels_inside_rect.append(pixel)\n    else:\n        pixels_outside_rect.append(pixel)\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":{"execution":{"iopub.status.busy":"2023-03-18T19:36:12.918857Z","iopub.execute_input":"2023-03-18T19:36:12.919757Z","iopub.status.idle":"2023-03-18T19:40:54.015554Z","shell.execute_reply.started":"2023-03-18T19:36:12.919707Z","shell.execute_reply":"2023-03-18T19:40:54.014048Z"},"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:","metadata":{}},{"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[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":{"execution":{"iopub.status.busy":"2023-03-18T19:40:54.018324Z","iopub.execute_input":"2023-03-18T19:40:54.019120Z"},"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.cpu(), cmap='gray')\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## **Submission**\n\nFinally, Kaggle expects a runlength-encoded submission.csv file, so let's output that.","metadata":{}},{"cell_type":"code","source":"def rle(output):\n    flat_img = np.where(output.flatten().cpu() > THRESHOLD, 1, 0).astype(np.uint8)\n    starts = np.array((flat_img[:-1] == 0) & (flat_img[1:] == 1))\n    ends = np.array((flat_img[:-1] == 1) & (flat_img[1:] == 0))\n    starts_ix = np.where(starts)[0] + 2\n    ends_ix = np.where(ends)[0] + 2\n    lengths = ends_ix - starts_ix\n    return \" \".join(map(str, sum(zip(starts_ix, lengths), ())))\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":{"trusted":true},"execution_count":null,"outputs":[]}]}