{"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 tensorflow as tf\nfrom tensorflow import keras\nfrom tensorflow.keras import layers\nimport numpy as np\n\nimport time\nimport glob\nimport PIL.Image as Image\nimport matplotlib.pyplot as plt\nimport matplotlib.patches as patches\nfrom tqdm import tqdm","metadata":{"execution":{"iopub.status.busy":"2023-04-09T10:18:55.514554Z","iopub.execute_input":"2023-04-09T10:18:55.514915Z","iopub.status.idle":"2023-04-09T10:19:00.981949Z","shell.execute_reply.started":"2023-04-09T10:18:55.514878Z","shell.execute_reply":"2023-04-09T10:19:00.980846Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Data config\nDATA_DIR = '/kaggle/input/vesuvius-challenge-ink-detection/'\nBUFFER = 32  \nZ_DIM = 16 \nZ_START = 25","metadata":{"execution":{"iopub.status.busy":"2023-04-09T10:19:00.983683Z","iopub.execute_input":"2023-04-09T10:19:00.984492Z","iopub.status.idle":"2023-04-09T10:19:00.991278Z","shell.execute_reply.started":"2023-04-09T10:19:00.984449Z","shell.execute_reply":"2023-04-09T10:19:00.990164Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Model config\nBATCH_SIZE = 64\nplt.imshow(Image.open(DATA_DIR + \"/train/1/ir.png\"), cmap=\"gray\")","metadata":{"execution":{"iopub.status.busy":"2023-04-09T10:19:00.992950Z","iopub.execute_input":"2023-04-09T10:19:00.993463Z","iopub.status.idle":"2023-04-09T10:19:04.225961Z","shell.execute_reply.started":"2023-04-09T10:19:00.993420Z","shell.execute_reply":"2023-04-09T10:19:04.224855Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Loadig the dataset","metadata":{}},{"cell_type":"code","source":"def load_mask(split, index):\n    img = Image.open(f\"{DATA_DIR}/{split}/{index}/mask.png\").convert('1')\n    return tf.convert_to_tensor(img, dtype=\"bool\")\n\ndef load_labels(split, index):\n    img = Image.open(f\"{DATA_DIR}/{split}/{index}/inklabels.png\")\n    return tf.convert_to_tensor(img, dtype=\"bool\")\n\nmask = load_mask(split=\"train\", index=1)\nlabels = load_labels(split=\"train\", index=1)\n\nfig, (ax1, ax2) = plt.subplots(1, 2)\nax1.set_title(\"mask.png\")\nax1.imshow(mask, cmap='gray')\nax2.set_title(\"inklabels.png\")\nax2.imshow(labels, cmap='gray')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-04-09T10:19:04.229672Z","iopub.execute_input":"2023-04-09T10:19:04.230445Z","iopub.status.idle":"2023-04-09T10:19:09.259928Z","shell.execute_reply.started":"2023-04-09T10:19:04.230401Z","shell.execute_reply":"2023-04-09T10:19:09.258621Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def load_volume(split, index):\n    # Load the 3d x-ray scan, one slice at a time\n    z_slices_fnames = sorted(glob.glob(f\"{DATA_DIR}/{split}/{index}/surface_volume/*.tif\"))[Z_START:Z_START + Z_DIM]\n    z_slices = []\n    for z, filename in  tqdm(enumerate(z_slices_fnames)):\n        img = Image.open(filename)\n        z_slice = np.array(img, dtype=np.int16)\n        z_slices.append(z_slice)\n\n    return tf.stack(z_slices, axis=-1)\n\nvolume = load_volume(split=\"train\", index=1)\n\nprint(f\"volume {volume.shape} - {volume.dtype}\")","metadata":{"execution":{"iopub.status.busy":"2023-04-09T10:19:09.261768Z","iopub.execute_input":"2023-04-09T10:19:09.262242Z","iopub.status.idle":"2023-04-09T10:19:19.844049Z","shell.execute_reply.started":"2023-04-09T10:19:09.262199Z","shell.execute_reply":"2023-04-09T10:19:19.842803Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, axes = plt.subplots(1, 5, figsize=(15, 3))\nfor z, ax in enumerate(axes):\n    ax.imshow(volume[:, :, z], cmap='gray')\n    ax.set_xticks([]); ax.set_yticks([])\nfig.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-04-09T10:19:19.845758Z","iopub.execute_input":"2023-04-09T10:19:19.846150Z","iopub.status.idle":"2023-04-09T10:19:26.838692Z","shell.execute_reply.started":"2023-04-09T10:19:19.846091Z","shell.execute_reply":"2023-04-09T10:19:26.837188Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"val_location = (3500, 1100)\nval_zone_size = (800, 800)\n\nfig, ax = plt.subplots()\nax.imshow(labels)\npatch = patches.Rectangle([val_location[1], val_location[0]], val_zone_size[0], val_zone_size[1], linewidth=2, edgecolor='g', facecolor='none')\nax.add_patch(patch)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-04-09T10:19:26.840228Z","iopub.execute_input":"2023-04-09T10:19:26.841321Z","iopub.status.idle":"2023-04-09T10:19:28.296207Z","shell.execute_reply.started":"2023-04-09T10:19:26.841282Z","shell.execute_reply":"2023-04-09T10:19:28.295077Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Create a dataset that randoml;y samplesw coordinates(Locations)","metadata":{}},{"cell_type":"code","source":"def sample_random_location(shape):\n    random_train_x = tf.random.uniform(shape=(), minval=BUFFER, maxval=shape[0] - BUFFER - 1, dtype=\"int32\")\n    random_train_y = tf.random.uniform(shape=(), minval=BUFFER, maxval=shape[1] - BUFFER - 1, dtype=\"int32\")\n    random_train_location = tf.stack([random_train_x, random_train_y])\n    return random_train_location\n\ndef is_in_masked_zone(location, mask):\n    return mask[location[0], location[1]]\n\nsample_random_location_train = lambda x: sample_random_location(mask.shape)\nis_in_mask_train = lambda x: is_in_masked_zone(x, mask)\n\ndef is_in_val_zone(location, val_location, val_zone_size):\n    x = location[0]\n    y = location[1]\n    x_match = val_location[0] - BUFFER <= x <= val_location[0] + val_zone_size[0] + BUFFER\n    y_match = val_location[1] - BUFFER <= y <= val_location[1] + val_zone_size[1] + BUFFER\n    return x_match and y_match\n\ndef is_proper_train_location(location):\n    return not is_in_val_zone(location, val_location, val_zone_size) and is_in_mask_train(location)\n","metadata":{"execution":{"iopub.status.busy":"2023-04-09T10:19:28.297622Z","iopub.execute_input":"2023-04-09T10:19:28.298356Z","iopub.status.idle":"2023-04-09T10:19:28.308079Z","shell.execute_reply.started":"2023-04-09T10:19:28.298313Z","shell.execute_reply":"2023-04-09T10:19:28.307010Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_locations_ds = tf.data.Dataset.from_tensor_slices([0]).repeat().map(sample_random_location_train, num_parallel_calls=tf.data.AUTOTUNE)\ntrain_locations_ds = train_locations_ds.filter(is_proper_train_location)","metadata":{"execution":{"iopub.status.busy":"2023-04-09T10:19:28.309537Z","iopub.execute_input":"2023-04-09T10:19:28.311325Z","iopub.status.idle":"2023-04-09T10:19:28.581563Z","shell.execute_reply.started":"2023-04-09T10:19:28.311286Z","shell.execute_reply":"2023-04-09T10:19:28.580552Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Visualize training patch location ","metadata":{}},{"cell_type":"code","source":"fig, ax = plt.subplots()\nax.imshow(labels)\n\nfor x, y in train_locations_ds.take(200):\n    patch = patches.Rectangle([y - BUFFER // 2, x - BUFFER // 2], 2 * BUFFER, 2 * BUFFER, linewidth=2, edgecolor='r', facecolor='none')\n    ax.add_patch(patch)\n\nval_patch = patches.Rectangle([val_location[1], val_location[0]], val_zone_size[0], val_zone_size[1], linewidth=2, edgecolor='g', facecolor='none')\nax.add_patch(val_patch)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-04-09T10:19:28.582918Z","iopub.execute_input":"2023-04-09T10:19:28.583283Z","iopub.status.idle":"2023-04-09T10:19:31.265088Z","shell.execute_reply.started":"2023-04-09T10:19:28.583246Z","shell.execute_reply":"2023-04-09T10:19:31.264085Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Creating traing dataset that yields volume the patches +  lables","metadata":{}},{"cell_type":"code","source":"def extract_subvolume(location, volume):\n    x = location[0]\n    y = location[1]\n    subvolume = volume[x-BUFFER:x+BUFFER, y-BUFFER:y+BUFFER, :]\n    subvolume = tf.cast(subvolume, dtype=\"float32\") / 65535.\n    return subvolume\n\ndef extract_labels(location, labels):\n    x = location[0]\n    y = location[1]\n    label = labels[x, y]\n    label = tf.cast(label, dtype=\"float32\")\n    return label\n\ndef extract_subvolume_and_label(location):\n    subvolume = extract_subvolume(location, volume)\n    label = extract_labels(location, labels)\n    return subvolume, label\n\nshuffle_buffer_size = BATCH_SIZE * 4\n\ntrain_ds = train_locations_ds.map(extract_subvolume_and_label, num_parallel_calls=tf.data.AUTOTUNE)\ntrain_ds = train_ds.prefetch(tf.data.AUTOTUNE).batch(BATCH_SIZE)","metadata":{"execution":{"iopub.status.busy":"2023-04-09T10:19:31.266889Z","iopub.execute_input":"2023-04-09T10:19:31.267667Z","iopub.status.idle":"2023-04-09T10:19:32.746014Z","shell.execute_reply.started":"2023-04-09T10:19:31.267624Z","shell.execute_reply":"2023-04-09T10:19:32.744897Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for subvolume_batch, label_batch in train_ds.take(1):\n    print(f\"subvolume shape: {subvolume_batch.shape[1:]}\")","metadata":{"execution":{"iopub.status.busy":"2023-04-09T10:19:32.747361Z","iopub.execute_input":"2023-04-09T10:19:32.747777Z","iopub.status.idle":"2023-04-09T10:19:34.161449Z","shell.execute_reply.started":"2023-04-09T10:19:32.747737Z","shell.execute_reply":"2023-04-09T10:19:34.150054Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Create a validation dataset that yields patches from the validation area and their labels","metadata":{}},{"cell_type":"code","source":"val_locations_stride = 3\nval_locations = []\nfor x in range(val_location[0], val_location[0] + val_zone_size[0], val_locations_stride):\n    for y in range(val_location[1], val_location[1] + val_zone_size[1], val_locations_stride):\n        val_locations.append((x, y))\n\nval_ds = tf.data.Dataset.from_tensor_slices(val_locations).map(extract_subvolume_and_label, num_parallel_calls=tf.data.AUTOTUNE)\nval_ds = val_ds.prefetch(tf.data.AUTOTUNE).batch(BATCH_SIZE)","metadata":{"execution":{"iopub.status.busy":"2023-04-09T10:19:34.166423Z","iopub.execute_input":"2023-04-09T10:19:34.166838Z","iopub.status.idle":"2023-04-09T10:19:34.578657Z","shell.execute_reply.started":"2023-04-09T10:19:34.166798Z","shell.execute_reply":"2023-04-09T10:19:34.577618Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"t0 = time.time()\nn = 200\nfor _ in train_ds.take(n):\n    pass\nprint(f\"Time per batch: {(time.time() - t0) / n:.4f}s\")","metadata":{"execution":{"iopub.status.busy":"2023-04-09T10:19:34.581024Z","iopub.execute_input":"2023-04-09T10:19:34.581899Z","iopub.status.idle":"2023-04-09T10:19:40.468106Z","shell.execute_reply.started":"2023-04-09T10:19:34.581854Z","shell.execute_reply":"2023-04-09T10:19:40.465906Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def trivial_baseline(dataset):\n    total = 0\n    matches = 0.\n    for _, batch_label in tqdm(dataset):\n        matches += tf.reduce_sum(tf.cast(batch_label, \"float32\"))\n        total += tf.reduce_prod(tf.shape(batch_label))\n    return 1. - matches / tf.cast(total, \"float32\")\n\nscore = trivial_baseline(val_ds).numpy()\nprint(f\"Best validation score achievable trivially: {score * 100:.2f}% accuracy\")","metadata":{"execution":{"iopub.status.busy":"2023-04-09T10:19:40.469690Z","iopub.execute_input":"2023-04-09T10:19:40.470972Z","iopub.status.idle":"2023-04-09T10:19:58.619066Z","shell.execute_reply.started":"2023-04-09T10:19:40.470925Z","shell.execute_reply":"2023-04-09T10:19:58.617806Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"augmenter = keras.Sequential([\n    layers.RandomFlip(\"horizontal\"),\n])\n\ndef augment_train_data(data, label):\n    data = augmenter(data)\n    return data, label\n\naugmented_train_ds = train_ds.map(augment_train_data, num_parallel_calls=tf.data.AUTOTUNE).prefetch(tf.data.AUTOTUNE)","metadata":{"execution":{"iopub.status.busy":"2023-04-09T10:19:58.620575Z","iopub.execute_input":"2023-04-09T10:19:58.621690Z","iopub.status.idle":"2023-04-09T10:19:58.740413Z","shell.execute_reply.started":"2023-04-09T10:19:58.621635Z","shell.execute_reply":"2023-04-09T10:19:58.739422Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_model(input_shape):\n    inputs = keras.Input(input_shape)\n    x = layers.Conv2D(kernel_size=3, filters=128, activation=\"relu\")(inputs)\n\n    for size in (128, 256, 512):\n        residual = x\n        x = layers.BatchNormalization()(x)\n        x = layers.SeparableConv2D(kernel_size=3, filters=size, activation=\"relu\", padding=\"same\")(x)\n        x = layers.SeparableConv2D(kernel_size=3, filters=size, activation=\"relu\", padding=\"same\")(x)\n        x = layers.MaxPooling2D(2, padding=\"same\")(x)\n        residual = layers.Conv2D(kernel_size=1, filters=size, strides=2, padding=\"same\")(residual)\n        x = layers.Add()([residual, x])\n\n    x = layers.GlobalAveragePooling2D()(x)\n    x = layers.Dropout(0.3)(x)\n    outputs = layers.Dense(1, activation=\"sigmoid\")(x)\n    model = keras.Model(inputs, outputs)\n    return model\n\nstrategy = tf.distribute.MirroredStrategy()\nwith strategy.scope():\n    model = get_model(input_shape=(BUFFER * 2, BUFFER * 2, Z_DIM))\n    model.compile(optimizer=\"adam\", loss=\"binary_crossentropy\", metrics=[\"accuracy\"])","metadata":{"execution":{"iopub.status.busy":"2023-04-09T10:19:58.743763Z","iopub.execute_input":"2023-04-09T10:19:58.744060Z","iopub.status.idle":"2023-04-09T10:19:59.471180Z","shell.execute_reply.started":"2023-04-09T10:19:58.744031Z","shell.execute_reply":"2023-04-09T10:19:59.470149Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model.summary()","metadata":{"execution":{"iopub.status.busy":"2023-04-09T10:19:59.473189Z","iopub.execute_input":"2023-04-09T10:19:59.473961Z","iopub.status.idle":"2023-04-09T10:19:59.541417Z","shell.execute_reply.started":"2023-04-09T10:19:59.473914Z","shell.execute_reply":"2023-04-09T10:19:59.540649Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Because train_ds is infinitely streaming, you need to specify steps_per_epoch\n# i.e. the dataset has no epoch boundary on its own\nmodel.fit(augmented_train_ds,\n          validation_data=val_ds,\n          epochs=10,\n          steps_per_epoch=300,)","metadata":{"execution":{"iopub.status.busy":"2023-04-09T10:19:59.542446Z","iopub.execute_input":"2023-04-09T10:19:59.542792Z","iopub.status.idle":"2023-04-09T10:32:32.367536Z","shell.execute_reply.started":"2023-04-09T10:19:59.542755Z","shell.execute_reply":"2023-04-09T10:32:32.366391Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model.save(\"model.keras\")","metadata":{"execution":{"iopub.status.busy":"2023-04-09T10:32:32.369399Z","iopub.execute_input":"2023-04-09T10:32:32.369807Z","iopub.status.idle":"2023-04-09T10:32:32.523269Z","shell.execute_reply.started":"2023-04-09T10:32:32.369769Z","shell.execute_reply":"2023-04-09T10:32:32.522168Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"del volume\ndel mask\ndel labels\ndel train_ds\ndel val_ds\n\n# Manually trigger garbage collection\nkeras.backend.clear_session()\nimport gc\ngc.collect()\n\nmodel = keras.models.load_model(\"model.keras\")","metadata":{"execution":{"iopub.status.busy":"2023-04-09T10:32:32.524962Z","iopub.execute_input":"2023-04-09T10:32:32.525345Z","iopub.status.idle":"2023-04-09T10:32:33.263692Z","shell.execute_reply.started":"2023-04-09T10:32:32.525307Z","shell.execute_reply":"2023-04-09T10:32:33.262527Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"INFER_BATCH_SIZE = 128\n\ndef compute_predictions_map(split, index):\n    print(f\"Load data for {split}/{index}\")\n\n    test_volume = load_volume(split=split, index=index)\n    test_mask = load_mask(split=split, index=index)\n\n    test_locations = []\n    stride = BUFFER // 2\n    for x in range(BUFFER, test_volume.shape[0] - BUFFER, stride):\n        for y in range(BUFFER, test_volume.shape[1] - BUFFER, stride):\n            test_locations.append((x, y))\n\n    print(f\"{len(test_locations)} test locations (before filtering by mask)\")\n\n    sample_random_location_test = lambda x: sample_random_location(test_mask.shape)\n    is_in_mask_test = lambda x: is_in_masked_zone(x, test_mask)\n    extract_subvolume_test = lambda x: extract_subvolume(x, test_volume)\n\n    test_locations_ds = tf.data.Dataset.from_tensor_slices(test_locations).filter(is_in_mask_test)\n    test_ds = test_locations_ds.map(extract_subvolume_test, num_parallel_calls=tf.data.AUTOTUNE)\n\n    predictions_map = np.zeros(test_volume.shape[:2], dtype=\"float16\")\n\n    print(f\"Compute predictions\")\n\n    for loc_batch, patch_batch in tqdm(zip(test_locations_ds.batch(INFER_BATCH_SIZE), test_ds.batch(INFER_BATCH_SIZE))):\n        predictions = model.predict_on_batch(patch_batch)\n        for (x, y), pred in zip(loc_batch, predictions):\n            predictions_map[x, y] = pred\n            \n    del test_volume\n    del test_mask\n    return predictions_map","metadata":{"execution":{"iopub.status.busy":"2023-04-09T11:01:15.321174Z","iopub.execute_input":"2023-04-09T11:01:15.322181Z","iopub.status.idle":"2023-04-09T11:01:15.332584Z","shell.execute_reply.started":"2023-04-09T11:01:15.322088Z","shell.execute_reply":"2023-04-09T11:01:15.331193Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"predictions_map_1 = compute_predictions_map(split=\"train\", index=1)","metadata":{"execution":{"iopub.status.busy":"2023-04-09T11:02:10.020413Z","iopub.execute_input":"2023-04-09T11:02:10.021243Z","iopub.status.idle":"2023-04-09T11:02:17.335320Z","shell.execute_reply.started":"2023-04-09T11:02:10.021201Z","shell.execute_reply":"2023-04-09T11:02:17.334112Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"predictions_map_a = compute_predictions_map(split=\"test\", index=\"a\")\npredictions_map_b = compute_predictions_map(split=\"test\", index=\"b\")","metadata":{"execution":{"iopub.status.busy":"2023-04-09T11:02:35.430312Z","iopub.execute_input":"2023-04-09T11:02:35.431489Z","iopub.status.idle":"2023-04-09T11:03:02.826389Z","shell.execute_reply.started":"2023-04-09T11:02:35.431439Z","shell.execute_reply":"2023-04-09T11:03:02.825178Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def stats(predictions_map):\n    size = predictions_map.shape[0] * predictions_map.shape[1]\n    for thr in (0.1, 0.25, 0.5, 0.75):\n        num_nonzero = np.count_nonzero(predictions_map_1 > thr)\n        percent_nonzero = 100. * num_nonzero / size\n        print(f\"Percent above {thr}: {percent_nonzero:.2f}%\")","metadata":{"execution":{"iopub.status.busy":"2023-04-09T11:10:28.387344Z","iopub.execute_input":"2023-04-09T11:10:28.387769Z","iopub.status.idle":"2023-04-09T11:10:28.394568Z","shell.execute_reply.started":"2023-04-09T11:10:28.387735Z","shell.execute_reply":"2023-04-09T11:10:28.392672Z"},"trusted":true},"execution_count":null,"outputs":[]}]}