{"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":"# Method to flatten the papyrus surface to have better training and testing data\nI looked at the fragments and noticed they are not completely flat.\nHere I want to show a way to flatten it more.\nMy hope is that this will lead to similar data on similar layers.\nIf similar data is on similar layers it should be easier to train a model on it.\nAlso, we can remove layers if most relevant data is on the same layer and by this reduce memory footprint.\n\nLet's have a look at the current data.\n","metadata":{}},{"cell_type":"code","source":"import gc\nfrom scipy.ndimage import gaussian_filter\nfrom scipy import ndimage\nimport numpy as np\nimport glob\nimport PIL.Image as Image\nimport matplotlib.pyplot as plt\nfrom tqdm import tqdm\nfrom os import path","metadata":{"collapsed":false,"ExecuteTime":{"start_time":"2023-04-12T07:51:04.574690Z","end_time":"2023-04-12T07:51:04.586807Z"},"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2023-04-12T06:50:02.669626Z","iopub.execute_input":"2023-04-12T06:50:02.670961Z","iopub.status.idle":"2023-04-12T06:50:02.677473Z","shell.execute_reply.started":"2023-04-12T06:50:02.670896Z","shell.execute_reply":"2023-04-12T06:50:02.675607Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"SAVE_IMAGE_STACK=True\nTRANSFORM_WHOLE_STACK=False\nPATH_TO_FRAGMENT = '/kaggle/input/vesuvius-challenge-ink-detection/train/1/surface_volume/'\nWORK_DIR='/kaggle/working/'\nGAUSSIAN_BLUR_TOPOGRAPHIC_MAP=False\n","metadata":{"collapsed":false,"ExecuteTime":{"start_time":"2023-04-12T07:51:04.831926Z","end_time":"2023-04-12T07:51:04.864574Z"},"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2023-04-12T06:50:03.090300Z","iopub.execute_input":"2023-04-12T06:50:03.090635Z","iopub.status.idle":"2023-04-12T06:50:03.096069Z","shell.execute_reply.started":"2023-04-12T06:50:03.090607Z","shell.execute_reply":"2023-04-12T06:50:03.095262Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Load images to numpy array. If you want you can save numpy array to load faster later.","metadata":{}},{"cell_type":"code","source":"image_stack=None\nif SAVE_IMAGE_STACK:\n    images = [np.array(Image.open(filename), dtype=np.float32)/65535.0 for filename in tqdm(sorted(glob.glob(PATH_TO_FRAGMENT+\"*.tif\")))]\n    image_stack = np.stack(images)\n    del images\n    gc.collect()\n    with open(path.join(WORK_DIR,\"image_stack.npy\"), 'wb') as f:\n\n        np.save(f, image_stack,allow_pickle=True)\nelse:\n    with open(path.join(WORK_DIR,\"image_stack.npy\"), 'rb') as f:\n        image_stack = np.load(f)\n","metadata":{"collapsed":false,"ExecuteTime":{"start_time":"2023-04-12T07:35:38.864644Z","end_time":"2023-04-12T07:35:53.817166Z"},"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2023-04-12T06:50:05.182112Z","iopub.execute_input":"2023-04-12T06:50:05.182466Z","iopub.status.idle":"2023-04-12T06:52:33.416955Z","shell.execute_reply.started":"2023-04-12T06:50:05.182438Z","shell.execute_reply":"2023-04-12T06:52:33.415827Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"image_stack=np.flip(image_stack,axis=0) #turn data upside down so papyrus is right side up","metadata":{"collapsed":false,"ExecuteTime":{"start_time":"2023-04-12T07:38:02.400370Z","end_time":"2023-04-12T07:38:02.432398Z"},"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2023-04-12T06:52:33.418884Z","iopub.execute_input":"2023-04-12T06:52:33.419222Z","iopub.status.idle":"2023-04-12T06:52:33.430705Z","shell.execute_reply.started":"2023-04-12T06:52:33.419188Z","shell.execute_reply":"2023-04-12T06:52:33.429324Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We have seen the surface volumes from the top. Let's have a look from the side.","metadata":{}},{"cell_type":"code","source":"slice_at=4000\nslice_length=slice(2000,2500)\nplt.imshow(image_stack[:,slice_length,slice_at], cmap='gray')\nplt.show()","metadata":{"collapsed":false,"ExecuteTime":{"start_time":"2023-04-12T07:51:13.419750Z","end_time":"2023-04-12T07:51:13.558961Z"},"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2023-04-12T06:52:44.799870Z","iopub.execute_input":"2023-04-12T06:52:44.800208Z","iopub.status.idle":"2023-04-12T06:52:45.045016Z","shell.execute_reply.started":"2023-04-12T06:52:44.800181Z","shell.execute_reply":"2023-04-12T06:52:45.043149Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The papyrus is the noisy stuff between layer  20 and 64. We can see that the surface of the papyrus is not completely flat. So the ink is not on the same layer all the time.\nLet's change that. Because my way takes a lot of memory we will first try on a smaller portion of the data. Then I will show a way to do it more memory friendly.","metadata":{}},{"cell_type":"code","source":"image_stack=image_stack[:,1000:3000,1000:3000] #smaller portion","metadata":{"collapsed":false,"ExecuteTime":{"start_time":"2023-04-12T07:51:14.914235Z","end_time":"2023-04-12T07:51:14.923330Z"},"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2023-04-12T06:52:47.670051Z","iopub.execute_input":"2023-04-12T06:52:47.670445Z","iopub.status.idle":"2023-04-12T06:52:47.676185Z","shell.execute_reply.started":"2023-04-12T06:52:47.670418Z","shell.execute_reply":"2023-04-12T06:52:47.674536Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"image_stack_=gaussian_filter(image_stack,sigma=1) #blur data a little bit","metadata":{"collapsed":false,"ExecuteTime":{"start_time":"2023-04-12T07:51:15.278338Z","end_time":"2023-04-12T07:51:59.675092Z"},"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2023-04-12T06:52:48.110384Z","iopub.execute_input":"2023-04-12T06:52:48.110732Z","iopub.status.idle":"2023-04-12T06:52:56.619732Z","shell.execute_reply.started":"2023-04-12T06:52:48.110704Z","shell.execute_reply":"2023-04-12T06:52:56.618383Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"image_stack_=ndimage.sobel(image_stack_,axis=0) # detect edges in top-down direction","metadata":{"collapsed":false,"ExecuteTime":{"start_time":"2023-04-12T07:51:59.677092Z","end_time":"2023-04-12T07:52:34.517604Z"},"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2023-04-12T06:52:56.622059Z","iopub.execute_input":"2023-04-12T06:52:56.622871Z","iopub.status.idle":"2023-04-12T06:53:02.083378Z","shell.execute_reply.started":"2023-04-12T06:52:56.622817Z","shell.execute_reply":"2023-04-12T06:53:02.082482Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"image_stack_ = gaussian_filter(image_stack_,sigma=1) #blur again","metadata":{"collapsed":false,"ExecuteTime":{"start_time":"2023-04-12T07:52:34.525694Z","end_time":"2023-04-12T07:53:01.175893Z"},"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2023-04-12T06:53:02.084277Z","iopub.execute_input":"2023-04-12T06:53:02.085748Z","iopub.status.idle":"2023-04-12T06:53:10.481056Z","shell.execute_reply.started":"2023-04-12T06:53:02.085682Z","shell.execute_reply":"2023-04-12T06:53:10.479620Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now for every pixel on the 2d-plane we find the first depth-index where a value is >=0.5. The 0.5 was chosen by trying out.","metadata":{}},{"cell_type":"code","source":"filtered_stack=np.where(image_stack_ >= 0.5, 1, 0)\ntopographic_map = np.argmax(filtered_stack, axis=0)\nif GAUSSIAN_BLUR_TOPOGRAPHIC_MAP:\n    topographic_map=gaussian_filter(topographic_map,sigma=1)","metadata":{"collapsed":false,"ExecuteTime":{"start_time":"2023-04-12T07:53:01.186862Z","end_time":"2023-04-12T07:53:06.145049Z"},"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2023-04-12T06:53:10.483734Z","iopub.execute_input":"2023-04-12T06:53:10.484006Z","iopub.status.idle":"2023-04-12T06:53:14.457176Z","shell.execute_reply.started":"2023-04-12T06:53:10.483980Z","shell.execute_reply":"2023-04-12T06:53:14.455430Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's have a look at the data","metadata":{}},{"cell_type":"code","source":"fig, (ax1, ax2,ax3,ax4) = plt.subplots(4, 1)\nslice_at=500\nslice_length=slice(500,1500)\nax1.imshow(image_stack[:,slice_length,slice_at], cmap='gray')\nax2.imshow(image_stack_[:,slice_length,slice_at], cmap='gray')\nax3.imshow(filtered_stack[:,slice_length,slice_at], cmap='gray')\nax4.plot(topographic_map[slice_length,slice_at]*-1)\nplt.show()","metadata":{"collapsed":false,"ExecuteTime":{"start_time":"2023-04-12T07:53:15.938177Z","end_time":"2023-04-12T07:53:19.020198Z"},"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2023-04-12T06:53:14.459575Z","iopub.execute_input":"2023-04-12T06:53:14.459981Z","iopub.status.idle":"2023-04-12T06:53:14.815983Z","shell.execute_reply.started":"2023-04-12T06:53:14.459953Z","shell.execute_reply":"2023-04-12T06:53:14.814429Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Looks like the topograhic_map is accurate enough for our slice. Some points are obviously wrong but overall it seems like a nice fit. Another round of gaussian filtering could help but i will leave that out for now.\nIf we look at the whole topograhic_map we can clearly see the surface of the papyrus.","metadata":{}},{"cell_type":"code","source":"plt.imshow(topographic_map, cmap='gray')","metadata":{"collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2023-04-12T06:53:14.817649Z","iopub.execute_input":"2023-04-12T06:53:14.817997Z","iopub.status.idle":"2023-04-12T06:53:15.379148Z","shell.execute_reply.started":"2023-04-12T06:53:14.817963Z","shell.execute_reply":"2023-04-12T06:53:15.377520Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now we can use the topographic map to flatten the data.","metadata":{}},{"cell_type":"code","source":"z_buffer=5 # how much air we will leave above the papyrus\nis_idx=np.indices(image_stack.shape)\nimage_stack_flattened=image_stack[(is_idx[0] + topographic_map-z_buffer) % image_stack.shape[0],is_idx[1],is_idx[2]]","metadata":{"collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2023-04-12T06:53:15.380822Z","iopub.execute_input":"2023-04-12T06:53:15.381111Z","iopub.status.idle":"2023-04-12T06:53:27.592872Z","shell.execute_reply.started":"2023-04-12T06:53:15.381084Z","shell.execute_reply":"2023-04-12T06:53:27.591769Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's have a look at the flattened data. We can see the surface is even, but we can also see that the layers under it are warped. So we probably should combine this in training with other data to not lose information. I will leave that to you for now.","metadata":{}},{"cell_type":"code","source":"plt.imshow(image_stack_flattened[:,slice_length,slice_at], cmap='gray')","metadata":{"collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2023-04-12T06:53:27.594408Z","iopub.execute_input":"2023-04-12T06:53:27.594687Z","iopub.status.idle":"2023-04-12T06:53:27.720020Z","shell.execute_reply.started":"2023-04-12T06:53:27.594662Z","shell.execute_reply":"2023-04-12T06:53:27.719299Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Looks very promising to me. Let's put it in a function we can use.","metadata":{}},{"cell_type":"code","source":"def flatten(arr, z_buffer, z_layers):\n    \"\"\"\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 = np.flip(arr, axis=0)\n    arr = gaussian_filter(arr, sigma=1)\n    arr = ndimage.sobel(arr, axis=0)\n    arr = gaussian_filter(arr, sigma=1)\n    topo = np.argmax(np.where(arr < 0.5, 0, 1), axis=0)\n    if GAUSSIAN_BLUR_TOPOGRAPHIC_MAP:\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]\n    return arr","metadata":{"collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2023-04-12T06:53:27.721373Z","iopub.execute_input":"2023-04-12T06:53:27.721864Z","iopub.status.idle":"2023-04-12T06:53:27.729306Z","shell.execute_reply.started":"2023-04-12T06:53:27.721828Z","shell.execute_reply":"2023-04-12T06:53:27.728571Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"If you don't have enough memory to flatten the whole image_stack in one go, here is a way to do it:","metadata":{}},{"cell_type":"code","source":"def do_in_stripes(arr, stripe_width,stripe_overlay,file_path, func,*args,**kwargs):\n    \"\"\"\n\n    :param arr: array to work on\n    :param stripe_width: width of stripes in pixels\n    :param stripe_overlay: overlay between stripes in pixels so we don't have artefacts at the borders\n    :param file_path: A file path to temporarily save data to\n    :param func: The function we want to apply to our array\n    :param args: args for the func\n    :param kwargs: kwargs for the func\n    :return: The number of stripes saved, so we can load them easily\n    \"\"\"\n    sc=0\n    with open(file_path, 'wb') as f:\n        for i in tqdm(range(0,arr.shape[1],stripe_width)):\n            sc=sc+1\n            start=max(i-stripe_overlay,0)\n            end=min(i+stripe_width+stripe_overlay,arr.shape[1])\n            stripe=func(arr[:,start:end,:],*args,**kwargs)\n            if end==arr.shape[1]:\n                np.save(f, stripe[:,stripe_overlay:stripe_width+stripe_overlay,:])\n                break\n            elif i==0:\n                np.save(f, stripe[:,0:stripe_width,:])\n            else:\n                np.save(f, stripe[:,stripe_overlay:stripe_width+stripe_overlay,:])\n    return sc","metadata":{"collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2023-04-12T06:53:27.732378Z","iopub.execute_input":"2023-04-12T06:53:27.732855Z","iopub.status.idle":"2023-04-12T06:53:27.744973Z","shell.execute_reply.started":"2023-04-12T06:53:27.732818Z","shell.execute_reply":"2023-04-12T06:53:27.743553Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if TRANSFORM_WHOLE_STACK:\n    #Load from disk. We saved it earlier\n    with open(path.join(PATH_TO_FRAGMENT,\"image_stack.npy\"), 'rb') as f:\n        image_stack = np.load(f)\n    stripe_count=do_in_stripes(image_stack,500,20,path.join(WORK_DIR,\"stripes\"),flatten,3,15) #3 is z_buffer and 15 is z_layers for flatten function\n    del image_stack\n    with open(path.join(WORK_DIR,\"stripes\"), 'rb') as f:\n        stripes=[]\n        for i in range(0,stripe_count):\n            stripes.append(np.load(f))\n    image_stack_flattened=np.concatenate(stripes,axis=1)","metadata":{"collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2023-04-12T06:53:27.746305Z","iopub.execute_input":"2023-04-12T06:53:27.746699Z","iopub.status.idle":"2023-04-12T06:53:27.755605Z","shell.execute_reply.started":"2023-04-12T06:53:27.746661Z","shell.execute_reply":"2023-04-12T06:53:27.754886Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"I hope you can use this in some way and understood the method. Sorry for my language and the short sentences but my written english is not the best.\nHave fun in the competition!","metadata":{}}]}