{"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":"**EDA & Visualization**","metadata":{}},{"cell_type":"code","source":"import seaborn as sns\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#from PIL import Image\nfrom pathlib import Path\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))","metadata":{"execution":{"iopub.status.busy":"2023-06-13T21:49:18.745625Z","iopub.execute_input":"2023-06-13T21:49:18.746362Z","iopub.status.idle":"2023-06-13T21:49:18.772218Z","shell.execute_reply.started":"2023-06-13T21:49:18.746312Z","shell.execute_reply":"2023-06-13T21:49:18.771176Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"N_SCAN_LAYERS = 65\n\n# Load the volume of the given papyrus fragment\ndef load_volume(set_name, folder_name):\n    path = f'/kaggle/input/vesuvius-challenge-ink-detection/{set_name}/{folder_name}/surface_volume'\n    file_format = \"{:02d}.tif\"\n    volume = []\n    \n    for i in tqdm(range(N_SCAN_LAYERS)):\n        file = f'{path}/{file_format.format(i)}'\n        img = cv2.imread(file, cv2.IMREAD_ANYDEPTH)\n        volume.append(img)\n        \n    return volume\n\n# Display the volume as a video\ndef display_volume(set_name, folder_name):\n    volume = load_volume(set_name, folder_name)\n    fig, ax = plt.subplots()\n    camera = Camera(fig)\n    \n    for i in range(N_SCAN_LAYERS):\n        ax.axis('off')\n        ax.text(0.5, 1.08, f\"Layer {i+1}/{N_SCAN_LAYERS}\", fontweight='bold', fontsize=18,\n                transform=ax.transAxes, horizontalalignment='center')\n        ax.imshow(volume[0], cmap='gray')\n        camera.snap()\n        del volume[0]\n        gc.collect()\n        \n    plt.close(fig)\n    \n    animation = camera.animate()\n    fix_video_adjust = \\\n    '<style> video {margin: 0px; padding: 0px; width:100%; height:auto;} </style>'\n    display(HTML(fix_video_adjust + animation.to_html5_video()))\n    \n    del camera\n    del animation\n    gc.collect()","metadata":{"execution":{"iopub.status.busy":"2023-06-13T21:17:26.390965Z","iopub.execute_input":"2023-06-13T21:17:26.391422Z","iopub.status.idle":"2023-06-13T21:17:26.402414Z","shell.execute_reply.started":"2023-06-13T21:17:26.391387Z","shell.execute_reply":"2023-06-13T21:17:26.401426Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_shape(set_name, folder_name):\n    file = f'/kaggle/input/vesuvius-challenge-ink-detection/{set_name}/{folder_name}/mask.png'\n    img = cv2.imread(file, cv2.IMREAD_ANYDEPTH)\n    return img.shape\n\ndef get_shape_info():\n    samples = []\n    shapes = []\n    px_nbs = []\n    fragments = {'train': ['1', '2', '3'], 'test': ['a', 'b']}\n    for set_name in fragments:\n        for folder_name in fragments[set_name]:\n            samples.append(f'{set_name} {folder_name}')\n            shape = get_shape(set_name, folder_name)\n            shapes.append(shape)\n            px_nbs.append(shape[0] * shape[1])\n    df = pd.DataFrame([shapes, px_nbs], columns=samples, index=['Shape', 'Nb of pixels']).T\n    \n    return df\n\ndf = get_shape_info()\n\nplt.figure(figsize=(10, 4))\nplt.subplot(1, 2, 1)\nsns.barplot(x=df.index, y=df['Nb of pixels'], )\nplt.title(\"Size in pixels for each scan\")\nplt.subplot(1, 2, 2)\nplt.pie(df['Nb of pixels'], labels=df.index, autopct='%0.1f%%')\nplt.title(\"Pixel size as percentage\")\nplt.show()\n\ndf","metadata":{"execution":{"iopub.status.busy":"2023-06-13T21:17:26.524668Z","iopub.execute_input":"2023-06-13T21:17:26.525855Z","iopub.status.idle":"2023-06-13T21:17:28.409851Z","shell.execute_reply.started":"2023-06-13T21:17:26.525809Z","shell.execute_reply":"2023-06-13T21:17:28.408960Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"TRAIN_PATH = Path(\"/kaggle/input/vesuvius-challenge-ink-detection/train\")","metadata":{"execution":{"iopub.status.busy":"2023-06-13T21:45:42.324404Z","iopub.execute_input":"2023-06-13T21:45:42.324778Z","iopub.status.idle":"2023-06-13T21:45:42.330196Z","shell.execute_reply.started":"2023-06-13T21:45:42.324747Z","shell.execute_reply":"2023-06-13T21:45:42.329050Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(15,6))\n\nimg = Image.open(TRAIN_PATH / \"1\" / \"surface_volume\" / \"00.tif\")\nprint(f\"Surface scanned image size: ({img.size[0]}, {img.size[1]})\")\n\nplt.subplot(1,3,1)\nplt.title(\"First layer\")\nplt.imshow(img)\n\nimg_labels=Image.open(TRAIN_PATH / \"1/inklabels.png\")\nprint(f\"Ink Labels Mask image size: ({img.size[0]}, {img.size[1]})\")\nplt.subplot(1,3,2)\nplt.imshow(img_labels)\nplt.title(\"Ink Labels Mask\")\n\nimg_ir=Image.open(TRAIN_PATH / \"1/ir.png\")\nprint(f\"IR image size: ({img.size[0]}, {img.size[1]})\")\nplt.subplot(1, 3, 3)\nplt.imshow(img_ir)\nplt.title(\"IR Image\")","metadata":{"execution":{"iopub.status.busy":"2023-06-13T21:45:44.386493Z","iopub.execute_input":"2023-06-13T21:45:44.387368Z","iopub.status.idle":"2023-06-13T21:45:58.343714Z","shell.execute_reply.started":"2023-06-13T21:45:44.387313Z","shell.execute_reply":"2023-06-13T21:45:58.342595Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, axs = plt.subplots(nrows=4, ncols=3, figsize=(12,16))\n\nfor i, ax in tqdm(enumerate(axs.flat)):\n    img = Image.open(TRAIN_PATH / \"1\" / \"surface_volume\" / f\"{(i*5):02d}.tif\")\n    ax.imshow(img)\n    ax.set_title(f\"{(i*5):02d}.tif\")\n    \nfig.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-13T21:17:40.449094Z","iopub.execute_input":"2023-06-13T21:17:40.449680Z","iopub.status.idle":"2023-06-13T21:18:01.800067Z","shell.execute_reply.started":"2023-06-13T21:17:40.449647Z","shell.execute_reply":"2023-06-13T21:18:01.799219Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**library**","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-ink-detection/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":{"execution":{"iopub.status.busy":"2023-06-13T21:18:01.802825Z","iopub.execute_input":"2023-06-13T21:18:01.803427Z","iopub.status.idle":"2023-06-13T21:18:04.629019Z","shell.execute_reply.started":"2023-06-13T21:18:01.803393Z","shell.execute_reply":"2023-06-13T21:18:04.628186Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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-06-13T21:18:04.630590Z","iopub.execute_input":"2023-06-13T21:18:04.631197Z","iopub.status.idle":"2023-06-13T21:18:08.310108Z","shell.execute_reply.started":"2023-06-13T21:18:04.631158Z","shell.execute_reply":"2023-06-13T21:18:08.309181Z"},"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)\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-06-13T21:18:08.311414Z","iopub.execute_input":"2023-06-13T21:18:08.312209Z","iopub.status.idle":"2023-06-13T21:18:32.033361Z","shell.execute_reply.started":"2023-06-13T21:18:08.312173Z","shell.execute_reply":"2023-06-13T21:18:32.032446Z"},"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-06-13T21:18:32.035667Z","iopub.execute_input":"2023-06-13T21:18:32.037235Z","iopub.status.idle":"2023-06-13T21:18:33.918710Z","shell.execute_reply.started":"2023-06-13T21:18:32.037197Z","shell.execute_reply":"2023-06-13T21:18:33.917841Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Dataset Class**","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    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)\n","metadata":{"execution":{"iopub.status.busy":"2023-06-13T21:18:33.920212Z","iopub.execute_input":"2023-06-13T21:18:33.920845Z","iopub.status.idle":"2023-06-13T21:18:33.942005Z","shell.execute_reply.started":"2023-06-13T21:18:33.920800Z","shell.execute_reply":"2023-06-13T21:18:33.940706Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"Generating pixel lists...\")\n\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\ninside_rect[rect[1]:rect[1]+rect[3]+1, rect[0]:rect[0]+rect[2]+1] = True\n\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","metadata":{"execution":{"iopub.status.busy":"2023-06-13T21:18:33.943745Z","iopub.execute_input":"2023-06-13T21:18:33.944124Z","iopub.status.idle":"2023-06-13T21:26:12.381864Z","shell.execute_reply.started":"2023-06-13T21:18:33.944091Z","shell.execute_reply":"2023-06-13T21:26:12.380723Z"},"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":{"execution":{"iopub.status.busy":"2023-06-13T21:26:12.385341Z","iopub.execute_input":"2023-06-13T21:26:12.385791Z","iopub.status.idle":"2023-06-13T21:28:08.408067Z","shell.execute_reply.started":"2023-06-13T21:26:12.385754Z","shell.execute_reply":"2023-06-13T21:28:08.407165Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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":{"execution":{"iopub.status.busy":"2023-06-13T21:28:08.409556Z","iopub.execute_input":"2023-06-13T21:28:08.410559Z","iopub.status.idle":"2023-06-13T21:28:11.074036Z","shell.execute_reply.started":"2023-06-13T21:28:08.410524Z","shell.execute_reply":"2023-06-13T21:28:11.073145Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"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)\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":{"execution":{"iopub.status.busy":"2023-06-13T21:28:11.075670Z","iopub.execute_input":"2023-06-13T21:28:11.076360Z","iopub.status.idle":"2023-06-13T21:28:11.549431Z","shell.execute_reply.started":"2023-06-13T21:28:11.076323Z","shell.execute_reply":"2023-06-13T21:28:11.548450Z"},"trusted":true},"execution_count":null,"outputs":[]}]}