{"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":"# Vesuvius EDA ([kaggle notebook](https://www.kaggle.com/code/thenoodleninja/exploratory-data-analysis))\n\nIn this notebook some basic EDA and preprocessing will be performed on the [vesuvius-challenge-ink-detection](https://www.kaggle.com/competitions/vesuvius-challenge-ink-detection) dataset.\n\nThis dataset consists of three Herculaneum roll fragments of differing sizes. If one looks into the folder of one of those fragments, one finds the following files:\n - **ir.png**: An IR image of the fragment on which the ink is clearly visible. This was used to generate the ink mask.\n - **inklabels.png**: The target mask, specifying the location of ink on the papyrus.\nfig, ax = plt.subplots(3, 2)\nfor i, fragment_path in enumerate(train_fragments):\n    plot_cross_section(ax[i, 0], ScanData(fragment_path))\n    ax[i, 0].set_title(f\"fragment {fragment_path.name}, y=4000, x=1500:1900\")\n    plot_cross_section(ax[i, 1], ScanData(fragment_path), y=slice(3800, 4200), x=1700)\n    ax[i, 1].set_title(f\"fragment {fragment_path.name}, y=4000, x=1500:1900\")\nfig.set_figwidth(20)\nfig.tight_layout()\nplt.show() \n - **mask.png**: A binary mask, seperating the papyrus from the background.\n - **surface_volume/\\<i\\>.tif**: The i-th layer of the x-ray volume. ","metadata":{}},{"cell_type":"code","source":"# imports\nimport cv2\nimport numpy as np\nimport os\nimport gc\nimport glob\nimport json\nfrom pathlib import Path\nfrom tqdm import tqdm\nfrom scipy import ndimage\nimport matplotlib.pyplot as plt\n\nfrom PIL import Image, ImageDraw, ImageFont\n\nMAX_UINT16 = int(2**16-1)","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-06-22T14:38:36.997676Z","iopub.execute_input":"2023-06-22T14:38:36.998207Z","iopub.status.idle":"2023-06-22T14:38:37.314184Z","shell.execute_reply.started":"2023-06-22T14:38:36.998143Z","shell.execute_reply":"2023-06-22T14:38:37.312615Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# get fragment paths\nbase_path = Path(\"/kaggle/input/vesuvius-challenge-ink-detection\")\ntrain_path = base_path / \"train\"\ntrain_fragments = sorted([train_path  / f.name for f in train_path.iterdir()])\ntest_path = base_path / \"test\"\ntest_fragments = sorted([test_path  / f.name for f in test_path.iterdir()])\n\nallFragments = train_fragments + test_fragments","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-06-22T14:38:37.316667Z","iopub.execute_input":"2023-06-22T14:38:37.317210Z","iopub.status.idle":"2023-06-22T14:38:37.340902Z","shell.execute_reply.started":"2023-06-22T14:38:37.317156Z","shell.execute_reply":"2023-06-22T14:38:37.339576Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Show individual fragments with corresponding masks and labels\n\nfig, axs = plt.subplots(3, 3)\nfig.set_figwidth(12)\nfig.set_figheight(10)\nfor i, fragment_path in enumerate(train_fragments):\n    mask = cv2.imread(str(fragment_path / \"mask.png\"))\n    ir = cv2.imread(str(fragment_path / \"ir.png\"))\n    label = cv2.imread(str(fragment_path / \"inklabels.png\"))\n    axs[0, i].imshow(mask)\n    axs[0, i].set_title(f\"Fragment {fragment_path.name}\")\n    axs[1, i].imshow(ir)\n    axs[2, i].imshow(label)\nfig.suptitle(\"Fragment masks, IR images and ink labels\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-16T11:46:54.678729Z","iopub.execute_input":"2023-06-16T11:46:54.679185Z","iopub.status.idle":"2023-06-16T11:48:17.056535Z","shell.execute_reply.started":"2023-06-16T11:46:54.679145Z","shell.execute_reply":"2023-06-16T11:48:17.055352Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Create gifs for individual fragment that loop through the 65 layers \n\noutput_folder = Path(\"/kaggle/working/visualizations\")\nfont = ImageFont.load_default()\n\nos.makedirs(output_folder, exist_ok=True)\nfor fragment_id in range(1, 4):\n    \n    volume_folder = base_path / 'train' / str(fragment_id) / \"surface_volume\"\n    frames = []\n    # load and resize layers\n    for i, image in tqdm(enumerate(sorted(glob.glob(str(volume_folder / \"*.tif\"))))):\n        img = Image.open(image)\n        width, height = img.size\n        img = Image.fromarray((np.array(img.resize((width//8, height//8)))//255).astype(np.uint8))\n        d = ImageDraw.Draw(img)\n        d.text( (30, 30), f\"{i:02}\", fill=(255))\n        frames.append(img)\n\n    # save individual layers as a looping gif\n    frames[0].save(output_folder / f\"x-ray_{fragment_id}.gif\", format=\"GIF\", append_images=frames[1:],\n               save_all=True, duration=100, loop=0)","metadata":{"execution":{"iopub.status.busy":"2023-06-16T11:48:17.058277Z","iopub.execute_input":"2023-06-16T11:48:17.059671Z","iopub.status.idle":"2023-06-16T11:59:46.932117Z","shell.execute_reply.started":"2023-06-16T11:48:17.059613Z","shell.execute_reply":"2023-06-16T11:59:46.930430Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def mem_efficient_index(X, mask):\n    masked_layers = []\n    for x in X:\n        masked_layers.append(x[mask])\n    return np.concatenate(masked_layers)\n\ndef flatten(arr, z_buffer=2, z_layers=7, blur_topo = True):\n    \"\"\"\n    :param arr: numpy array with the surface_volume_data\n    :param z_buffer: how much air we will leave above the papyrus\n    :param z_layers: how much layers we want to keep after transsforming\n    :return:\n    \"\"\"\n    arr = arr.astype(float) / (2**16-1) # convert to float\n    arr = np.flip(arr, axis=0)\n    arr = gaussian_filter(arr, sigma=1)\n    arr_sob = ndimage.sobel(arr, axis=0)\n    arr_sob = gaussian_filter(arr, sigma=1)\n    topo = np.argmax(np.where(arr_sob < 0.5, 0, 1), axis=0)\n    del arr_sob\n    if blur_topo:\n        topo=gaussian_filter(topo,sigma=1)\n    arr_idx = np.indices(arr.shape)\n    arr = arr[\n        (arr_idx[0] + topo - z_buffer) % arr.shape[0], arr_idx[1], arr_idx[2]]\n    arr = (arr[0:z_layers]*(2**16-1)).astype(np.uint16)\n    return arr, topo.astype(np.uint8)\n\nclass ScanData:\n    def __init__(self, baseDir, cache_folder=\"/tmp\"):\n        self.id = hash(str(baseDir))\n        self.dir = baseDir.resolve()\n        maskName = str(baseDir / \"mask.png\")\n        mask = np.array(Image.open(maskName).convert(\"1\"))        \n        self.mask = mask\n        (self.h,self.w) = mask.shape\n\n        labelName = str(baseDir / \"inklabels.png\" )\n        try:\n            self.label = np.array(Image.open(labelName).convert(\"1\"))        \n        except:\n            self.label = None\n\n        self.sliceNames = sorted( (baseDir / \"surface_volume\").rglob(\"*.tif\") )\n        dim = len(self.sliceNames)\n        self.dim = dim\n        self.cache_folder=cache_folder\n    \n    def get_img_uint16(self, z_slc: slice=slice(None, None), y_slc:slice=slice(None, None), x_slc:slice=slice(None, None)):\n        if not os.path.isfile(f\"{self.cache_folder}/{self.id}_uint16.npy\"):\n            print(\"caching img\")\n            os.makedirs(f\"{self.cache_folder}\", exist_ok= True)\n            \n            img = np.zeros((self.dim, self.h, self.w), dtype=np.uint16)\n            for idx, filename in enumerate(tqdm(self.sliceNames, leave=False)):\n                fname = str(filename)\n                img[idx, :, :] = np.array(Image.open(fname))\n            np.save(f\"{self.cache_folder}/{self.id}_uint16.npy\", img)\n            \n        return np.load(f'{self.cache_folder}/{self.id}_uint16.npy', mmap_mode='r')[z_slc, y_slc, x_slc].copy()\n    \n    def get_img_uint8(self, z_slc: slice=slice(None, None), y_slc:slice=slice(None, None), x_slc:slice=slice(None, None)):\n        if not os.path.isfile(f\"{self.cache_folder}/{self.id}_uint8.npy\"):\n            print(\"caching img\")\n            os.makedirs(f\"{self.cache_folder}\", exist_ok= True)\n\n            img = np.zeros((self.dim, self.h, self.w), dtype=np.uint8)\n            for idx, filename in enumerate(tqdm(self.sliceNames, leave=False)):\n                fname = str(filename)\n                img[idx, :, :] = (np.array(Image.open(fname))//256).astype(np.uint8)\n            np.save(f\"{self.cache_folder}/{self.id}_uint8.npy\", img)\n\n        return np.load(f'{self.cache_folder}/{self.id}_uint8.npy', mmap_mode='r')[z_slc, y_slc, x_slc].copy()\n    \n    def flatten(self, z_layers=7, stripe_width=500, stripe_overlay = 20):\n        flattened = np.zeros((z_layers, self.h, self.w), dtype=np.uint16)\n        topo = np.zeros((self.h, self.w), dtype=np.uint8)\n        for i in tqdm(range(0,self.h, stripe_width)):\n            y0=max(i-stripe_overlay,0)\n            y1=min(i+stripe_width+stripe_overlay,self.h)\n\n            img_stripe = self.get_img_uint16(y0=y0, y1=y1)\n            stripe_flat, stripe_topo = flatten(img_stripe, 2, z_layers)\n            if y0 == 0:\n                flattened[:, i:min(i+stripe_width, self.h), :] = stripe_flat[:, 0:min(stripe_width, stripe_flat.shape[1]), :]\n                topo[i:min(i+stripe_width, self.h), :] = stripe_topo[0:min(stripe_width, stripe_topo.shape[0]), :]\n            else:\n                flattened[:, i:min(i+stripe_width, self.h), :] = stripe_flat[:, stripe_overlay:min(stripe_width+stripe_overlay, stripe_flat.shape[1]), :] \n                topo[i:min(i+stripe_width, self.h), :] = stripe_topo[stripe_overlay:min(stripe_width+stripe_overlay, stripe_topo.shape[0]), :] \n        return flattened, topo\n    \n    def save_flattened(self, folder):\n        os.makedirs(f\"{folder}/surface_volume\", exist_ok= True)\n        \n        flattened, topo = self.flatten()\n        # saving layers\n        for i in tqdm(range(flattened.shape[0])):\n            layer = Image.fromarray(flattened[i])\n            layer.save(f\"{folder}/surface_volume/{i:02}.tif\")\n        layer = Image.fromarray(topo)\n        layer.save(f\"{folder}/topo.png\")","metadata":{"execution":{"iopub.status.busy":"2023-06-22T14:38:59.665114Z","iopub.execute_input":"2023-06-22T14:38:59.665549Z","iopub.status.idle":"2023-06-22T14:38:59.708768Z","shell.execute_reply.started":"2023-06-22T14:38:59.665508Z","shell.execute_reply":"2023-06-22T14:38:59.707548Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# visualize intensity distributions for ink (red), papyrus (blue) and everything (grey)\nfig, axs = plt.subplots(1, 3)\nfig.set_figwidth(20)\nfor i, fragment_path in enumerate(train_fragments):\n    # read data\n    scan = ScanData(fragment_path)\n    full = scan.get_img_uint8()      \n    label_name = fragment_path / \"inklabels.png\"\n    label = np.array(Image.open(label_name), dtype=bool)\n\n    # plot papyrus and ink hist\n    axs[i].hist(mem_efficient_index(full, scan.mask), bins=256, alpha=0.3, color=\"gray\")\n    \n    ink = mem_efficient_index(full, label)\n    ink_mean = ink.mean()\n    axs[i].hist(ink, bins=256, alpha=0.3, color=\"red\")\n    del ink\n    gc.collect()\n    \n    pap = mem_efficient_index(full, scan.mask & np.logical_not(label))\n    pap_mean = pap.mean()\n    axs[i].hist(pap, bins=256, alpha=0.3, color=\"blue\")\n    del pap\n    gc.collect()\n    \n    axs[i].set_title(f\"ink_mean:{ink_mean:.2f}, pap_mean:{pap_mean:.2f}\")\n    del scan, full, label\n    gc.collect()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-16T11:59:46.988212Z","iopub.execute_input":"2023-06-16T11:59:46.989777Z","iopub.status.idle":"2023-06-16T12:17:11.395613Z","shell.execute_reply.started":"2023-06-16T11:59:46.989703Z","shell.execute_reply":"2023-06-16T12:17:11.394344Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"N = 5\n\nfragment_data = []\nfor i, fragment_path in enumerate(train_fragments):\n    scan = ScanData(fragment_path)\n    label_name = fragment_path / \"inklabels.png\"\n    label = np.array(Image.open(label_name), dtype=bool)\n    \n    paps, inks = [], []\n    for n in tqdm(range(65-N)):\n        z0, z1 = n, n+N\n        slc = scan.get_img_uint8(z_slc=slice(z0, z1))\n        \n        paps.append(mem_efficient_index(slc, scan.mask & np.logical_not(label)).mean())\n        inks.append(mem_efficient_index(slc, label).mean())\n    fragment_data.append((paps, inks))","metadata":{"execution":{"iopub.status.busy":"2023-06-16T12:17:11.399595Z","iopub.execute_input":"2023-06-16T12:17:11.400532Z","iopub.status.idle":"2023-06-16T12:26:52.164496Z","shell.execute_reply.started":"2023-06-16T12:17:11.400480Z","shell.execute_reply":"2023-06-16T12:26:52.162607Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, axs = plt.subplots(2, 3)\nfor i, fragment_path in enumerate(train_fragments):\n    inks, paps = fragment_data[i]\n    ax0, ax1 = axs[:, i]\n    ink, = ax0.plot(inks, color=\"red\", label=\"ink\")\n    pap, = ax0.plot(paps, color = \"blue\", label=\"papyrus\")\n    ax0.set_title(\"mean intensities\")\n    ax0.set_xlabel(\"sliding window starting height\")\n    ax0.legend(handles=[ink, pap])\n    ax1.plot(abs(np.array(inks)-np.array(paps)), color=\"black\", label=\"|ink-pap|\")\n    ax1.set_title(\"absolute diff mean intensity\")\n    ax1.set_xlabel(\"sliding window starting height\")\n    gc.collect()\nfig.suptitle(f\"mean intensities for sliding window of size {N}\")\nfig.set_figwidth(20)\nfig.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-16T12:26:52.166552Z","iopub.execute_input":"2023-06-16T12:26:52.167013Z","iopub.status.idle":"2023-06-16T12:26:53.926641Z","shell.execute_reply.started":"2023-06-16T12:26:52.166967Z","shell.execute_reply":"2023-06-16T12:26:53.925066Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from skimage.filters import gaussian\n\n#visualize cross_sections\ndef plot_cross_section(ax, scan, y=4000, x=slice(1500, 1900), denoise=False):\n    slc = scan.get_img_uint8(x_slc=x, y_slc=y)\n    if denoise:\n        slc = gaussian(slc)\n    ax.imshow(slc, cmap='gray', origin='lower')\n\nfig, ax = plt.subplots(3, 2)\nfor i, fragment_path in enumerate(train_fragments):\n    plot_cross_section(ax[i, 0], ScanData(fragment_path))\n    ax[i, 0].set_title(f\"fragment {fragment_path.name}, y=4000, x=1500:1900\")\n    plot_cross_section(ax[i, 1], ScanData(fragment_path), y=slice(3800, 4200), x=1700)\n    ax[i, 1].set_title(f\"fragment {fragment_path.name}, y=3800:4200, x=1700\")\nfig.set_figwidth(20)\nfig.tight_layout()\nplt.savefig(\"/kaggle/working/fig.png\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-22T14:58:35.679781Z","iopub.execute_input":"2023-06-22T14:58:35.680246Z","iopub.status.idle":"2023-06-22T14:58:45.675209Z","shell.execute_reply.started":"2023-06-22T14:58:35.680205Z","shell.execute_reply":"2023-06-22T14:58:45.673823Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# visiualize denoised\nfig, ax = plt.subplots(3, 2)\nfor i, fragment_path in enumerate(train_fragments):\n    plot_cross_section(ax[i, 0], ScanData(fragment_path), denoise=True)\n    ax[i, 0].set_title(f\"fragment {fragment_path.name}, y=4000, x=1500:1900\")\n    plot_cross_section(ax[i, 1], ScanData(fragment_path), y=slice(3800, 4200), x=1700, denoise=True)\n    ax[i, 1].set_title(f\"fragment {fragment_path.name}, y=3800:4200, x=1700\")\nfig.set_figwidth(20)\nfig.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-16T12:27:00.364831Z","iopub.execute_input":"2023-06-16T12:27:00.366428Z","iopub.status.idle":"2023-06-16T12:27:06.328199Z","shell.execute_reply.started":"2023-06-16T12:27:00.366360Z","shell.execute_reply":"2023-06-16T12:27:06.326474Z"},"trusted":true},"execution_count":null,"outputs":[]}]}