{"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 libraries\nimport 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","metadata":{"execution":{"iopub.status.busy":"2023-04-20T11:20:39.476157Z","iopub.execute_input":"2023-04-20T11:20:39.476474Z","iopub.status.idle":"2023-04-20T11:20:42.208231Z","shell.execute_reply.started":"2023-04-20T11:20:39.476444Z","shell.execute_reply":"2023-04-20T11:20:42.207072Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Goal:\n\nAccurately predict the presence of ink in a carbonized fragment of a vesuvius scroll using x-ray scans of the fragment as input to the model.\nGround truth to use for training are inklabels retrieved from IR scans of the same fragment where ink is visible on the IR-scan.\n\nIn the end this model should be used on scrolls where IR-scanning is not possible/feasible.","metadata":{}},{"cell_type":"markdown","source":"### Load mask and ink labels","metadata":{}},{"cell_type":"code","source":"# define initital hyperparameters\nPREFIX = '/kaggle/input/vesuvius-challenge-ink-detection/train/1/'\nBUFFER = 70  # Buffer size in x and y direction\nZ_START = 27 # First slice in the z direction to use\nZ_DIM = 11  # Number of slices in the z direction\nTRAINING_STEPS = 30000\nLEARNING_RATE = 0.01\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":{"execution":{"iopub.status.busy":"2023-04-20T11:20:46.411854Z","iopub.execute_input":"2023-04-20T11:20:46.412962Z","iopub.status.idle":"2023-04-20T11:20:50.171674Z","shell.execute_reply.started":"2023-04-20T11:20:46.412912Z","shell.execute_reply":"2023-04-20T11:20:50.170672Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# plot mask and y-label images\nmask = np.array(Image.open(PREFIX+\"mask.png\").convert('1')) # .convert(1) means 1-bit black and white pixel\nlabel = torch.from_numpy(np.array(Image.open(PREFIX+\"inklabels.png\"))).gt(0).float().to(DEVICE) #.gt(0).float() checks if pixel value is greater than 0, gives True if greater than zero, False if not and converts to float (1. or 0.) \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-04-20T11:20:55.167703Z","iopub.execute_input":"2023-04-20T11:20:55.168158Z","iopub.status.idle":"2023-04-20T11:21:02.399242Z","shell.execute_reply.started":"2023-04-20T11:20:55.168115Z","shell.execute_reply":"2023-04-20T11:21:02.398266Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Display x-ray images","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":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-20T11:21:05.758640Z","iopub.execute_input":"2023-04-20T11:21:05.759600Z","iopub.status.idle":"2023-04-20T11:21:28.238767Z","shell.execute_reply.started":"2023-04-20T11:21:05.759542Z","shell.execute_reply":"2023-04-20T11:21:28.237603Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Define dataset and model","metadata":{}},{"cell_type":"code","source":"# define Dataset\nclass 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# define initial model\nmodel = nn.Sequential(\n    nn.Conv3d(1, 16, 3, 1, padding=\"same\"),\n    nn.BatchNorm3d(16),\n    nn.ReLU(),\n    nn.MaxPool3d((1,2,2), (1,2,2)),\n    nn.Conv3d(16, 32, 3, 1, padding=\"same\"),\n    nn.BatchNorm3d(32),\n    nn.ReLU(),\n    nn.MaxPool3d(2, 2),\n    nn.Conv3d(32, 64, 3, 1, padding=\"same\"),\n    nn.BatchNorm3d(64),\n    nn.ReLU(),\n    nn.MaxPool3d(2, 2),\n    nn.Conv3d(64, 256, 1, 1, padding=\"same\"),\n    nn.BatchNorm3d(256),\n    nn.ReLU(),\n    nn.Conv3d(256, 512, 3, 1, padding=\"same\"),\n    nn.BatchNorm3d(512),\n    nn.ReLU(),\n    nn.MaxPool3d(2, 2),\n    nn.Flatten(start_dim=1),\n    nn.Linear(32768, 256),\n    nn.ReLU(),\n    nn.Linear(256, 128),\n    nn.ReLU(),\n    nn.Linear(128, 1),\n    nn.Sigmoid()\n).to(DEVICE)","metadata":{"execution":{"iopub.status.busy":"2023-04-12T14:36:46.484454Z","iopub.execute_input":"2023-04-12T14:36:46.485209Z","iopub.status.idle":"2023-04-12T14:36:46.619089Z","shell.execute_reply.started":"2023-04-12T14:36:46.485166Z","shell.execute_reply":"2023-04-12T14:36:46.618053Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Training","metadata":{}},{"cell_type":"code","source":"# to speed up training we initially predict on a small rectangle which we define below\nrect = (1100, 3500, 2500, 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-04-12T14:36:46.851481Z","iopub.execute_input":"2023-04-12T14:36:46.851835Z","iopub.status.idle":"2023-04-12T14:36:48.806963Z","shell.execute_reply.started":"2023-04-12T14:36:46.851802Z","shell.execute_reply":"2023-04-12T14:36:48.805873Z"},"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...\")\n# load train dataset\ntrain_dataset = SubvolumeDataset(image_stack, label, pixels_outside_rect)\n# load dataset in DataLoader\ntrain_loader = data.DataLoader(train_dataset, batch_size=BATCH_SIZE, shuffle=True)\n# create Binary Cross Entropy Loss function\ncriterion = nn.BCELoss()\n# SGD optimizer\noptimizer = optim.SGD(model.parameters(), lr=LEARNING_RATE)\n# Create scheduler\nscheduler = torch.optim.lr_scheduler.OneCycleLR(optimizer, max_lr=LEARNING_RATE, total_steps=TRAINING_STEPS)\nmodel.train()\nrunning_loss = 0.0\n# run epochs and print losses\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-04-12T14:36:52.295703Z","iopub.execute_input":"2023-04-12T14:36:52.296086Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Evaluation","metadata":{}},{"cell_type":"code","source":"# evaluation of model by predicting the unseen x-ray data in the rectangle defined above\neval_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()\n# we don't want to calculate gradients\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# show ink prediction\nfig, (ax1, ax2) = plt.subplots(1, 2)\nax1.imshow(output.cpu(), cmap='gray')\nax2.imshow(label.cpu(), cmap='gray')\nplt.show()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Save modeling","metadata":{}},{"cell_type":"code","source":"# Pickling model\nimport pickle\n# open a file, where you ant to store the data\nfile = open('/kaggle/working/model.pickle', 'wb')\n\n# dump information to that file\npickle.dump(model, file)\n\n# close the file\nfile.close()","metadata":{"execution":{"iopub.status.busy":"2023-04-12T16:22:25.130234Z","iopub.execute_input":"2023-04-12T16:22:25.131270Z","iopub.status.idle":"2023-04-12T16:22:25.142845Z","shell.execute_reply.started":"2023-04-12T16:22:25.131224Z","shell.execute_reply":"2023-04-12T16:22:25.141436Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# define threshold to predict as ink\nTHRESHOLD = 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":{"execution":{"iopub.status.busy":"2023-04-12T13:23:06.612678Z","iopub.status.idle":"2023-04-12T13:23:06.613203Z","shell.execute_reply.started":"2023-04-12T13:23:06.612935Z","shell.execute_reply":"2023-04-12T13:23:06.612963Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Save submission file as .csv","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\n# def 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# rle_output = rle(output)\n\n# print(\"Id,Predicted\\na,\" + rle_output + \"\\nb,\" + rle_output, file=open('submission.csv', 'w'))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}