doc/jupyter_notebooks/pylibCZIrw_4_1_0.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/zeiss/pylibczirw/blob/main/doc/jupyter_notebooks/pylibCZIrw_4_1_0.ipynb)\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Install pylibCZIrw and dependencies\n",
    "! pip install --upgrade pip\n",
    "! pip install \"pylibCZIrw~=4.1\" \"cztile>=0.0,<1.0\" matplotlib tqdm pooch"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "id": "ad16447f",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## pylibCZIrw - Tutorial\n",
    "\n",
    "In this notebook one can see how to read and write CZI (Carl Zeiss Image) files using from Python using **[pylibCZIrw](https://pypi.org/project/pylibczirw/)**. This package is wrapper for the **libCZIrw** (C++) library and can be used both in Windows and Linux.\n",
    "\n",
    "### Remarks\n",
    "\n",
    "* At the moment, [pylibCZIrw](https://pypi.org/project/pylibczirw/) completely abstracts away the subblock concept, both in the reading and in the writing APIs.\n",
    "* If [pylibCZIrw](https://pypi.org/project/pylibczirw/) is extended in the future to support subblock-based access (e.g. accessing acquisition tiles), this API must not be altered.\n",
    "\n",
    "### Disclaimer\n",
    "\n",
    "This notebook is free to use for everybody. Carl Zeiss Microscopy GmbH's ZEN software undertakes no warranty concerning the use of those scripts, image analysis settings and ZEN experiments. Use them on your own risk.\n",
    "\n",
    "Additionally Carl Zeiss Microscopy GmbH's ZEN software allows connection and usage to the third party software packages.\n",
    "Therefore Carl Zeiss Microscopy GmbH undertakes no warranty concerning those software packages, makes no representation that they will work on your system and/or hardware and will not be liable for any damages caused by the use of this extension. **By using any of those examples you agree to this disclaimer.**\n",
    "\n",
    "Version: 2024.01.22\n",
    "\n",
    "Copyright (c) 2024 Carl Zeiss AG, Germany. All Rights Reserved."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "sufficient-scheme",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### Required Module Imports\n",
    "\n",
    "Both reading and writing APIs are contained within the **czi** object in pylibCZIrw, so let's start by importing it.\n",
    "\n",
    "Additional modules are imported alongside the **pylibCZIrw** package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "martial-craps",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "from pylibCZIrw import czi as pyczi\n",
    "from pylibCZIrw.czi import ReaderFileInputTypes, CacheOptions, CacheType\n",
    "import json\n",
    "import timeit\n",
    "from matplotlib import pyplot as plt\n",
    "import matplotlib.cm as cm\n",
    "import numpy as np\n",
    "import os, sys\n",
    "from tqdm import tqdm\n",
    "from tqdm.contrib import itertools as it\n",
    "from matplotlib.patches import Rectangle\n",
    "from IPython.display import display, JSON\n",
    "\n",
    "# show the used python env\n",
    "print(\"Using:\", sys.executable)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cde4f2a6",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Define parameters for data loading"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6fd4232a",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# Define CZI locations\n",
    "\n",
    "czifile_DAPI_PGC = r\"https://github.com/zeiss-microscopy/OAD/raw/master/jupyter_notebooks/pylibCZIrw/data/DAPI_PGC.czi\"\n",
    "\n",
    "czifile_5dstack = r\"https://github.com/zeiss-microscopy/OAD/raw/master/jupyter_notebooks/pylibCZIrw/data/T%3D3_Z%3D5_CH%3D2_X%3D240_Y%3D170.czi\"\n",
    "\n",
    "czifile_scenes = r\"https://github.com/zeiss-microscopy/OAD/raw/master/jupyter_notebooks/pylibCZIrw/data/w96_A1%2BA2.czi\"\n",
    "\n",
    "czifile_timeseries = r\"https://github.com/zeiss-microscopy/OAD/blob/master/jupyter_notebooks/pylibCZIrw/data/CellDivision_T%3D16_CH%3D2_ZSTD.czi\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b25068ab",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### Opening a CZI (read-only)\n",
    "\n",
    "A CZI file can be opened from disk in a context manager or alternatively it can also be opened directly from a seekable http url.\n",
    "\n",
    "* `pyczi.ReaderFileInputTypes.Curl` reads from http\n",
    "* * `pyczi.ReaderFileInputTypes.Standard` reads from file directly (default)\n",
    "\n",
    "```python\n",
    "from pylibCZIrw import pyczi\n",
    "\n",
    "with czi.open_czi(file_path) as czi:\n",
    "\n",
    "    # add your code here\n",
    "\n",
    "with czi.open_czi(file_url, pyczi.ReaderFileInputTypes.Curl) as czi:\n",
    "\n",
    "    # add your code here\n",
    "```\n",
    "\n",
    "Remark: Inside this notebook only the option to read from http is used which removes the need to download the data when using Colab notebooks."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6699726f",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### Open a CZI and read metadata\n",
    "\n",
    "We can open a czi using a context manager and then read the full metadata either as a raw XML string or as a dictionary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "urban-consumption",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# open the CZI for reading using a context manager (preferred way to do it)\n",
    "with pyczi.open_czi(czifile_5dstack, pyczi.ReaderFileInputTypes.Curl) as czidoc:\n",
    "    # get the raw metadata as XML\n",
    "    md_xml = czidoc.raw_metadata\n",
    "    print(md_xml[:1000])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d357f57a",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# open the CZI for reading\n",
    "with pyczi.open_czi(czifile_scenes, pyczi.ReaderFileInputTypes.Curl) as czidoc:\n",
    "    # get the raw metadata as a dictionary\n",
    "    md_dict = czidoc.metadata\n",
    "\n",
    "    # show some parts of it\n",
    "    # print(json.dumps(md_dict[\"ImageDocument\"][\"Metadata\"][\"Information\"][\"Image\"], sort_keys=False, indent=4))\n",
    "    display(JSON(md_dict[\"ImageDocument\"][\"Metadata\"][\"Information\"][\"Image\"], expanded=False))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "freelance-virtue",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Get Image dimensions\n",
    "There are different properties that allow us to retrieve information about dimensions.\n",
    "The **`total_bounding_box`** gives us all the dimensions of all **orthogonal planes** of the CZI image file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "desirable-layer",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "with pyczi.open_czi(czifile_5dstack, pyczi.ReaderFileInputTypes.Curl) as czidoc:\n",
    "    # get the image dimensions as an dictionary, where the key identifies the dimension\n",
    "    total_bounding_box = czidoc.total_bounding_box\n",
    "\n",
    "print(total_bounding_box)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "rapid-savings",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "The **`total_bounding_rectangle`** gives us the X and Y dimensions of the CZI, i.e. the (X, Y) of the **`total_bounding_box`**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b8177a12",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "with pyczi.open_czi(czifile_scenes, pyczi.ReaderFileInputTypes.Curl) as czidoc:\n",
    "    # get the total bounding box for all scenes\n",
    "    total_bounding_rectangle = czidoc.total_bounding_rectangle\n",
    "\n",
    "print(total_bounding_rectangle)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bce02b4e",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "**Scenes are not orthogonal to the other dimensions.**\n",
    "They are contained within the 2D planes and can be seen simply as tags for certain regions.\n",
    "\n",
    "The **`scene_bounding_rectangle`** give us the bounding rectangle for each scene."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18a7c877",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "with pyczi.open_czi(czifile_scenes, pyczi.ReaderFileInputTypes.Curl) as czidoc:\n",
    "    # get the bounding boxes for each individual scene\n",
    "    scenes_bounding_rectangle = czidoc.scenes_bounding_rectangle\n",
    "\n",
    "print(scenes_bounding_rectangle)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "381aa4b3",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Get the pixel type\n",
    "A channel's pixel type can be discovered with using: **`get_channel_pixel_type`**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a55bc247",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "with pyczi.open_czi(czifile_scenes, pyczi.ReaderFileInputTypes.Curl) as czidoc:\n",
    "    # get the pixel type for the 1st channel\n",
    "    c0_pixel_type = czidoc.get_channel_pixel_type(0)\n",
    "\n",
    "print(c0_pixel_type)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e105838",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Or we can simply get all pixel types using: **`pixel_types`**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d5d8f6c",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "with pyczi.open_czi(czifile_scenes, pyczi.ReaderFileInputTypes.Curl) as czidoc:\n",
    "    # get all pixel types as a dictionary, where the key is the channel index\n",
    "    pixel_type = czidoc.pixel_types\n",
    "\n",
    "print(pixel_type)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "143fbcb0",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Read pixel data\n",
    "Pixel data is read and returned as a numpy arry. The shape of the returned array and its data type will depend on the specified `pixel_type`.\n",
    "- BGR pixel type -> [y, x, 3]\n",
    "- Gray pixel type -> [y, x, 1]\n",
    "\n",
    "The background_pixel determines the pixel value for the background (i.e. regions where no pixel data exists),\n",
    " and must be consistent with the `pixel_type`\n",
    "\n",
    "|pixel_type | background_pixel type | Default value  | \n",
    "--- | --- | ---\n",
    "|BGR|3-value tuple |(0, 0, 0)|\n",
    "|Gray|Scalar value bounded by the gray scale|0|\n",
    "\n",
    "### Read whole 2D plane\n",
    "\n",
    "In order to **`read`** pixel data from a CZI image file it is important to understand the concept of \"planewise reading\".\n",
    "\n",
    "Reading data planewise means that every read operation yields a 2D array with the shape (Y, X, 1) or (Y, X, 3) depending on the requested pixel type.\n",
    "We can specify which 2D plane to read using the **`plane`** argument (dictionary).\n",
    "\n",
    "* **`plane = {'C': 0, 'Z': 2, 'T': 1}`** - this will read a 2D plane from C=0, Z=2 and T=1\n",
    "* if **`plane`** is not specified at all the default is to use the 1st index for all dimensions\n",
    "\n",
    "In order to read a complete Z-Stack, one would have to read each of the Z-slices using the plane argument. This can be easily created by \"wrapping\" those calls into some utility functions.\n",
    "To get started let's read a complete 2D plane."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "114e59e5",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "with pyczi.open_czi(czifile_5dstack, pyczi.ReaderFileInputTypes.Curl) as czidoc:\n",
    "    # define some plane coordinates\n",
    "    plane_1 = {'C': 0, 'Z': 2, 'T': 1}\n",
    "    plane_2 = {'C': 1, 'Z': 3, 'T': 2}\n",
    "\n",
    "    # equivalent to reading {'C': 0, 'Z': 0, 'T': 0}\n",
    "    frame_0 = czidoc.read()\n",
    "\n",
    "    # get the shape of the 2d plane - the last dime indicates the pixel type\n",
    "    # 3 = BGR and 1 = Gray\n",
    "    print(\"Array Shape: \", frame_0.shape)\n",
    "\n",
    "    # get specific planes \n",
    "    frame_1 = czidoc.read(plane=plane_1)\n",
    "    frame_2 = czidoc.read(plane=plane_2)\n",
    "\n",
    "fig, ax = plt.subplots(1, 3, figsize=(15, 5))\n",
    "ax[0].imshow(frame_0[..., 0], cmap=cm.inferno)\n",
    "ax[0].set_title(\"Frame_0\")\n",
    "ax[1].imshow(frame_1[..., 0], cmap=cm.inferno)\n",
    "ax[1].set_title(\"Frame_1\")\n",
    "ax[2].imshow(frame_2[..., 0], cmap=cm.Greens_r)\n",
    "ax[2].set_title(\"Frame_2\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cdd64f81",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### Read ROI inside a specific 2D plane\n",
    "\n",
    "We can also read an ROI within the 2D plane by using the **`roi`** argument.\n",
    "This is particularly useful when dealing with large images whose planes would not fit into memory.\n",
    "\n",
    "The **`roi`** is a tuple defined as a (axis-aligned) rectangle in **(x, y, w, h)** form, with:\n",
    "\n",
    "* x = x coordinate of upper-left point\n",
    "* y = y coordinate of upper-left point\n",
    "* w = rectangle width\n",
    "* h = rectangle height"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "71c47e67",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# define the ROI using upper left, widht and height as coordinates\n",
    "my_roi = (200, 400, 800, 600)\n",
    "\n",
    "with pyczi.open_czi(czifile_scenes, pyczi.ReaderFileInputTypes.Curl) as czidoc:\n",
    "    ch0 = czidoc.read(roi=my_roi, plane={'C': 0})\n",
    "    ch1 = czidoc.read(roi=my_roi, plane={'C': 1})\n",
    "\n",
    "fig, ax = plt.subplots(1, 2, figsize=(10, 5))\n",
    "ax[0].imshow(ch0[..., 0], cmap=cm.inferno, vmin=100, vmax=4000)\n",
    "ax[0].set_title(\"ch0\")\n",
    "ax[1].imshow(ch1[..., 0], cmap=cm.Greens_r, vmin=100, vmax=4000)\n",
    "ax[1].set_title(\"ch1\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc3927ef",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### Read ROI only for one scene\n",
    "We can also get pixel data exclusively from one scene using the **`scene`** argument. This argument acts as\n",
    "a filter returning only data from the specified scene. If not set, data from all scenes is returned."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "776290df",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "with pyczi.open_czi(czifile_scenes, pyczi.ReaderFileInputTypes.Curl) as czidoc:\n",
    "    # read a 2D image from a specific channel and scene\n",
    "    c0_s0 = czidoc.read(plane={'C': 0}, scene=0)\n",
    "    c1_s0 = czidoc.read(plane={'C': 1}, scene=0)\n",
    "    c0_s1 = czidoc.read(plane={'C': 0}, scene=1)\n",
    "    c1_s1 = czidoc.read(plane={'C': 1}, scene=1)\n",
    "\n",
    "fig, ax = plt.subplots(2, 2, figsize=(16, 16))\n",
    "ax[0, 0].imshow(c0_s0[..., 0], cmap=cm.inferno, vmin=100, vmax=4000)\n",
    "ax[0, 0].set_title(\"scene:0 ch:0\")\n",
    "ax[0, 1].imshow(c1_s0[..., 0], cmap=cm.Greens_r, vmin=100, vmax=4000)\n",
    "ax[0, 1].set_title(\"scene:0 ch:1\")\n",
    "ax[1, 0].imshow(c0_s1[..., 0], cmap=cm.inferno, vmin=100, vmax=4000)\n",
    "ax[1, 0].set_title(\"scene:1 ch:0\")\n",
    "ax[1, 1].imshow(c1_s1[..., 0], cmap=cm.Greens_r, vmin=100, vmax=4000)\n",
    "ax[1, 1].set_title(\"scene:1 ch:0\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb589b3d",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### Additional parameters\n",
    "\n",
    "We can specify the **`pixel_type`** in which we want the data to be read. If the desired pixel type is not specified, the data is returned in its original pixel format."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb09392f",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "with pyczi.open_czi(czifile_5dstack, pyczi.ReaderFileInputTypes.Curl) as czidoc:\n",
    "    \n",
    "    # get image pixel and specify the desired pixel type\n",
    "    frame_gray8 = czidoc.read(pixel_type='Gray8')\n",
    "    frame_gray32float = czidoc.read(pixel_type='Gray32Float')\n",
    "    frame_bgr48 = czidoc.read(pixel_type='Bgr48')\n",
    "\n",
    "print('shape {0}, data type {1}'.format(frame_gray8.shape, frame_gray8.dtype))\n",
    "print('shape {0}, data type {1}'.format(frame_gray32float.shape, frame_gray32float.dtype))\n",
    "print('shape {0}, data type {1}'.format(frame_bgr48.shape, frame_bgr48.dtype))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7ff5d1d",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "We can also specify a **`zoom`** factor between 0 and 1. This can be very useful to read images at a lower resolution,\n",
    "because not every application requires the full resolution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "78245625",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# define the ROI using upper left, width and height as coordinates\n",
    "my_roi = (300, 400, 300, 300)\n",
    "\n",
    "# define two different zoom factors\n",
    "zl1 = 1.0\n",
    "zl2 = 0.25\n",
    "\n",
    "with pyczi.open_czi(czifile_scenes, pyczi.ReaderFileInputTypes.Curl) as czidoc:\n",
    "    \n",
    "    z1 = czidoc.read(plane={'C': 0}, roi=my_roi, scene=0, zoom=zl1, pixel_type='Gray8')\n",
    "    z2 = czidoc.read(plane={'C': 0}, roi=my_roi, scene=0, zoom=zl2, pixel_type='Gray8')\n",
    "\n",
    "fig, ax = plt.subplots(1, 2, figsize=(16, 8))\n",
    "ax[0].imshow(z1[..., 0], cmap=cm.inferno)\n",
    "ax[1].imshow(z2[..., 0], cmap=cm.inferno)\n",
    "\n",
    "print(\"Shape Zoomlevel = \", zl1, z1.shape)\n",
    "print(\"Shape Zoomlevel = \", zl2, z2.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a0b8077d",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Writing a CZI file\n",
    "\n",
    "* **The writing API only allows creating new CZI files.**\n",
    "* editing existing czi file is not supported.\n",
    "* writing API is follows the same principles as the reading API.\n",
    "\n",
    "### Create a new CZI\n",
    "\n",
    "Similarly to reading, before writing, we must create a CZI instance.\n",
    "Ideally this is done using a context manager again.\n",
    "\n",
    "### Write pixel data\n",
    "\n",
    "Writing pixel data is very similar to reading albeit with fewer arguments:\n",
    "\n",
    "* the **`plane`** to which the pixel data should be written.\n",
    "* the data to be written, as a numpy array.\n",
    "\n",
    "The pixel type in which the data is written is automatically determined from the shape and data type of the provided array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7422d863",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "array = np.random.randint(low=1, high=10, size=(32, 64, 64)).astype(np.uint8)\n",
    "print(array.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bbc578c6",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "In order to create a new CZI file we must provide the file name.\n",
    "The optional flag (**`exist_ok`**) allows overwriting if set to **`True`** (default is **`False`**)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44ba6c63",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# create the filename for the new CZI image file\n",
    "newczi_4dstack = os.path.join(os.getcwd(), \"newCZI_stack.czi\")\n",
    "numZ = 32\n",
    "\n",
    "# open a new CZI and allow overwrite (!!!) to play around ...\n",
    "with pyczi.create_czi(newczi_4dstack, exist_ok=True) as czidoc_w:\n",
    "    # loop over all z-planes and channels\n",
    "    for frame in tqdm(range(numZ)):\n",
    "        # get the 2d array for the current plane and add axis to get (Y, X, 1) as shape\n",
    "        # write the plane with shape (Y, X, 1) to the new CZI file\n",
    "        czidoc_w.write(\n",
    "            data=array[frame, ...][..., np.newaxis],\n",
    "            plane={\"Z\": frame}\n",
    "        )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7df77b11",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "We can also write data to a specific **`location`** within the plane by providing the upper left pixel coordinates.\n",
    "\n",
    "To illustrate this idea let's write all z-planes from the image \"side-by-side\" into a new CZI image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af9d0d03",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# create the filename for the new CZI image file\n",
    "newczi_zloc = os.path.join(os.getcwd(), \"newCZI_zloc.czi\")\n",
    "xstart = 0\n",
    "\n",
    "with pyczi.create_czi(newczi_zloc, exist_ok=True) as czidoc_w:\n",
    "    # loop over all z-planes\n",
    "    for frame in tqdm(range(4)):\n",
    "        # get the 2d array for the current plane and add axis to get (Y, X, 1) as shape\n",
    "        # for fun - write the z-planes to different locations as a single channel\n",
    "        czidoc_w.write(\n",
    "            data=array[frame, ...][..., np.newaxis],\n",
    "            plane={\"C\": 0},\n",
    "            location=(xstart, 0)\n",
    "        )\n",
    "\n",
    "        # change the x-position for the next round to write \"side-by-side\"\n",
    "        xstart = xstart + 80"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74a8b9a3",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "![newCZI_stack.czi - shown in ZEN blue](screenshots/new_czi_stack.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c2e0fce",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### Write scenes\n",
    "\n",
    "Writing a scene is in fact \"just\" writing pixel data to a specific region within a plane and tagging it with a **`scene`** index.\n",
    "To illustrate this we now will write the z-planes to different scenes at different locations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5943b3a9",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# first step is to create some kind of grid of locations\n",
    "locx = []\n",
    "locy = []\n",
    "xystart = 0\n",
    "offset = 128\n",
    "\n",
    "# create the list for the 4x4 grid locations for the 16 planes \n",
    "for x, y in it.product(range(4), range(4)):\n",
    "    locx.append(xystart + offset * x)\n",
    "    locy.append(xystart + offset * y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "286a855a",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# create the filename for the new CZI image file\n",
    "newczi_zscenes = os.path.join(os.getcwd(), \"newCZI_scenes.czi\")\n",
    "\n",
    "with pyczi.create_czi(newczi_zscenes, exist_ok=True) as czidoc_w:\n",
    "    # loop over all z-planes\n",
    "    for frame in tqdm(range(16)):\n",
    "        # get the 2d array for the current plane and add axis to get (Y, X, 1) as shape\n",
    "        # for \"fun\" - write the z-planes to different scenes using the locations we just created\n",
    "        czidoc_w.write(\n",
    "            data=array[frame, ...][..., np.newaxis],\n",
    "            plane={\"C\": 0},\n",
    "            scene=frame,\n",
    "            location=(locx[frame], locy[frame])\n",
    "        )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63f19967",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "![new_czi_scenes - shown in ZEN blue](screenshots/new_czi_scenes.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b08dfc02",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### ZSTD Compression\n",
    "\n",
    "#### Why we have used ZSTD Compression?\n",
    "It is a fast compression technique, which provides the high compression ratios. We are handling these three modes:\n",
    "\n",
    "1. uncompressed: - Default setting, the image would not be compressed.\n",
    "2. zstd0 - You can specify the parameters of it like \"explicitlevel\" between -131072 and 22. <b>(e.g. \"zstd0:ExplicitLevel=10\")</b>\n",
    "3. zstd1 - You can specify the parameters of it like \"explicitlevel\" between -131072 and 22 or \"preprocess\" - only have HiLoByteUnpack. This is only done for zstd1 and for pixel types gray16 and Bgr48. <b>(e.g. \"zstd1:ExplicitLevel=10\" or \"zstd1:Preprocess=HiLoByteUnpack\")</b>\n",
    "\n",
    "To read more: https://facebook.github.io/zstd/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18a1fd30",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# create some array to write to CZI\n",
    "array_big = np.random.randint(low=1, high=50, size=(4096, 4096)).astype(np.uint16)\n",
    "\n",
    "# create the filename for the new CZI image file\n",
    "newczi_uncompressed = os.path.join(os.getcwd(), \"newCZI_uncompressed.czi\")\n",
    "\n",
    "# create the filename for the new CZI image file\n",
    "newczi_compressed = os.path.join(os.getcwd(), \"newCZI_compressed.czi\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9df34f48",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "with pyczi.create_czi(newczi_compressed, exist_ok=True) as czidoc_w:\n",
    "    \n",
    "    # The CZIWriter class has a compression parameter that allows to define a default compression for each write call.\n",
    "    # The write call will the allow you to change the compression for a specific case by setting its compression parameter explicitly as here changing to the zstd0: \n",
    "    # writing the image with \"zstd0\" compression.\n",
    "    czidoc_w.write(\n",
    "        data=array_big,\n",
    "        plane={\"C\": 0, \"T\": 0, \"Z\": 0},\n",
    "        scene=0,\n",
    "        compression_options = \"zstd0:ExplicitLevel=10\",\n",
    "        location=(0,0)\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "881320d9",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# The compression option could be defined while creating an czi file and could still be overwritten \n",
    "# with each individual call to the write function.\n",
    "with pyczi.create_czi(newczi_uncompressed, exist_ok=True, compression_options = \"uncompressed:\") as czidoc_w:\n",
    "                \n",
    "        czidoc_w.write(\n",
    "            data=array_big,\n",
    "            plane={\"C\": 0, \"T\": 0, \"Z\": 0},\n",
    "            scene=0,\n",
    "            location=(0,0)\n",
    "        )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "09e7d48b",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "size = os.path.getsize(newczi_uncompressed)/1024.0**2\n",
    "print(f\"Size of the uncompressed CZI file : {np.round(size, 2)} MB \")\n",
    "\n",
    "size_comp = os.path.getsize(newczi_compressed)/1024.0**2\n",
    "print(f\"Size of the compressed CZI file 1 : {np.round(size_comp, 2)} MB \")\n",
    "\n",
    "print(f\"ZSTD Compression Ratio 1: {np.round(size/size_comp, 2)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca687707",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### Write Metadata\n",
    "\n",
    "**Metadata is automatically written**, however, we can write some additional metadata by explicitly calling\n",
    " **`write_metadata`** with the following **optional** parameters:\n",
    "\n",
    "- **`document_name`** (in ZEN: \"Title\")\n",
    "- **`channel_names`** (in ZEN: \"Channel Names\")\n",
    "- **`custom_attributes`** (Not visible in ZEN)\n",
    "- **`scale_x`**, **`scale_y`**, **`scale_z`** in **meters (!)** (in ZEN: \"Scaling (per Pixel)\")\n",
    "\n",
    "The parameter **`custom_attributes`** allows customized attributes with type int, float, boolean or string. If an attribute already exists, writing the same attribute again will replace the existing value.\n",
    "\n",
    "This list of parameters will be extended on a case-by-case basis in the future.\n",
    "\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "654ebfd3",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# create the filename for the new CZI image file\n",
    "newczi_mdata = os.path.join(os.getcwd(), \"newczi_mdata.czi\")\n",
    "numZ = 32\n",
    "\n",
    "# open a new CZI and allow overwrite (!!!) to play around ...\n",
    "with pyczi.create_czi(newczi_mdata, compression_options = \"zstd0:ExplicitLevel=10\", exist_ok=True) as czidoc_w:\n",
    "    # loop over all z-planes\n",
    "    for z in tqdm(range(numZ)):\n",
    "        # get the 2d array for the current plane and add axis to get (Y, X, 1) as shape\n",
    "        # write the plane with shape (Y, X, 1) to the new CZI file\n",
    "        czidoc_w.write(\n",
    "            data=array[z, ...][..., np.newaxis],\n",
    "            plane={\"Z\": z}\n",
    "        )\n",
    "\n",
    "    # write the document title, channel names, custom attributes and XYZ scaling to the CZI file\n",
    "    czidoc_w.write_metadata(\n",
    "        document_name='new_czi_z=32',\n",
    "        channel_names={0: \"Red\"},\n",
    "        custom_attributes={\"key1\": \"value1\", \"key2\": \"value2\"},\n",
    "        scale_x=0.1 * 10 ** -6,\n",
    "        scale_y=0.1 * 10 ** -6,\n",
    "        scale_z=0.2 * 10 ** -6\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c54dc25",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "![new_czi_z=32_metadata - shown in ZEN blue](screenshots/new_czi_z32_metadata.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e35b693",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Subblock caching\n",
    "\n",
    "In certain scenarios (e.g. reading from a stream or reading compressed images) caching uncompressed subblocks might yield a higher reading performance. A CZI reader allows for specifying caching options that lead to subblocks being cached after their first use."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b04210d",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "num_executions = 100\n",
    "\n",
    "for use_cache in [True, False]:\n",
    "    cache_options = CacheOptions(CacheType.Standard) if use_cache else None\n",
    "\n",
    "    with pyczi.open_czi(newczi_compressed, ReaderFileInputTypes.Standard, cache_options=cache_options) as czi_document:\n",
    "        total_time = timeit.timeit(czi_document.read, number=num_executions)\n",
    "        average_time = total_time / num_executions\n",
    "\n",
    "        print(\n",
    "            f\"Average time for read_czi (use_cache={use_cache}) over \"\n",
    "            f\"{num_executions} executions: {average_time:.6f} seconds\"\n",
    "        )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca1b952f",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Tiling a large image\n",
    "\n",
    "The read and write API allow the user to work with large images by allowing breaking down the access using the **`roi`** argument.\n",
    "However, defining these rois (aka tiles) can be a cumbersome process.\n",
    "\n",
    "The **czitile** package was created to answer this need.\n",
    "\n",
    "**The package contains pure logic for tiling a rectangle. This means it is not file or format bound in any way.**\n",
    "\n",
    "Let's read the original file using cztile and write a new with the same dimensions and scenes but with\n",
    "random data instead."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee1f8b52",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### Define a 2D processing function\n",
    "\n",
    "The benefits of tiling are evident when the data needs to be broken down for processing.\n",
    "Segmentation using Deep Neural Networks is a great example of a task where images must be cut to fit.\n",
    "We can use **czitile** to break down the image into tiles with some overlap.\n",
    "\n",
    "We can illustrate this with a simpler 2D processing function first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03987642",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "from skimage.filters import gaussian\n",
    "\n",
    "\n",
    "# define a generic 2D processing function - could by a DL segmentation as well\n",
    "def process2d(image2d: np.ndarray, **kwargs: int) -> np.ndarray:\n",
    "    # insert or modify the desired processing function here\n",
    "    image2d = gaussian(\n",
    "        image2d, sigma=kwargs[\"sigma\"],\n",
    "        preserve_range=True,\n",
    "        mode='nearest'\n",
    "    ).astype(image2d.dtype)\n",
    "\n",
    "    return image2d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "179a48c9",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# import the \"tiling strategy\" from the cztile package\n",
    "from cztile.fixed_total_area_strategy import AlmostEqualBorderFixedTotalAreaStrategy2D\n",
    "\n",
    "# create the filename for the new CZI image file\n",
    "newczi_tile = os.path.join(os.getcwd(), \"newczi_tilewise.czi\")\n",
    "\n",
    "# create a \"tile\" by specifying the desired tile dimension and the\n",
    "# minimum required overlap between tiles (depends on the processing)\n",
    "tiler = AlmostEqualBorderFixedTotalAreaStrategy2D(\n",
    "    total_tile_width=1600,\n",
    "    total_tile_height=1400,\n",
    "    min_border_width=128\n",
    ")\n",
    "\n",
    "# create CZI instance to read some metadata \n",
    "with pyczi.open_czi(czifile_DAPI_PGC, ReaderFileInputTypes.Curl) as czidoc_r:\n",
    "    # get the size of the bounding rectange for the scence\n",
    "    tiles = tiler.tile_rectangle(czidoc_r.scenes_bounding_rectangle[0])\n",
    "\n",
    "# show the created tile locations\n",
    "for tile in tiles:\n",
    "    print(tile.roi.x, tile.roi.y, tile.roi.w, tile.roi.h)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2ffd9c10",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### Opening a CZI file for reading and write to a new CZI in parallel\n",
    "\n",
    "The idea of tile-wise processing can enable the following approach:\n",
    "\n",
    "* open the \"original CZI image file for reading\n",
    "* tile the image using **czitile**\n",
    "* use the computed tiles as the **roi** parameters to read a \"tile\" of data\n",
    "* apply the desired processing to the read tile\n",
    "* write the processed tile into the new CZI at the correct location using the tile properties"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34d9eadf",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# open an existing CZI file for reading planes and a new CZI to insert\n",
    "# the newly processed chunks into the new CZI file\n",
    "\n",
    "# create to figures\n",
    "fig1, ax1 = plt.subplots(1, 4, figsize=(16, 8))\n",
    "fig2, ax2 = plt.subplots(1, 4, figsize=(16, 8))\n",
    "\n",
    "# counter figures\n",
    "a = 0\n",
    "\n",
    "# open a CZI instance to read and in parallel one to write\n",
    "with pyczi.open_czi(czifile_DAPI_PGC, ReaderFileInputTypes.Curl) as czidoc_r:\n",
    "    with pyczi.create_czi(newczi_tile, exist_ok=True) as czidoc_w:\n",
    "        # loop over all tiles created by the \"tiler\"\n",
    "        for tile in tqdm(tiles):\n",
    "            # read a specific tile from the CZI using the roi parameter\n",
    "            tile2d = czidoc_r.read(\n",
    "                plane={\"C\": 0},\n",
    "                roi=(\n",
    "                    tile.roi.x,\n",
    "                    tile.roi.y,\n",
    "                    tile.roi.w,\n",
    "                    tile.roi.h\n",
    "                )\n",
    "            )\n",
    "\n",
    "            # process the current tile using a function\n",
    "            tile2d_processed = process2d(tile2d, sigma=11)\n",
    "\n",
    "            # show sthe tiles to illustrate the idea\n",
    "            ax1[a].imshow(tile2d[..., 0], interpolation=\"nearest\", cmap=\"gray\", vmin=100, vmax=4000)\n",
    "            ax1[a].set_title(\"Tile: \" + str(a))\n",
    "            ax2[a].imshow(tile2d_processed[..., 0], interpolation=\"nearest\", cmap=\"gray\", vmin=100, vmax=4000)\n",
    "            ax2[a].set_title(\"Tile: \" + str(a))\n",
    "\n",
    "            # write the new CZI file using the processed data\n",
    "            czidoc_w.write(\n",
    "                tile2d_processed,\n",
    "                plane={\"C\": 0},\n",
    "                location=(\n",
    "                    tile.roi.x,\n",
    "                    tile.roi.y\n",
    "                )\n",
    "            )\n",
    "\n",
    "            a += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1cfbc6f3",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import random\n",
    "\n",
    "# read the CZI pixel data\n",
    "with pyczi.open_czi(czifile_DAPI_PGC, ReaderFileInputTypes.Curl) as czidoc_r:\n",
    "    data_orig = czidoc_r.read(plane={\"C\": 0})\n",
    "\n",
    "with pyczi.open_czi(newczi_tile) as czidoc_r:\n",
    "    data_processed = czidoc_r.read(plane={\"C\": 0})\n",
    "\n",
    "# show the results\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 8))\n",
    "ax1.imshow(data_orig[..., 0], interpolation=\"nearest\", cmap=\"gray\")\n",
    "ax1.set_title('Original CZI')\n",
    "ax2.imshow(data_processed[..., 0], interpolation=\"nearest\", cmap=\"gray\")\n",
    "ax2.set_title('Processed CZI - TileWise')\n",
    "\n",
    "# define the tile to be shown\n",
    "tn = 0\n",
    "\n",
    "# add rectangle to plot to visualize the tile\n",
    "ax1.add_patch(\n",
    "    Rectangle((tiles[tn].roi.x, tiles[tn].roi.y), tiles[tn].roi.w, tiles[tn].roi.h,\n",
    "              edgecolor=\"red\",\n",
    "              fill=False,\n",
    "              lw=3,\n",
    "              ls=\"--\")\n",
    ")\n",
    "\n",
    "ax2.add_patch(\n",
    "    Rectangle((tiles[tn].roi.x, tiles[tn].roi.y), tiles[tn].roi.w, tiles[tn].roi.h,\n",
    "              edgecolor=\"red\",\n",
    "              fill=False,\n",
    "              lw=3,\n",
    "              ls=\"--\")\n",
    ")\n",
    "\n",
    "# show the results\n",
    "fig, (ax3, ax4) = plt.subplots(1, 2, figsize=(16, 8))\n",
    "ax3.imshow(data_orig[..., 0], interpolation=\"nearest\", cmap=\"gray\")\n",
    "ax3.set_title('Original CZI')\n",
    "ax4.imshow(data_processed[..., 0], interpolation=\"nearest\", cmap=\"gray\")\n",
    "ax4.set_title('Processed CZI- TileWise')\n",
    "\n",
    "edge_color_opacity = 1.0\n",
    "face_color_opacity = 0.6\n",
    "\n",
    "# visualize a single tile\n",
    "for tile in tiles:\n",
    "    # use random colors to visualize the tiles\n",
    "    r = random.random()\n",
    "    b = random.random()\n",
    "    g = random.random()\n",
    "    edgecolor = (r, g, b, face_color_opacity)\n",
    "    facecolor = (r, g, b, face_color_opacity)\n",
    "\n",
    "    # add rectangle to plot\n",
    "    ax3.add_patch(\n",
    "        Rectangle((tile.roi.x, tile.roi.y), tile.roi.w, tile.roi.h,\n",
    "                  edgecolor=edgecolor,\n",
    "                  fill=False,\n",
    "                  lw=3,\n",
    "                  ls=\"--\")\n",
    "    )\n",
    "\n",
    "    ax4.add_patch(\n",
    "        Rectangle((tile.roi.x, tile.roi.y), tile.roi.w, tile.roi.h,\n",
    "                  facecolor=facecolor,\n",
    "                  fill=True,\n",
    "                  lw=3,\n",
    "                  ls=\"--\")\n",
    "    )\n",
    "\n",
    "    c = + 1"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Raw Cell Format",
  "kernelspec": {
   "display_name": "zen",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}