{"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":"## Fast & efficient image storing and manipulation\n### Brett Olsen, April 2023\n\nThis is an extension of my previous work that now has code stored in github that we can use for fast and efficient storage and manipulation of the papyrus image data, both here in Kaggle and in the larger challenge.\n\nThere are two classes of interest:  `PapyrusImage`, which is for storing flat images as chunked 2D arrays, and `PapyrusVolume`, which is for storing image stacks as chunked 3D arrays.\nBy default, the chunk size in 1024 pixels in the XY dimension and only 1 pixel in the Z dimension.\nTo use these stores, you will need to convert (once) the input file(s) to the new format and then can rapidly load subsections of the\ndata into memory.\nWe'll demonstrate some simple ways to use this capability.\nWhen generating the new format, you may also limit the range of the input data to convert in X, Y, or Z, if, e.g., some of the images are not useful for the work you are trying to do.\nBy default, the conversion additionally adds downsampled copies of the image into its internal store.\nThis allows for faster visualization of zoomed-out data at a slight cost of increased storage requirements.","metadata":{}},{"cell_type":"markdown","source":"First we're going to download my `vesuvius-image` package from github and install its requirements into our local conda environment.\nThis can take a couple of minutes for mamba to download all the necessary packages.","metadata":{}},{"cell_type":"code","source":"GITHUB_REPO_URL = \"https://github.com/caethan/vesuvius_image.git\"\n!git clone {GITHUB_REPO_URL}\n!git -C /kaggle/working/vesuvius_image/ pull","metadata":{"execution":{"iopub.status.busy":"2023-04-02T04:19:49.255364Z","iopub.execute_input":"2023-04-02T04:19:49.256476Z","iopub.status.idle":"2023-04-02T04:19:51.896813Z","shell.execute_reply.started":"2023-04-02T04:19:49.256433Z","shell.execute_reply":"2023-04-02T04:19:51.895395Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!mamba env update --name base --file vesuvius_image/conda_env.yaml","metadata":{"execution":{"iopub.status.busy":"2023-04-02T04:19:55.083485Z","iopub.execute_input":"2023-04-02T04:19:55.083933Z","iopub.status.idle":"2023-04-02T04:24:25.800071Z","shell.execute_reply.started":"2023-04-02T04:19:55.083891Z","shell.execute_reply":"2023-04-02T04:24:25.798211Z"},"collapsed":true,"jupyter":{"outputs_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Great.  Now let's add our new package into `sys.path`, import the `PapyrusImage` and `PapyrusVolume` classes,\nas well as all the other stuff we'll need.","metadata":{}},{"cell_type":"code","source":"import sys\nif \"/kaggle/working/vesuvius_image/\" not in sys.path:\n    sys.path.append(\"/kaggle/working/vesuvius_image/\")\n\nfrom vesuvius import PapyrusImage, PapyrusVolume\nfrom vesuvius.utils import Timer","metadata":{"execution":{"iopub.status.busy":"2023-04-02T04:38:57.820129Z","iopub.execute_input":"2023-04-02T04:38:57.821162Z","iopub.status.idle":"2023-04-02T04:38:57.826554Z","shell.execute_reply.started":"2023-04-02T04:38:57.821117Z","shell.execute_reply":"2023-04-02T04:38:57.825425Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport plotly\nimport plotly.graph_objects as go\nimport dash\nimport resource\nfrom dash import dcc\nfrom dash import html\nfrom jupyter_dash import JupyterDash\nfrom dash.dependencies import Input, Output\n\nINPUT_FOLDER = \"/kaggle/input/vesuvius-challenge-ink-detection\"\nWORKING_FOLDER = \"/kaggle/working/\"\nTEMP_FOLDER = \"/kaggle/temp/\"","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-04-02T04:33:25.709983Z","iopub.execute_input":"2023-04-02T04:33:25.711375Z","iopub.status.idle":"2023-04-02T04:33:25.721663Z","shell.execute_reply.started":"2023-04-02T04:33:25.711288Z","shell.execute_reply":"2023-04-02T04:33:25.720436Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"OK.  The first thing we're going to do is convert some of the input data into the new format in our working directory.\nWe'll go ahead and do the conversion on all of the image files for the first training set.\nThis will take a couple of minutes, especially for the surface volume.","metadata":{}},{"cell_type":"code","source":"input_dir = os.path.join(INPUT_FOLDER, \"train/1/\")\noutput_dir = os.path.join(WORKING_FOLDER, \"train_1\")\nos.makedirs(output_dir)\n\nir_input = os.path.join(input_dir, \"ir.png\")\nir_output = os.path.join(output_dir, \"ir.imagezarr\")\nif not os.path.exists(ir_output):\n    PapyrusImage.build_from_image(ir_input, ir_output)\n    \nink_input = os.path.join(input_dir, \"inklabels.png\")\nink_output = os.path.join(output_dir, \"inklabels.imagezarr\")\nif not os.path.exists(ink_output):\n    PapyrusImage.build_from_image(ink_input, ink_output)\n\nmask_input = os.path.join(input_dir, \"mask.png\")\nmask_output = os.path.join(output_dir, \"mask.imagezarr\")\nif not os.path.exists(mask_output):\n    PapyrusImage.build_from_image(mask_input, mask_output)\n    \nfrag_input = os.path.join(input_dir, \"surface_volume\")\nfrag_output = os.path.join(output_dir, \"frag.volzarr\")\nif not os.path.exists(frag_output):\n    PapyrusVolume.build_from_tiffdir(frag_input, frag_output)","metadata":{"execution":{"iopub.status.busy":"2023-04-02T04:24:26.419130Z","iopub.execute_input":"2023-04-02T04:24:26.419483Z","iopub.status.idle":"2023-04-02T04:29:27.631854Z","shell.execute_reply.started":"2023-04-02T04:24:26.419448Z","shell.execute_reply":"2023-04-02T04:29:27.629078Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Most of the time spent was in downsampling the surface volume images for the lower-resolution arrays.\nIf you're not planning on using these, then you can save time by passing the argument `multiscale=False` and only\nthe full-resolution arrays will be stored.\nLet's take a look at how we can use these files, but first let's look at the sizes.\nWe're currently using an `SQLiteStore` to store all the array data inside a single file.\nThis is a fair bit slower to create the zarr file than just using the default `DirectoryStore` that generates new files for each chunk, but it's much more manageable to work with, especially for large arrays.","metadata":{}},{"cell_type":"code","source":"!du -h {frag_input}\n!ls -l -h {input_dir}\n!ls -l -h {output_dir}","metadata":{"execution":{"iopub.status.busy":"2023-04-02T04:33:29.072816Z","iopub.execute_input":"2023-04-02T04:33:29.073223Z","iopub.status.idle":"2023-04-02T04:33:32.205189Z","shell.execute_reply.started":"2023-04-02T04:33:29.073188Z","shell.execute_reply":"2023-04-02T04:33:32.203617Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can see that for the small single images, we're not saving space:  the combination of the additional structure to chunk and store the data along with the multiscaling giving us extra data is significantly more than the savings we get from compressing and chunking the data.\nFor the boolean masks, we go from tens of KB to hundreds of KB, while for the more complex infrared image, we go from about 8MB to about 30MB.\nHowever, for the much larger surface volume array, we start actually saving space:  the full directory of `.tif` files takes about 6.26GB, but the single zarr file containing not just the full array but the downsampled ones as well fits into only 4.2GB.\nSo for larger datasets, this file format can not only speed access but also cut down on space used in both memory and on disk.\n\nLet's go ahead and load these volumes in and take a look at the memory usage.","metadata":{}},{"cell_type":"code","source":"def show_memuse():\n    rusage = resource.getrusage(resource.RUSAGE_SELF)\n    # Convert from kB to GB\n    rss_gb = rusage.ru_maxrss / 1024 ** 2\n    print(f\"Memory used: {rss_gb:.2f}GB\")\n    \nshow_memuse()\nir_image = PapyrusImage(ir_output)\nshow_memuse()\nink_image = PapyrusImage(ink_output)\nshow_memuse()\nmask_image = PapyrusImage(mask_output)\nshow_memuse()\nsurface_volume = PapyrusVolume(frag_output)\nshow_memuse()","metadata":{"execution":{"iopub.status.busy":"2023-04-02T04:36:29.736207Z","iopub.execute_input":"2023-04-02T04:36:29.737251Z","iopub.status.idle":"2023-04-02T04:36:29.829092Z","shell.execute_reply.started":"2023-04-02T04:36:29.737207Z","shell.execute_reply":"2023-04-02T04:36:29.828153Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We now have access to all this data, but the memory is barely touched!\nThat's because the data is not loaded into memory until we actually need it, and when we do, we only load the chunks we need.\nLet's do a quick view of the full infrared image with matplotlib.","metadata":{}},{"cell_type":"code","source":"with Timer():\n    plt.figure(figsize=(6,6))\n    plt.imshow(ir_image[:,:], cmap=\"gray\")\n    show_memuse()","metadata":{"execution":{"iopub.status.busy":"2023-04-02T04:39:01.915574Z","iopub.execute_input":"2023-04-02T04:39:01.916298Z","iopub.status.idle":"2023-04-02T04:39:02.423618Z","shell.execute_reply.started":"2023-04-02T04:39:01.916254Z","shell.execute_reply":"2023-04-02T04:39:02.422389Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"See how the image is actually much smaller, only about 1K pixels?\nThat's because it's not actually loading the whole figure.\nBy default, it will try to load a downsampled image, the one that is about 1K pixels wide.\nIf we want, we can force a load of the full-scale image by adding a third index.\n0 is the full image, and larger indices will provide more heavily downsampled images.\nLet's compare.","metadata":{}},{"cell_type":"code","source":"for ds_level in range(ir_image.scale_levels + 1):\n    with Timer():\n        plt.figure(figsize=(3,3))\n        plt.imshow(ir_image[:,:,ds_level], cmap=\"gray\")\n        show_memuse()","metadata":{"execution":{"iopub.status.busy":"2023-04-02T04:41:34.782176Z","iopub.execute_input":"2023-04-02T04:41:34.782609Z","iopub.status.idle":"2023-04-02T04:41:38.880181Z","shell.execute_reply.started":"2023-04-02T04:41:34.782570Z","shell.execute_reply":"2023-04-02T04:41:38.878835Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Looks like this fragment has 4 levels of downsampled images.\nLet's try a different plotting library that will let us do a little bit more.\nThe team has confirmed that the voxel size of the fragments is 4 microns, so we can show size directly on the axis.\nNow, for large-scale visualization we'll probably want to use the lowest-resolution data.","metadata":{}},{"cell_type":"code","source":"def plotly_image(img_array, size=600, pixel_size=4):\n    fig = go.Figure(layout={\n        \"autosize\": True,\n        \"xaxis\": {\"mirror\": True, \"side\": \"top\", \"autorange\": True, \"title\": \"mm\"},\n        \"yaxis\": {\"mirror\": True, \"side\": \"top\", \"autorange\": \"reversed\", \"scaleanchor\": \"x\", \"title\": \"mm\"},\n        \"height\": size,\n        \"width\": size,\n        \"dragmode\": \"pan\",\n    })\n    fig.add_trace(\n        go.Heatmap(\n            z = img_array[:,:],\n            colorscale=\"gray\",\n            x0 = 0,\n            dx = pixel_size,\n            y0 = 0,\n            dy = pixel_size,\n            transpose=False,\n            showscale=False,\n            hoverinfo=\"skip\",\n        )\n    )\n    fig.show(config={'scrollZoom': True})\n    \n# in mm\nDEFAULT_PIXEL_SIZE = 4 / 1000\nscale_level = 3\nplotly_image(ir_image[:,:,scale_level], pixel_size=DEFAULT_PIXEL_SIZE * 2 ** scale_level)","metadata":{"execution":{"iopub.status.busy":"2023-04-02T05:02:42.160582Z","iopub.execute_input":"2023-04-02T05:02:42.161026Z","iopub.status.idle":"2023-04-02T05:02:42.225566Z","shell.execute_reply.started":"2023-04-02T05:02:42.160988Z","shell.execute_reply":"2023-04-02T05:02:42.224576Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's take a look at the surface volume as well.\nFor volumes, we'll need four indices in order to specify the resolution rather than just 3.\nNote as well that the downsampling happens in z as well, so not all z-frames are stored.","metadata":{}},{"cell_type":"code","source":"scale_level = 3\nz_level = 31\n# TODO: fix the bug with indexing by ints\nplotly_image(surface_volume[:,:,slice(z_level, z_level + 1, None),scale_level][:,:,0], pixel_size=DEFAULT_PIXEL_SIZE * 2 ** scale_level)","metadata":{"execution":{"iopub.status.busy":"2023-04-02T05:10:07.878216Z","iopub.execute_input":"2023-04-02T05:10:07.879215Z","iopub.status.idle":"2023-04-02T05:10:08.013037Z","shell.execute_reply.started":"2023-04-02T05:10:07.879167Z","shell.execute_reply":"2023-04-02T05:10:08.011432Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can also do some non-plotting manipulation.  \nLet's look at the non-zero elements of each z-stack and plot the mean and variance.\nThat should let us see where the stack is empty.","metadata":{}},{"cell_type":"code","source":"with Timer():\n    zindices = np.arange(surface_volume.shape[2])\n    means = np.zeros_like(zindices)\n    stdevs = np.zeros_like(zindices)\n    mask = mask_image[:,:,0]\n    for z in zindices:\n        print(z, end=\"\\r\")\n        array = surface_volume.zarr.volume[:,:,z][mask]\n        means[z] = array.mean()\n        stdevs[z] = np.std(array)","metadata":{"execution":{"iopub.status.busy":"2023-04-02T05:20:34.250769Z","iopub.execute_input":"2023-04-02T05:20:34.251704Z","iopub.status.idle":"2023-04-02T05:20:55.628404Z","shell.execute_reply.started":"2023-04-02T05:20:34.251659Z","shell.execute_reply":"2023-04-02T05:20:55.627123Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure()\nplt.grid()\nplt.errorbar(zindices, means, stdevs, marker='o')\nplt.plot(zindices, stdevs)","metadata":{"execution":{"iopub.status.busy":"2023-04-02T05:20:59.857759Z","iopub.execute_input":"2023-04-02T05:20:59.858962Z","iopub.status.idle":"2023-04-02T05:21:00.095053Z","shell.execute_reply.started":"2023-04-02T05:20:59.858916Z","shell.execute_reply":"2023-04-02T05:21:00.093884Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Interesting!  We can see an increase in signal between z 20 and 28 or so, followed by a drop-off and recovery.\nThe variance drops substantially past z 30 as well, presumably indicating that we're out of frame of the papyrus.","metadata":{}},{"cell_type":"code","source":"scale_level = 2\nz_level = 39\n# TODO: fix the bug with indexing by ints\nplotly_image(surface_volume[:,:,slice(z_level, z_level + 1, None),scale_level][:,:,0], pixel_size=DEFAULT_PIXEL_SIZE * 2 ** scale_level)","metadata":{"execution":{"iopub.status.busy":"2023-04-02T05:28:08.061923Z","iopub.execute_input":"2023-04-02T05:28:08.062388Z","iopub.status.idle":"2023-04-02T05:28:08.631797Z","shell.execute_reply.started":"2023-04-02T05:28:08.062347Z","shell.execute_reply":"2023-04-02T05:28:08.630501Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}