{"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":"# Fragment flattening ([kaggle notebook](https://www.kaggle.com/code/thenoodleninja/fragment-flattening))\n\nScript to flatten the [vesuvius scrolls](https://www.kaggle.com/competitions/vesuvius-challenge-ink-detection).","metadata":{}},{"cell_type":"code","source":"import gc\nfrom scipy.ndimage import gaussian_filter\nfrom scipy import ndimage\nfrom pathlib import Path\nimport numpy as np\nimport glob\nimport PIL.Image as Image\nimport matplotlib.pyplot as plt\nimport seaborn as sn\nfrom tqdm.auto import tqdm \nimport os","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def 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_sob, 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=\"/kaggle/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, x_slc=slice(None, None), y_slc=slice(None, None), z_slc=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, x_slc=slice(None, None), y_slc=slice(None, None), z_slc=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(y_slc=slice(y0, 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-23T07:58:31.637645Z","iopub.execute_input":"2023-06-23T07:58:31.638203Z","iopub.status.idle":"2023-06-23T07:58:31.682011Z","shell.execute_reply.started":"2023-06-23T07:58:31.638148Z","shell.execute_reply":"2023-06-23T07:58:31.680737Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# get fragment paths\n\noutput_folder = Path(\"/kaggle/working/flat\")\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\n\nfor fragment in allFragments:\n    scan = ScanData(fragment)\n    print(fragment.name)\n    scan.save_flattened( output_folder / fragment.name)","metadata":{"execution":{"iopub.status.busy":"2023-06-22T14:47:07.871803Z","iopub.execute_input":"2023-06-22T14:47:07.873252Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Comparison of segmentation approaches","metadata":{}},{"cell_type":"code","source":"output_folder = Path(\"/kaggle/working/flat\")\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()])\nallFragments = train_fragments + test_fragments\nscan = ScanData(allFragments[3])","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:58:45.509344Z","iopub.execute_input":"2023-06-23T07:58:45.509773Z","iopub.status.idle":"2023-06-23T07:58:45.685327Z","shell.execute_reply.started":"2023-06-23T07:58:45.509735Z","shell.execute_reply":"2023-06-23T07:58:45.684242Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Clustering\n\nIn this section several clustering algorithms are compared. All of these use the scaled intensity and z-position of pixels features.","metadata":{}},{"cell_type":"code","source":"slc = scan.get_img_uint8()\ntest = gaussian_filter(slc[:, 1000, 2400:6000], sigma=1)\ntest_z = np.stack((test, np.indices(test.shape)[0]), axis=2).astype(float)\ntest_z[:, :, 0]/=255.0*0.5\ntest_z[:, :, 1]/=64.0","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:38:17.580703Z","iopub.execute_input":"2023-06-23T07:38:17.581880Z","iopub.status.idle":"2023-06-23T07:38:18.413818Z","shell.execute_reply.started":"2023-06-23T07:38:17.581824Z","shell.execute_reply":"2023-06-23T07:38:18.411373Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.scatter(x=test_z[:, :, 0].reshape(-1, 1).squeeze(), y=test_z[:, :, 1].reshape(-1, 1).squeeze(), marker=\".\", alpha=0.005)","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:38:29.300684Z","iopub.execute_input":"2023-06-23T07:38:29.301554Z","iopub.status.idle":"2023-06-23T07:38:30.014665Z","shell.execute_reply.started":"2023-06-23T07:38:29.301511Z","shell.execute_reply":"2023-06-23T07:38:30.013435Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.cluster import KMeans\nN=3\nkmeans = KMeans(n_clusters=N).fit(test_z.reshape((-1, 2)))\na = kmeans.predict(test_z.reshape(-1, 2)).reshape(test.shape)\nplt.imshow(test[:, 1500:2000])\nplt.show()\nplt.imshow(a[:, 1500:2000])\nplt.show()\nsn.scatterplot(x=test_z[:, : ,0].flatten(), y=test_z[:, :, 1].flatten(), hue=a.flatten(), alpha=0.005)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:41:27.897293Z","iopub.execute_input":"2023-06-23T07:41:27.897693Z","iopub.status.idle":"2023-06-23T07:41:42.243643Z","shell.execute_reply.started":"2023-06-23T07:41:27.897658Z","shell.execute_reply":"2023-06-23T07:41:42.242345Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.cluster import AgglomerativeClustering\nlinkage=\"ward\"\nwrd = AgglomerativeClustering(linkage=linkage, n_clusters = N).fit(test_z[:, 1500:2000, :].reshape((-1, 2)))\n#a = wrd.predict(test_z.reshape(-1, 2)).reshape(test.shape)\nplt.imshow(test[:, 1500:2000])\nplt.show()\nplt.imshow(wrd.labels_.reshape((65, -1)))\nplt.show()\nsn.scatterplot(x=test_z[:, 1500:2000, 0].flatten(), y=test_z[:, 1500:2000, 1].flatten(), hue=wrd.labels_, alpha=0.05)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:41:42.246185Z","iopub.execute_input":"2023-06-23T07:41:42.246600Z","iopub.status.idle":"2023-06-23T07:42:06.718947Z","shell.execute_reply.started":"2023-06-23T07:41:42.246557Z","shell.execute_reply":"2023-06-23T07:42:06.717555Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.cluster import AgglomerativeClustering\nlinkage=\"average\"\nwrd = AgglomerativeClustering(linkage=linkage, n_clusters = N).fit(test_z[:, 1500:2000, :].reshape((-1, 2)))\n#a = wrd.predict(test_z.reshape(-1, 2)).reshape(test.shape)\nplt.imshow(test[:, 1500:2000])\nplt.show()\nplt.imshow(wrd.labels_.reshape((65, -1)))\nplt.show()\nsn.scatterplot(x=test_z[:, 1500:2000, 0].flatten(), y=test_z[:, 1500:2000, 1].flatten(), hue=wrd.labels_, alpha=0.05)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:42:06.720660Z","iopub.execute_input":"2023-06-23T07:42:06.721482Z","iopub.status.idle":"2023-06-23T07:42:27.073113Z","shell.execute_reply.started":"2023-06-23T07:42:06.721430Z","shell.execute_reply":"2023-06-23T07:42:27.071833Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.cluster import AgglomerativeClustering\nlinkage=\"complete\"\nwrd = AgglomerativeClustering(linkage=linkage, n_clusters = N).fit(test_z[:, 1500:2000, :].reshape((-1, 2)))\n#a = wrd.predict(test_z.reshape(-1, 2)).reshape(test.shape)\nplt.imshow(test[:, 1500:2000])\nplt.show()\nplt.imshow(wrd.labels_.reshape((65, -1)))\nplt.show()\nsn.scatterplot(x=test_z[:, 1500:2000, 0].flatten(), y=test_z[:, 1500:2000, 1].flatten(), hue=wrd.labels_, alpha=0.05)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:42:27.075191Z","iopub.execute_input":"2023-06-23T07:42:27.075532Z","iopub.status.idle":"2023-06-23T07:42:47.002991Z","shell.execute_reply.started":"2023-06-23T07:42:27.075499Z","shell.execute_reply":"2023-06-23T07:42:47.001948Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Edge Detection","metadata":{}},{"cell_type":"code","source":"arr = scan.get_img_uint16(x_slc=slice(2400, 6000), y_slc=slice(500, 1500))\narr = arr.astype(float) / (2**16-1) # convert to float\narr = np.flip(arr, axis=0)\narr_filtered = gaussian_filter(arr, sigma=1)\narr_sob = ndimage.sobel(arr_filtered, axis=0)\narr_sob = gaussian_filter(arr_sob, sigma=1)\ntopo = np.argmax(np.where(arr_sob < 0.5, 0, 1), axis=0)","metadata":{"execution":{"iopub.status.busy":"2023-06-23T08:11:08.931110Z","iopub.execute_input":"2023-06-23T08:11:08.931581Z","iopub.status.idle":"2023-06-23T08:11:35.400278Z","shell.execute_reply.started":"2023-06-23T08:11:08.931540Z","shell.execute_reply":"2023-06-23T08:11:35.398988Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, axs = plt.subplots(4, 1)\naxs[0].imshow(arr[:, 500, 1500:2000])\naxs[1].imshow(arr_sob[:, 500, 1500:2000])\naxs[2].imshow(np.where(arr_sob < 0.5, 0, 1)[:, 500, 1500:2000])\naxs[3].set_ylim(-64, 0)\naxs[3].plot(-topo[500, 1500:2000], )","metadata":{"execution":{"iopub.status.busy":"2023-06-23T08:13:28.947898Z","iopub.execute_input":"2023-06-23T08:13:28.948415Z","iopub.status.idle":"2023-06-23T08:13:30.467708Z","shell.execute_reply.started":"2023-06-23T08:13:28.948369Z","shell.execute_reply":"2023-06-23T08:13:30.465591Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Thresholding","metadata":{}},{"cell_type":"code","source":"from skimage.filters import threshold_otsu\n\nslc = scan.get_img_uint8()\ntest = gaussian_filter(slc[:, 1000, 2400:6000], sigma=1)\nthresh = threshold_otsu(test)","metadata":{"execution":{"iopub.status.busy":"2023-06-23T08:22:41.731749Z","iopub.execute_input":"2023-06-23T08:22:41.732247Z","iopub.status.idle":"2023-06-23T08:22:42.619990Z","shell.execute_reply.started":"2023-06-23T08:22:41.732207Z","shell.execute_reply":"2023-06-23T08:22:42.618642Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.hist(test.flatten(), bins=256, range=(0,256))\nplt.axvline(thresh, color='r')\nplt.show()\nthreshed = np.zeros_like(test)\nthreshed[test > thresh] = 1\nplt.imshow(threshed[:, 1500:2000])","metadata":{"execution":{"iopub.status.busy":"2023-06-23T08:24:49.199137Z","iopub.execute_input":"2023-06-23T08:24:49.199595Z","iopub.status.idle":"2023-06-23T08:24:50.106837Z","shell.execute_reply.started":"2023-06-23T08:24:49.199552Z","shell.execute_reply":"2023-06-23T08:24:50.105612Z"},"trusted":true},"execution_count":null,"outputs":[]}]}