posts/2024-10-02-primitive-blocks.ipynb

{
 "cells": [
  {
   "cell_type": "raw",
   "id": "43b47afe",
   "metadata": {},
   "source": [
    "---\n",
    "author: \n",
    "    - name: Tobias Pietzsch \n",
    "    - name: John Bogovic\n",
    "    - name: Caleb Hulbert\n",
    "image: 2024-04-02-n5-attribute-paths-thumbnail.png\n",
    "badges: true\n",
    "branch: master\n",
    "categories:\n",
    "- imglib2\n",
    "- algorithm\n",
    "date: '2024-10-02'\n",
    "description: Learn to use algorithms implemented using PrimitiveBlocks\n",
    "output-file: 2024-10-02-primitive-blocks.html\n",
    "title: Primitive Blocks\n",
    "toc: true\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "7e0c32f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "//| code-fold: true\n",
    "\n",
    "%mavenRepo scijava.public https://maven.scijava.org/content/groups/public\n",
    "%maven org.scijava:scijava-common:2.97.0\n",
    "%maven net.imglib2:imglib2:7.1.2\n",
    "%maven net.imglib2:imglib2-algorithm:0.16.0\n",
    "%maven net.imglib2:imglib2-cache:1.0.0-beta-18\n",
    "\n",
    "import net.imglib2.*;\n",
    "import net.imglib2.converter.*;\n",
    "import net.imglib2.img.array.*;\n",
    "import net.imglib2.img.cell.*;\n",
    "import net.imglib2.type.numeric.real.*;\n",
    "import net.imglib2.type.numeric.integer.*;\n",
    "import net.imglib2.view.*;\n",
    "import net.imglib2.util.*;\n",
    "\n",
    "import net.imglib2.blocks.*;\n",
    "import net.imglib2.algorithm.blocks.*;\n",
    "import net.imglib2.algorithm.blocks.downsample.*;\n",
    "import net.imglib2.algorithm.gauss3.*;\n",
    "\n",
    "import net.imglib2.cache.img.*;\n",
    "\n",
    "public static long timeIt( Runnable r ) {\n",
    "    long start = System.currentTimeMillis();\n",
    "    r.run();\n",
    "    return System.currentTimeMillis() - start;\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc9cbd3e",
   "metadata": {},
   "source": [
    "This post summarizes content from the following pull requests:\n",
    "\n",
    "* [PrimitiveBlocks copies RandomAccessible blocks to primitive arrays](https://github.com/imglib/imglib2/pull/330)\n",
    "* [Block algorithms: type conversion and downsampling](https://github.com/imglib/imglib2-algorithm/pull/97)\n",
    "* [BlockSupplier API and block algorithm improvements](https://github.com/imglib/imglib2-algorithm/pull/104)\n",
    "\n",
    "## 1 PrimitiveBlocks\n",
    "\n",
    "Imglib2's `PrimitiveBlocks` adds functionality to extract blocks from RandomAccessible into flat primitive arrays.\n",
    "\n",
    "This will be useful for example to copy data to [Tensor](https://github.com/bioimage-io/model-runner-java/blob/main/src/main/java/io/bioimage/modelrunner/tensor/Tensor.java), to extract blocks for storing into [N5](https://github.com/saalfeldlab/n5), for interfacing CLIJ, for running small algorithm kernels directly on primitive arrays, etc..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "d5c89867",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6\n"
     ]
    }
   ],
   "source": [
    "// Some Img\n",
    "CellImg< UnsignedByteType, ? > cellImg3D = new CellImgFactory(\n",
    "    new UnsignedByteType()).create(60, 60, 60);\n",
    "\n",
    "// A potentially complicated view of an Img\n",
    "RandomAccessible< FloatType > view = Converters.convert(\n",
    "    Views.extendBorder(\n",
    "        Views.hyperSlice(\n",
    "            Views.zeroMin(\n",
    "                Views.rotate( cellImg3D, 1, 0 )\n",
    "            ),\n",
    "            2, 80 )\n",
    "    ),\n",
    "    new RealFloatConverter<>(),\n",
    "    new FloatType()\n",
    ");\n",
    "\n",
    "// Copy the data into a PrimitiveBlock (float[]) in order to optimize processing\n",
    "\n",
    "// allocate destination array\n",
    "float[] data = new float[ 40 * 50 ];\n",
    "\n",
    "// starting position and size for copy operation\n",
    "int[] position = new int[]{ 10, 20 };\n",
    "int[] size = new int[]{ 40, 50 };\n",
    "\n",
    "// copy from 'view' into 'data'\n",
    "PrimitiveBlocks< FloatType > blocks = PrimitiveBlocks.of( view );\n",
    "\n",
    "long timeBlockCopy = timeIt( () -> {\n",
    "    blocks.copy( position, data, size );\n",
    "});\n",
    "System.out.println(\"Block \" + timeBlockCopy + \"ms\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c0feabb",
   "metadata": {},
   "source": [
    "The idea of the optimized copier is:\n",
    "\n",
    "Instead of using `RandomAccess` that checks for every pixel whether it enters a new `Cell`, whether it is out-of-bounds, etc., all these checks are precomputed and then relevant data from each `Cell` is copied in one go. The speedup can be dramatic, in particular if the underlying source data is in a `CellImg`.\n",
    "\n",
    "Some benchmarks included, here is for example results of [`CopyBenchmarkViewPrimitiveBlocks`](https://github.com/imglib/imglib2/blob/1fd92e5385127f25d27695771f001810c6e11109/src/test/java/net/imglib2/blocks/CopyBenchmarkViewPrimitiveBlocks.java):\n",
    "\n",
    "```\n",
    "# JMH version: 1.35\n",
    "# VM version: JDK 17.0.3, OpenJDK 64-Bit Server VM, 17.0.3+7-LTS\n",
    "...\n",
    "\n",
    "Benchmark                                                  (oob)  (permute)  Mode  Cnt   Score   Error  Units\n",
    "CopyBenchmarkViewPrimitiveBlocks.benchmarkLoopBuilder       true       true  avgt    5  12,789 ± 0,285  ms/op\n",
    "CopyBenchmarkViewPrimitiveBlocks.benchmarkLoopBuilder       true      false  avgt    5   9,682 ± 0,152  ms/op\n",
    "CopyBenchmarkViewPrimitiveBlocks.benchmarkLoopBuilder      false       true  avgt    5  14,333 ± 0,099  ms/op\n",
    "CopyBenchmarkViewPrimitiveBlocks.benchmarkLoopBuilder      false      false  avgt    5  12,721 ± 0,123  ms/op\n",
    "CopyBenchmarkViewPrimitiveBlocks.benchmarkPrimitiveBlocks   true       true  avgt    5   0,541 ± 0,010  ms/op\n",
    "CopyBenchmarkViewPrimitiveBlocks.benchmarkPrimitiveBlocks   true      false  avgt    5   0,315 ± 0,024  ms/op\n",
    "CopyBenchmarkViewPrimitiveBlocks.benchmarkPrimitiveBlocks  false       true  avgt    5   0,570 ± 0,013  ms/op\n",
    "CopyBenchmarkViewPrimitiveBlocks.benchmarkPrimitiveBlocks  false      false  avgt    5   0,322 ± 0,008  ms/op\n",
    "```\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5f4dac6",
   "metadata": {},
   "source": [
    "If a source `RandomAccessible` cannot be understood, `PrimitiveBlocks.of(...)` will return\n",
    "a fall-back implementation (based on `LoopBuilder`). With the optional `OnFallback` argument\n",
    "of `PrimitiveBlocks.of(...)` it can be configured whether fall-back should be\n",
    "\n",
    "* silently accepted (`ACCEPT`),\n",
    "* a warning should be printed (`WARN`) -- the default,\n",
    "* or an IllegalArgumentException thrown (`FAIL`).\n",
    "    * The warning/exception message explains why the source `RandomAccessible` requires fall-back.\n",
    "\n",
    "The only really un-supported case is if the pixel type `T` does not map one-to-one to a \n",
    "primitive type. For example, `ComplexDoubleType` or `Unsigned4BitType` are not supported \n",
    "(at least not yet).\n",
    "\n",
    "`PrimitiveBlocks.copy` is single-threaded, the idea being to parallelize over blocks instead\n",
    "of  the copying within a block.  `PrimitiveBlocks` is not thread-safe in general, but has\n",
    "a method `threadSafe()` to obtain a thread-safe instance (implemented using `ThreadLocal` copies). \n",
    "For example,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "95abff23",
   "metadata": {},
   "outputs": [],
   "source": [
    "PrimitiveBlocks< FloatType > blocks = PrimitiveBlocks.of( view ).threadSafe();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e05adc37",
   "metadata": {},
   "source": [
    "## 2 Block algorithms"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39a2f21e",
   "metadata": {},
   "source": [
    "1. A framework for implementing block-processing algorithms that build off the `PrimitiveBlocks` utility introduced in imglib2 6.2.0. Package `net.imglib2.algorithm.block` contains helpers to define and chain block operators. In particular,  the `BlockProcessor<I,O>` interface specifies an algorithm that computes values in a flattened primitive output array (type `O`) from values in a flattened primitive input array (type `I`). To avoid mistakes in unchecked casting of primitive array types, `BlockProcessor` should be typically wrapped in `UnaryBlockOperator<S,T>` where `S` and `T` are source and target ImgLib2 `RealType` corresponding to the respective `I` and `O` array types. `UnaryBlockOperator` can be chained using `UnaryBlockOperator.andThen(...)`.\n",
    "2. An operator to *convert* between common `RealTypes` (signed and unsigned, with and without clamping). This code is partially auto-generated using the `bin/generate.groovy` script lifted from scijava-ops.\n",
    "3. An operator to *downsample*, by factor 2, blocks of common `RealType`. It can be specified whether intermediate calculations should happen in `FLOAT` or `DOUBLE` precision (or `AUTO`matically determined). There are versions for downsampling 2x2x2... blocks with half-pixel offset and 3x3x3... blocks on pixel centers.\n",
    "\n",
    "Together with `PrimitiveBlocks` from imglib2 core, the downsampling algorithm achieves speedup of factor ~4 over the downsampling used in bigdataviewer-core, and speedup of factor ~10 to ~30 over the downsampling used in BigStitcher (depending on Java version, block size, etc). (This is on `ArrayImg` inputs. Speedups are probably even larger for `CellImg` or `PlanarImg` inputs, where `PrimitveBlocks` has a larger impact.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "957e0f81",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Gaussian downsampling took 216ms\n",
      "Blocks downsampling 12ms\n"
     ]
    }
   ],
   "source": [
    "public static RandomAccessibleInterval<FloatType> downsampleGauss(RandomAccessibleInterval<FloatType> img) { \n",
    "    \n",
    "    final boolean[] downsampleInDim = { true, true, true };\n",
    "    final long[] downsampledDimensions = Downsample.getDownsampledDimensions( \n",
    "        img.dimensionsAsLongArray(), downsampleInDim );\n",
    "        \n",
    "    final CellImg<FloatType, ?> tgt = new CellImgFactory<FloatType>(new FloatType(), cellDimensions)\n",
    "        .create(downsampledDimensions);\n",
    "\n",
    "    Gauss3.gauss(1.0, Views.extendBorder(img), tgt);\n",
    "    return tgt;\n",
    "}\n",
    "\n",
    "public static RandomAccessibleInterval<FloatType> downsampleBlocks(RandomAccessibleInterval<FloatType> img) { \n",
    "    \n",
    "    final boolean[] downsampleInDim = { true, true, true };\n",
    "    final long[] downsampledDimensions = Downsample.getDownsampledDimensions( \n",
    "        img.dimensionsAsLongArray(), downsampleInDim );\n",
    "        \n",
    "    final BlockSupplier< FloatType > blocks = BlockSupplier.of( Views.extendBorder( img ) );\n",
    "    final CachedCellImg< FloatType, ? > downsampled = BlockAlgoUtils.cellImg(\n",
    "                    blocks.andThen( Downsample.downsample(\n",
    "                            ComputationType.DOUBLE,\n",
    "                            Downsample.Offset.CENTERED )\n",
    "                    ),\n",
    "                    downsampledDimensions, cellDimensions );\n",
    "\n",
    "    return downsampled;\n",
    "}\n",
    "\n",
    "final int[] cellDimensions = { 64, 64, 64 };\n",
    "final CellImg<FloatType, ?> img = new CellImgFactory<FloatType>(new FloatType(), cellDimensions)\n",
    "                .create(256, 256, 256);\n",
    "img.forEach(x -> x.setReal(Math.random()));\n",
    "\n",
    "long timeGauss = timeIt( () -> {\n",
    "    downsampleGauss( img );\n",
    "});\n",
    "System.out.println(\"Gaussian downsampling took \" + timeGauss + \"ms\");\n",
    "\n",
    "long timeBlocks = timeIt( () -> {\n",
    "    downsampleBlocks( img );\n",
    "});\n",
    "System.out.println(\"Blocks downsampling \" + timeBlocks + \"ms\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "584d88c7",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Java",
   "language": "java",
   "name": "java"
  },
  "language_info": {
   "codemirror_mode": "java",
   "file_extension": ".jshell",
   "mimetype": "text/x-java-source",
   "name": "Java",
   "pygments_lexer": "java",
   "version": "17.0.12+7-Ubuntu-1ubuntu222.04"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}