docs/examples/examples_h5writing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " ## Writing workflow\n",
    " The library accepts any 2d or 3d numpy arrays: (z,y,x) axis order as `stack`, or (y,x) as `plane`. \n",
    " The input array is converted to `uint16` automatically.\n",
    " \n",
    " 1. A writer object is created: it opens a new h5 file and requires info about setups and saving options: \n",
    " number of setup attributes (e.g. channels, angles), compression, subsampling. \n",
    " 2. Numpy stacks (3d arrays) can be appended to the h5 file \n",
    " as views by `BdvWriter.append_view(stack, ...)`. \n",
    " Stacks can be normal (whole stack appended at once), or virtual, appended by plane or by substack.\n",
    " The time point `time` and attributes such as `channel`, `angle` etc must be specified \n",
    " for each view.\n",
    " \n",
    " 3. An XML file with all meta information is created in the end by `BdvWriter.write_xml(...)`.\n",
    "  \n",
    " 4. Writing is finalized by calling `BdvWriter.close()`.\n",
    " \n",
    "**Todo:** \n",
    "- examples of reading and cropping views in existing file using `npy2bdv.BdvEditor()` class.\n",
    "- examples of appending and reading transforms to/from XML file. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "import sys\n",
    "import os\n",
    "import numpy as np\n",
    "import npy2bdv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# optional: upgrade to latest version\n",
    "! pip install -U npy2bdv "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# optional: check installed version\n",
    "! pip show npy2bdv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_test_image(dim_yx, iz, nz):\n",
    "    \"\"\"Gaussian blob spanning the whole range of uint16 type\"\"\"\n",
    "    x = np.linspace(-3, 3, dim_yx[1])\n",
    "    y = np.linspace(-3, 3, dim_yx[0])\n",
    "    sigma = 1.0 - abs(iz - nz/2) / nz\n",
    "    x, y = np.meshgrid(x, y)\n",
    "    return (65535 * np.exp(- ((x ** 2) + (y ** 2)) / (2 * sigma**2) )).astype('uint16')\n",
    "\n",
    "def generate_test_stack(nz, ny, nx):\n",
    "    stack = np.empty((nz, ny, nx))\n",
    "    for z in range(nz):\n",
    "        stack[z, :, :] = generate_test_image((ny, nx), z, nz)\n",
    "    return stack.astype('uint16')\n",
    "    \n",
    "examples_dir = \"./example_files/\"\n",
    "if not os.path.exists(examples_dir):\n",
    "    os.mkdir(examples_dir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Basic stack writing\n",
    "\n",
    "### Stacks (views) are written as raw data.  \n",
    "No compression, no downsampling. This mode ensures highest speed, but may be sub-optimal for post-processing and data storage that requires compression and downsampling. See Example 5 for more options."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dataset in ./example_files/ex1_t2_ch2_illum2_angle2_raw.h5\n",
      "Wall time: 9.03 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "nz, ny, nx = 50, 1024, 2048\n",
    "stack = generate_test_stack(nz, ny, nx)\n",
    "fname = examples_dir + \"ex1_t2_ch2_illum2_angle2_raw.h5\"\n",
    "bdv_writer = npy2bdv.BdvWriter(fname, nchannels=2, nilluminations=2, nangles=2)\n",
    "\n",
    "bdv_writer.set_attribute_labels('channel', ('488', '561'))\n",
    "bdv_writer.set_attribute_labels('illumination', ('L', 'R'))\n",
    "bdv_writer.set_attribute_labels('angle', ('90', '180'))\n",
    "\n",
    "for t in range(2):\n",
    "    for i_ch in range(2):\n",
    "        for i_illum in range(2):\n",
    "            for i_angle in range(2):\n",
    "                bdv_writer.append_view(stack, time=t, channel=i_ch, illumination=i_illum, angle=i_angle)\n",
    "\n",
    "bdv_writer.write_xml()\n",
    "bdv_writer.close()\n",
    "print(f\"dataset in {fname}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Writing speed test\n",
    "\n",
    "Speed test for raw data writing, 20 time points and 2 channels. File size is 17 GB."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "clock on: 1721819333\n",
      "clock off: 1721819398\n",
      "H5 mean writing time per stack: 0.046 sec.\n",
      "H5 mean writing speed: 9042 MB/s\n",
      "dataset in ./example_files/ex2_t20_chan2.h5\n"
     ]
    }
   ],
   "source": [
    "nt, nc, nz, ny, nx = 20, 2, 50, 2048, 2048\n",
    "stack = generate_test_stack(nz, ny, nx)\n",
    "\n",
    "fname = examples_dir + \"ex2_t20_chan2.h5\"\n",
    "bdv_writer = npy2bdv.BdvWriter(fname, nchannels=nc)\n",
    "\n",
    "start_time = time.time()\n",
    "print('clock on:',int(start_time))\n",
    "\n",
    "for ichannel in range(nc):\n",
    "    for itime in range(nt):\n",
    "        start_time = time.time()\n",
    "        bdv_writer.append_view(stack, time=itime, channel=ichannel)\n",
    "\n",
    "bdv_writer.write_xml()\n",
    "bdv_writer.close()\n",
    "\n",
    "stop_time = time.time()\n",
    "print('clock off:', int(stop_time))\n",
    "time_per_stack = (stop_time - start_time) / (nt*nc)\n",
    "\n",
    "print(f\"H5 mean writing time per stack: {time_per_stack:1.3f} sec.\")\n",
    "print(f\"H5 mean writing speed: {int(sys.getsizeof(stack) / time_per_stack / 1e6)} MB/s\")\n",
    "print(f\"dataset in {fname}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Warning:** This speed should not be trusted: time variables `start_time` and `stop_time` change during cell execution, see below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "clock on: 1721819396\n",
      "clock off: 1721819398\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "1.855459213256836"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print('clock on:',int(start_time))\n",
    "print('clock off:', int(stop_time))\n",
    "stop_time - start_time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Writing with affine transformations defined in XML file\n",
    "\n",
    "Affine transformations define translation, rotation, scaling, and shear.\n",
    "\n",
    "In this example we write 1 time point and 1 channel with 10-px shear transformation along X axis. With non-isotropic voxel size calibration.\n",
    "      \n",
    "The affine transformation defined in XML file is automatically applied every time you open the dataset in BigDataViewer/BigStitcher."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sheared dataset in ./example_files/ex3_t1_ch1_shear.h5\n"
     ]
    }
   ],
   "source": [
    "shear_x_px = 10\n",
    "affine_matrix = np.array(((1.0, 0.0, -shear_x_px, 0.0),\n",
    "                          (0.0, 1.0, 0.0, 0.0),\n",
    "                          (0.0, 0.0, 1.0, 0.0)))\n",
    "\n",
    "fname = examples_dir + \"ex3_t1_ch1_shear.h5\"\n",
    "bdv_writer = npy2bdv.BdvWriter(fname)\n",
    "bdv_writer.append_view(stack, time=0, channel=0,\n",
    "                       m_affine=affine_matrix,\n",
    "                       name_affine=\"shearing transformation\",\n",
    "                       calibration=(1, 1, 1))\n",
    "bdv_writer.write_xml()\n",
    "bdv_writer.close()\n",
    "print(f\"sheared dataset in {fname}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.1 Writing multiple tiles on a grid\n",
    "\n",
    "Here, affine transformation matrix defines translation of individual tiles in the global coordinate system.\n",
    "Tile position is the position of its lower left corner. The last (4-th) column of the affine matrix defines the translation terms (x,y,z).\n",
    "\n",
    "Optionally, make some tiles missing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      " File with multiple tiles: ./example_files/ex3p1_tiles.h5\n"
     ]
    }
   ],
   "source": [
    "DEBUG_MODE = False\n",
    "MISSING_TILES = False\n",
    "\n",
    "n_tiles_x = 4\n",
    "n_tiles_y = 2\n",
    "\n",
    "tile_w, tile_h = stack.shape[2] // n_tiles_x, stack.shape[1] // n_tiles_y\n",
    "\n",
    "tile_position_yx = ((0,0), (0, tile_w), (tile_h, 0), (tile_h, tile_w))\n",
    "\n",
    "unit_matrix = np.array(((1.0, 0.0, 0.0, 0.0), # change the 4. value for x_translation (px)\n",
    "                        (0.0, 1.0, 0.0, 0.0), # change the 4. value for y_translation (px)\n",
    "                        (0.0, 0.0, 1.0, 0.0)))# change the 4. value for z_translation (px)\n",
    "\n",
    "fname = examples_dir + \"ex3p1_tiles.h5\"\n",
    "bdv_writer = npy2bdv.BdvWriter(fname, ntiles=n_tiles_x*n_tiles_y)\n",
    "for j in range(n_tiles_y):\n",
    "    for i in range(n_tiles_x):\n",
    "        tile_stack = stack[:, j*tile_h:(j+1)*tile_h, i*tile_w:(i+1)*tile_w]\n",
    "        tile_index = i + j*n_tiles_x\n",
    "        affine_matrix = unit_matrix\n",
    "        affine_matrix[0,3] = i * tile_w # x-translation\n",
    "        affine_matrix[1,3] = j * tile_h # y-translation\n",
    "        if DEBUG_MODE:\n",
    "            print(f\"Tile index: {tile_index}\")\n",
    "            print(f\"Transform matrix: \\n {affine_matrix}\")\n",
    "        if MISSING_TILES and tile_index in (1,2): # make two tiles missing \n",
    "            pass\n",
    "        else:\n",
    "            bdv_writer.append_view(tile_stack, time=0, \n",
    "                                   tile=tile_index, \n",
    "                                   m_affine=affine_matrix, \n",
    "                                   name_affine=f\"tile {tile_index} translation\")\n",
    "bdv_writer.write_xml()\n",
    "bdv_writer.close()\n",
    "print(f\"\\n File with multiple tiles: {fname}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Writing with experiment metadata\n",
    "Dataset with 1 time point and 1 channel with voxel size, exposure, camera and microscope properties stored in XML."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dataset is in ./example_files/ex4_cam_props.h5\n"
     ]
    }
   ],
   "source": [
    "fname = examples_dir + \"ex4_cam_props.h5\"\n",
    "bdv_writer = npy2bdv.BdvWriter(fname)\n",
    "bdv_writer.append_view(stack, time=0, \n",
    "                       voxel_size_xyz=(1, 1, 5), voxel_units='um',\n",
    "                       exposure_time=10, exposure_units='ms')\n",
    "bdv_writer.write_xml(camera_name=\"Hamamatsu OrcaFlash100\",\n",
    "                          microscope_name='Superscope',\n",
    "                          user_name='nvladimus')\n",
    "bdv_writer.close()\n",
    "print(f\"dataset is in {fname}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5 Writing with downsampling and compression\n",
    "### 5.1 Stacks are written with downsampling and compression on-the-fly. \n",
    "This mode is usually too slow for image writing during experiment, but offers user convenience for offline processing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dataset in ./example_files/ex5_t2_ch2_illum2_angle2_realtime_downsamp.h5\n",
      "Wall time: 24.9 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "fname = examples_dir + \"ex5_t2_ch2_illum2_angle2_realtime_downsamp.h5\"\n",
    "bdv_writer = npy2bdv.BdvWriter(fname, nchannels=2, nilluminations=2, nangles=2, \n",
    "                               subsamp=((1, 4, 4),    (2, 8, 8),    (4, 16, 16)),\n",
    "                               blockdim=((8, 32, 32), (16, 32, 32), (32, 32, 32)),\n",
    "                               compression='gzip')\n",
    "\n",
    "bdv_writer.set_attribute_labels('channel', ('488', '561'))\n",
    "bdv_writer.set_attribute_labels('illumination', ('L', 'R'))\n",
    "bdv_writer.set_attribute_labels('angle', ('90', '180'))\n",
    "\n",
    "for t in range(2):\n",
    "    for i_ch in range(2):\n",
    "        for i_illum in range(2):\n",
    "            for i_angle in range(2):\n",
    "                bdv_writer.append_view(stack, time=t, channel=i_ch, illumination=i_illum, angle=i_angle)\n",
    "\n",
    "bdv_writer.write_xml()\n",
    "bdv_writer.close()\n",
    "print(f\"dataset in {fname}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.2 Stacks are written as raw data followed by lazy downsampling and compression. \n",
    "This mode combines highest speed during acqusition (raw data writing) with delayed (lazy) downsampling and compression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "time points:   0%|                                                                               | 0/2 [00:00<?, ?it/s]\n",
      "views:   0%|                                                                                     | 0/8 [00:00<?, ?it/s]\u001b[A"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Raw data writing time, total: 4.94 sec.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n",
      "views:  12%|█████████▋                                                                   | 1/8 [00:01<00:12,  1.72s/it]\u001b[A\n",
      "views:  25%|███████████████████▎                                                         | 2/8 [00:03<00:09,  1.64s/it]\u001b[A\n",
      "views:  38%|████████████████████████████▉                                                | 3/8 [00:04<00:07,  1.59s/it]\u001b[A\n",
      "views:  50%|██████████████████████████████████████▌                                      | 4/8 [00:06<00:06,  1.56s/it]\u001b[A\n",
      "views:  62%|████████████████████████████████████████████████▏                            | 5/8 [00:07<00:04,  1.54s/it]\u001b[A\n",
      "views:  75%|█████████████████████████████████████████████████████████▊                   | 6/8 [00:09<00:03,  1.54s/it]\u001b[A\n",
      "views:  88%|███████████████████████████████████████████████████████████████████▍         | 7/8 [00:11<00:01,  1.56s/it]\u001b[A\n",
      "views: 100%|█████████████████████████████████████████████████████████████████████████████| 8/8 [00:12<00:00,  1.57s/it]\u001b[A\n",
      "time points:  50%|███████████████████████████████████▌                                   | 1/2 [00:12<00:12, 12.57s/it]\n",
      "views:   0%|                                                                                     | 0/8 [00:00<?, ?it/s]\u001b[A\n",
      "views:  12%|█████████▋                                                                   | 1/8 [00:01<00:10,  1.52s/it]\u001b[A\n",
      "views:  25%|███████████████████▎                                                         | 2/8 [00:03<00:09,  1.52s/it]\u001b[A\n",
      "views:  38%|████████████████████████████▉                                                | 3/8 [00:04<00:07,  1.52s/it]\u001b[A\n",
      "views:  50%|██████████████████████████████████████▌                                      | 4/8 [00:06<00:06,  1.52s/it]\u001b[A\n",
      "views:  62%|████████████████████████████████████████████████▏                            | 5/8 [00:07<00:04,  1.52s/it]\u001b[A\n",
      "views:  75%|█████████████████████████████████████████████████████████▊                   | 6/8 [00:09<00:03,  1.52s/it]\u001b[A\n",
      "views:  88%|███████████████████████████████████████████████████████████████████▍         | 7/8 [00:10<00:01,  1.52s/it]\u001b[A\n",
      "views: 100%|█████████████████████████████████████████████████████████████████████████████| 8/8 [00:12<00:00,  1.52s/it]\u001b[A\n",
      "time points: 100%|███████████████████████████████████████████████████████████████████████| 2/2 [00:24<00:00, 12.36s/it]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Downsampling & compression & pyramid writing time, total: 24.75 sec.\n",
      "dataset in ./example_files/ex5_t2_ch2_illum2_angle2_lazy_downsamp.h5\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "fname = examples_dir + \"ex5_t2_ch2_illum2_angle2_lazy_downsamp.h5\"\n",
    "bdv_writer = npy2bdv.BdvWriter(fname, nchannels=2, nilluminations=2, nangles=2)\n",
    "\n",
    "bdv_writer.set_attribute_labels('channel', ('488', '561'))\n",
    "bdv_writer.set_attribute_labels('illumination', ('L', 'R'))\n",
    "bdv_writer.set_attribute_labels('angle', ('90', '180'))\n",
    "\n",
    "start_time = time.time()\n",
    "for t in range(2):\n",
    "    for i_ch in range(2):\n",
    "        for i_illum in range(2):\n",
    "            for i_angle in range(2):\n",
    "                bdv_writer.append_view(stack, time=t, channel=i_ch, illumination=i_illum, angle=i_angle)\n",
    "\n",
    "print(f\"Raw data writing time, total: {(time.time() - start_time):2.2f} sec.\")\n",
    "\n",
    "start_time = time.time()\n",
    "bdv_writer.create_pyramids(subsamp=((1, 4, 4),    (2, 8, 8),    (4, 16, 16)), \n",
    "                           blockdim=((8, 32, 32), (16, 32, 32), (32, 32, 32)),\n",
    "                          compression='gzip')    \n",
    "print(f\"Downsampling & compression & pyramid writing time, total: {(time.time() - start_time):2.2f} sec.\")\n",
    "\n",
    "bdv_writer.write_xml()\n",
    "bdv_writer.close()\n",
    "print(f\"dataset in {fname}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Writing virtual stacks that are too big to fit RAM\n",
    "### 6.1 Writing by plane\n",
    "Dataset has 1 time point, 2 channels, and large virtual stack, written by plane."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "H5 mean writing speed: 216 MB/s\n",
      "H5 mean writing time per plane: 0.039 sec.\n",
      "virtual stack in ./example_files/ex6_t1_ch2_virtual_by_plane.h5\n",
      "Wall time: 31.1 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "nc, nz, ny, nx = 2, 400, 2048, 2048\n",
    "test_image = generate_test_image((ny, nx), nz/2, nz)\n",
    "fname = examples_dir + \"ex6_t1_ch2_virtual_by_plane.h5\"\n",
    "bdv_writer = npy2bdv.BdvWriter(fname, nchannels=2, blockdim=((1, 256, 256),))\n",
    "\n",
    "start_time = time.time()\n",
    "for i_ch in range(nc):\n",
    "    bdv_writer.append_view(stack=None, virtual_stack_dim=(nz,ny,nx), time=0, channel=i_ch)\n",
    "    for z in range(nz):\n",
    "        bdv_writer.append_plane(plane=test_image, z=z, time=0, channel=i_ch)\n",
    "\n",
    "bdv_writer.write_xml()\n",
    "bdv_writer.close()\n",
    "\n",
    "time_total = time.time() - start_time\n",
    "ave_time_per_plane = time_total/(nc*nz)\n",
    "print(f\"H5 mean writing speed: {int(sys.getsizeof(test_image)*nc*nz / time_total / 1e6)} MB/s\")\n",
    "print(f\"H5 mean writing time per plane: {ave_time_per_plane:1.3f} sec.\")\n",
    "print(f\"virtual stack in {fname}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.2 Virtual stack, writing by substack"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "clock started\n",
      "H5 mean writing speed: 325 MB/s\n",
      "H5 mean writing time per plane: 0.026 sec.\n",
      "virtual stack in ./example_files/ex6_t1_ch2_virtual_by_substack.h5\n",
      "Wall time: 2min 3s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "nc, nz, ny, nx, n_substacks = 2, 400, 2048, 2048, 4\n",
    "stack = generate_test_stack(nz, ny, nx)\n",
    "fname = examples_dir + \"ex6_t1_ch2_virtual_by_substack.h5\"\n",
    "bdv_writer = npy2bdv.BdvWriter(fname, nchannels=2, blockdim=((32, 32, 32),))\n",
    "\n",
    "print('clock started')\n",
    "start_time = time.time()\n",
    "#initialize virtual stacks\n",
    "for i_ch in range(nc):\n",
    "    bdv_writer.append_view(stack=None, virtual_stack_dim=(nz,ny,nx), time=0, channel=i_ch)\n",
    "\n",
    "# populate virtual stacks \n",
    "for i_ch in range(nc):\n",
    "    for isub in range(n_substacks):\n",
    "        zslice = slice(isub*(nz//n_substacks), (isub+1)*(nz//n_substacks))\n",
    "        bdv_writer.append_substack(substack=stack[zslice, :, :],\n",
    "                                   z_start=zslice.start,\n",
    "                                   channel=i_ch)\n",
    "\n",
    "bdv_writer.write_xml()\n",
    "bdv_writer.close()\n",
    "\n",
    "time_total = time.time() - start_time\n",
    "ave_time_per_plane = time_total/(nc*nz)\n",
    "\n",
    "print(f\"H5 mean writing speed: {int(sys.getsizeof(stack) * nc / time_total / 1e6)} MB/s\")\n",
    "print(f\"H5 mean writing time per plane: {ave_time_per_plane:1.3f} sec.\")\n",
    "print(f\"virtual stack in {fname}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.a Virtual stacks by plane, with downsampling\n",
    "\n",
    "Subsampling allows faster viewing in BDV, because only the necessary level of details is loaded into RAM, instead of full resolution image.\n",
    "\n",
    "When writing by plane, the subsampling level Z must be 1, so there is no way to subsample in Z on the fly. For example, `subsamp=((1, 1, 1), (1, 4, 4), (1, 8, 8))`.\n",
    "\n",
    "The newly generated file is slightly bigger (by ca. 8%), but the processing time increases 5x.\n",
    "\n",
    "Compression is supported."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "INFO: blockdim levels (1) < subsamp levels (3): First-level block size (1, 128, 128) will be used for all levels\n",
      "clock started\n",
      "H5 mean writing speed: 161 MB/s\n",
      "H5 mean writing time per plane: 0.052 sec.\n",
      "virtual stack in ./example_files/ex6a_t1_ch1_virtual_by_plane_subsamp.h5\n",
      "Wall time: 41.8 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "nc, nz, ny, nx = 2, 400, 2048, 2048\n",
    "test_image = generate_test_image((ny, nx), nz/2, nz)\n",
    "fname = examples_dir + \"ex6a_t1_ch1_virtual_by_plane_subsamp.h5\"\n",
    "bdv_writer = npy2bdv.BdvWriter(fname, nchannels=2,\n",
    "                               blockdim=((1, 128, 128),),\n",
    "                               subsamp=((1, 1, 1), (1, 4, 4), (1, 8, 8)))\n",
    "print('clock started')\n",
    "start_time = time.time()\n",
    "for i_ch in range(nc):\n",
    "    bdv_writer.append_view(stack=None, virtual_stack_dim=(nz,ny,nx), time=0, channel=i_ch)\n",
    "    for z in range(nz):\n",
    "        bdv_writer.append_plane(plane=test_image, z=z, time=0, channel=i_ch)\n",
    "\n",
    "bdv_writer.write_xml()\n",
    "bdv_writer.close()\n",
    "\n",
    "time_total = time.time() - start_time\n",
    "ave_time_per_plane = time_total/(nc*nz)\n",
    "\n",
    "print(f\"H5 mean writing speed: {int(sys.getsizeof(test_image) * nc * nz / time_total / 1e6)} MB/s\")\n",
    "print(f\"H5 mean writing time per plane: {ave_time_per_plane:1.3f} sec.\")\n",
    "print(f\"virtual stack in {fname}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7. Missing views, normal stack\n",
    "Missing views are detected automatically and marked as such in XML file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dataset with missing views in ./example_files/ex7_missing_views.h5\n"
     ]
    }
   ],
   "source": [
    "nz, ny, nx = 50, 1024, 2048\n",
    "stack = generate_test_stack(nz, ny, nx)\n",
    "fname = examples_dir + \"ex7_missing_views.h5\"\n",
    "bdv_writer = npy2bdv.BdvWriter(fname, nchannels=2)\n",
    "bdv_writer.append_view(stack, time=0, channel=0)\n",
    "bdv_writer.append_view(stack, time=1, channel=1)\n",
    "bdv_writer.write_xml()\n",
    "bdv_writer.close()\n",
    "print(f\"dataset with missing views in {fname}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8. Missing views, virtual stack by plane"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dataset with missing views in ./example_files/ex8_virtual_stack_missing_views.h5\n"
     ]
    }
   ],
   "source": [
    "nz, ny, nx = 50, 1024, 2048\n",
    "stack = generate_test_stack(nz, ny, nx)\n",
    "fname = examples_dir + \"ex8_virtual_stack_missing_views.h5\"\n",
    "bdv_writer = npy2bdv.BdvWriter(fname, nchannels=2)\n",
    "bdv_writer.append_view(stack=None, virtual_stack_dim=(nz, ny, nx), time=0, channel=0)\n",
    "bdv_writer.append_view(stack=None, virtual_stack_dim=(nz, ny, nx), time=1, channel=1)\n",
    "\n",
    "for z in range(nz):\n",
    "    bdv_writer.append_plane(plane=test_image, z=z, time=0, channel=0)\n",
    "    bdv_writer.append_plane(plane=test_image, z=z, time=1, channel=1)\n",
    "\n",
    "bdv_writer.write_xml()\n",
    "bdv_writer.close()\n",
    "print(f\"dataset with missing views in {fname}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Clean up the data files generated"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "import shutil\n",
    "if os.path.exists(examples_dir):\n",
    "    shutil.rmtree(examples_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}