Merge remote-tracking branch 'origin/main'

- add pytest workflow

.github/workflows/publish_crates.yml

@@ -1,6 +1,6 @@
-name: Publish
+name: Publish_rs
 
-on: [workflow_dispatch]
+on: [ workflow_dispatch ]
 
 permissions:
   contents: read
@@ -10,23 +10,27 @@ jobs:
     name: Publish (crates.io)
     runs-on: ubuntu-latest
     timeout-minutes: 25
+    env:
+      CARGO_HOME: ${{ github.workspace }}/.cargo
     steps:
       - uses: actions/checkout@v6
+        with:
+          submodules: recursive
       - name: Restore cache
         uses: actions/cache/restore@v4
         with:
-          path: |
-            ~/.cargo
+          path: ${{ env.CARGO_HOME }}
           key: cache-ubuntu-cargo-publish
 
       - name: Install Rust
         run: |-
-          export PATH="$HOME/.cargo/bin:$PATH"
+          export PATH="$CARGO_HOME/bin:$PATH"
           if ! command -v rustc >/dev/null 2>&1; then
             curl --proto '=https' --tlsv1.2 -sSf https://sh.rustup.rs | sh -s -- -y
           else
             rustup update
           fi
+          echo "$CARGO_HOME/bin" >> "$GITHUB_PATH"
           cargo install cargo-release
         shell: bash
 
@@ -46,9 +50,9 @@ jobs:
       # which is reasonable, since this job only runs after the Linux, Windows, and WASM builds
       # have passed.
       - name: "cargo release publish"
+        env:
+          CARGO_REGISTRY_TOKEN: ${{ secrets.CRATES_IO_API_TOKEN }}
         run: |-
-          export PATH="$HOME/.osxcross/bin:$PATH"
-          cargo login ${{ secrets.CRATES_IO_API_TOKEN }}
           cargo release \
             publish \
             --workspace \
@@ -61,6 +65,5 @@ jobs:
       - name: Store cache
         uses: actions/cache/save@v4
         with:
-          path: |
-            ~/.cargo
+          path: ${{ env.CARGO_HOME }}
           key: cache-ubuntu-cargo-publish

.github/workflows/publish_pypi.yml

@@ -1,6 +1,6 @@
-name: CI
+name: Publish_py
 
-on: [workflow_dispatch]
+on: [ workflow_dispatch ]
 
 permissions:
   contents: read
@@ -9,8 +9,14 @@ jobs:
   pypi_publish:
     name: Publish (pypi.org)
     runs-on: ubuntu-latest
+    env:
+      CARGO_HOME: ${{ github.workspace }}/.cargo
+      RUSTUP_HOME: ${{ github.workspace }}/.rustup
+      RUSTUP_TOOLCHAIN: stable
     steps:
       - uses: actions/checkout@v6
+        with:
+          submodules: recursive
       - uses: actions/setup-python@v5
         with:
           python-version: 3.x
@@ -23,33 +29,62 @@ jobs:
             ~/.cache/pip-wheel
             ~/.cache/sccache
             ~/.cache/cargo-xwin
-            ~/.cargo
+            ${{ env.CARGO_HOME }}
             ~/.osxcross
           key: cache-ubuntu-maturin-cross-compile
 
-      - name: Install llvm
+      - name: Install build dependencies
         run: |
-          if ! command -v llvm-dlltool >/dev/null 2>&1; then
-            sudo apt update
-            sudo apt install -y llvm
-          fi
+          # the stock clang on noble (18) is too old for the 2025+ MSVC STL that
+          # cargo-xwin downloads (yvals_core.h: "error STL1000: expected Clang
+          # 19.0.0 or newer"), so take a newer clang toolchain from apt.llvm.org.
+          # apt.llvm.org sometimes half-publishes a new version, in which case the
+          # unversioned meta packages become uninstallable ("clang: Depends:
+          # clang-NN ... not installable"); so pick the newest versioned set that
+          # actually resolves (a dry-run checks the whole dependency closure) -
+          # nothing is hardcoded, this tracks the newest working release.
+          sudo apt update
+          echo "deb [trusted=yes] https://apt.llvm.org/noble/ llvm-toolchain-noble main" | sudo tee /etc/apt/sources.list.d/llvm-toolchain-noble.list >/dev/null
+          sudo apt update
+          sudo apt install -y build-essential cmake
+          NEWEST=""
+          for name in $(apt-cache search '^clang-[0-9]+$' | awk '{print $1}' | sort -rV); do
+            v=${name#clang-}
+            if sudo apt-get install --dry-run -y "clang-$v" "llvm-$v" "lld-$v" "libclang-$v-dev" >/dev/null 2>&1; then
+              NEWEST=$v
+              break
+            fi
+          done
+          [ -n "$NEWEST" ] || { echo "no installable clang toolchain found in apt.llvm.org" >&2; exit 1; }
+          echo "using clang/llvm toolchain $NEWEST"
+          sudo apt install -y "clang-$NEWEST" "llvm-$NEWEST" "lld-$NEWEST" "libclang-$NEWEST-dev"
+          # expose the toolchain's own bin/lib dirs (clang, clang-cl, lld-link,
+          # llvm-lib, libclang) to the later steps
+          echo "LLVM_BIN=/usr/lib/llvm-$NEWEST/bin" >> "$GITHUB_ENV"
+          echo "LIBCLANG_PATH=/usr/lib/llvm-$NEWEST/lib" >> "$GITHUB_ENV"
+          export PATH="/usr/lib/llvm-$NEWEST/bin:$PATH"
+          clang --version | head -1
+          command -v clang clang-cl lld-link llvm-lib
         shell: bash
 
       - name: Install Rust
         run: |
-          export PATH="$HOME/.cargo/bin:$PATH"
-          if ! command -v rustc >/dev/null 2>&1; then
+          export PATH="$CARGO_HOME/bin:$PATH"
+          if ! command -v rustup >/dev/null 2>&1; then
             curl --proto '=https' --tlsv1.2 -sSf https://sh.rustup.rs | sh -s -- -y
-          else
-            rustup update
           fi
+          rustup toolchain install stable --profile minimal
+          rustup default stable
+          echo "$CARGO_HOME/bin" >> "$GITHUB_PATH"
+          cargo --version
         shell: bash
 
       - name: Install sccache and maturin
         run: |
-          export PATH="$HOME/.cargo/bin:$PATH"
-          python -m pip install --upgrade pip
-          pip install maturin ziglang
+          python -m pip install --upgrade pip maturin ziglang cargo-zigbuild
+          # in case an older copy was cargo-installed into the CARGO_HOME cache,
+          # prevent it from shadowing the pip-installed one
+          rm -f "${CARGO_HOME:?}/bin/cargo-zigbuild" || true
           if ! command -v sccache >/dev/null 2>&1; then
             cargo install sccache || pip install sccache
           fi
@@ -57,7 +92,6 @@ jobs:
 
       - name: Install xwin
         run: |
-          export PATH="$HOME/.cargo/bin:$PATH"
           if ! command -v cargo-xwin >/dev/null 2>&1; then
             cargo install cargo-xwin || pip install cargo-xwin
             cargo xwin cache xwin
@@ -67,32 +101,87 @@ jobs:
       - name: Install osxcross
         run: |
           export PATH="$HOME/.osxcross/bin:$PATH"
-          if ! command -v osxcross >/dev/null 2>&1; then
+          if [ ! -d "$HOME/.osxcross/bin" ]; then
             wget ${{ secrets.OSXCROSS_LINK }} -O osxcross.tar.gz
             tar -xzf osxcross.tar.gz -C ~/
             mv ~/osxcross ~/.osxcross
           fi
+        shell: bash
+
+      - name: Shim clang for ring's aarch64-windows build
+        run: |
+          # ring 0.17.14 forces a plain `clang` for its C build on
+          # aarch64-pc-windows-msvc (its "must use Clang on Win AArch64"
+          # workaround), while the CFLAGS injected by cargo-xwin carry
+          # clang-cl-style `/imsvc <dir>` pairs that plain clang rejects.
+          # This shim rewrites them to `-I <dir>` and execs the real clang
+          # found further down PATH; all other invocations pass through
+          # untouched, and tools with other names (clang-cl, lld-link) are not
+          # affected.
+          mkdir -p /root/ci-shims
+          cat > /root/ci-shims/clang <<'EOF'
+          #!/usr/bin/env bash
+          args=()
+          prev_imsvc=0
+          for a in "$@"; do
+            if [ "$prev_imsvc" = "1" ]; then
+              args+=("-I$a")
+              prev_imsvc=0
+              continue
+            fi
+            case "$a" in
+              /imsvc)   prev_imsvc=1 ;;
+              /imsvc:*) args+=("-I${a#/imsvc:}") ;;
+              *)        args+=("$a") ;;
+            esac
+          done
+          self_dir=$(CDPATH= cd -- "$(dirname -- "$0")" && pwd)
+          real=""
+          IFS=":" read -r -a _dirs <<< "$PATH"
+          IFS=""
+          for d in "${_dirs[@]}"; do
+            [ "$d" = "$self_dir" ] && continue
+            if [ -x "$d/clang" ]; then real="$d/clang"; break; fi
+          done
+          [ -n "$real" ] || real=/usr/bin/clang
+          exec "$real" "${args[@]}"
+          EOF
+          chmod +x /root/ci-shims/clang
+        shell: bash
 
       - name: Build wheels
+        shell: bash
         run: |
-          export PATH="$HOME/.cargo/bin:$HOME/.osxcross/bin:$PATH"
-          maturin sdist --out dist
-          rustup default nightly
+          export PATH="$HOME/.osxcross/bin:$PATH"
+          # LLVM_BIN and LIBCLANG_PATH were set by the "Install build
+          # dependencies" step; put the shim (for ring's forced plain `clang`)
+          # in front of the toolchain so clang-cl / lld-link (distinct names)
+          # still resolve to the real toolchain
+          export PATH="/root/ci-shims:$LLVM_BIN:$PATH"
+          rm -rf dist
+          maturin sdist --out dist -m ndbioimage/Cargo.toml
 
-          rustup target add x86_64-unknown-linux-gnu --toolchain nightly
-          maturin build --release --out dist --target x86_64-unknown-linux-gnu
-          rustup target add aarch64-unknown-linux-gnu --toolchain nightly
-          maturin build --release --out dist --target aarch64-unknown-linux-gnu --zig
+          rustup target add x86_64-apple-darwin
+          maturin build --release --out dist --target x86_64-apple-darwin --zig --auditwheel repair -m ndbioimage/Cargo.toml
+          rustup target add aarch64-apple-darwin
+          maturin build --release --out dist --target aarch64-apple-darwin --zig --auditwheel repair -m ndbioimage/Cargo.toml
 
-          rustup target add x86_64-pc-windows-msvc --toolchain nightly
-          maturin build --release --out dist --target x86_64-pc-windows-msvc
-          rustup target add aarch64-pc-windows-msvc --toolchain nightly
-          maturin build --release --out dist --target aarch64-pc-windows-msvc
+          rustup target add x86_64-unknown-linux-gnu
+          # the x86_64-gnu build links the runner's host liblzma (lzma-sys),
+          # which needs glibc symbols up to 2.34 - higher than the
+          # manylinux_2_17 baseline
+          maturin build --release --out dist --target x86_64-unknown-linux-gnu --zig --auditwheel repair --compatibility manylinux_2_34 -m ndbioimage/Cargo.toml
+          rustup target add aarch64-unknown-linux-gnu
+          maturin build --release --out dist --target aarch64-unknown-linux-gnu --zig --auditwheel repair --compatibility manylinux_2_34 -m ndbioimage/Cargo.toml
+          rustup target add x86_64-unknown-linux-musl
+          maturin build --release --out dist --target x86_64-unknown-linux-musl --zig --auditwheel repair -m ndbioimage/Cargo.toml
+          rustup target add aarch64-unknown-linux-musl
+          maturin build --release --out dist --target aarch64-unknown-linux-musl --zig --auditwheel repair -m ndbioimage/Cargo.toml
 
-          rustup target add x86_64-apple-darwin --toolchain nightly
-          maturin build --release --out dist --target x86_64-apple-darwin --zig
-          rustup target add aarch64-apple-darwin --toolchain nightly
-          maturin build --release --out dist --target aarch64-apple-darwin --zig
+          rustup target add x86_64-pc-windows-msvc
+          maturin build --release --out dist --target x86_64-pc-windows-msvc --auditwheel repair -m ndbioimage/Cargo.toml
+          rustup target add aarch64-pc-windows-msvc
+          maturin build --release --out dist --target aarch64-pc-windows-msvc --auditwheel repair -m ndbioimage/Cargo.toml
 
       - name: Store cache
         uses: actions/cache/save@v4
@@ -102,7 +191,7 @@ jobs:
             ~/.cache/pip-wheel
             ~/.cache/sccache
             ~/.cache/cargo-xwin
-            ~/.cargo
+            ${{ env.CARGO_HOME }}
             ~/.osxcross
           key: cache-ubuntu-maturin-cross-compile
 
@@ -110,5 +199,6 @@ jobs:
         uses: pypa/gh-action-pypi-publish@release/v1
         with:
           password: ${{ secrets.PYPI_API_TOKEN }}
+          skip-existing: true
         env:
           GITHUB_WORKFLOW_REF: 1.10.1

.github/workflows/pytest.yml

@@ -0,0 +1,79 @@
+name: pytest
+
+on: [ workflow_dispatch, workflow_call, push, pull_request ]
+
+permissions:
+  contents: read
+
+jobs:
+  pypi_publish:
+    name: Publish (pypi.org)
+    runs-on: ubuntu-latest
+    env:
+      CARGO_HOME: ${{ github.workspace }}/.cargo
+      RUSTUP_HOME: ${{ github.workspace }}/.rustup
+      RUSTUP_TOOLCHAIN: stable
+    steps:
+      - uses: actions/checkout@v6
+        with:
+          submodules: recursive
+      - uses: actions/setup-python@v5
+        with:
+          python-version: 3.x
+
+      - name: Restore cache
+        uses: actions/cache/restore@v4
+        with:
+          path: |
+            ~/.cache/pip
+            ~/.cache/pip-wheel
+            ~/.cache/sccache
+            ${{ env.CARGO_HOME }}
+          key: cache-ubuntu-maturin-cross-compile
+
+      - name: Install build dependencies
+        run: |
+          sudo apt update
+          sudo apt install -y llvm cmake build-essential libclang-dev clang
+        shell: bash
+
+      - name: Install Rust
+        run: |
+          export PATH="$CARGO_HOME/bin:$PATH"
+          if ! command -v rustup >/dev/null 2>&1; then
+            curl --proto '=https' --tlsv1.2 -sSf https://sh.rustup.rs | sh -s -- -y
+          fi
+          rustup toolchain install stable --profile minimal
+          rustup default stable
+          echo "$CARGO_HOME/bin" >> "$GITHUB_PATH"
+          cargo --version
+        shell: bash
+
+      - name: Install sccache and maturin
+        run: |
+          python -m pip install --upgrade pip maturin ziglang cargo-zigbuild
+          # in case an older copy was cargo-installed into the CARGO_HOME cache,
+          # prevent it from shadowing the pip-installed one
+          rm -f "${CARGO_HOME:?}/bin/cargo-zigbuild" || true
+          if ! command -v sccache >/dev/null 2>&1; then
+            cargo install sccache || pip install sccache
+          fi
+        shell: bash
+
+      - name: PyTest
+        run: |
+          cd ndbioimage
+          python -m pip install .[test]
+          python -m pytest 
+
+      - name: Store cache
+        uses: actions/cache/save@v4
+        with:
+          path: |
+            ~/.cache/pip
+            ~/.cache/pip-wheel
+            ~/.cache/sccache
+            ~/.cache/cargo-xwin
+            ${{ env.CARGO_HOME }}
+            ~/.osxcross
+          key: cache-ubuntu-maturin-cross-compile

.gitignore

@@ -29,6 +29,7 @@ venv/
 *.egg-info/
 .installed.cfg
 *.egg
+.cargo
 
 # Installer logs
 pip-log.txt
@@ -72,11 +73,11 @@ docs/_build/
 # Pyenv
 .python-version
 
-/tests/files/*
 AGENTS.md
 .agentbridge
 .agent-work
 
-py/ndbioimage/jassets
-py/ndbioimage/deps
+/ndbioimage/tests/files/*
+/ndbioimage/py/ndbioimage/jassets
+/ndbioimage/py/ndbioimage/deps
 /build_wheels.sh

.gitmodules

@@ -0,0 +1,6 @@
+[submodule "image-tiff"]
+	path = image-tiff
+	url = https://git.pomppervova.nl/wim/image-tiff.git
+[submodule "image-registration"]
+	path = image-registration
+	url = https://git.pomppervova.nl/wim/image-registration.git

Cargo.toml

@@ -1,90 +1,26 @@
-[package]
-name = "ndbioimage"
-version = "0.2.0"
+[workspace]
+resolver = "3"
+members = ["ndbioimage", "image-tiff", "image-registration"]
+
+[workspace.package]
+version = "0.3.2"
 edition = "2024"
-rust-version = "1.94.0"
+rust-version = "1.94.1"
 authors = ["Wim Pomp <w.pomp@nki.nl>"]
 license = "MIT OR Apache-2.0"
 description = "Read bio image formats using the bio-formats java package."
-homepage = "https://git.pomppervova.nl/wim/ndbioimage/src/branch/rs"
-repository = "https://git.pomppervova.nl/wim/ndbioimage/src/branch/rs"
+repository = "https://git.pomppervova.nl/wim/ndbioimage"
 documentation = "https://docs.rs/ndbioimage"
 readme = "README.md"
 keywords = ["bioformats", "imread", "ndarray", "metadata"]
 categories = ["multimedia::images", "science"]
-exclude = ["/tests"]
-
-[lib]
-name = "ndbioimage"
-crate-type = ["cdylib", "rlib"]
-
-[dependencies]
-bioformats = { version = "0.1", optional = true, default-features = false, features = ["zarr", "tissuefaxs"] }
-clap = { version = "4", features = ["derive"] }
-console = { version = "0.16", optional = true }
-downloader = { version = "0.2", optional = true, default-features = false, features = ["rustls-tls"] }
-ffmpeg-sidecar = { version = "2", optional = true }
-itertools = "0.15"
-image-registration = { path = "../image-registration", optional = true }
-indexmap = { version = "2", features = ["serde"] }
-indicatif = { version = "0.18", features = ["rayon"], optional = true }
-j4rs = { version = "0.25", optional = true }
-libczirw-sys = { version = "0.5", optional = true }
-ndarray = { version = "0.17", features = ["serde"] }
-num = "0.4"
-numpy = { version = "0.29", optional = true }
-ome-metadata = "0.5"
-ordered-float = { version = "5", optional = true }
-phf = { version = "0.14", features = ["macros"] }
-postcard = { version = "1", features = ["use-std"], optional = true }
-pyo3 = { version = "0.29", features = ["abi3-py310", "anyhow", "generate-import-lib"], optional = true }
-pyo3-stub-gen = { version = "0.23", optional = true }
-rayon = { version = "1", optional = true }
-regex = "1"
-serde = { version = "1", features = ["rc", "derive"] }
-serde_with = "3"
-strum = { version = "0.28", features = ["derive"] }
-thiserror = "2"
-tiff = { version = "0.11", features = ["zstd"], optional = true }
-tiffwrite = { version = "2026.6.0", optional = true }
-tokio = { version = "1", features = ["rt", "rt-multi-thread"], optional = true }
-thread_local = { version = "1", optional = true }
-xmltree = { version = "0.12", optional = true }
-yaml_serde = { version = "0.10", optional = true }
-
-[dev-dependencies]
-rayon = "1"
-regex = "1"
-
-[build-dependencies]
-j4rs = "0.25"
-ffmpeg-sidecar = "2"
-retry = "2"
-toml = "1"
-
-[features]
-default = ["bioformats_java", "gpl-formats", "czi", "tiff", "tiffseq", "movie", "tiffwrite"]
-all = ["bioformats_java", "bioformats_rust", "czi", "gpl-formats", "movie", "tiffseq", "tiffwrite", "tiff", "transforms"]
-gpl-formats = []
-python = ["dep:pyo3", "dep:numpy", "dep:pyo3-stub-gen", "dep:postcard", "ome-metadata/python"]
-czi = ["dep:libczirw-sys", "dep:xmltree", "dep:thread_local"]
-bioformats_rust = ["dep:bioformats", "dep:thread_local"]
-bioformats_java = ["dep:j4rs", "dep:thread_local", "dep:downloader"]
-tiffwrite = ["dep:tiffwrite", "dep:indicatif", "dep:console", "dep:rayon"]
-tiffseq = ["dep:tiff", "dep:yaml_serde"]
-tiff = ["dep:tiff", "dep:thread_local"]
-movie = ["dep:ffmpeg-sidecar", "dep:tokio", "dep:ordered-float", "dep:indicatif", "dep:console"]
-transforms = ["dep:image-registration", "dep:yaml_serde"]
-
-[package.metadata.docs.rs]
-no-default-features = true
-features = ["bioformats_java", "bioformats_rust", "czi", "movie", "tiffseq", "tiffwrite", "tiff", "transforms"]
 
 [profile.test]
 inherits = "release"
+debug = true
 
 [profile.release]
-debug = true
+#debug = true
 
 [profile.dev.package.backtrace]
 opt-level = 3

README.md

@@ -1,11 +1,10 @@
-# ndbioimage
-
-[![Pytest](https://github.com/pomppervova/ndbioimage/actions/workflows/pytest.yml/badge.svg)](https://github.com/pomppervova/ndbioimage/actions/workflows/pytest.yml)
+[![pytest](https://git.pomppervova.nl/wim/ndbioimage/actions/workflows/pytest.yml/badge.svg)](https://git.pomppervova.nl/wim/ndbioimage/actions?workflow=pytest.yml)
 
+# NDBioImage
 ## Work in progress
 
-Rust rewrite of python version. Read bio image formats using the bio-formats java package.
-[https://www.openmicroscopy.org/bio-formats/](https://www.openmicroscopy.org/bio-formats/)
+Rust rewrite of python version. Read bio image formats using some Rust code, the libczi api, and if needed the
+bio-formats java package.
 
 Exposes (bio) images as a numpy ndarray-like object (Python) or a struct that can be sliced like an ndarray Array
 (Rust), but without loading the whole image into memory, reading from the file only when needed. Some metadata is read
@@ -18,20 +17,6 @@ Currently, it supports imagej tif files, czi files, micromanager tif sequences a
 
 To transition to semver, versions before 0.1.0 were yanked from crates.io.
 
-## Installation
-
-```sh
-pip install ndbioimage
-```
-
-### Installation with option to write mp4 or mkv:
-
-Work in progress! Make sure ffmpeg is installed.
-
-```sh
-pip install ndbioimage[write]
-```
-
 ## Usage
 
 ### Python

ndbioimage/Cargo.toml

@@ -0,0 +1,84 @@
+[package]
+name = "ndbioimage"
+version.workspace = true
+edition.workspace = true
+rust-version.workspace = true
+authors.workspace = true
+license.workspace = true
+description.workspace = true
+repository.workspace = true
+documentation.workspace = true
+keywords.workspace = true
+categories.workspace = true
+readme = "../README.md"
+exclude = ["/tests"]
+
+[lib]
+name = "ndbioimage"
+crate-type = ["cdylib", "rlib"]
+
+[dependencies]
+bioformats = { version = "0.1", optional = true, default-features = false, features = ["zarr", "tissuefaxs"] }
+clap = { version = "4", features = ["derive"] }
+console = { version = "0.16", optional = true }
+downloader = { version = "0.2", optional = true, default-features = false, features = ["rustls-tls"] }
+ffmpeg-sidecar = { version = "2", optional = true }
+itertools = "0.15"
+image-registration = { path = "../image-registration", version = "0.1", optional = true }
+indexmap = { version = "2", features = ["serde"] }
+indicatif = { version = "0.18", features = ["rayon"], optional = true }
+j4rs = { version = "0.25", optional = true }
+libczirw-sys = { version = "0.5.2", optional = true }
+ndarray = { version = "0.17", features = ["serde"] }
+num = "0.4"
+numpy = { version = "0.29", optional = true }
+ome-metadata = "0.6"
+ordered-float = { version = "5", optional = true }
+phf = { version = "0.14", features = ["macros"] }
+postcard = { version = "1", features = ["use-std"], optional = true }
+pyo3 = { version = "0.29", features = ["abi3-py310", "anyhow", "generate-import-lib", "multiple-pymethods"], optional = true }
+pyo3-stub-gen = { version = "0.23", optional = true }
+rayon = { version = "1", optional = true }
+regex = "1"
+serde = { version = "1", features = ["rc", "derive"] }
+serde_with = "3"
+strum = { version = "0.28", features = ["derive"] }
+thiserror = "2"
+tiff = { package = "ndbioimage-tiff", path = "../image-tiff", version = "0.11", features = ["zstd"], optional = true }
+tiffwrite = { version = "2026.6.0", optional = true }
+tokio = { version = "1", features = ["rt", "rt-multi-thread"], optional = true }
+thread_local = { version = "1", optional = true }
+xmltree = { version = "0.12", optional = true }
+yaml_serde = { version = "0.10", optional = true }
+# "bioformats 0.1.8" is written for "zip 2.4.x", where zip::read::ZipFile has a single lifetime parameter.
+# "zip >= 2.5" added a second "R: Read" generic (2.6.1 breaks bioformats), so the build fails.
+zip = { version = "=2.4.2", optional = true, default-features = false }
+
+[dev-dependencies]
+rayon = "1"
+regex = "1"
+postcard = { version = "1", features = ["use-std"] }
+
+[build-dependencies]
+j4rs = "0.25"
+ffmpeg-sidecar = "2"
+retry = "2"
+toml = "1"
+
+[features]
+default = ["bioformats_java", "gpl-formats", "czi", "tiff", "tiffseq", "movie", "tiffwrite"]
+all = ["bioformats_java", "bioformats_rust", "czi", "gpl-formats", "movie", "tiffseq", "tiffwrite", "tiff", "transforms"]
+gpl-formats = []
+python = ["dep:pyo3", "dep:numpy", "dep:pyo3-stub-gen", "dep:postcard", "ome-metadata/python"]
+czi = ["dep:libczirw-sys", "dep:xmltree", "dep:thread_local"]
+bioformats_rust = ["dep:bioformats", "dep:thread_local", "dep:zip"]
+bioformats_java = ["dep:j4rs", "dep:thread_local", "dep:downloader"]
+tiffwrite = ["dep:tiffwrite", "dep:indicatif", "dep:console", "dep:rayon"]
+tiffseq = ["dep:tiff", "dep:yaml_serde"]
+tiff = ["dep:tiff", "dep:thread_local"]
+movie = ["dep:ffmpeg-sidecar", "dep:tokio", "dep:ordered-float", "dep:indicatif", "dep:console"]
+transforms = ["dep:image-registration", "dep:yaml_serde"]
+
+[package.metadata.docs.rs]
+no-default-features = true
+features = ["bioformats_java", "bioformats_rust", "czi", "movie", "tiffseq", "tiffwrite", "tiff", "transforms"]

ndbioimage/build.rs

@@ -9,21 +9,41 @@ mod bioformats {
     #[derive(Clone, Debug)]
     pub(crate) enum BuildError {
         J4rsVersionNotFound,
+        LockFileNotFound,
     }
 
     impl Display for BuildError {
         fn fmt(&self, fmt: &mut Formatter) -> Result<(), std::fmt::Error> {
             match self {
                 Self::J4rsVersionNotFound => write!(fmt, "J4rsVersion not found in Cargo.lock"),
+                Self::LockFileNotFound => write!(
+                    fmt,
+                    "Cargo.lock not found in any parent of CARGO_MANIFEST_DIR ({})",
+                    std::env::var("CARGO_MANIFEST_DIR").unwrap_or_default()
+                ),
             }
         }
     }
 
     impl Error for BuildError {}
 
+    /// the `Cargo.lock` that cargo resolved for this build: next to the manifest for a standalone
+    /// package, or at the workspace root when this crate is a workspace member
+    fn lock_file() -> Result<std::path::PathBuf, BuildError> {
+        let manifest_dir = std::path::PathBuf::from(
+            std::env::var("CARGO_MANIFEST_DIR").expect("CARGO_MANIFEST_DIR not set"),
+        );
+        manifest_dir
+            .ancestors()
+            .map(|dir| dir.join("Cargo.lock"))
+            .find(|lock| lock.is_file())
+            .ok_or(BuildError::LockFileNotFound)
+    }
+
     fn get_j4rs_version() -> Result<String, Box<dyn Error>> {
-        let manifest_dir = std::env::var("CARGO_MANIFEST_DIR").expect("CARGO_MANIFEST_DIR not set");
-        let lock_path = std::path::Path::new(&manifest_dir).join("Cargo.lock");
+        let lock_path = lock_file()?;
+        // the resolved j4rs version is baked into the generated constants, so rerun when it changes
+        println!("cargo::rerun-if-changed={}", lock_path.display());
         let lock_toml = std::fs::read_to_string(&lock_path)?;
         let value: toml::Value = toml::from_str(&lock_toml)?;
         if let Some(packages) = value.get("package").and_then(|v| v.as_array()) {

ndbioimage/py/ndbioimage/__init__.py

@@ -5,7 +5,7 @@ import sys
 import warnings
 from importlib.metadata import version
 from pathlib import Path
-from typing import TypeVar
+from typing import TYPE_CHECKING
 
 import numpy as np
 from numpy.typing import ArrayLike
@@ -13,14 +13,18 @@ from numpy.typing import ArrayLike
 os.environ["RUST_BACKTRACE"] = "full"
 os.environ["COLORBT_SHOW_HIDDEN"] = "1"
 
-from .ndbioimage_rs import Imread, generate_stub
+from .ndbioimage_rs import Imread, Shape, generate_stub
 
-try:
-    from .ndbioimage_rs import transforms
+if TYPE_CHECKING:
+    from .transforms import Transform as Transform
+else:
+    try:
+        from .ndbioimage_rs import transforms
 
-    sys.modules["ndbioimage.transforms"] = transforms
-except (ImportError, ModuleNotFoundError):
-    pass
+        sys.modules["ndbioimage.transforms"] = transforms
+        Transform = transforms.Transform
+    except (ImportError, ModuleNotFoundError):
+        pass
 
 try:
     from .ndbioimage_rs import batch_to_tiff
@@ -50,7 +54,16 @@ except Exception:  # noqa
     __git_commit_hash__ = "unknown"
 
 warnings.filterwarnings("ignore", "Reference to unknown ID")
-Number = int | float | np.integer | np.floating
+
+__all__ = [
+    "Imread",
+    "JVMException",
+    "ReaderNotFoundError",
+    "Shape",
+    "Transform",
+    "get_positions",
+    "is_noise",
+]
 
 
 # backwards compatibility
@@ -90,9 +103,6 @@ def ndbioimage_generate_stub():
             (nd / "ndbioimage.pyi").rename(nd / "__init__.pyi")
 
 
-R = TypeVar("R")
-
-
 def get_positions(path: str | Path) -> set[int]:
     return Imread.get_available_series(path)
 

ndbioimage/py/ndbioimage/ndbioimage_rs/__init__.pyi

@@ -2,6 +2,8 @@
 # ruff: noqa: E501, F401, F403, F405
 
 import builtins
+import collections
+import collections.abc
 import os
 import pathlib
 import typing
@@ -52,6 +54,8 @@ class Flags:
     def __new__(cls) -> Flags: ...
     def __getnewargs__(self) -> None: ...
     def __getitem__(self, item: builtins.str) -> builtins.bool: ...
+    def __str__(self) -> builtins.str: ...
+    def __repr__(self) -> builtins.str: ...
 
 class Imread:
     r"""
@@ -96,16 +100,6 @@ class Imread:
         the name of the reader used to open the file
         """
     @property
-    def transforms(self) -> transforms.Transforms:
-        r"""
-        get the transformation
-        """
-    @transforms.setter
-    def transforms(self, value: transforms.Transforms) -> None:
-        r"""
-        set the transformation
-        """
-    @property
     def path(self) -> pathlib.Path:
         r"""
         the file path
@@ -220,6 +214,38 @@ class Imread:
         r"""
         the name of the tube lens
         """
+    @property
+    def channel_names(self) -> builtins.dict[builtins.int, builtins.str]:
+        r"""
+        names of the channels in order
+        """
+    @property
+    def summary(self) -> builtins.str:
+        r"""
+        gives a helpful summary of the recorded experiment
+        """
+    @property
+    def interpolation_order(self) -> typing.Optional[builtins.int]:
+        r"""
+        the bspline order used for interpolation when using transforms
+        None -> 0 (nearest neighbor) for int dtypes, 3 (cubic) for float dtypes
+        """
+    @interpolation_order.setter
+    def interpolation_order(self, value: typing.Optional[builtins.int]) -> None:
+        r"""
+        the bspline order used for interpolation when using transforms
+        None -> 0 (nearest neighbor) for int dtypes, 3 (cubic) for float dtypes
+        """
+    @property
+    def transforms(self) -> transforms.Transforms:
+        r"""
+        get the transformation
+        """
+    @transforms.setter
+    def transforms(self, value: transforms.Transforms) -> None:
+        r"""
+        set the transformation
+        """
     def __new__(
         cls,
         path: str | pathlib.Path | Imread | bytes,
@@ -228,7 +254,33 @@ class Imread:
         reader: typing.Optional[builtins.str] = None,
     ) -> Imread:
         r"""
-        new view on a file at path, open series #, open as dtype: (u)int(8/16/32) or float(32/64)
+        new view on an image file
+        path: path to image file to open, can include /Pos\d
+        dtype: when not None; will cast image to this dtype,
+            choices: (u)int(8/16/32) or float(32/64)
+        axes: the resulting object will represent the image with the axes in this order
+        reader: use a specific reader to open the file,
+            view available choices with Imread.available_readers
+        """
+    def debug(self) -> None:
+        r"""
+        some debug info
+        """
+    @staticmethod
+    def warm_jvm() -> None:
+        r"""
+        creates the jvm (if it does not exist yet) so that hotspot installs its signal handlers before
+        anything else does: see `crate::readers::bioformats_java::warm_jvm`
+
+        in a python process, call this before `ray.init()`, which installs CPython's faulthandler
+        """
+    @staticmethod
+    def jvm_is_running() -> builtins.bool:
+        r"""
+        returns whether the jvm was already created in this process
+
+        this is false in a process that inherited a jvm across a fork, since such a jvm cannot be
+        used; reading a bio-formats file there raises instead
         """
     def set_flags(
         self,
@@ -251,47 +303,19 @@ class Imread:
         r"""
         only remains for backwards compatibility
         """
-    def reshape(self, order: builtins.str, copy: builtins.bool) -> typing.Any:
-        r"""
-        reshape the view with a new axis order
-        """
-    def with_transform(
+    def reshape(
         self,
-        channels: builtins.bool = True,
-        drift: builtins.bool = False,
-        file: typing.Optional[typing.Any] = None,
-        bead_files: typing.Optional[typing.Any] = None,
-        main_channel: typing.Optional[builtins.int] = None,
-        default_transform: typing.Optional[typing.Sequence[builtins.float]] = None,
-    ) -> Imread:
+        *shape: typing.Any,
+        order: builtins.str = "C",
+        copy: typing.Optional[builtins.bool] = None,
+    ) -> Imread | numpy.NDArray:
         r"""
-        return a new view with transformations applied (channel alignment, drift correction)
-        backward compatibility
+        reshape the view with a new axis order, arguments order and copy cannot be used and will
+        raise an exception if order != "C" or copy is False
         """
-    def with_transforms(
-        self,
-        channels: builtins.bool = True,
-        drift: builtins.bool = False,
-        file: typing.Optional[typing.Any] = None,
-        bead_files: typing.Optional[typing.Any] = None,
-        main_channel: typing.Optional[builtins.int] = None,
-        default_transform: typing.Optional[typing.Sequence[builtins.float]] = None,
-    ) -> Imread:
-        r"""
-        return a new view with transformations applied (channel alignment, drift correction)
-        """
-    def load_transforms_from_yaml(
-        self, path: builtins.str | os.PathLike | pathlib.Path
-    ) -> None: ...
-    def calculate_channel_transforms_2d(
-        self, main_channel: builtins.int
-    ) -> builtins.dict[builtins.str, transforms.Transform]: ...
-    def calculate_channel_transforms_3d(
-        self, main_channel: builtins.int
-    ) -> builtins.dict[builtins.str, transforms.Transform]: ...
-    def calculate_drift_transforms_2d(self) -> builtins.list[transforms.Transform]: ...
-    def calculate_drift_transforms_3d(self) -> builtins.list[transforms.Transform]: ...
-    def squeeze(self) -> numpy.ndarray | int | float: ...
+    def squeeze(
+        self, axes: str | collections.abc.Sequence[int | str] = None
+    ) -> Imread | float | int: ...
     def close(self) -> None:
         r"""
         close the file: does nothing as this is handled automatically
@@ -304,7 +328,19 @@ class Imread:
         r"""
         change the data type of the view: (u)int(8/16/32) or float(32/64)
         """
-    def __getitem__(self, n: typing.Any) -> typing.Any:
+    def __call__(
+        self,
+        c: typing.Optional[int | slice | numpy.ArrayLike] = None,
+        z: typing.Optional[int | slice | numpy.ArrayLike] = None,
+        t: typing.Optional[int | slice | numpy.ArrayLike] = None,
+        y: typing.Optional[int | slice | numpy.ArrayLike] = None,
+        x: typing.Optional[int | slice | numpy.ArrayLike] = None,
+    ) -> typing.Optional[int | Imread]:
+        r"""
+        same as im[] but allowing keyword axes, but slices need to made with slice()
+        also, no ellipsis or new axes
+        """
+    def __getitem__(self, n: typing.Any) -> typing.Optional[int | Imread]:
         r"""
         slice the view and return a new view or a single number
         """
@@ -312,7 +348,7 @@ class Imread:
         self,
         dtype: typing.Optional[numpy.typing.DTypeLike] = None,
         copy: typing.Optional[builtins.bool] = None,
-    ) -> typing.Any:
+    ) -> numpy.typing.NDArray:
         r"""
         convert to a numpy array, optionally with a different dtype
         """
@@ -320,139 +356,199 @@ class Imread:
         r"""
         check if an item is contained in the view (not implemented)
         """
-    def __lt__(self, other: typing.Any) -> typing.Any:
+    def __lt__(
+        self, other: float | numpy.ArrayLike | Imread
+    ) -> numpy.typing.NDArray | bool:
         r"""
         element-wise less than comparison
         """
-    def __le__(self, other: typing.Any) -> typing.Any:
+    def __le__(
+        self, other: float | numpy.ArrayLike | Imread
+    ) -> numpy.typing.NDArray | bool:
         r"""
         element-wise less than or equal comparison
         """
-    def __eq__(self, other: typing.Any) -> typing.Any:
+    def __eq__(
+        self, other: float | numpy.ArrayLike | Imread
+    ) -> numpy.typing.NDArray | bool:
         r"""
         element-wise equality comparison
         """
-    def __ne__(self, other: typing.Any) -> typing.Any:
+    def __ne__(
+        self, other: float | numpy.ArrayLike | Imread
+    ) -> numpy.typing.NDArray | bool:
         r"""
         element-wise not equal comparison
         """
-    def __gt__(self, other: typing.Any) -> typing.Any:
+    def __gt__(
+        self, other: float | numpy.ArrayLike | Imread
+    ) -> numpy.typing.NDArray | bool:
         r"""
         element-wise greater than comparison
         """
-    def __ge__(self, other: typing.Any) -> typing.Any:
+    def __ge__(
+        self, other: float | numpy.ArrayLike | Imread
+    ) -> numpy.typing.NDArray | bool:
         r"""
         element-wise greater than or equal comparison
         """
-    def __add__(self, other: typing.Any) -> typing.Any:
+    def __add__(
+        self, other: float | numpy.ArrayLike | Imread
+    ) -> numpy.typing.NDArray | float:
         r"""
         element-wise addition
         """
-    def __radd__(self, other: typing.Any) -> typing.Any:
+    def __radd__(
+        self, other: float | numpy.ArrayLike | Imread
+    ) -> numpy.typing.NDArray | float:
         r"""
         element-wise addition (reflected)
         """
-    def __sub__(self, other: typing.Any) -> typing.Any:
+    def __sub__(
+        self, other: float | numpy.ArrayLike | Imread
+    ) -> numpy.typing.NDArray | float:
         r"""
         element-wise subtraction
         """
-    def __rsub__(self, other: typing.Any) -> typing.Any:
+    def __rsub__(
+        self, other: float | numpy.ArrayLike | Imread
+    ) -> numpy.typing.NDArray | float:
         r"""
         element-wise subtraction (reflected)
         """
-    def __mul__(self, other: typing.Any) -> typing.Any:
+    def __mul__(
+        self, other: float | numpy.ArrayLike | Imread
+    ) -> numpy.typing.NDArray | float:
         r"""
         element-wise multiplication
         """
-    def __rmul__(self, other: typing.Any) -> typing.Any:
+    def __rmul__(
+        self, other: float | numpy.ArrayLike | Imread
+    ) -> numpy.typing.NDArray | float:
         r"""
         element-wise multiplication (reflected)
         """
-    def __truediv__(self, other: typing.Any) -> typing.Any:
+    def __truediv__(
+        self, other: float | numpy.ArrayLike | Imread
+    ) -> numpy.typing.NDArray | float:
         r"""
         element-wise true division
         """
-    def __rtruediv__(self, other: typing.Any) -> typing.Any:
+    def __rtruediv__(
+        self, other: float | numpy.ArrayLike | Imread
+    ) -> numpy.typing.NDArray | float:
         r"""
         element-wise true division (reflected)
         """
-    def __floordiv__(self, other: typing.Any) -> typing.Any:
+    def __floordiv__(
+        self, other: float | numpy.ArrayLike | Imread
+    ) -> numpy.typing.NDArray | float:
         r"""
         element-wise floor division
         """
-    def __rfloordiv__(self, other: typing.Any) -> typing.Any:
+    def __rfloordiv__(
+        self, other: float | numpy.ArrayLike | Imread
+    ) -> numpy.typing.NDArray | float:
         r"""
         element-wise floor division (reflected)
         """
-    def __mod__(self, other: typing.Any) -> typing.Any:
+    def __mod__(
+        self, other: float | numpy.ArrayLike | Imread
+    ) -> numpy.typing.NDArray | float:
         r"""
         element-wise modulo
         """
-    def __rmod__(self, other: typing.Any) -> typing.Any:
+    def __rmod__(
+        self, other: float | numpy.ArrayLike | Imread
+    ) -> numpy.typing.NDArray | float:
         r"""
         element-wise modulo (reflected)
         """
-    def __matmul__(self, other: typing.Any) -> typing.Any:
+    def __matmul__(
+        self, other: float | numpy.ArrayLike | Imread
+    ) -> numpy.typing.NDArray | float:
         r"""
         element-wise matrix multiplication
         """
-    def __rmatmul__(self, other: typing.Any) -> typing.Any:
+    def __rmatmul__(
+        self, other: float | numpy.ArrayLike | Imread
+    ) -> numpy.typing.NDArray | float:
         r"""
         element-wise matrix multiplication (reflected)
         """
-    def __and__(self, other: typing.Any) -> typing.Any:
+    def __and__(
+        self, other: float | numpy.ArrayLike | Imread
+    ) -> numpy.typing.NDArray | bool:
         r"""
         element-wise bitwise AND
         """
-    def __rand__(self, other: typing.Any) -> typing.Any:
+    def __rand__(
+        self, other: float | numpy.ArrayLike | Imread
+    ) -> numpy.typing.NDArray | bool:
         r"""
         element-wise bitwise AND (reflected)
         """
-    def __or__(self, other: typing.Any) -> typing.Any:
+    def __or__(
+        self, other: float | numpy.ArrayLike | Imread
+    ) -> numpy.typing.NDArray | bool:
         r"""
         element-wise bitwise OR
         """
-    def __ror__(self, other: typing.Any) -> typing.Any:
+    def __ror__(
+        self, other: float | numpy.ArrayLike | Imread
+    ) -> numpy.typing.NDArray | bool:
         r"""
         element-wise bitwise OR (reflected)
         """
-    def __xor__(self, other: typing.Any) -> typing.Any:
+    def __xor__(
+        self, other: float | numpy.ArrayLike | Imread
+    ) -> numpy.typing.NDArray | bool:
         r"""
         element-wise bitwise XOR
         """
-    def __rxor__(self, other: typing.Any) -> typing.Any:
+    def __rxor__(
+        self, other: float | numpy.ArrayLike | Imread
+    ) -> numpy.typing.NDArray | bool:
         r"""
         element-wise bitwise XOR (reflected)
         """
-    def __lshift__(self, other: typing.Any) -> typing.Any:
+    def __lshift__(
+        self, other: float | numpy.ArrayLike | Imread
+    ) -> numpy.typing.NDArray | float:
         r"""
         element-wise left shift
         """
-    def __rlshift__(self, other: typing.Any) -> typing.Any:
+    def __rlshift__(
+        self, other: float | numpy.ArrayLike | Imread
+    ) -> numpy.typing.NDArray | float:
         r"""
         element-wise left shift (reflected)
         """
-    def __rshift__(self, other: typing.Any) -> typing.Any:
+    def __rshift__(
+        self, other: float | numpy.ArrayLike | Imread
+    ) -> numpy.typing.NDArray | float:
         r"""
         element-wise right shift
         """
-    def __rrshift__(self, other: typing.Any) -> typing.Any:
+    def __rrshift__(
+        self, other: float | numpy.ArrayLike | Imread
+    ) -> numpy.typing.NDArray | float:
         r"""
         element-wise right shift (reflected)
         """
-    def __neg__(self) -> typing.Any:
+    def __neg__(self) -> numpy.typing.NDArray | float:
         r"""
         element-wise negation
         """
-    def __pos__(self) -> typing.Any:
+    def __pos__(self) -> numpy.typing.NDArray | float:
         r"""
         element-wise positive
         """
-    def __abs__(self) -> typing.Any:
+    def __abs__(self) -> numpy.typing.NDArray | float:
         r"""
         element-wise absolute value
         """
-    def __invert__(self) -> typing.Any:
+    def __invert__(self) -> numpy.typing.NDArray | float:
         r"""
         element-wise bitwise inversion
         """
@@ -505,13 +601,19 @@ class Imread:
         r"""
         the file path as a string
         """
+    def frame(
+        self, c: builtins.int, z: builtins.int, t: builtins.int
+    ) -> numpy.typing.NDArray | float:
+        r"""
+        retrieve a single frame at czt, sliced accordingly, alias of get_frame
+        """
     def get_frame(
         self, c: builtins.int, z: builtins.int, t: builtins.int
-    ) -> typing.Any:
+    ) -> numpy.typing.NDArray | float:
         r"""
         retrieve a single frame at czt, sliced accordingly
         """
-    def flatten(self) -> typing.Any:
+    def flatten(self) -> numpy.typing.NDArray | float:
         r"""
         flatten the view into a 1D numpy array
         """
@@ -527,17 +629,21 @@ class Imread:
         r"""
         find the position of an axis
         """
+    def swapaxes(self, ax0: int | str, ax1: int | str) -> Imread:
+        r"""
+        swap two axes
+        """
     def swap_axes(self, ax0: int | str, ax1: int | str) -> Imread:
         r"""
         swap two axes
         """
     def transpose(
-        self, axes: typing.Optional[typing.Sequence[int | str]] = None
+        self, axes: typing.Optional[collections.abc.Sequence[int | str] | str] = None
     ) -> Imread:
         r"""
         permute the order of the axes
         """
-    def as_array(self) -> typing.Any:
+    def as_array(self) -> numpy.typing.NDArray | float:
         r"""
         collect data into a numpy array
         """
@@ -567,9 +673,9 @@ class Imread:
         r"""
         the detector gain for a given channel
         """
-    def summary(self) -> builtins.str:
+    def channel_name(self, channel: builtins.int) -> builtins.str:
         r"""
-        gives a helpful summary of the recorded experiment
+        the name of a particular channel
         """
     @staticmethod
     def get_available_series(
@@ -673,6 +779,73 @@ class Imread:
         r"""
         Return the sum along a given axis. Arguments beyond axis are not implemented
         """
+    def with_interpolation_order(self, order: typing.Optional[builtins.int]) -> Imread:
+        r"""
+        the bspline order used for interpolation when using transforms
+        None -> 0 (nearest neighbor) for int dtypes, 3 (cubic) for float dtypes
+        """
+    def get_default_interpolation_value(
+        self, channel: builtins.int
+    ) -> transforms.DefaultValue:
+        r"""
+        the value of the interpolated image outside the data
+        """
+    def set_default_interpolation_value(
+        self,
+        channel: builtins.int,
+        value: float | transforms.DefaultValue | collections.abc.Sequence[float],
+    ) -> None:
+        r"""
+        the value of the interpolated image outside the data
+        """
+    def with_default_interpolation_value(
+        self,
+        channel: builtins.int,
+        value: float | transforms.DefaultValue | collections.abc.Sequence[float],
+    ) -> Imread:
+        r"""
+        the value of the interpolated image outside the data
+        """
+    def with_transform(
+        self,
+        channels: builtins.bool = True,
+        drift: builtins.bool = False,
+        file: typing.Optional[builtins.str | os.PathLike | pathlib.Path] = None,
+        bead_files: typing.Optional[
+            typing.Sequence[builtins.str | os.PathLike | pathlib.Path]
+        ] = None,
+        main_channel: typing.Optional[builtins.int] = None,
+        default_transform: typing.Optional[typing.Sequence[builtins.float]] = None,
+    ) -> Imread:
+        r"""
+        return a new view with transformations applied (channel alignment, drift correction)
+        backward compatibility
+        """
+    def with_transforms(
+        self,
+        channels: builtins.bool = True,
+        drift: builtins.bool = False,
+        file: typing.Optional[builtins.str | os.PathLike | pathlib.Path] = None,
+        bead_files: typing.Optional[
+            typing.Sequence[builtins.str | os.PathLike | pathlib.Path]
+        ] = None,
+        main_channel: typing.Optional[builtins.int] = None,
+        default_transform: typing.Optional[typing.Sequence[builtins.float]] = None,
+    ) -> Imread:
+        r"""
+        returns a view where channels and/or frames are registered with an affine transformation
+        """
+    def load_transforms_from_yaml(
+        self, path: builtins.str | os.PathLike | pathlib.Path
+    ) -> None: ...
+    def calculate_channel_transforms_2d(
+        self, main_channel: builtins.int
+    ) -> builtins.dict[builtins.str, transforms.Transform]: ...
+    def calculate_channel_transforms_3d(
+        self, main_channel: builtins.int
+    ) -> builtins.dict[builtins.str, transforms.Transform]: ...
+    def calculate_drift_transforms_2d(self) -> builtins.list[transforms.Transform]: ...
+    def calculate_drift_transforms_3d(self) -> builtins.list[transforms.Transform]: ...
 
 class Shape:
     r"""
@@ -742,7 +915,7 @@ class Shape:
         arguments for pickling
         """
     def __getitem__(
-        self, idx: str | int | None | Ellipsis | slice | list[int] | tuple[int]
+        self, idx: str | int | None | Ellipsis | slice | typing.Sequence[int]
     ) -> typing.Optional[int | tuple[int]]:
         r"""
         get dimension size by index or axis name
@@ -773,7 +946,7 @@ def batch_to_tiff(
 
 def generate_stub(dest_path: builtins.str) -> None:
     r"""
-    generates ndbioimage/__init__.pyi
+    generates .pyi files
     """
 
 def main() -> None:

ndbioimage/py/ndbioimage/ndbioimage_rs/transforms.pyi

@@ -2,6 +2,7 @@
 # ruff: noqa: E501, F401, F403, F405
 
 import builtins
+import collections.abc
 import os
 import pathlib
 import typing
@@ -10,6 +11,7 @@ import numpy
 import numpy.typing
 
 __all__ = [
+    "DefaultValue",
     "Optimizer",
     "RegistrationStep",
     "SamplingArg",
@@ -18,6 +20,21 @@ __all__ = [
     "Transforms",
 ]
 
+class DefaultValue:
+    @property
+    def values(self) -> builtins.list[builtins.float]: ...
+    @property
+    def weights(self) -> builtins.list[builtins.float]: ...
+    def __eq__(self, other: builtins.object, /) -> builtins.bool: ...
+    def __new__(
+        cls,
+        values: typing.Sequence[builtins.float],
+        weights: typing.Optional[typing.Sequence[builtins.float]] = None,
+    ) -> DefaultValue: ...
+    def __getnewargs__(
+        self,
+    ) -> tuple[builtins.list[builtins.float], builtins.list[builtins.float]]: ...
+
 class Optimizer:
     def __eq__(self, other: builtins.object, /) -> builtins.bool: ...
     def __new__(cls, optimizer: builtins.str) -> Optimizer: ...
@@ -116,14 +133,6 @@ class Transform:
     @dparameters.setter
     def dparameters(self, value: typing.Sequence[builtins.float]) -> None: ...
     @property
-    def center(self) -> builtins.list[builtins.float]: ...
-    @center.setter
-    def center(self, value: typing.Sequence[builtins.float]) -> None: ...
-    @property
-    def shape(self) -> builtins.list[builtins.int]: ...
-    @shape.setter
-    def shape(self, value: typing.Sequence[builtins.int]) -> None: ...
-    @property
     def ndim(self) -> builtins.int: ...
     @property
     def matrix(self) -> numpy.typing.NDArray[numpy.float64]: ...
@@ -139,18 +148,12 @@ class Transform:
     def __new__(
         cls,
         parameters: typing.Optional[typing.Sequence[builtins.float]] = None,
-        shape: typing.Optional[typing.Sequence[builtins.int]] = None,
+        dparameters: typing.Optional[typing.Sequence[builtins.float]] = None,
         center: typing.Optional[typing.Sequence[builtins.float]] = None,
     ) -> Transform: ...
     def __getnewargs__(
         self,
-    ) -> tuple[
-        builtins.list[builtins.float],
-        builtins.list[builtins.int],
-        typing.Optional[builtins.list[builtins.float]],
-    ]: ...
-    def __getstate__(self) -> builtins.list[builtins.float]: ...
-    def __setstate__(self, state: typing.Sequence[builtins.float]) -> None: ...
+    ) -> tuple[builtins.list[builtins.float], builtins.list[builtins.float]]: ...
     def __add__(self, other: Transform) -> Transform: ...
     def __radd__(self, other: Transform) -> Transform: ...
     def __sub__(self, other: Transform) -> Transform: ...
@@ -158,11 +161,6 @@ class Transform:
     def __mul__(self, other: Transform | float) -> Transform: ...
     def __rmul__(self, other: Transform | float) -> Transform: ...
     def __truediv__(self, other: builtins.float) -> Transform: ...
-    def adapt(
-        self,
-        center: typing.Sequence[builtins.float],
-        shape: typing.Sequence[builtins.int],
-    ) -> None: ...
     @staticmethod
     def from_scaling(scaling: typing.Sequence[builtins.float]) -> Transform: ...
     @staticmethod
@@ -179,17 +177,32 @@ class Transform:
         self, theta: builtins.float, center: typing.Sequence[builtins.float]
     ) -> Transform: ...
     def interpolate(
-        self, order: builtins.int, image: numpy.typing.ArrayLike
+        self,
+        order: builtins.int,
+        image: numpy.typing.ArrayLike,
+        default_value: typing.Optional[DefaultValue] = None,
     ) -> numpy.typing.NDArray[numpy.float64]: ...
     def interpolate_par(
-        self, order: builtins.int, image: numpy.typing.ArrayLike
+        self,
+        order: builtins.int,
+        image: numpy.typing.ArrayLike,
+        default_value: typing.Optional[DefaultValue] = None,
     ) -> numpy.typing.NDArray[numpy.float64]: ...
+    def frame(
+        self,
+        image: numpy.typing.ArrayLike,
+        default: typing.Optional[typing.Any] = None,
+        default_value: typing.Optional[DefaultValue] = None,
+    ) -> numpy.ndarray:
+        r"""
+        backwards compatibility, default parameter is not implemented
+        """
     def is_unity(self) -> builtins.bool: ...
     def transform_point(
-        self, point: numpy.typing.ArrayLike
+        self, point: numpy.typing.ArrayLike, center: numpy.typing.ArrayLike
     ) -> numpy.typing.NDArray[numpy.float64]: ...
     def transform_points(
-        self, points: numpy.typing.ArrayLike
+        self, points: numpy.typing.ArrayLike, center: numpy.typing.ArrayLike
     ) -> numpy.typing.NDArray[numpy.float64]: ...
     @staticmethod
     def register(
@@ -206,9 +219,25 @@ class Transform:
     def register_translation(
         fixed: numpy.typing.ArrayLike, moving: numpy.typing.ArrayLike
     ) -> Transform: ...
+    @staticmethod
+    def register_translation_fft(
+        fixed: numpy.typing.ArrayLike, moving: numpy.typing.ArrayLike
+    ) -> Transform: ...
 
 class Transforms:
     @property
+    def interpolation_order(self) -> typing.Optional[builtins.int]:
+        r"""
+        the bspline order used for interpolation when using transforms,
+        None -> 0 (nearest neighbor) for int dtypes, 3 (cubic) for float dtypes
+        """
+    @interpolation_order.setter
+    def interpolation_order(self, value: typing.Optional[builtins.int]) -> None:
+        r"""
+        the bspline order used for interpolation when using transforms
+        None -> 0 (nearest neighbor) for int dtypes, 3 (cubic) for float dtypes
+        """
+    @property
     def channel_transforms(self) -> builtins.dict[builtins.str, Transform]: ...
     @channel_transforms.setter
     def channel_transforms(
@@ -223,10 +252,48 @@ class Transforms:
         cls,
         channel: typing.Mapping[builtins.str, Transform] = {},
         drift: typing.Sequence[Transform] = [],
+        interpolation_order: typing.Optional[builtins.int] = None,
+        default_value: typing.Optional[
+            typing.Mapping[builtins.str, DefaultValue]
+        ] = None,
     ) -> Transforms: ...
     def __getnewargs__(
         self,
-    ) -> tuple[builtins.dict[builtins.str, Transform], builtins.list[Transform]]: ...
+    ) -> tuple[
+        builtins.dict[builtins.str, Transform],
+        builtins.list[Transform],
+        typing.Optional[builtins.int],
+        builtins.dict[builtins.str, DefaultValue],
+    ]: ...
+    def with_interpolation_order(
+        self, order: typing.Optional[builtins.int]
+    ) -> Transforms:
+        r"""
+        the bspline order used for interpolation when using transforms
+        None -> 0 (nearest neighbor) for int dtypes, 3 (cubic) for float dtypes
+        """
+    def get_default_interpolation_value(
+        self, channel_name: builtins.str
+    ) -> DefaultValue:
+        r"""
+        the value of the interpolated image outside the data
+        """
+    def set_default_interpolation_value(
+        self,
+        channel_name: builtins.str,
+        value: float | DefaultValue | collections.abc.Sequence[float],
+    ) -> None:
+        r"""
+        the value of the interpolated image outside the data
+        """
+    def with_default_interpolation_value(
+        self,
+        channel_name: builtins.str,
+        value: float | DefaultValue | collections.abc.Sequence[float],
+    ) -> Transforms:
+        r"""
+        the value of the interpolated image outside the data
+        """
     @staticmethod
     def load(path: builtins.str | os.PathLike | pathlib.Path) -> Transforms: ...
     def save(self, path: builtins.str | os.PathLike | pathlib.Path) -> None: ...

ndbioimage/py/ndbioimage/transforms.pyi

@@ -1,4 +1,5 @@
 from .ndbioimage_rs.transforms import (
+    DefaultValue,
     Optimizer,
     RegistrationStep,
     SamplingArg,
@@ -8,6 +9,7 @@ from .ndbioimage_rs.transforms import (
 )
 
 __all__ = [
+    "DefaultValue",
     "Optimizer",
     "RegistrationStep",
     "SamplingArg",

ndbioimage/pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "maturin"
 
 [project]
 name = "ndbioimage"
-version = "2027.0.5"
+version = "2027.0.8"
 requires-python = ">=3.10"
 classifiers = [
     "License :: OSI Approved :: MIT License",
@@ -35,7 +35,7 @@ ndbioimage_generate_stub = "ndbioimage:ndbioimage_generate_stub"
 
 [tool.maturin]
 python-source = "py"
-features = ["python", "bioformats_java", "gpl-formats", "czi", "tiff", "tiffseq", "tiffwrite", "movie", "transforms"]
+features = ["python", "bioformats_java", "gpl-formats", "bioformats_rust", "czi", "tiff", "tiffseq", "tiffwrite", "movie", "transforms"]
 no-default-features = true
 module-name = "ndbioimage.ndbioimage_rs"
 include = ["py/ndbioimage/jassets/j4rs*", "py/ndbioimage/deps/libj4rs*"]

ndbioimage/src/axes.rs

@@ -196,11 +196,11 @@ impl<'de> DeserializeAs<'de, SliceInfoElem> for SliceInfoElemDef {
     }
 }
 
-#[derive(Clone, Debug)]
+#[derive(Clone, Debug, Copy)]
 pub(crate) struct Slice {
-    start: isize,
-    end: isize,
-    step: isize,
+    pub(crate) start: isize,
+    pub(crate) end: isize,
+    pub(crate) step: isize,
 }
 
 impl Slice {
@@ -236,7 +236,7 @@ impl IntoIterator for &Slice {
     type IntoIter = Slice;
 
     fn into_iter(self) -> Self::IntoIter {
-        self.clone()
+        *self
     }
 }
 
@@ -314,16 +314,16 @@ impl Index<&Axis> for Shape {
 impl Index<usize> for Shape {
     type Output = usize;
 
-    fn index(&self, dim: usize) -> &Self::Output {
-        &self[self.order[dim % self.order.len()]]
+    fn index(&self, index: usize) -> &Self::Output {
+        &self[self.order[index]]
     }
 }
 
 impl Index<&usize> for Shape {
     type Output = usize;
 
-    fn index(&self, dim: &usize) -> &Self::Output {
-        &self[self.order[dim % self.order.len()]]
+    fn index(&self, index: &usize) -> &Self::Output {
+        &self[self.order[*index]]
     }
 }
 
@@ -406,6 +406,10 @@ impl Shape {
         }
     }
 
+    pub fn get(&self, index: &usize) -> Option<&usize> {
+        self.order.get(*index).map(|axis| &self[axis])
+    }
+
     /// the number of axes in the shape
     pub fn len(&self) -> usize {
         self.order.len()
@@ -421,6 +425,11 @@ impl Shape {
         self.order.iter().map(|axis| self[axis]).collect()
     }
 
+    /// the sizes in cztyx order
+    pub fn cztyx(&self) -> [usize; 5] {
+        [self.c, self.z, self.t, self.y, self.x]
+    }
+
     /// the sizes along the axes as a map from axis to size
     pub fn to_hashmap(&self) -> HashMap<Axis, usize> {
         let mut map = HashMap::new();

ndbioimage/src/cache.rs

@@ -1,13 +1,22 @@
 use crate::axes::{Axis, Operation};
+use crate::readers::Frame;
+#[cfg(feature = "transforms")]
+use crate::transforms::{DefaultValues, Transforms};
+#[cfg(feature = "transforms")]
+use image_registration::transform::Transform;
 use indexmap::{Equivalent, IndexMap};
+#[cfg(feature = "transforms")]
+use ndarray::IxDyn;
 use ndarray::{ArrayD, SliceInfoElem};
 use std::any::Any;
+#[cfg(feature = "transforms")]
+use std::borrow::Cow;
+use std::cmp::PartialEq;
 use std::fmt::{Debug, Formatter};
 use std::hash::{Hash, Hasher};
 use std::path::{Path, PathBuf};
 use std::sync::{Arc, Mutex, OnceLock};
-
-use crate::readers::Frame;
+use std::time::SystemTime;
 
 /// maximum number of frames held in the cache
 const DEFAULT_FRAME_CACHE_SIZE: usize = 128;
@@ -15,6 +24,68 @@ const DEFAULT_FRAME_CACHE_SIZE: usize = 128;
 /// maximum number of materialized arrays held in the cache
 const DEFAULT_ARRAY_CACHE_SIZE: usize = 32;
 
+/// fingerprint of the file(s) a reader reads from.
+///
+/// The cache is keyed partly on this stamp, so a file that is rewritten (for
+/// example an output file that is written and then read back) gets a new cache
+/// key instead of returning stale pixels. The stamp is read from the
+/// filesystem on every lookup, so a long-lived view also notices a change.
+#[derive(Debug, Clone, Copy, Default, PartialEq, Eq, Hash)]
+pub(crate) struct FileStamp {
+    /// size of the file in bytes
+    len: u64,
+    /// last modification time
+    modified: Option<SystemTime>,
+    /// last status change time; on unix this catches same-size rewrites that
+    /// happen within the timestamp resolution of `modified`
+    #[cfg(unix)]
+    changed: Option<(i64, u32)>,
+    /// inode number; on unix this catches an atomic replace that preserves the
+    /// modification time and size
+    #[cfg(unix)]
+    inode: u64,
+    /// creation time; windows has no stable `ctime` or file index, but this
+    /// changes when a file is deleted/recreated or atomically replaced, which is
+    /// the case the unix inode covers. in-place rewrites are caught by
+    /// `modified`, which has 100 ns resolution on windows.
+    #[cfg(windows)]
+    created: Option<SystemTime>,
+}
+
+impl FileStamp {
+    /// read the current stamp for `path`, returning the default (all-zero)
+    /// stamp when the metadata cannot be read
+    pub(crate) fn from_path(path: &Path) -> Self {
+        let Ok(meta) = std::fs::metadata(path) else {
+            return Self::default();
+        };
+        let len = meta.len();
+        let modified = meta.modified().ok();
+        #[cfg(unix)]
+        {
+            use std::os::unix::fs::MetadataExt;
+            Self {
+                len,
+                modified,
+                changed: Some((meta.ctime(), meta.ctime_nsec() as u32)),
+                inode: meta.ino(),
+            }
+        }
+        #[cfg(windows)]
+        {
+            Self {
+                len,
+                modified,
+                created: meta.created().ok(),
+            }
+        }
+        #[cfg(not(any(unix, windows)))]
+        {
+            Self { len, modified }
+        }
+    }
+}
+
 /// identity of the reader a frame was read from
 #[derive(Debug, Clone, PartialEq, Eq, Hash)]
 pub(crate) struct ReaderKey {
@@ -22,6 +93,7 @@ pub(crate) struct ReaderKey {
     pub(crate) path: PathBuf,
     pub(crate) series: usize,
     pub(crate) position: usize,
+    pub(crate) stamp: FileStamp,
 }
 
 /// borrowed view of [`ReaderKey`] for cache lookups without allocation.
@@ -33,6 +105,7 @@ pub(crate) struct ReaderKeyRef<'a> {
     pub(crate) path: &'a Path,
     pub(crate) series: usize,
     pub(crate) position: usize,
+    pub(crate) stamp: FileStamp,
 }
 
 impl Hash for ReaderKeyRef<'_> {
@@ -41,6 +114,7 @@ impl Hash for ReaderKeyRef<'_> {
         self.path.hash(state);
         self.series.hash(state);
         self.position.hash(state);
+        self.stamp.hash(state);
     }
 }
 
@@ -50,6 +124,7 @@ impl Equivalent<ReaderKey> for ReaderKeyRef<'_> {
             && self.path == key.path.as_path()
             && self.series == key.series
             && self.position == key.position
+            && self.stamp == key.stamp
     }
 }
 
@@ -60,11 +135,38 @@ impl ReaderKeyRef<'_> {
             path: self.path.to_path_buf(),
             series: self.series,
             position: self.position,
+            stamp: self.stamp,
         }
     }
 }
 
-pub(crate) type FrameKey = (ReaderKey, usize, usize, usize);
+#[derive(Debug, PartialEq, Eq)]
+pub(crate) struct FrameKey {
+    pub(crate) reader: ReaderKey,
+    pub(crate) c: usize,
+    pub(crate) z: usize,
+    pub(crate) t: usize,
+    #[cfg(feature = "transforms")]
+    pub(crate) transform: Option<(Transform<IxDyn>, usize, DefaultValues)>,
+}
+
+impl Hash for FrameKey {
+    fn hash<H: Hasher>(&self, state: &mut H) {
+        self.reader.hash(state);
+        self.c.hash(state);
+        self.z.hash(state);
+        self.t.hash(state);
+        #[cfg(feature = "transforms")]
+        if let Some(transform) = self.transform.as_ref() {
+            transform.0.hash(state);
+            transform.1.hash(state);
+            transform.2.iter().for_each(|(k, v)| {
+                k.hash(state);
+                v.hash(state);
+            })
+        }
+    }
+}
 
 /// borrowed view of [`FrameKey`] for cache lookups without allocation.
 #[derive(Debug, PartialEq, Eq)]
@@ -73,6 +175,8 @@ pub(crate) struct FrameKeyRef<'a> {
     pub(crate) c: usize,
     pub(crate) z: usize,
     pub(crate) t: usize,
+    #[cfg(feature = "transforms")]
+    pub(crate) transform: Option<(Cow<'a, Transform<IxDyn>>, usize, Cow<'a, DefaultValues>)>,
 }
 
 impl Hash for FrameKeyRef<'_> {
@@ -81,19 +185,60 @@ impl Hash for FrameKeyRef<'_> {
         self.c.hash(state);
         self.z.hash(state);
         self.t.hash(state);
+        #[cfg(feature = "transforms")]
+        if let Some(transform) = self.transform.as_ref() {
+            transform.0.hash(state);
+            transform.1.hash(state);
+            transform.2.iter().for_each(|(k, v)| {
+                k.hash(state);
+                v.hash(state);
+            })
+        }
     }
 }
 
 impl Equivalent<FrameKey> for FrameKeyRef<'_> {
     fn equivalent(&self, key: &FrameKey) -> bool {
-        let (rk, c, z, t) = key;
-        self.reader.equivalent(rk) && self.c == *c && self.z == *z && self.t == *t
+        let FrameKey {
+            reader: reader_key,
+            c,
+            z,
+            t,
+            #[cfg(feature = "transforms")]
+            transform,
+        } = key;
+        #[cfg(feature = "transforms")]
+        {
+            self.reader.equivalent(reader_key)
+                && self.c == *c
+                && self.z == *z
+                && self.t == *t
+                && self
+                    .transform
+                    .as_ref()
+                    .map(|(t, i, d)| (t.as_ref(), i, d.as_ref()))
+                    == transform.as_ref().map(|(t, i, d)| (t, i, d))
+        }
+        #[cfg(not(feature = "transforms"))]
+        {
+            self.reader.equivalent(reader_key) && self.c == *c && self.z == *z && self.t == *t
+        }
     }
 }
 
 impl FrameKeyRef<'_> {
     pub(crate) fn to_owned(&self) -> FrameKey {
-        (self.reader.to_owned(), self.c, self.z, self.t)
+        FrameKey {
+            reader: self.reader.to_owned(),
+            c: self.c,
+            z: self.z,
+            t: self.t,
+            #[cfg(feature = "transforms")]
+            transform: self
+                .transform
+                .as_ref()
+                .map(|(t, i, d)| (t.clone().into_owned(), *i, d.clone().into_owned())),
+        }
     }
 }
 
@@ -105,6 +250,8 @@ pub(crate) struct ArrayKey {
     pub(crate) slice: Vec<SliceInfoElem>,
     pub(crate) axes: Vec<Axis>,
     pub(crate) operations: Vec<(Axis, Operation)>,
+    #[cfg(feature = "transforms")]
+    pub(crate) transforms: Transforms,
 }
 
 /// borrowed view of [`ArrayKey`] for cache lookups without allocation.
@@ -115,6 +262,8 @@ pub(crate) struct ArrayKeyRef<'a> {
     pub(crate) slice: &'a [SliceInfoElem],
     pub(crate) axes: &'a [Axis],
     pub(crate) operations: &'a IndexMap<Axis, Operation>,
+    #[cfg(feature = "transforms")]
+    pub(crate) transforms: &'a Transforms,
 }
 
 impl Hash for ArrayKeyRef<'_> {
@@ -128,21 +277,40 @@ impl Hash for ArrayKeyRef<'_> {
             ax.hash(state);
             op.hash(state);
         }
+        #[cfg(feature = "transforms")]
+        self.transforms.hash(state);
     }
 }
 
 impl Equivalent<ArrayKey> for ArrayKeyRef<'_> {
     fn equivalent(&self, key: &ArrayKey) -> bool {
-        self.reader.equivalent(&key.reader)
-            && self.dtype == key.dtype
-            && self.slice == key.slice
-            && self.axes == key.axes
-            && self.operations.len() == key.operations.len()
-            && self
-                .operations
-                .iter()
-                .zip(&key.operations)
-                .all(|((ax, op), (key_ax, key_op))| ax == key_ax && op == key_op)
+        #[cfg(feature = "transforms")]
+        {
+            self.reader.equivalent(&key.reader)
+                && self.dtype == key.dtype
+                && self.slice == key.slice
+                && self.axes == key.axes
+                && self.operations.len() == key.operations.len()
+                && self
+                    .operations
+                    .iter()
+                    .zip(&key.operations)
+                    .all(|((ax, op), (key_ax, key_op))| ax == key_ax && op == key_op)
+                && self.transforms == key.transforms
+        }
+        #[cfg(not(feature = "transforms"))]
+        {
+            self.reader.equivalent(&key.reader)
+                && self.dtype == key.dtype
+                && self.slice == key.slice
+                && self.axes == key.axes
+                && self.operations.len() == key.operations.len()
+                && self
+                    .operations
+                    .iter()
+                    .zip(&key.operations)
+                    .all(|((ax, op), (key_ax, key_op))| ax == key_ax && op == key_op)
+        }
     }
 }
 
@@ -158,6 +326,8 @@ impl ArrayKeyRef<'_> {
                 .iter()
                 .map(|(ax, op)| (*ax, op.clone()))
                 .collect(),
+            #[cfg(feature = "transforms")]
+            transforms: self.transforms.to_owned(),
         }
     }
 }
@@ -320,3 +490,42 @@ impl Debug for ArrayCache {
         f.debug_struct("ArrayCache").finish_non_exhaustive()
     }
 }
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    fn temp_path(name: &str) -> PathBuf {
+        let dir = std::env::temp_dir().join(format!("ndbioimage-cache-{}", std::process::id()));
+        std::fs::create_dir_all(&dir).unwrap();
+        dir.join(name)
+    }
+
+    #[test]
+    fn file_stamp_detects_rewrite() {
+        let path = temp_path("stamp.bin");
+        std::fs::write(&path, b"aaaa").unwrap();
+        let before = FileStamp::from_path(&path);
+        std::fs::write(&path, b"bbbbbb").unwrap();
+        let after = FileStamp::from_path(&path);
+        assert_ne!(before, after);
+        let _ = std::fs::remove_file(&path);
+    }
+
+    #[test]
+    fn reader_key_includes_file_stamp() {
+        let key = |stamp| ReaderKey {
+            name: "test".to_string(),
+            path: PathBuf::from("test.tif"),
+            series: 0,
+            position: 0,
+            stamp,
+        };
+        let a = key(FileStamp::from_path(Path::new("does-not-exist.tif")));
+        let b = key(FileStamp {
+            len: 1,
+            ..FileStamp::default()
+        });
+        assert_ne!(a, b);
+    }
+}

ndbioimage/src/error.rs

@@ -47,7 +47,7 @@ pub enum Error {
     #[error(transparent)]
     Tiff(#[from] tiff::TiffError),
     /// a postcard (de)serialization error
-    #[cfg(feature = "python")]
+    #[cfg(any(test, feature = "python"))]
     #[error(transparent)]
     PostCard(#[from] postcard::Error),
     /// an error from the libczi binding
@@ -65,10 +65,6 @@ pub enum Error {
     #[cfg(feature = "czi")]
     #[error(transparent)]
     XmlTreeParse(#[from] xmltree::ParseError),
-    /// a czi-specific error
-    #[cfg(feature = "czi")]
-    #[error(transparent)]
-    Czi(#[from] crate::readers::czi::CziError),
     /// an error joining a tokio task
     #[cfg(feature = "movie")]
     #[error(transparent)]
@@ -121,6 +117,9 @@ pub enum Error {
     /// the view is empty
     #[error("empty view")]
     EmptyView,
+    /// the view cannot be parsed as a scalar
+    #[error("not a scalar")]
+    NotScalar,
     /// the color string could not be parsed
     #[error("invalid color: {0}")]
     InvalidColor(String),
@@ -172,6 +171,24 @@ pub enum Error {
     /// file mismatch
     #[error("{0} do not match in {1} and {2}")]
     FileMismatch(String, String, String),
+    /// a jvm created in another process was inherited across a fork
+    #[error(
+        "the jvm was created in process {0} but this is process {1}: a forked child inherits a jvm \
+         it cannot use. read the file before forking, or start the child with spawn/exec instead."
+    )]
+    /// jvm was forked, so it cannot run in this process
+    ForkedJvm(u32, u32),
+    /// only 2d and 3d transformations are implemented
+    #[error("{0}D transformation is not implemented")]
+    TransformationNotImplemented(usize),
+    /// cannot find the name for some channel because it does not exist
+    #[error("no channel {0} present")]
+    NoChannel(usize),
+    /// no valid blocks in czi file
+    #[error("czi file has no valid blocks")]
+    CziNoValidBlocks,
+    #[error("f128 is not supported")]
+    F128NotSupported,
 }
 
 impl Error {

ndbioimage/src/lib.rs

@@ -6,8 +6,8 @@
 //! structure. Additionally, it can automatically calculate an affine transform that corrects for
 //! chromatic aberrations etc. and apply it on the fly to the image.
 //!
-//! Currently, it supports imagej tif files, czi files, micromanager tif sequences and anything
-//! [bioformats](https://www.openmicroscopy.org/bio-formats/) can handle.
+//! Currently, it supports imagej and ome tif files, czi files, micromanager tif sequences and
+//! anything [bioformats](https://www.openmicroscopy.org/bio-formats/) can handle.
 //!
 //! ```rust,no_run
 //! use ndarray::Array2;
@@ -16,7 +16,7 @@
 //! # fn main() -> Result<(), ndbioimage::error::Error> {
 //! let path = "/path/to/file";
 //! let reader = DynReader::new(&path, 0, 0)?;
-//! println!("shape: {}", reader.shape());
+//! println!("shape: {}", reader.raw_shape());
 //! let frame = reader.get_frame(0, 0, 0)?;
 //! if let Ok(arr) = <Frame as TryInto<Array2<i8>>>::try_into(frame) {
 //!     println!("{:?}", arr);
@@ -70,7 +70,7 @@ pub mod readers;
 pub mod tiffwrite;
 #[cfg(feature = "transforms")]
 /// image registration and warping images using affine transforms
-mod transforms;
+pub mod transforms;
 mod utils;
 
 /// main entry point for the application
@@ -338,7 +338,7 @@ mod tests {
     fn read_ser() -> Result<(), Error> {
         let file = "czi/Experiment-2029.czi";
         let reader = open(file)?;
-        println!("shape: {}", reader.shape());
+        println!("shape: {}", reader.raw_shape());
         let frame = reader.get_frame(0, 0, 0)?;
         if let Ok(arr) = <Frame as TryInto<Array2<i8>>>::try_into(frame) {
             println!("{:?}", arr);

ndbioimage/src/py.rs

@@ -40,7 +40,7 @@ pub(crate) fn ndbioimage_file() -> PathBuf {
     PathBuf::from(file)
 }
 
-/// generates ndbioimage/__init__.pyi
+/// generates .pyi files
 #[gen_stub_pyfunction]
 #[pyfunction]
 pub fn generate_stub(dest_path: String) -> PyResult<()> {
@@ -67,6 +67,7 @@ fn main() -> PyResult<()> {
 #[pymodule]
 #[pyo3(name = "ndbioimage_rs")]
 mod ndbioimage_rs {
+    #[cfg(feature = "transforms")]
     use pyo3::prelude::*;
 
     #[pymodule_export]
@@ -100,6 +101,9 @@ mod ndbioimage_rs {
         #[pymodule_export]
         use super::super::transforms::PyTransforms;
 
+        #[pymodule_export]
+        use super::super::transforms::PyDefaultValue;
+
         #[pymodule_export]
         use super::super::transforms::PyRegistrationStep;
 

ndbioimage/src/py/transforms.rs

@@ -1,13 +1,17 @@
 use crate::error::Error;
+use crate::py::imread::PyView;
 use crate::transforms::Transforms;
+use image_registration::bspline::DefaultValue;
 use image_registration::metric::{SamplingArg, Sigma};
 use image_registration::register::{Optimizer, RegistrationStep};
 use image_registration::transform::Transform;
 use ndarray::{Ix0, Ix1, Ix2, Ix3, Ix4, Ix5, IxDyn};
 use numpy::{
-    AllowTypeChange, IntoPyArray, PyArray, PyArray1, PyArray2, PyArrayDyn, PyArrayLike1,
-    PyArrayLike2, PyArrayLikeDyn,
+    AllowTypeChange, IntoPyArray, PyArray, PyArray1, PyArray2, PyArrayDescrMethods, PyArrayDyn,
+    PyArrayLike1, PyArrayLike2, PyArrayLikeDyn, PyArrayMethods, PyUntypedArray,
+    PyUntypedArrayMethods, dtype,
 };
+use pyo3::IntoPyObjectExt;
 use pyo3::exceptions::{PyNotImplementedError, PyValueError};
 use pyo3::prelude::*;
 use pyo3_stub_gen::derive::*;
@@ -15,6 +19,45 @@ use serde::{Deserialize, Serialize};
 use std::collections::HashMap;
 use std::path::PathBuf;
 
+#[gen_stub_pyclass]
+#[pyclass(
+    subclass,
+    from_py_object,
+    eq,
+    name = "DefaultValue",
+    module = "ndbioimage.transforms"
+)]
+#[derive(Clone, Debug, Serialize, Deserialize, PartialEq)]
+pub(crate) struct PyDefaultValue {
+    inner: DefaultValue,
+}
+
+#[gen_stub_pymethods]
+#[pymethods]
+impl PyDefaultValue {
+    #[new]
+    #[pyo3(signature = (values, weights = None))]
+    fn new(values: Vec<f64>, weights: Option<Vec<f64>>) -> Self {
+        Self {
+            inner: DefaultValue::new(values, weights.unwrap_or_default()),
+        }
+    }
+
+    #[getter]
+    fn values(&self) -> Vec<f64> {
+        self.inner.values().to_vec()
+    }
+
+    #[getter]
+    fn weights(&self) -> Vec<f64> {
+        self.inner.weights().to_vec()
+    }
+
+    fn __getnewargs__(&self) -> (Vec<f64>, Vec<f64>) {
+        (self.inner.values().to_vec(), self.inner.weights().to_vec())
+    }
+}
+
 #[gen_stub_pyclass]
 #[pyclass(
     subclass,
@@ -32,41 +75,31 @@ pub(crate) struct PyTransform {
 #[pymethods]
 impl PyTransform {
     #[new]
-    #[pyo3(signature = (parameters = None, shape = None, center = None))]
+    #[pyo3(signature = (parameters = None, dparameters = None, center = None))]
     fn new(
         parameters: Option<Vec<f64>>,
-        shape: Option<Vec<usize>>,
+        dparameters: Option<Vec<f64>>,
         center: Option<Vec<f64>>,
     ) -> Self {
         let parameters = parameters.unwrap_or_else(|| vec![1.0, 0.0, 0.0, 1.0, 0.0, 0.0]);
-        let shape = shape.unwrap_or_else(|| vec![512, 512]);
-        if let Some(center) = center {
-            Self {
-                inner: Transform::new_with_center(parameters, center, shape),
-            }
+        let mut transform = if let Some(center) = center {
+            Transform::new_around_center(parameters, &center)
         } else {
-            Self {
-                inner: Transform::new(parameters, shape),
-            }
+            Transform::new(parameters)
+        };
+        if let Some(dparameters) = dparameters {
+            transform.dparameters = dparameters;
         }
+        Self { inner: transform }
     }
 
-    fn __getnewargs__(&self) -> (Vec<f64>, Vec<usize>, Option<Vec<f64>>) {
+    fn __getnewargs__(&self) -> (Vec<f64>, Vec<f64>) {
         (
             self.inner.parameters.clone(),
-            self.inner.shape.clone(),
-            Some(self.inner.center.clone()),
+            self.inner.dparameters.clone(),
         )
     }
 
-    fn __getstate__(&self) -> Vec<f64> {
-        self.inner.dparameters.clone()
-    }
-
-    fn __setstate__(&mut self, state: Vec<f64>) {
-        self.inner.dparameters = state;
-    }
-
     fn __add__(&self, other: &PyTransform) -> PyTransform {
         PyTransform {
             inner: &self.inner + &other.inner,
@@ -96,9 +129,9 @@ impl PyTransform {
         py: Python,
         #[gen_stub(override_type(type_repr = "Transform | float"))] other: &Bound<PyAny>,
     ) -> PyResult<PyTransform> {
-        if other.is_instance_of::<PyTransform>() {
+        if let Ok(other) = other.extract::<PyTransform>() {
             Ok(PyTransform {
-                inner: &self.inner * &other.extract::<PyTransform>()?.inner,
+                inner: &self.inner * other.inner,
             })
         } else {
             let builtins = PyModule::import(py, "builtins")?;
@@ -114,9 +147,9 @@ impl PyTransform {
         py: Python,
         #[gen_stub(override_type(type_repr = "Transform | float"))] other: &Bound<PyAny>,
     ) -> PyResult<PyTransform> {
-        if other.is_instance_of::<PyTransform>() {
+        if let Ok(other) = other.extract::<PyTransform>() {
             Ok(PyTransform {
-                inner: &other.extract::<PyTransform>()?.inner * &self.inner,
+                inner: &other.inner * &self.inner,
             })
         } else {
             let builtins = PyModule::import(py, "builtins")?;
@@ -155,26 +188,6 @@ impl PyTransform {
         self.inner.dparameters = dparameters;
     }
 
-    #[getter]
-    fn get_center(&self) -> Vec<f64> {
-        self.inner.center.clone()
-    }
-
-    #[setter]
-    fn set_center(&mut self, center: Vec<f64>) {
-        self.inner.center = center;
-    }
-
-    #[getter]
-    fn get_shape(&self) -> Vec<usize> {
-        self.inner.shape.clone()
-    }
-
-    #[setter]
-    fn set_shape(&mut self, shape: Vec<usize>) {
-        self.inner.shape = shape;
-    }
-
     #[getter]
     fn get_ndim(&self) -> usize {
         self.inner.ndim()
@@ -217,10 +230,6 @@ impl PyTransform {
         })
     }
 
-    fn adapt(&mut self, center: Vec<f64>, shape: Vec<usize>) {
-        self.inner.adapt(center.as_slice(), shape.as_slice());
-    }
-
     #[staticmethod]
     fn from_scaling(scaling: Vec<f64>) -> PyTransform {
         PyTransform {
@@ -261,74 +270,178 @@ impl PyTransform {
         }
     }
 
+    #[pyo3(signature = (order, image, default_value = None))]
     fn interpolate<'py>(
         &self,
         py: Python<'py>,
         order: usize,
         #[gen_stub(override_type(type_repr="numpy.typing.ArrayLike", imports=("numpy.typing")))]
         image: PyArrayLikeDyn<f64, AllowTypeChange>,
+        default_value: Option<PyDefaultValue>,
     ) -> PyResult<Bound<'py, PyArrayDyn<f64>>> {
         let image = image.as_array();
+        let default_value = default_value.map(|i| i.inner).unwrap_or_default();
         Ok(py
             .detach(|| match order {
-                0 => Ok(self
-                    .inner
-                    .interpolate::<0, _, _>(image.into_dimensionality().map_err(Error::from)?)),
-                1 => Ok(self
-                    .inner
-                    .interpolate::<1, _, _>(image.into_dimensionality().map_err(Error::from)?)),
-                2 => Ok(self
-                    .inner
-                    .interpolate::<2, _, _>(image.into_dimensionality().map_err(Error::from)?)),
-                3 => Ok(self
-                    .inner
-                    .interpolate::<3, _, _>(image.into_dimensionality().map_err(Error::from)?)),
-                4 => Ok(self
-                    .inner
-                    .interpolate::<4, _, _>(image.into_dimensionality().map_err(Error::from)?)),
-                5 => Ok(self
-                    .inner
-                    .interpolate::<5, _, _>(image.into_dimensionality().map_err(Error::from)?)),
+                0 => Ok(self.inner.interpolate::<0, _, _>(
+                    image.into_dimensionality().map_err(Error::from)?,
+                    &default_value,
+                )),
+                1 => Ok(self.inner.interpolate::<1, _, _>(
+                    image.into_dimensionality().map_err(Error::from)?,
+                    &default_value,
+                )),
+                2 => Ok(self.inner.interpolate::<2, _, _>(
+                    image.into_dimensionality().map_err(Error::from)?,
+                    &default_value,
+                )),
+                3 => Ok(self.inner.interpolate::<3, _, _>(
+                    image.into_dimensionality().map_err(Error::from)?,
+                    &default_value,
+                )),
+                4 => Ok(self.inner.interpolate::<4, _, _>(
+                    image.into_dimensionality().map_err(Error::from)?,
+                    &default_value,
+                )),
+                5 => Ok(self.inner.interpolate::<5, _, _>(
+                    image.into_dimensionality().map_err(Error::from)?,
+                    &default_value,
+                )),
                 _ => Err(PyValueError::new_err("order must be 0 <= order < 6")),
             })?
             .map_err(Error::from)?
             .into_pyarray(py))
     }
 
+    #[pyo3(signature = (order, image, default_value = None))]
     fn interpolate_par<'py>(
         &self,
         py: Python<'py>,
         order: usize,
         #[gen_stub(override_type(type_repr="numpy.typing.ArrayLike", imports=("numpy.typing")))]
         image: PyArrayLikeDyn<f64, AllowTypeChange>,
+        default_value: Option<PyDefaultValue>,
     ) -> PyResult<Bound<'py, PyArrayDyn<f64>>> {
         let image = image.as_array();
+        let default_value = default_value.map(|i| i.inner).unwrap_or_default();
         Ok(py
             .detach(|| match order {
-                0 => Ok(self
-                    .inner
-                    .interpolate_par::<0, _, _>(image.into_dimensionality().map_err(Error::from)?)),
-                1 => Ok(self
-                    .inner
-                    .interpolate_par::<1, _, _>(image.into_dimensionality().map_err(Error::from)?)),
-                2 => Ok(self
-                    .inner
-                    .interpolate_par::<2, _, _>(image.into_dimensionality().map_err(Error::from)?)),
-                3 => Ok(self
-                    .inner
-                    .interpolate_par::<3, _, _>(image.into_dimensionality().map_err(Error::from)?)),
-                4 => Ok(self
-                    .inner
-                    .interpolate_par::<4, _, _>(image.into_dimensionality().map_err(Error::from)?)),
-                5 => Ok(self
-                    .inner
-                    .interpolate_par::<5, _, _>(image.into_dimensionality().map_err(Error::from)?)),
+                0 => Ok(self.inner.interpolate_par::<0, _, _>(
+                    image.into_dimensionality().map_err(Error::from)?,
+                    &default_value,
+                )),
+                1 => Ok(self.inner.interpolate_par::<1, _, _>(
+                    image.into_dimensionality().map_err(Error::from)?,
+                    &default_value,
+                )),
+                2 => Ok(self.inner.interpolate_par::<2, _, _>(
+                    image.into_dimensionality().map_err(Error::from)?,
+                    &default_value,
+                )),
+                3 => Ok(self.inner.interpolate_par::<3, _, _>(
+                    image.into_dimensionality().map_err(Error::from)?,
+                    &default_value,
+                )),
+                4 => Ok(self.inner.interpolate_par::<4, _, _>(
+                    image.into_dimensionality().map_err(Error::from)?,
+                    &default_value,
+                )),
+                5 => Ok(self.inner.interpolate_par::<5, _, _>(
+                    image.into_dimensionality().map_err(Error::from)?,
+                    &default_value,
+                )),
                 _ => Err(PyValueError::new_err("order must be 0 <= order < 6")),
             })?
             .map_err(Error::from)?
             .into_pyarray(py))
     }
 
+    /// backwards compatibility, default parameter is not implemented
+    #[gen_stub(override_return_type(type_repr="numpy.ndarray", imports=("numpy")))]
+    #[pyo3(signature = (image, default = None, default_value = None))]
+    fn frame<'py>(
+        &self,
+        py: Python<'py>,
+        #[gen_stub(override_type(type_repr="numpy.typing.ArrayLike", imports=("numpy.typing")))]
+        image: Bound<'py, PyUntypedArray>,
+        default: Option<Bound<'py, PyAny>>,
+        default_value: Option<PyDefaultValue>,
+    ) -> PyResult<Bound<'py, PyAny>> {
+        if default.is_some() {
+            return Err(PyNotImplementedError::new_err(
+                "default parameter is not implemented",
+            ));
+        }
+        if image.ndim() != self.inner.ndim() {
+            return Err(PyValueError::new_err(format!(
+                "dimensionality of image ({}) and transform ({}) do not correspond",
+                image.ndim(),
+                self.inner.ndim()
+            )));
+        }
+        let default_value = default_value.map(|i| i.inner).unwrap_or_default();
+
+        macro_rules! interp {
+            ($order:expr, $t:tt) => {
+                self.inner
+                    .interpolate_par::<$order, _, _>(
+                        unsafe { image.cast::<PyArrayDyn<$t>>()?.as_array() }
+                            .mapv(|i| i as f64)
+                            .into_dimensionality()
+                            .map_err(Error::from)?
+                            .view(),
+                        &default_value,
+                    )
+                    .map_err(Error::from)?
+                    .mapv(|i| i as $t)
+                    .into_pyarray(py)
+                    .into_bound_py_any(py)
+            };
+        }
+
+        let dt = image.dtype();
+        if dt.is_equiv_to(&dtype::<bool>(py)) {
+            self.inner
+                .interpolate_par::<0, _, _>(
+                    unsafe { image.cast::<PyArrayDyn<bool>>()?.as_array() }
+                        .mapv(|i| if i { 1.0 } else { 0.0 })
+                        .into_dimensionality()
+                        .map_err(Error::from)?
+                        .view(),
+                    &default_value,
+                )
+                .map_err(Error::from)?
+                .mapv(|i| i > 0.5)
+                .into_pyarray(py)
+                .into_bound_py_any(py)
+        } else if dt.is_equiv_to(&dtype::<u8>(py)) {
+            interp!(0, u8)
+        } else if dt.is_equiv_to(&dtype::<i8>(py)) {
+            interp!(0, i8)
+        } else if dt.is_equiv_to(&dtype::<u16>(py)) {
+            interp!(0, u16)
+        } else if dt.is_equiv_to(&dtype::<i16>(py)) {
+            interp!(0, i16)
+        } else if dt.is_equiv_to(&dtype::<u32>(py)) {
+            interp!(0, u32)
+        } else if dt.is_equiv_to(&dtype::<i32>(py)) {
+            interp!(0, i32)
+        } else if dt.is_equiv_to(&dtype::<u64>(py)) {
+            interp!(0, u64)
+        } else if dt.is_equiv_to(&dtype::<i64>(py)) {
+            interp!(0, i64)
+        } else if dt.is_equiv_to(&dtype::<f32>(py)) {
+            interp!(3, f32)
+        } else if dt.is_equiv_to(&dtype::<f64>(py)) {
+            interp!(3, f64)
+        } else {
+            Err(PyValueError::new_err(format!(
+                "invalid dtype: {}; must be bool, (u)int(8/16/32/64) or float(32/64)",
+                dt
+            )))
+        }
+    }
+
     fn is_unity(&self) -> bool {
         self.inner.is_unity()
     }
@@ -338,13 +451,19 @@ impl PyTransform {
         py: Python<'py>,
         #[gen_stub(override_type(type_repr="numpy.typing.ArrayLike", imports=("numpy.typing")))]
         point: PyArrayLike1<f64, AllowTypeChange>,
+        #[gen_stub(override_type(type_repr="numpy.typing.ArrayLike", imports=("numpy.typing")))]
+        center: PyArrayLike1<f64, AllowTypeChange>,
     ) -> Bound<'py, PyArray1<f64>> {
         let point = point.as_array();
-        if let Some(slice) = point.as_slice() {
-            self.inner.transform_point(slice)
-        } else {
-            let point = point.to_vec();
-            self.inner.transform_point(&point)
+        let center = center.as_array();
+
+        match (point.as_slice(), center.as_slice()) {
+            (Some(point), Some(center)) => self.inner.transform_point(point, center),
+            (Some(point), None) => self.inner.transform_point(point, &center.to_vec()),
+            (None, Some(center)) => self.inner.transform_point(&point.to_vec(), center),
+            (None, None) => self
+                .inner
+                .transform_point(&point.to_vec(), &center.to_vec()),
         }
         .into_pyarray(py)
     }
@@ -354,10 +473,19 @@ impl PyTransform {
         py: Python<'py>,
         #[gen_stub(override_type(type_repr="numpy.typing.ArrayLike", imports=("numpy.typing")))]
         points: PyArrayLike2<f64, AllowTypeChange>,
+        #[gen_stub(override_type(type_repr="numpy.typing.ArrayLike", imports=("numpy.typing")))]
+        center: PyArrayLike1<f64, AllowTypeChange>,
     ) -> PyResult<Bound<'py, PyArray2<f64>>> {
         let points = points.as_array();
+        let center = center.as_array();
         Ok(py
-            .detach(|| self.inner.transform_points(points))
+            .detach(|| {
+                if let Some(center) = center.as_slice() {
+                    self.inner.transform_points(points, center)
+                } else {
+                    self.inner.transform_points(points, &center.to_vec())
+                }
+            })
             .map_err(Error::from)?
             .into_pyarray(py))
     }
@@ -598,6 +726,29 @@ impl PyTransform {
             ))),
         })
     }
+
+    #[staticmethod]
+    fn register_translation_fft<'py>(
+        py: Python<'py>,
+        #[gen_stub(override_type(type_repr="numpy.typing.ArrayLike", imports=("numpy.typing")))]
+        fixed: PyArrayLikeDyn<f64, AllowTypeChange>,
+        #[gen_stub(override_type(type_repr="numpy.typing.ArrayLike", imports=("numpy.typing")))]
+        moving: PyArrayLikeDyn<f64, AllowTypeChange>,
+    ) -> PyResult<PyTransform> {
+        let fixed = fixed.as_array();
+        let moving = moving.as_array();
+        if fixed.shape() != moving.shape() {
+            return Err(PyErr::from(Error::ShapeMismatch(
+                fixed.shape().to_vec(),
+                moving.shape().to_vec(),
+            )));
+        }
+        py.detach(|| {
+            Ok(PyTransform {
+                inner: Transform::register_translation_fft(fixed, moving).map_err(Error::from)?,
+            })
+        })
+    }
 }
 
 #[gen_stub_pyclass]
@@ -613,22 +764,375 @@ pub(crate) struct PyTransforms {
     pub(crate) inner: Transforms,
 }
 
+#[gen_stub_pymethods]
+#[pymethods]
+impl PyView {
+    /// the bspline order used for interpolation when using transforms
+    /// None -> 0 (nearest neighbor) for int dtypes, 3 (cubic) for float dtypes
+    #[getter]
+    fn get_interpolation_order(&self) -> Option<usize> {
+        self.view.transforms.interpolation_order
+    }
+
+    /// the bspline order used for interpolation when using transforms
+    /// None -> 0 (nearest neighbor) for int dtypes, 3 (cubic) for float dtypes
+    #[setter]
+    fn set_interpolation_order(&mut self, order: Option<usize>) {
+        self.view.transforms.interpolation_order = order;
+    }
+
+    /// the bspline order used for interpolation when using transforms
+    /// None -> 0 (nearest neighbor) for int dtypes, 3 (cubic) for float dtypes
+    fn with_interpolation_order(&self, order: Option<usize>) -> Self {
+        let mut new = self.clone();
+        new.view.transforms.interpolation_order = order;
+        new
+    }
+
+    /// the value of the interpolated image outside the data
+    fn get_default_interpolation_value(&self, channel: usize) -> PyResult<PyDefaultValue> {
+        let channel_name = self.view.channel_name(channel)?.to_string();
+        Ok(PyDefaultValue {
+            inner: self
+                .view
+                .transforms
+                .default_value
+                .get(&channel_name)
+                .cloned()
+                .unwrap_or_default(),
+        })
+    }
+
+    /// the value of the interpolated image outside the data
+    fn set_default_interpolation_value(
+        &mut self,
+        channel: usize,
+        #[gen_stub(override_type(type_repr="float | transforms.DefaultValue | collections.abc.Sequence[float]", imports=("collections.abc")))]
+        value: Bound<PyAny>,
+    ) -> PyResult<()> {
+        let channel_name = self.view.channel_name(channel)?.to_string();
+        if let Ok(value) = value.extract::<PyDefaultValue>() {
+            self.view
+                .transforms
+                .default_value
+                .insert(channel_name, value.inner);
+        } else if let Ok(value) = value.extract::<f64>() {
+            self.view
+                .transforms
+                .default_value
+                .insert(channel_name, value.into());
+        } else if let Ok(value) = value.try_iter() {
+            self.view.transforms.default_value.insert(
+                channel_name,
+                value
+                    .map(|i| i?.extract::<f64>())
+                    .collect::<PyResult<Vec<_>>>()?
+                    .into(),
+            );
+        } else {
+            return Err(PyValueError::new_err(format!(
+                "cannot parse value: {:?} into DefaultValue",
+                value
+            )));
+        }
+        Ok(())
+    }
+
+    /// the value of the interpolated image outside the data
+    fn with_default_interpolation_value(
+        &self,
+        channel: usize,
+        #[gen_stub(override_type(type_repr="float | transforms.DefaultValue | collections.abc.Sequence[float]", imports=("collections.abc")))]
+        value: Bound<PyAny>,
+    ) -> PyResult<Self> {
+        let mut new = self.clone();
+        new.set_default_interpolation_value(channel, value)?;
+        Ok(new)
+    }
+
+    /// return a new view with transformations applied (channel alignment, drift correction)
+    /// backward compatibility
+    #[cfg(feature = "transforms")]
+    #[pyo3(signature = (channels = true, drift = false, file = None, bead_files = None, main_channel = None, default_transform = None))]
+    fn with_transform(
+        &self,
+        channels: bool,
+        drift: bool,
+        file: Option<PathBuf>,
+        bead_files: Option<Vec<PathBuf>>,
+        main_channel: Option<usize>,
+        default_transform: Option<Vec<f64>>,
+    ) -> PyResult<PyView> {
+        self.with_transforms(
+            channels,
+            drift,
+            file,
+            bead_files,
+            main_channel,
+            default_transform,
+        )
+    }
+
+    /// returns a view where channels and/or frames are registered with an affine transformation
+    //         channels: True/False register channels using bead_files
+    //         drift: True/False register frames to correct drift
+    //         file: load registration from file with name file, default: transform.yml in self.path.parent
+    //         bead_files: files used to register channels, default: files in self.path.parent,
+    //             with names starting with 'beads'
+    //         main channel: which channel to transform with default_transform
+    //         default_transform: parameters of the default transform: [a, b, c, d, e, f] where the affine transform matrix is
+    //             | a b e |
+    //             | c d f |
+    //             | 0 0 1 |
+    #[cfg(feature = "transforms")]
+    #[pyo3(signature = (channels = true, drift = false, file = None, bead_files = None, main_channel = None, default_transform = None))]
+    fn with_transforms(
+        &self,
+        channels: bool,
+        drift: bool,
+        file: Option<PathBuf>,
+        bead_files: Option<Vec<PathBuf>>,
+        main_channel: Option<usize>,
+        default_transform: Option<Vec<f64>>,
+    ) -> PyResult<PyView> {
+        let mut new = self.clone();
+        let bead_files = bead_files
+            .as_ref()
+            .map(|i| i.iter().map(|j| j.as_path()).collect::<Vec<_>>());
+        new.view.transforms = Transforms::calculate_transforms_from(
+            file.as_deref(),
+            if channels {
+                bead_files.as_deref()
+            } else {
+                None
+            },
+            if drift { Some(self.view.path()) } else { None },
+            main_channel,
+            default_transform,
+        )?;
+        Ok(new)
+    }
+
+    /// get the transformation
+    #[getter]
+    fn get_transforms(&self) -> PyResult<PyTransforms> {
+        Ok(PyTransforms {
+            inner: self.view.get_transforms().clone(),
+        })
+    }
+
+    /// set the transformation
+    #[setter]
+    fn set_transforms(&mut self, transform: PyTransforms) {
+        self.view.set_transforms(transform.inner);
+    }
+
+    fn load_transforms_from_yaml(&mut self, path: PathBuf) -> PyResult<()> {
+        self.view.load_transforms_from_yaml(path.as_path())?;
+        Ok(())
+    }
+
+    fn calculate_channel_transforms_2d(
+        &self,
+        py: Python,
+        main_channel: usize,
+    ) -> PyResult<HashMap<String, PyTransform>> {
+        py.detach(|| {
+            Ok(self
+                .view
+                .calculate_channel_transforms_2d(main_channel)?
+                .into_iter()
+                .map(|(n, t)| {
+                    (
+                        n,
+                        PyTransform {
+                            inner: t.into_dyn(),
+                        },
+                    )
+                })
+                .collect())
+        })
+    }
+
+    fn calculate_channel_transforms_3d(
+        &self,
+        py: Python,
+        main_channel: usize,
+    ) -> PyResult<HashMap<String, PyTransform>> {
+        py.detach(|| {
+            Ok(self
+                .view
+                .calculate_channel_transforms_3d(main_channel)?
+                .into_iter()
+                .map(|(n, t)| {
+                    (
+                        n,
+                        PyTransform {
+                            inner: t.into_dyn(),
+                        },
+                    )
+                })
+                .collect())
+        })
+    }
+
+    fn calculate_drift_transforms_2d(&self, py: Python) -> PyResult<Vec<PyTransform>> {
+        py.detach(|| {
+            Ok(self
+                .view
+                .calculate_drift_transforms_2d()?
+                .into_iter()
+                .map(|t| PyTransform {
+                    inner: t.into_dyn(),
+                })
+                .collect())
+        })
+    }
+
+    fn calculate_drift_transforms_3d(&self, py: Python) -> PyResult<Vec<PyTransform>> {
+        py.detach(|| {
+            Ok(self
+                .view
+                .calculate_drift_transforms_3d()?
+                .into_iter()
+                .map(|t| PyTransform {
+                    inner: t.into_dyn(),
+                })
+                .collect())
+        })
+    }
+}
+
 #[gen_stub_pymethods]
 #[pymethods]
 impl PyTransforms {
     #[new]
-    #[pyo3(signature = (channel = HashMap::new(), drift = Vec::new()))]
+    #[pyo3(signature = (channel = HashMap::new(), drift = Vec::new(), interpolation_order = None, default_value = None))]
     fn new(
         channel: HashMap<String, PyTransform>,
         drift: Vec<PyTransform>,
+        interpolation_order: Option<usize>,
+        default_value: Option<HashMap<String, PyDefaultValue>>,
     ) -> PyResult<PyTransforms> {
         let channel = channel.into_iter().map(|(n, t)| (n, t.inner)).collect();
         let drift = drift.into_iter().map(|t| t.inner).collect();
         Ok(PyTransforms {
-            inner: Transforms { channel, drift },
+            inner: Transforms {
+                default_value: default_value
+                    .into_iter()
+                    .flatten()
+                    .map(|(k, v)| (k, v.inner))
+                    .collect::<HashMap<_, _>>(),
+                interpolation_order,
+                channel,
+                drift,
+            },
         })
     }
 
+    #[allow(clippy::type_complexity)]
+    fn __getnewargs__(
+        &self,
+    ) -> (
+        HashMap<String, PyTransform>,
+        Vec<PyTransform>,
+        Option<usize>,
+        HashMap<String, PyDefaultValue>,
+    ) {
+        (
+            self.inner
+                .channel
+                .iter()
+                .map(|(k, v)| (k.to_string(), PyTransform { inner: v.clone() }))
+                .collect(),
+            self.inner
+                .drift
+                .iter()
+                .map(|i| PyTransform { inner: i.clone() })
+                .collect(),
+            self.inner.interpolation_order,
+            self.inner
+                .default_value
+                .iter()
+                .map(|(k, v)| (k.to_string(), PyDefaultValue { inner: v.clone() }))
+                .collect(),
+        )
+    }
+
+    /// the bspline order used for interpolation when using transforms,
+    /// None -> 0 (nearest neighbor) for int dtypes, 3 (cubic) for float dtypes
+    #[getter]
+    fn get_interpolation_order(&self) -> Option<usize> {
+        self.inner.interpolation_order
+    }
+
+    /// the bspline order used for interpolation when using transforms
+    /// None -> 0 (nearest neighbor) for int dtypes, 3 (cubic) for float dtypes
+    #[setter]
+    fn set_interpolation_order(&mut self, order: Option<usize>) {
+        self.inner.interpolation_order = order;
+    }
+
+    /// the bspline order used for interpolation when using transforms
+    /// None -> 0 (nearest neighbor) for int dtypes, 3 (cubic) for float dtypes
+    fn with_interpolation_order(&self, order: Option<usize>) -> Self {
+        let mut new = self.clone();
+        new.inner.interpolation_order = order;
+        new
+    }
+
+    /// the value of the interpolated image outside the data
+    fn get_default_interpolation_value(&self, channel_name: &str) -> PyResult<PyDefaultValue> {
+        Ok(PyDefaultValue {
+            inner: self
+                .inner
+                .default_value
+                .get(channel_name)
+                .cloned()
+                .unwrap_or_default(),
+        })
+    }
+
+    /// the value of the interpolated image outside the data
+    fn set_default_interpolation_value(
+        &mut self,
+        channel_name: String,
+        #[gen_stub(override_type(type_repr="float | DefaultValue | collections.abc.Sequence[float]", imports=("collections.abc")))]
+        value: Bound<PyAny>,
+    ) -> PyResult<()> {
+        if let Ok(value) = value.extract::<PyDefaultValue>() {
+            self.inner.default_value.insert(channel_name, value.inner);
+        } else if let Ok(value) = value.extract::<f64>() {
+            self.inner.default_value.insert(channel_name, value.into());
+        } else if let Ok(value) = value.try_iter() {
+            self.inner.default_value.insert(
+                channel_name,
+                value
+                    .map(|i| i?.extract::<f64>())
+                    .collect::<PyResult<Vec<_>>>()?
+                    .into(),
+            );
+        } else {
+            return Err(PyValueError::new_err(format!(
+                "cannot parse value: {:?} into DefaultValue",
+                value
+            )));
+        }
+        Ok(())
+    }
+
+    /// the value of the interpolated image outside the data
+    fn with_default_interpolation_value(
+        &self,
+        channel_name: String,
+        #[gen_stub(override_type(type_repr="float | DefaultValue | collections.abc.Sequence[float]", imports=("collections.abc")))]
+        value: Bound<PyAny>,
+    ) -> PyResult<Self> {
+        let mut new = self.clone();
+        new.set_default_interpolation_value(channel_name, value)?;
+        Ok(new)
+    }
+
     #[getter]
     fn get_channel_transforms(&self) -> HashMap<String, PyTransform> {
         self.inner
@@ -659,10 +1163,6 @@ impl PyTransforms {
         Ok(())
     }
 
-    fn __getnewargs__(&self) -> (HashMap<String, PyTransform>, Vec<PyTransform>) {
-        (self.get_channel_transforms(), self.get_drift_transforms())
-    }
-
     #[staticmethod]
     pub(crate) fn load(path: PathBuf) -> PyResult<PyTransforms> {
         Ok(PyTransforms {

ndbioimage/src/readers.rs

@@ -2,7 +2,7 @@ use crate::axes::{Axis, Shape};
 use crate::error::Error;
 use crate::view::View;
 use ndarray::{Array, Dimension, Ix2, Ix5, s};
-use num::{FromPrimitive, Zero};
+use num::{FromPrimitive, One, Zero};
 use ome_metadata::Ome;
 use regex::Regex;
 use serde::{Deserialize, Serialize};
@@ -128,7 +128,7 @@ pub enum PixelType {
     I128,
     /// unsigned 128-bit integer
     U128,
-    /// 128-bit float (emulated)
+    /// 128-bit float (not supported)
     F128,
 }
 
@@ -148,6 +148,8 @@ impl PixelType {
 /// array data with a specific pixel type
 #[derive(Clone, Debug)]
 pub enum ArrayT<D: Dimension> {
+    /// boolean array
+    Bool(Array<bool, D>),
     /// signed 8-bit integer array
     I8(Array<i8, D>),
     /// unsigned 8-bit integer array
@@ -168,24 +170,28 @@ pub enum ArrayT<D: Dimension> {
     I64(Array<i64, D>),
     /// unsigned 64-bit integer array
     U64(Array<u64, D>),
-    /// signed 128-bit integer array
+    /// signed 128-bit integer array (not supported by numpy)
     I128(Array<i128, D>),
-    /// unsigned 128-bit integer array
+    /// unsigned 128-bit integer array (not supported by numpy)
     U128(Array<u128, D>),
-    /// 128-bit float array (emulated as f64)
-    F128(Array<f64, D>),
 }
 
 /// type alias for a single frame (2D image)
 pub type Frame = ArrayT<Ix2>;
 
 /// trait for reading image files
-pub trait Reader: Clone + Sized + Debug + Send + Hash + Into<DynReader> {
+pub trait Reader: Clone + Sized + Debug + Sync + Send + Hash + Into<DynReader> {
     /// create a new reader for the given path, series, and position
     fn new<P>(path: P, series: usize, position: usize) -> Result<Self, Error>
     where
         P: AsRef<Path>;
 
+    /// a simple fast test to see if the reader can open a path without parsing the whole file
+    /// not a guarantee of success, but rather an early failure
+    fn can_open<P>(path: P) -> Result<bool, Error>
+    where
+        P: AsRef<Path>;
+
     /// the name of the reader type
     fn reader_name(&self) -> &'static str {
         type_name::<Self>()
@@ -197,9 +203,9 @@ pub trait Reader: Clone + Sized + Debug + Send + Hash + Into<DynReader> {
 
     /// get a sliceable view on the image file
     fn view(&self) -> View<Ix5, Self> {
-        let shape = self.shape();
+        let shape = self.raw_shape();
         let slice = s![0..shape.c, 0..shape.z, 0..shape.t, 0..shape.y, 0..shape.x,];
-        View::new(
+        View::new_with_slice_and_axes(
             self.clone(),
             slice.as_ref().to_vec(),
             vec![Axis::C, Axis::Z, Axis::T, Axis::Y, Axis::X],
@@ -216,7 +222,7 @@ pub trait Reader: Clone + Sized + Debug + Send + Hash + Into<DynReader> {
     /// the position index
     fn position(&self) -> usize;
     /// the shape of the image data
-    fn shape(&self) -> &Shape;
+    fn raw_shape(&self) -> &Shape;
     /// the pixel type of the image data
     fn pixel_type(&self) -> &PixelType;
     /// get all available positions for a given series
@@ -314,6 +320,7 @@ macro_rules! try_into_array_body {
     ($self:expr, $map:ident) => {{
         let mut err = Ok(());
         let arr = match $self {
+            ArrayT::Bool(v) => v.$map(|x| if x { T::one() } else { T::zero() }),
             ArrayT::I8(v) => v.$map(|x| {
                 T::from_i8(x).unwrap_or_else(|| {
                     err = Err(Error::Cast(x.to_string(), type_name::<T>().to_string()));
@@ -356,7 +363,7 @@ macro_rules! try_into_array_body {
                     T::zero()
                 })
             }),
-            ArrayT::F64(v) | ArrayT::F128(v) => v.$map(|x| {
+            ArrayT::F64(v) => v.$map(|x| {
                 T::from_f64(x).unwrap_or_else(|| {
                     err = Err(Error::Cast(x.to_string(), type_name::<T>().to_string()));
                     T::zero()
@@ -397,7 +404,7 @@ macro_rules! try_into_array_body {
 impl<D, T> TryInto<Array<T, D>> for ArrayT<D>
 where
     D: Dimension,
-    T: FromPrimitive + Zero + 'static,
+    T: FromPrimitive + Zero + One + 'static,
 {
     type Error = Error;
 
@@ -409,7 +416,7 @@ where
 impl<D, T> TryInto<Array<T, D>> for &ArrayT<D>
 where
     D: Dimension,
-    T: Clone + FromPrimitive + Zero + 'static,
+    T: Clone + FromPrimitive + Zero + One + 'static,
 {
     type Error = Error;
 
@@ -475,6 +482,34 @@ impl Reader for DynReader {
         ))
     }
 
+    fn can_open<P>(path: P) -> Result<bool, Error>
+    where
+        P: AsRef<Path>,
+    {
+        let mut can_open = false;
+        #[cfg(feature = "tiff")]
+        {
+            can_open |= tiff::TiffReader::can_open(&path)?;
+        }
+        #[cfg(feature = "tiffseq")]
+        {
+            can_open |= tiffseq::TiffSeqReader::can_open(&path)?;
+        }
+        #[cfg(feature = "czi")]
+        {
+            can_open |= czi::CziReader::can_open(&path)?;
+        }
+        #[cfg(feature = "bioformats_rust")]
+        {
+            can_open |= bioformats_rust::BioFormatsRustReader::can_open(&path)?;
+        }
+        #[cfg(feature = "bioformats_java")]
+        {
+            can_open |= bioformats_java::BioFormatsJavaReader::can_open(&path)?;
+        }
+        Ok(can_open)
+    }
+
     fn reader_name(&self) -> &'static str {
         match self {
             #[cfg(feature = "tiff")]
@@ -577,18 +612,18 @@ impl Reader for DynReader {
         }
     }
 
-    fn shape(&self) -> &Shape {
+    fn raw_shape(&self) -> &Shape {
         match self {
             #[cfg(feature = "tiff")]
-            DynReader::Tiff(r) => r.shape(),
+            DynReader::Tiff(r) => r.raw_shape(),
             #[cfg(feature = "tiffseq")]
-            DynReader::TiffSeq(r) => r.shape(),
+            DynReader::TiffSeq(r) => r.raw_shape(),
             #[cfg(feature = "czi")]
-            DynReader::Czi(r) => r.shape(),
+            DynReader::Czi(r) => r.raw_shape(),
             #[cfg(feature = "bioformats_rust")]
-            DynReader::BioFormatsRust(r) => r.shape(),
+            DynReader::BioFormatsRust(r) => r.raw_shape(),
             #[cfg(feature = "bioformats_java")]
-            DynReader::BioFormatsJava(r) => r.shape(),
+            DynReader::BioFormatsJava(r) => r.raw_shape(),
             #[expect(unreachable_patterns)]
             _ => unreachable!(),
         }

ndbioimage/src/readers/bioformats_java.rs

@@ -14,6 +14,7 @@ use std::hash::{Hash, Hasher};
 use std::ops::Deref;
 use std::rc::Rc;
 use std::sync::Mutex;
+use std::sync::atomic::{AtomicU32, Ordering};
 use thread_local::ThreadLocal;
 
 include!(concat!(env!("OUT_DIR"), "/constants.rs"));
@@ -24,10 +25,33 @@ thread_local! {
 
 static DOWNLOAD_LOCK: Mutex<()> = Mutex::new(());
 static JVM_BUILT: Mutex<bool> = Mutex::new(false);
+/// the process that called [JvmBuilder::build], 0 if no jvm has been built yet
+static JVM_PID: AtomicU32 = AtomicU32::new(0);
+
+/// Returns `Ok(())` if the jvm in this process can still be used, and an error if it was inherited
+/// across a fork.
+///
+/// A fork copies the address space, so the child gets the parent's jvm wholesale: its threads do
+/// not exist and any lock another thread held at the moment of the fork stays held forever. There
+/// is no way to repair that, and no way to start a replacement either, because HotSpot's
+/// "a jvm already exists" flag lives in libjvm's private data and is copied along with everything
+/// else -- verified by `JNI_GetCreatedJavaVMs`, which still reports 1 in the child, so a second
+/// `JNI_CreateJavaVM` (which is what j4rs' `Jvm::attach_thread` calls) takes the already-exists
+/// path. An `exec` replaces the address space, and that does reset the count to 0.
+pub fn check_not_forked() -> Result<(), Error> {
+    let built_in = JVM_PID.load(Ordering::Relaxed);
+    if built_in != 0 && built_in != std::process::id() {
+        return Err(Error::ForkedJvm(built_in, std::process::id()));
+    }
+    Ok(())
+}
 
 /// Ensure 1 jvm per thread
-fn jvm() -> Rc<Jvm> {
-    JVM.with(|cell| {
+fn jvm() -> Result<Rc<Jvm>, Error> {
+    // before touching JVM: in a forked child that thread's cell may already hold an Rc<Jvm> whose
+    // vm is gone, and the guard in check_not_forked is the only thing that can catch that.
+    check_not_forked()?;
+    Ok(JVM.with(|cell| {
         cell.get_or_init(move || {
             #[cfg(feature = "python")]
             let path = crate::py::ndbioimage_file();
@@ -113,6 +137,7 @@ fn jvm() -> Rc<Jvm> {
                         .java_opt(JavaOpt::new("--sun-misc-unsafe-memory-access=allow"))
                         .build()
                         .expect("Failed to build JVM");
+                    JVM_PID.store(std::process::id(), Ordering::Relaxed);
                     if let Ok(e) = InvocationArg::try_from("ERROR") {
                         let _ = j.invoke_static(
                             "loci.common.DebugTools",
@@ -125,9 +150,103 @@ fn jvm() -> Rc<Jvm> {
             }
         })
             .clone()
-    })
+    }))
 }
 
+/// Returns `true` if the jvm was already created on the current thread.
+///
+/// The jvm is created lazily, on the first bio-formats read, and HotSpot installs its own signal
+/// handlers while starting up. So this tells whether those handlers are in place yet.
+///
+/// This is process-wide, not thread-local: a forked child inherited the parent's jvm but cannot
+/// use it, so reporting `true` there would only invite the deadlock that [check_not_forked] turns
+/// into an error.
+pub fn jvm_is_running() -> bool {
+    JVM_PID.load(Ordering::Relaxed) == std::process::id() && JVM.with(|cell| cell.get().is_some())
+}
+
+/// Creates the jvm on the current thread if it does not exist yet (no-op if it already does),
+/// and initializes bio-formats on it.
+///
+/// This exists because of how HotSpot probes the cpu features during VM init: on an avx-512
+/// capable cpu it deliberately reads through a NULL pointer (`VM_Version::_cpuinfo_segv_addr`, set
+/// to the address of the generated cpu-feature stub) and resumes from its own SIGSEGV handler after
+/// recognizing that address in `PosixSignals::pd_hotspot_signal_handler`.
+///
+/// Any SIGSEGV handler that was installed *before* the VM started is not that handler, and
+/// HotSpot chains to it instead. Then the deliberate fault is not recognized, the chained handler
+/// re-raises SIGSEGV, and the re-entrant pass no longer has a recognizable pc and ends in
+/// `VMError::report_and_die` -> `os::die()` -> `abort()`.
+///
+/// In a python process this bites when `ray.init()` has already installed CPython's faulthandler
+/// (ray/_private/worker.py), so call this before that.
+///
+/// Note that creating the VM is *not* enough: it has to be initialized far enough that the
+/// deliberate fault has been taken before the foreign handler shows up. Measured on an AMD EPYC
+/// 9124 (avx-512, XCR0 = 0x2e7) with ray: merely creating the VM leaves the crash rate unchanged,
+/// and so does loading the reader classes without opening a file. What does work is initializing a
+/// reader on an actual (tiny) image, so that is what this does.
+pub fn warm_jvm() -> Result<(), Error> {
+    let jvm = jvm()?;
+    let no_args: &[InvocationArg] = &[];
+    // jvm() only *creates* the VM; run one static method so that this only returns Ok once the VM
+    // can actually run java bytecode.
+    let version = jvm.invoke_static(
+        "java.lang.System",
+        "getProperty",
+        &[InvocationArg::try_from("java.version")?],
+    )?;
+    let _: String = jvm.to_rust(version)?;
+
+    // loads (and runs the static initializers of) every reader class in the registry
+    let _ = jvm.invoke_static(
+        "loci.formats.ImageReader",
+        "getDefaultReaderClasses",
+        no_args,
+    )?;
+
+    // then open a 1x1 image so that the reader and its OMEXML metadata store get initialized.
+    // ImageReader.setId needs a real file, so write the tiff out and clean it up again.
+    let path = std::env::temp_dir().join(format!("ndbioimage-warm-jvm-{}.tif", std::process::id()));
+    std::fs::write(&path, WARM_TIFF)?;
+    let result = (|| -> Result<(), Error> {
+        let reader = jvm.create_instance("loci.formats.ImageReader", no_args)?;
+        jvm.invoke(
+            &reader,
+            "setId",
+            &[InvocationArg::try_from(path.to_string_lossy().as_ref())?],
+        )?;
+        let _: i32 = jvm.to_rust(jvm.invoke(&reader, "getSizeY", no_args)?)?;
+        let _ = jvm.invoke(&reader, "getMetadataStore", no_args)?;
+        Ok(())
+    })();
+    let _ = std::fs::remove_file(&path);
+    result
+}
+
+/// A 1x1 8-bit grayscale tiff (written by tifffile), used by [warm_jvm] to give bio-formats
+/// something it can open.
+const WARM_TIFF: &[u8] = &[
+    0x49, 0x49, 0x2a, 0x00, 0x08, 0x00, 0x00, 0x00, 0x0e, 0x00, 0x00, 0x01, 0x04, 0x00, 0x01, 0x00,
+    0x00, 0x00, 0x01, 0x00, 0x00, 0x00, 0x01, 0x01, 0x04, 0x00, 0x01, 0x00, 0x00, 0x00, 0x01, 0x00,
+    0x00, 0x00, 0x02, 0x01, 0x03, 0x00, 0x01, 0x00, 0x00, 0x00, 0x08, 0x00, 0x00, 0x00, 0x03, 0x01,
+    0x03, 0x00, 0x01, 0x00, 0x00, 0x00, 0x01, 0x00, 0x00, 0x00, 0x06, 0x01, 0x03, 0x00, 0x01, 0x00,
+    0x00, 0x00, 0x01, 0x00, 0x00, 0x00, 0x0e, 0x01, 0x02, 0x00, 0x12, 0x00, 0x00, 0x00, 0xb6, 0x00,
+    0x00, 0x00, 0x11, 0x01, 0x04, 0x00, 0x01, 0x00, 0x00, 0x00, 0x00, 0x01, 0x00, 0x00, 0x15, 0x01,
+    0x03, 0x00, 0x01, 0x00, 0x00, 0x00, 0x01, 0x00, 0x00, 0x00, 0x16, 0x01, 0x04, 0x00, 0x01, 0x00,
+    0x00, 0x00, 0x01, 0x00, 0x00, 0x00, 0x17, 0x01, 0x04, 0x00, 0x01, 0x00, 0x00, 0x00, 0x01, 0x00,
+    0x00, 0x00, 0x1a, 0x01, 0x05, 0x00, 0x01, 0x00, 0x00, 0x00, 0xd8, 0x00, 0x00, 0x00, 0x1b, 0x01,
+    0x05, 0x00, 0x01, 0x00, 0x00, 0x00, 0xe0, 0x00, 0x00, 0x00, 0x28, 0x01, 0x03, 0x00, 0x01, 0x00,
+    0x00, 0x00, 0x01, 0x00, 0x00, 0x00, 0x31, 0x01, 0x02, 0x00, 0x0c, 0x00, 0x00, 0x00, 0xe8, 0x00,
+    0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x7b, 0x22, 0x73, 0x68, 0x61, 0x70, 0x65, 0x22, 0x3a,
+    0x20, 0x5b, 0x31, 0x2c, 0x20, 0x31, 0x5d, 0x7d, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
+    0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x01, 0x00,
+    0x00, 0x00, 0x01, 0x00, 0x00, 0x00, 0x01, 0x00, 0x00, 0x00, 0x01, 0x00, 0x00, 0x00, 0x01, 0x00,
+    0x00, 0x00, 0x00, 0x74, 0x69, 0x66, 0x66, 0x66, 0x69, 0x6c, 0x65, 0x2e, 0x70, 0x79, 0x00, 0x00,
+    0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
+    0x00, 0x00,
+];
+
 macro_rules! method_return {
     ($R:ty$(|c)?) => { Result<$R, Error> };
     () => { Result<(), Error> };
@@ -147,14 +266,14 @@ macro_rules! method {
         #[allow(dead_code)]
         pub(crate) fn $name(&self, $($($n: $t),*)?) -> method_return!($($tt)?) {
             let args: Vec<InvocationArg> = vec![$($( method_arg!($n:$t$(|$p)?) ),*)?];
-            let _result = jvm().invoke(&self.0, $method, &args)?;
+            let _result = jvm()?.invoke(&self.0, $method, &args)?;
 
             macro_rules! method_result {
                 ($R:ty|c) => {
-                    Ok(jvm().to_rust(_result)?)
+                    Ok(jvm()?.to_rust(_result)?)
                 };
                 ($R:ty|d) => {
-                    Ok(jvm().to_rust_deserialized(_result)?)
+                    Ok(jvm()?.to_rust_deserialized(_result)?)
                 };
                 ($R:ty) => {
                     Ok(_result)
@@ -187,7 +306,7 @@ pub struct DebugTools;
 impl DebugTools {
     /// set debug root level: ERROR, DEBUG, TRACE, INFO, OFF
     pub fn set_root_level(level: &str) -> Result<(), Error> {
-        jvm().invoke_static(
+        jvm()?.invoke_static(
             "loci.common.DebugTools",
             "setRootLevel",
             &[InvocationArg::try_from(level)?],
@@ -201,7 +320,7 @@ pub(crate) struct ChannelSeparator(Instance);
 
 impl ChannelSeparator {
     pub(crate) fn new(image_reader: &ImageReader) -> Result<Self, Error> {
-        let jvm = jvm();
+        let jvm = jvm()?;
         let channel_separator = jvm.create_instance(
             "loci.formats.ChannelSeparator",
             &[InvocationArg::from(jvm.clone_instance(&image_reader.0)?)],
@@ -228,7 +347,7 @@ impl Drop for ImageReader {
 
 impl ImageReader {
     pub(crate) fn new() -> Result<Self, Error> {
-        let reader = jvm().create_instance("loci.formats.ImageReader", InvocationArg::empty())?;
+        let reader = jvm()?.create_instance("loci.formats.ImageReader", InvocationArg::empty())?;
         Ok(ImageReader(reader))
     }
 
@@ -238,7 +357,7 @@ impl ImageReader {
 
     pub(crate) fn ome_xml(&self) -> Result<String, Error> {
         let mds = self.get_metadata_store()?;
-        Ok(jvm()
+        Ok(jvm()?
             .chain(&mds)?
             .cast("loci.formats.ome.OMEPyramidStore")?
             .invoke("dumpXML", InvocationArg::empty())?
@@ -276,7 +395,7 @@ pub(crate) struct MetadataTools(Instance);
 impl MetadataTools {
     pub(crate) fn new() -> Result<Self, Error> {
         let meta_data_tools =
-            jvm().create_instance("loci.formats.MetadataTools", InvocationArg::empty())?;
+            jvm()?.create_instance("loci.formats.MetadataTools", InvocationArg::empty())?;
         Ok(MetadataTools(meta_data_tools))
     }
 
@@ -402,7 +521,7 @@ impl BioFormatsJavaReader {
             PixelType::U64 => 8,
             PixelType::I128 => 16,
             PixelType::U128 => 16,
-            PixelType::F128 => 8,
+            PixelType::F128 => 16,
         };
         Ok(bytes
             .chunks(chunk_size)
@@ -468,7 +587,7 @@ impl BioFormatsJavaReader {
             (PixelType::U64, true) => get_frame!(u64, <8),
             (PixelType::I128, true) => get_frame!(i128, <16),
             (PixelType::U128, true) => get_frame!(u128, <16),
-            (PixelType::F128, true) => get_frame!(f64, <8),
+            (PixelType::F128, true) => Err(Error::F128NotSupported),
             (PixelType::Bool, false) => get_frame!(bool, >1),
             (PixelType::I8, false) => get_frame!(i8, >1),
             (PixelType::U8, false) => get_frame!(u8, >1),
@@ -482,7 +601,7 @@ impl BioFormatsJavaReader {
             (PixelType::U64, false) => get_frame!(u64, >8),
             (PixelType::I128, false) => get_frame!(i128, >16),
             (PixelType::U128, false) => get_frame!(u128, >16),
-            (PixelType::F128, false) => get_frame!(f64, >8),
+            (PixelType::F128, false) => Err(Error::F128NotSupported),
         }
     }
 }
@@ -517,6 +636,10 @@ impl Reader for BioFormatsJavaReader {
     where
         P: AsRef<Path>,
     {
+        // checked here as well as in jvm() so that this entry point fails with the error rather
+        // than reaching get_size_x() below, which goes through deref() -> get_reader().unwrap()
+        // and would turn it into a panic.
+        check_not_forked()?;
         let mut path = path.as_ref().to_path_buf();
         if path.is_dir() {
             let orig = path.clone();
@@ -543,6 +666,13 @@ impl Reader for BioFormatsJavaReader {
         Ok(new)
     }
 
+    fn can_open<P>(_path: P) -> Result<bool, Error>
+    where
+        P: AsRef<Path>,
+    {
+        Ok(true)
+    }
+
     fn metadata(&self) -> Result<Ome, Error> {
         self.get_ome()
     }
@@ -575,7 +705,7 @@ impl Reader for BioFormatsJavaReader {
         0
     }
 
-    fn shape(&self) -> &Shape {
+    fn raw_shape(&self) -> &Shape {
         &self.shape
     }
 
@@ -648,6 +778,7 @@ mod tests {
         metadata_g: "tiffseq/4-Pos_001_002/img_000000000_Cy3-Cy3_filter_000.tif",
         metadata_h: "tiffseq/20-Pos_005_005/img_000000000_Cy3-Cy3_filter_000.tif",
         metadata_i: "tiffseq/YTL1841B2-2-1_1hr_DMSO_galinduction_1",
+        metadata_j: "tiff/R1_NRF1_Gtf3a500nM_REAL_20260731_25pct_500ms_001_MMStack_PosR1_NRF1_Gtf3a500nM_REAL_20260731_25pct_500ms_001.ome.tif",
     }
 
     #[test]

ndbioimage/src/readers/bioformats_rust.rs

@@ -140,7 +140,7 @@ impl BioFormatsRustReader {
             PixelType::U64 => 8,
             PixelType::I128 => 16,
             PixelType::U128 => 16,
-            PixelType::F128 => 8,
+            PixelType::F128 => 16,
         };
         Ok(bytes
             .chunks(chunk_size)
@@ -206,7 +206,7 @@ impl BioFormatsRustReader {
             (PixelType::U64, true) => get_frame!(u64, <8),
             (PixelType::I128, true) => get_frame!(i128, <16),
             (PixelType::U128, true) => get_frame!(u128, <16),
-            (PixelType::F128, true) => get_frame!(f64, <8),
+            (PixelType::F128, true) => Err(Error::F128NotSupported),
             (PixelType::Bool, false) => get_frame!(bool, >1),
             (PixelType::I8, false) => get_frame!(i8, >1),
             (PixelType::U8, false) => get_frame!(u8, >1),
@@ -220,7 +220,7 @@ impl BioFormatsRustReader {
             (PixelType::U64, false) => get_frame!(u64, >8),
             (PixelType::I128, false) => get_frame!(i128, >16),
             (PixelType::U128, false) => get_frame!(u128, >16),
-            (PixelType::F128, false) => get_frame!(f64, >8),
+            (PixelType::F128, false) => Err(Error::F128NotSupported),
         }
     }
 }
@@ -274,6 +274,13 @@ impl Reader for BioFormatsRustReader {
         Ok(new)
     }
 
+    fn can_open<P>(path: P) -> Result<bool, Error>
+    where
+        P: AsRef<Path>,
+    {
+        Ok(path.as_ref().extension() != Some(std::ffi::OsStr::new("czi")))
+    }
+
     fn metadata(&self) -> Result<Ome, Error> {
         let reader = self.get_reader()?;
         let metadata = reader.metadata();
@@ -309,7 +316,7 @@ impl Reader for BioFormatsRustReader {
         0
     }
 
-    fn shape(&self) -> &Shape {
+    fn raw_shape(&self) -> &Shape {
         &self.shape
     }
 

ndbioimage/src/readers/czi.rs

@@ -122,14 +122,6 @@ impl From<Option<String>> for Version {
     }
 }
 
-/// errors specific to czi file reading
-#[derive(Debug, thiserror::Error)]
-pub enum CziError {
-    /// czi file has no valid blocks
-    #[error("czi file has no valid blocks")]
-    NoValidBlocks,
-}
-
 impl Reader for CziReader {
     fn new<P>(path: P, series: usize, position: usize) -> Result<Self, Error>
     where
@@ -195,7 +187,7 @@ impl Reader for CziReader {
         let mut max_y = i32::MIN;
         for &i in Ok(block_map.get(&(0, 0, 0)))
             .transpose()
-            .unwrap_or_else(|| block_map.values().next().ok_or(CziError::NoValidBlocks))?
+            .unwrap_or_else(|| block_map.values().next().ok_or(Error::CziNoValidBlocks))?
         {
             if let Ok(info) = reader.try_get_sub_block_info_for_index(i) {
                 let rect = info.get_logical_rect();
@@ -219,10 +211,17 @@ impl Reader for CziReader {
         new.shape.x = (max_y - min_y) as usize;
         new.shape.y = (max_x - min_x) as usize;
         new.block_map = block_map;
-        new.pixel_type = pixel_type.ok_or(CziError::NoValidBlocks)?;
+        new.pixel_type = pixel_type.ok_or(Error::CziNoValidBlocks)?;
         Ok(new)
     }
 
+    fn can_open<P>(path: P) -> Result<bool, Error>
+    where
+        P: AsRef<Path>,
+    {
+        Ok(path.as_ref().extension() == Some(std::ffi::OsStr::new("czi")))
+    }
+
     fn metadata(&self) -> Result<Ome, Error> {
         let mut ome = Ome::default();
         let xml = xmltree::Element::parse(self.metadata_xml()?.as_bytes())?;
@@ -1133,7 +1132,7 @@ impl Reader for CziReader {
             PixelType::U64 => get_frame!(u64, 8),
             PixelType::I128 => get_frame!(i128, 16),
             PixelType::U128 => get_frame!(u128, 16),
-            PixelType::F128 => get_frame!(f64, 8),
+            PixelType::F128 => Err(Error::F128NotSupported),
         }
     }
 
@@ -1149,7 +1148,7 @@ impl Reader for CziReader {
         self.position
     }
 
-    fn shape(&self) -> &Shape {
+    fn raw_shape(&self) -> &Shape {
         &self.shape
     }
 

ndbioimage/src/readers/tiff.rs

@@ -0,0 +1,585 @@
+use crate::error::Error;
+use crate::readers::{ArrayT, DynReader, Frame, PixelType, Reader, Shape};
+use ndarray::Array2;
+use ome_metadata::{Ome, ome};
+use serde::{Deserialize, Serialize};
+use std::cell::{RefCell, RefMut};
+use std::collections::{HashMap, HashSet};
+use std::fs::File;
+use std::hash::{Hash, Hasher};
+use std::io::{Read, Seek};
+use std::path::{Path, PathBuf};
+use thread_local::ThreadLocal;
+use tiff::decoder::ifd::Value;
+use tiff::decoder::{Decoder, DecodingSampleBuffer};
+use tiff::tags::Tag;
+
+/// reader for single tiff image files
+#[derive(Debug, Serialize, Deserialize)]
+pub struct TiffReader {
+    #[serde(skip)]
+    reader: ThreadLocal<RefCell<Decoder<File>>>,
+    path: PathBuf,
+    series: usize,
+    position: usize,
+    shape: Shape,
+    pixel_type: PixelType,
+    n_samples: usize,
+    p_ndim: u8,
+    planar: bool,
+}
+
+impl From<TiffReader> for DynReader {
+    fn from(value: TiffReader) -> Self {
+        DynReader::Tiff(value)
+    }
+}
+
+impl Hash for TiffReader {
+    fn hash<H: Hasher>(&self, state: &mut H) {
+        self.path.hash(state);
+        self.series.hash(state);
+        self.position.hash(state);
+    }
+}
+
+impl Clone for TiffReader {
+    fn clone(&self) -> Self {
+        TiffReader {
+            reader: ThreadLocal::default(),
+            path: self.path.clone(),
+            series: self.series,
+            position: self.position,
+            shape: self.shape.clone(),
+            pixel_type: self.pixel_type,
+            n_samples: self.n_samples,
+            p_ndim: self.p_ndim,
+            planar: self.planar,
+        }
+    }
+}
+
+impl PartialEq for TiffReader {
+    fn eq(&self, other: &Self) -> bool {
+        self.path == other.path && self.series == other.series && self.position == other.position
+    }
+}
+
+impl Eq for TiffReader {}
+
+impl TiffReader {
+    fn get_reader(&self) -> Result<RefMut<'_, Decoder<File>>, Error> {
+        self.reader
+            .get_or_try(|| {
+                let mut decoder = Decoder::open(File::open(&self.path)?)?;
+                decoder.next_image()?;
+                Ok(RefCell::new(decoder))
+            })
+            .map(|i| i.borrow_mut())
+    }
+}
+
+fn pixel_type_from_bits_sample(bits: u16, fmt: u16) -> Result<PixelType, Error> {
+    match (fmt, bits) {
+        (1, 8) => Ok(PixelType::U8),
+        (1, 16) => Ok(PixelType::U16),
+        (1, 32) => Ok(PixelType::U32),
+        (2, 8) => Ok(PixelType::I8),
+        (2, 16) => Ok(PixelType::I16),
+        (2, 32) => Ok(PixelType::I32),
+        (3, 32) => Ok(PixelType::F32),
+        (3, 64) => Ok(PixelType::F64),
+        _ => Err(Error::Parse(format!(
+            "unsupported pixel type: SampleFormat={}, BitsPerSample={}",
+            fmt, bits
+        ))),
+    }
+}
+
+fn parse_imagej_metadata(s: &str) -> HashMap<String, String> {
+    let mut map = HashMap::new();
+    for line in s.lines() {
+        let line = line.trim();
+        if line.is_empty() || line.starts_with('#') {
+            continue;
+        }
+        if let Some(pos) = line.find('=') {
+            let key = line[..pos].trim().to_string();
+            let val = line[pos + 1..].trim().to_string();
+            map.insert(key, val);
+        }
+    }
+    map
+}
+
+impl Reader for TiffReader {
+    fn new<P>(path: P, series: usize, position: usize) -> Result<Self, Error>
+    where
+        P: AsRef<Path>,
+    {
+        let mut new = Self {
+            reader: ThreadLocal::default(),
+            path: path.as_ref().to_path_buf(),
+            series,
+            position,
+            shape: Shape::default(),
+            pixel_type: PixelType::U16,
+            n_samples: 0,
+            p_ndim: 0,
+            planar: true,
+        };
+
+        let mut reader = new.get_reader()?;
+        let width = reader.current_ifd().get_tag(Tag::ImageWidth)?.into_u32()? as usize;
+        let height = reader.current_ifd().get_tag(Tag::ImageLength)?.into_u32()? as usize;
+
+        let n_samples = reader
+            .current_ifd()
+            .get_tag(Tag::SamplesPerPixel)
+            .ok()
+            .and_then(|t| t.into_u16().ok())
+            .unwrap_or(1) as usize;
+
+        let bits = tag_first_u16(&mut reader, Tag::BitsPerSample).unwrap_or(8);
+        let sfmt = tag_first_u16(&mut reader, Tag::SampleFormat).unwrap_or(1);
+
+        // indexed color pages (PhotometricInterpretation::RGBPalette) are read as plain
+        // intensities: the pixel values are indices into the ColorMap palette, which is display
+        // information that ndbioimage does not use.
+        let pixtype = pixel_type_from_bits_sample(bits, sfmt)?;
+
+        let planar = reader
+            .current_ifd()
+            .get_tag(Tag::PlanarConfiguration)
+            .ok()
+            .and_then(|t| t.into_u16().ok())
+            .unwrap_or(1)
+            == 2;
+
+        let image_description = reader.current_ifd().get_tag(Tag::ImageDescription);
+        let p_ndim = if n_samples > 1 { 3u8 } else { 2u8 };
+        let (size_c, size_t, size_z) = if let Ok(Value::Ascii(xml)) = image_description
+            && xml.starts_with("<?xml")
+        {
+            let ome = Ome::from_xml(xml)?;
+            (
+                ome.image[0].pixels.size_c as usize,
+                ome.image[0].pixels.size_t as usize,
+                ome.image[0].pixels.size_z as usize,
+            )
+        } else {
+            let imagej_bytes = reader
+                .current_ifd()
+                .get_tag(Tag::Unknown(50839))
+                .or_else(|_| reader.current_ifd().get_tag(Tag::ImageDescription))?;
+            let imagej_bytes = match imagej_bytes {
+                Value::Ascii(s) => s.into_bytes(),
+                other => other.into_u8_vec()?,
+            };
+            let metadata_str = String::from_utf8_lossy(&imagej_bytes);
+            let metadata_map = parse_imagej_metadata(&metadata_str);
+
+            let size_c = if p_ndim == 3 {
+                n_samples
+            } else {
+                metadata_map
+                    .get("channels")
+                    .and_then(|v| v.parse::<usize>().ok())
+                    .unwrap_or(1)
+            };
+            let size_z = metadata_map
+                .get("slices")
+                .and_then(|v| v.parse::<usize>().ok())
+                .unwrap_or(1);
+            let size_t = metadata_map
+                .get("frames")
+                .and_then(|v| v.parse::<usize>().ok())
+                .unwrap_or(1);
+
+            let mult = if p_ndim == 3 { n_samples } else { 1 };
+            let expected = size_c * size_z * size_t / mult;
+
+            let mut n_ifds = 1;
+            while reader.more_images() {
+                reader.next_directory()?;
+                n_ifds += 1;
+            }
+
+            let size_z = if n_ifds != expected && n_ifds > 0 {
+                n_ifds / (size_c * size_t / mult).max(1)
+            } else {
+                size_z
+            };
+            (size_c, size_t, size_z)
+        };
+        drop(reader);
+
+        new.shape = Shape {
+            c: size_c,
+            z: size_z,
+            t: size_t,
+            y: height,
+            x: width,
+            ..Default::default()
+        };
+        new.pixel_type = pixtype;
+        new.n_samples = n_samples;
+        new.p_ndim = p_ndim;
+        new.planar = planar;
+
+        // test reading first frame, if error another reader (bioformats) can take over
+        new.get_frame(0, 0, 0)?;
+        Ok(new)
+    }
+
+    fn can_open<P>(path: P) -> Result<bool, Error>
+    where
+        P: AsRef<Path>,
+    {
+        let path = path.as_ref();
+        let mut file = File::open(path)?;
+        let mut sig = [0; 2];
+        file.read_exact(&mut sig)?;
+        if sig == [73, 73] {
+            file.read_exact(&mut sig)?;
+            Ok((sig[0] == 42) || (sig[0] == 43))
+        } else if sig == [77, 77] {
+            Ok((sig[1] == 42) || (sig[1] == 43))
+        } else {
+            Ok(false)
+        }
+    }
+
+    fn metadata(&self) -> Result<Ome, Error> {
+        let mut reader = self.get_reader()?;
+        let image_description = reader.current_ifd().get_tag(Tag::ImageDescription);
+        if let Ok(Value::Ascii(xml)) = image_description
+            && xml.starts_with("<?xml")
+        {
+            Ok(Ome::from_xml(xml)?)
+        } else {
+            let imagej_bytes = reader
+                .current_ifd()
+                .get_tag(Tag::Unknown(50839))
+                .or_else(|_| reader.current_ifd().get_tag(Tag::ImageDescription))?;
+            let imagej_bytes = match imagej_bytes {
+                Value::Ascii(s) => s.into_bytes(),
+                other => other.into_u8_vec()?,
+            };
+            let metadata_str = String::from_utf8_lossy(&imagej_bytes);
+            let metadata_map = parse_imagej_metadata(&metadata_str);
+
+            let interval_t = metadata_map
+                .get("interval")
+                .and_then(|v| v.parse::<f64>().ok());
+            let ru = reader
+                .current_ifd()
+                .get_tag(Tag::ResolutionUnit)
+                .ok()
+                .and_then(|i| i.into_u16().ok())
+                .unwrap_or(2);
+            let pixel_size = reader
+                .current_ifd()
+                .get_tag(Tag::XResolution)
+                .ok()
+                .and_then(|i| i.into_u64_vec().ok())
+                .and_then(|i| {
+                    i.first()
+                        .and_then(|&n| i.get(1).map(|&d| (n as f64) / (d as f64)))
+                })
+                .and_then(|i| {
+                    if i <= 0.0 {
+                        None
+                    } else if ru == 3 {
+                        Some(1e4 / i)
+                    } else if ru == 2 {
+                        Some(25400.0 / i)
+                    } else {
+                        Some(1.0 / i)
+                    }
+                });
+            let interval_z = metadata_map
+                .get("spacing")
+                .and_then(|v| v.parse::<f64>().ok());
+
+            let ome_pixel_type = match self.pixel_type {
+                PixelType::I8 => ome::PixelType::Int8,
+                PixelType::U8 => ome::PixelType::Uint8,
+                PixelType::I16 => ome::PixelType::Int16,
+                PixelType::U16 => ome::PixelType::Uint16,
+                PixelType::I32 => ome::PixelType::Int32,
+                PixelType::U32 => ome::PixelType::Uint32,
+                PixelType::F32 => ome::PixelType::Float,
+                PixelType::F64 => ome::PixelType::Double,
+                _ => ome::PixelType::Uint16,
+            };
+
+            let mut pixels = ome::Pixels {
+                id: "Pixels:0".to_string(),
+                size_x: self.shape.x as i32,
+                size_y: self.shape.y as i32,
+                size_z: self.shape.z as i32,
+                size_c: self.shape.c as i32,
+                size_t: self.shape.t as i32,
+                dimension_order: ome::PixelsDimensionOrderType::Xyczt,
+                r#type: ome_pixel_type,
+                physical_size_x: pixel_size.map(|v| v as f32),
+                physical_size_x_unit: ome::UnitsLength::um,
+                physical_size_y: pixel_size.map(|v| v as f32),
+                physical_size_y_unit: ome::UnitsLength::um,
+                physical_size_z: interval_z.map(|v| v as f32),
+                physical_size_z_unit: ome::UnitsLength::um,
+                time_increment: interval_t.map(|v| v as f32),
+                time_increment_unit: ome::UnitsTime::s,
+                significant_bits: None,
+                interleaved: None,
+                big_endian: None,
+                channel: Vec::new(),
+                bin_data: Vec::new(),
+                tiff_data: Vec::new(),
+                metadata_only: None,
+                plane: Vec::new(),
+            };
+
+            for c in 0..self.shape.c {
+                pixels.channel.push(ome::Channel {
+                    id: format!("Channel:{}", c),
+                    name: Some(format!("Channel {}", c)),
+                    ..Default::default()
+                });
+            }
+
+            for c in 0..self.shape.c {
+                for z in 0..self.shape.z {
+                    for t in 0..self.shape.t {
+                        pixels.plane.push(ome::Plane {
+                            the_c: Some(c as i32),
+                            the_z: Some(z as i32),
+                            the_t: Some(t as i32),
+                            delta_t: interval_t.and_then(|v| {
+                                if v > 0.0 {
+                                    Some((t as f64 * v) as f32)
+                                } else {
+                                    None
+                                }
+                            }),
+                            ..Default::default()
+                        });
+                    }
+                }
+            }
+
+            let mut ome = Ome::default();
+            ome.image.push(ome::Image {
+                id: "Image:0".to_string(),
+                name: self
+                    .path
+                    .file_name()
+                    .map(|n| n.to_string_lossy().to_string()),
+                pixels,
+                instrument_ref: None,
+                objective_settings: None,
+                acquisition_date: None,
+                description: None,
+                experimenter_ref: None,
+                experiment_ref: None,
+                experimenter_group_ref: None,
+                imaging_environment: None,
+                stage_label: None,
+                roi_ref: Vec::new(),
+                microbeam_manipulation_ref: Vec::new(),
+                annotation_ref: Vec::new(),
+            });
+            Ok(ome)
+        }
+    }
+
+    fn get_frame(&self, c: usize, z: usize, t: usize) -> Result<Frame, Error> {
+        let (page_idx, offset, stride) = if self.p_ndim == 3 {
+            (t * self.shape.z + z, c, self.n_samples)
+        } else {
+            (c + z * self.shape.c + t * self.shape.c * self.shape.z, 0, 1)
+        };
+
+        let mut reader = self.get_reader()?;
+        reader.seek_to_image(page_idx)?;
+
+        let data: DecodingSampleBuffer = if self.planar {
+            let mut result = DecodingSampleBuffer::U16(Vec::new());
+            reader.read_image_to_buffer(&mut result)?;
+            result
+        } else {
+            reader.read_image()?
+        };
+
+        macro_rules! to_frame {
+            ($var:ident, $ty:ty) => {{
+                let v: Vec<$ty> = match data {
+                    DecodingSampleBuffer::$var(v) => v,
+                    _ => {
+                        return Err(Error::Parse(format!(
+                            "type mismatch: expected {}",
+                            stringify!($var)
+                        )));
+                    }
+                };
+                let yx = self.shape.y * self.shape.x;
+                if self.planar {
+                    let chan: Vec<$ty> = v.into_iter().skip(c * yx).take(yx).collect();
+                    Ok(ArrayT::$var(Array2::from_shape_vec(
+                        (self.shape.y, self.shape.x),
+                        chan,
+                    )?))
+                } else if stride > 1 {
+                    let chan: Vec<$ty> = v
+                        .into_iter()
+                        .skip(offset)
+                        .step_by(stride)
+                        .take(yx)
+                        .collect();
+                    Ok(ArrayT::$var(Array2::from_shape_vec(
+                        (self.shape.y, self.shape.x),
+                        chan,
+                    )?))
+                } else {
+                    Ok(ArrayT::$var(Array2::from_shape_vec(
+                        (self.shape.y, self.shape.x),
+                        v,
+                    )?))
+                }
+            }};
+        }
+
+        match self.pixel_type {
+            PixelType::U8 => to_frame!(U8, u8),
+            PixelType::U16 => to_frame!(U16, u16),
+            PixelType::U32 => to_frame!(U32, u32),
+            PixelType::U64 => to_frame!(U64, u64),
+            PixelType::I8 => to_frame!(I8, i8),
+            PixelType::I16 => to_frame!(I16, i16),
+            PixelType::I32 => to_frame!(I32, i32),
+            PixelType::I64 => to_frame!(I64, i64),
+            PixelType::F32 => to_frame!(F32, f32),
+            PixelType::F64 => to_frame!(F64, f64),
+            _ => Err(Error::NotImplemented(format!(
+                "unsupported TIFF pixel type: {:?}",
+                self.pixel_type
+            ))),
+        }
+    }
+
+    fn path(&self) -> &Path {
+        &self.path
+    }
+    fn series(&self) -> usize {
+        self.series
+    }
+    fn position(&self) -> usize {
+        self.position
+    }
+    fn raw_shape(&self) -> &Shape {
+        &self.shape
+    }
+    fn pixel_type(&self) -> &PixelType {
+        &self.pixel_type
+    }
+
+    fn get_available_positions<P>(path: P, _series: usize) -> Result<HashSet<usize>, Error>
+    where
+        P: AsRef<Path>,
+    {
+        let path = path.as_ref();
+        if Self::can_open(path)? {
+            Ok(HashSet::from([0]))
+        } else {
+            Err(Error::InvalidReader(
+                "TiffReader".to_string(),
+                path.to_string_lossy().to_string(),
+                "it is not a tiff".to_string(),
+            ))
+        }
+    }
+
+    fn get_available_series<P>(path: P) -> Result<HashSet<usize>, Error>
+    where
+        P: AsRef<Path>,
+    {
+        let path = path.as_ref();
+        if Self::can_open(path)? {
+            Ok(HashSet::from([0]))
+        } else {
+            Err(Error::InvalidReader(
+                "TiffReader".to_string(),
+                path.to_string_lossy().to_string(),
+                "it is not a tiff".to_string(),
+            ))
+        }
+    }
+}
+
+fn tag_first_u16<R: Read + Seek>(reader: &mut RefMut<Decoder<R>>, tag: Tag) -> Option<u16> {
+    let val = reader.current_ifd().get_tag(tag).ok()?;
+    match val {
+        Value::Short(v) => Some(v),
+        Value::List(list) => list.first().and_then(|v| match v {
+            Value::Short(v) => Some(*v),
+            _ => None,
+        }),
+        _ => None,
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use ndarray::s;
+
+    fn open(file: &str) -> Result<TiffReader, Error> {
+        let path = std::env::current_dir()?
+            .join("tests")
+            .join("files")
+            .join(file);
+        TiffReader::new(&path, 0, 0)
+    }
+
+    macro_rules! test_metadata {
+        ($($name:ident: $file:expr $(,)?)*) => {
+            $(
+                #[test]
+                fn $name() -> Result<(), Error> {
+                    let ts = open($file)?;
+                    println!("{}", ts.view().squeeze()?.summary()?);
+
+                    let view = ts.view();
+                    let first_frame = view.slice_cztyx(s![0, 0, 0, .., ..])?;
+                    let array = first_frame.as_array::<u16>()?;
+                    println!("{:?}", array);
+                    assert_eq!(ts.shape.y, array.shape()[0]);
+                    assert_eq!(ts.shape.x, array.shape()[1]);
+
+                    if view.shape.t > 1 {
+                        let second_frame = view.slice_cztyx(s![0, 0, 1, .., ..])?;
+                        let array = second_frame.as_array::<u16>()?;
+                        println!("{:?}", array);
+                        assert_eq!(ts.shape.y, array.shape()[0]);
+                        assert_eq!(ts.shape.x, array.shape()[1]);
+                    }
+
+                    Ok(())
+                }
+            )*
+        };
+    }
+
+    test_metadata! {
+        metadata_a: "tiff/20251014_20-Pos_000_000_mask.tif",
+        metadata_b: "tiff/20251014_20-Pos_000_000_max.tif",
+        metadata_c: "tiff/20251014_20-Pos_000_000_loc_results_Cy3.tif",
+        metadata_e: "tiff/1xp53-01-AP1.tiff",
+        metadata_f: "tiff/test.tif",
+        metadata_g: "tiff/YTL1849A111_2023_05_04__14_46_19_cellnr_1_track.tif",
+        metadata_h: "tiff/Wapl-AID_C46_13_tagRFPt_MCP_well1_video1_1_max.tif",
+        metadata_i: "tiff/R1_NRF1_Gtf3a500nM_REAL_20260731_25pct_500ms_001_MMStack_PosR1_NRF1_Gtf3a500nM_REAL_20260731_25pct_500ms_001.ome.tif"
+        metadata_j: "tiff/20260922_Pos0_max.tif"
+    }
+}

ndbioimage/src/readers/tiffseq.rs

@@ -8,7 +8,7 @@ use std::collections::{HashMap, HashSet};
 use std::hash::{Hash, Hasher};
 use std::path::{Path, PathBuf};
 use std::str::FromStr;
-use tiff::decoder::{Decoder, DecodingResult};
+use tiff::decoder::{Decoder, DecodingSampleBuffer};
 use tiff::tags::Tag;
 
 /// reader for sequences of tiff files (one file per z/time plane)
@@ -84,9 +84,13 @@ impl TiffSeqReader {
 
     fn read_tiff_dimensions(path: &Path) -> Result<(usize, usize), Error> {
         let file = std::fs::File::open(path)?;
-        let mut decoder = Decoder::new(file)?;
-        let width = decoder.get_tag(Tag::ImageWidth)?.into_u32()? as usize;
-        let height = decoder.get_tag(Tag::ImageLength)?.into_u32()? as usize;
+        let mut decoder = Decoder::open(file)?;
+        decoder.next_image()?;
+        let width = decoder.current_ifd().get_tag(Tag::ImageWidth)?.into_u32()? as usize;
+        let height = decoder
+            .current_ifd()
+            .get_tag(Tag::ImageLength)?
+            .into_u32()? as usize;
         Ok((width, height))
     }
 }
@@ -235,6 +239,13 @@ impl Reader for TiffSeqReader {
         })
     }
 
+    fn can_open<P>(path: P) -> Result<bool, Error>
+    where
+        P: AsRef<Path>,
+    {
+        Ok(path.as_ref().is_dir())
+    }
+
     fn metadata(&self) -> Result<Ome, Error> {
         let mut ome = Ome::default();
         let info = self.metadata_map.get("Info").and_then(|v| v.as_mapping());
@@ -514,38 +525,39 @@ impl Reader for TiffSeqReader {
             .ok_or_else(|| Error::OutOfBounds(c.max(z).max(t) as isize, 0))?;
 
         let rdr = std::fs::File::open(file)?;
-        let mut decoder = Decoder::new(rdr)?;
+        let mut decoder = Decoder::open(rdr)?;
+        decoder.next_image()?;
 
         match decoder.read_image()? {
-            DecodingResult::U8(d) => Ok(ArrayT::U8(Array2::from_shape_vec(
+            DecodingSampleBuffer::U8(d) => Ok(ArrayT::U8(Array2::from_shape_vec(
                 (self.shape.y, self.shape.x),
                 d,
             )?)),
-            DecodingResult::U16(d) => Ok(ArrayT::U16(Array2::from_shape_vec(
+            DecodingSampleBuffer::U16(d) => Ok(ArrayT::U16(Array2::from_shape_vec(
                 (self.shape.y, self.shape.x),
                 d,
             )?)),
-            DecodingResult::U32(d) => Ok(ArrayT::U32(Array2::from_shape_vec(
+            DecodingSampleBuffer::U32(d) => Ok(ArrayT::U32(Array2::from_shape_vec(
                 (self.shape.y, self.shape.x),
                 d,
             )?)),
-            DecodingResult::I8(d) => Ok(ArrayT::I8(Array2::from_shape_vec(
+            DecodingSampleBuffer::I8(d) => Ok(ArrayT::I8(Array2::from_shape_vec(
                 (self.shape.y, self.shape.x),
                 d,
             )?)),
-            DecodingResult::I16(d) => Ok(ArrayT::I16(Array2::from_shape_vec(
+            DecodingSampleBuffer::I16(d) => Ok(ArrayT::I16(Array2::from_shape_vec(
                 (self.shape.y, self.shape.x),
                 d,
             )?)),
-            DecodingResult::I32(d) => Ok(ArrayT::I32(Array2::from_shape_vec(
+            DecodingSampleBuffer::I32(d) => Ok(ArrayT::I32(Array2::from_shape_vec(
                 (self.shape.y, self.shape.x),
                 d,
             )?)),
-            DecodingResult::F32(d) => Ok(ArrayT::F32(Array2::from_shape_vec(
+            DecodingSampleBuffer::F32(d) => Ok(ArrayT::F32(Array2::from_shape_vec(
                 (self.shape.y, self.shape.x),
                 d,
             )?)),
-            DecodingResult::F64(d) => Ok(ArrayT::F64(Array2::from_shape_vec(
+            DecodingSampleBuffer::F64(d) => Ok(ArrayT::F64(Array2::from_shape_vec(
                 (self.shape.y, self.shape.x),
                 d,
             )?)),
@@ -567,7 +579,7 @@ impl Reader for TiffSeqReader {
         self.position
     }
 
-    fn shape(&self) -> &Shape {
+    fn raw_shape(&self) -> &Shape {
         &self.shape
     }
 
@@ -575,11 +587,20 @@ impl Reader for TiffSeqReader {
         &self.pixel_type
     }
 
-    fn get_available_positions<P>(_path: P, _series: usize) -> Result<HashSet<usize>, Error>
+    fn get_available_positions<P>(path: P, _series: usize) -> Result<HashSet<usize>, Error>
     where
         P: AsRef<Path>,
     {
-        Ok(HashSet::from([0]))
+        let path = path.as_ref();
+        if Self::can_open(path)? {
+            Ok(HashSet::from([0]))
+        } else {
+            Err(Error::InvalidReader(
+                "TiffSeqReader".to_string(),
+                path.to_string_lossy().to_string(),
+                "not a tiff sequence".to_string(),
+            ))
+        }
     }
 
     fn get_available_series<P>(path: P) -> Result<HashSet<usize>, Error>

ndbioimage/src/tiffwrite.rs

@@ -178,7 +178,7 @@ where
             PixelType::U64 => self.save_as_tiff_with_type::<u64, P>(path, options)?,
             PixelType::I128 => self.save_as_tiff_with_type::<i64, P>(path, options)?,
             PixelType::U128 => self.save_as_tiff_with_type::<u64, P>(path, options)?,
-            PixelType::F128 => self.save_as_tiff_with_type::<f64, P>(path, options)?,
+            PixelType::F128 => Err(Error::F128NotSupported)?,
         }
 
         Ok(())
@@ -321,4 +321,44 @@ mod tests {
         batch_to_tiff(&files_in, &files_out, None, None, true, true, None)?;
         Ok(())
     }
+
+    /// a rewritten file must not be served from the process-wide cache
+    #[cfg(feature = "tiff")]
+    #[test]
+    fn cache_invalidated_after_rewrite() -> Result<(), Error> {
+        use crate::view::View;
+        use ndarray::Array2;
+        use std::time::{Duration, SystemTime};
+        use tiffwrite::IJTiffFile;
+
+        let dir = std::env::temp_dir().join(format!("ndbioimage-rewrite-{}", std::process::id()));
+        create_dir_all(&dir)?;
+        let path = dir.join("cache_rewrite.tif");
+
+        let write = |value: u16| -> Result<(), Error> {
+            let mut tiff = IJTiffFile::new(&path)?;
+            tiff.save(Array2::<u16>::from_elem((4, 4), value), 0, 0, 0)?;
+            Ok(())
+        };
+
+        write(1)?;
+        let view: View<_, DynReader> = View::from_path(&path)?;
+        let frame: Array2<u16> = view.get_frame(0, 0, 0)?;
+        assert_eq!(frame[(0, 0)], 1);
+
+        write(2)?;
+        // make the modification time unambiguous even on filesystems with a
+        // coarse timestamp resolution
+        let file = std::fs::OpenOptions::new().write(true).open(&path)?;
+        file.set_modified(SystemTime::UNIX_EPOCH + Duration::from_secs(1_000_000_000))?;
+        drop(file);
+
+        let view: View<_, DynReader> = View::from_path(&path)?;
+        let frame: Array2<u16> = view.get_frame(0, 0, 0)?;
+        assert_eq!(frame[(0, 0)], 2);
+
+        let _ = std::fs::remove_file(&path);
+        let _ = std::fs::remove_dir(&dir);
+        Ok(())
+    }
 }

ndbioimage/src/transforms.rs

@@ -0,0 +1,877 @@
+use crate::axes::Axis;
+use crate::error::Error;
+use crate::readers::{DynReader, Frame, PixelType, Reader};
+use crate::view::View;
+use image_registration::{bspline::DefaultValue, filter, transform::Transform};
+use ndarray::{Array2, Dimension, Ix2, Ix3, IxDyn, array, concatenate, s};
+use rayon::iter::{IntoParallelIterator, ParallelIterator};
+use serde::{Deserialize, Serialize};
+use std::borrow::Cow;
+use std::collections::{HashMap, HashSet};
+use std::fmt::Debug;
+use std::hash::{Hash, Hasher};
+use std::io::Seek;
+use std::path::Path;
+use std::sync::Arc;
+use tiffwrite::IJTiffFile;
+
+#[derive(Clone, Debug, Default, PartialEq, Serialize, Deserialize)]
+struct TransformOld {
+    #[serde(alias = "CenterOfRotationPoint")]
+    center_of_rotation_point: Vec<f64>,
+    #[serde(alias = "Size")]
+    size: Vec<f64>,
+    #[serde(alias = "TransformParameters")]
+    transform_parameters: Vec<f64>,
+    #[serde(alias = "dTransformParameters")]
+    dtransform_parameters: Vec<f64>,
+}
+
+impl TryFrom<HashMap<String, TransformOld>> for Transforms {
+    type Error = Error;
+
+    fn try_from(t: HashMap<String, TransformOld>) -> Result<Self, Error> {
+        Ok(Self {
+            default_value: Default::default(),
+            interpolation_order: None,
+            channel: t
+                .into_iter()
+                .map(|(n, i)| {
+                    let center = i
+                        .size
+                        .iter()
+                        .zip(i.center_of_rotation_point.iter())
+                        .map(|(s, c)| c - s / 2.0 + 0.5)
+                        .collect::<Vec<_>>();
+                    Ok((
+                        n,
+                        Transform::new_around_center(i.transform_parameters, &center)
+                            .with_dparameters(i.dtransform_parameters)
+                            .inverse()?,
+                    ))
+                })
+                .collect::<Result<_, Error>>()?,
+            drift: Vec::new(),
+        })
+    }
+}
+
+pub type DefaultValues = HashMap<String, DefaultValue>;
+
+#[derive(Clone, Debug, Default, Eq, PartialEq, Deserialize)]
+pub struct Transforms {
+    #[serde(default)]
+    pub(crate) interpolation_order: Option<usize>,
+    #[serde(default)]
+    pub(crate) default_value: DefaultValues,
+    /// a transform for each channel by channel name
+    #[serde(default)]
+    pub channel: HashMap<String, Transform<IxDyn>>,
+    /// a transform for each time point
+    #[serde(default)]
+    pub drift: Vec<Transform<IxDyn>>,
+}
+
+impl Hash for Transforms {
+    fn hash<H: Hasher>(&self, state: &mut H) {
+        let mut keys: Vec<_> = self.channel.keys().collect();
+        keys.sort();
+        for key in keys {
+            key.hash(state);
+            self.channel[key].hash(state);
+        }
+        self.drift.hash(state);
+    }
+}
+
+impl PartialEq<Transforms> for &Transforms {
+    fn eq(&self, other: &Transforms) -> bool {
+        **self == *other
+    }
+}
+
+impl Serialize for Transforms {
+    fn serialize<S>(&self, serializer: S) -> Result<S::Ok, S::Error>
+    where
+        S: serde::Serializer,
+    {
+        #[derive(Serialize)]
+        struct HumanReadable<'a> {
+            #[serde(skip_serializing_if = "HashMap::is_empty")]
+            channel: &'a HashMap<String, Transform<IxDyn>>,
+            #[serde(skip_serializing_if = "Vec::is_empty")]
+            drift: &'a Vec<Transform<IxDyn>>,
+        }
+
+        #[derive(Serialize)]
+        struct Binary<'a> {
+            interpolation_order: Option<usize>,
+            default_value: &'a DefaultValues,
+            channel: &'a HashMap<String, Transform<IxDyn>>,
+            drift: &'a Vec<Transform<IxDyn>>,
+        }
+
+        if serializer.is_human_readable() {
+            HumanReadable {
+                channel: &self.channel,
+                drift: &self.drift,
+            }
+            .serialize(serializer)
+        } else {
+            Binary {
+                interpolation_order: self.interpolation_order,
+                default_value: &self.default_value,
+                channel: &self.channel,
+                drift: &self.drift,
+            }
+            .serialize(serializer)
+        }
+    }
+}
+
+/// the kind of corrections the transforms are meant to correct for
+#[derive(Clone, Debug, PartialEq, Copy)]
+pub enum TransformKind {
+    /// no tranforms
+    None,
+    /// correct warp between channels, ex: camera/filter misalignment, chromatic aberrations
+    Channel,
+    /// drift in the miscrocope between time points
+    Drift,
+    /// both warp between channels and drift
+    Both,
+}
+
+impl Transforms {
+    /// the kind of corrections the transforms are meant to correct for
+    pub fn kind(&self) -> TransformKind {
+        match (self.channel.is_empty(), self.drift.is_empty()) {
+            (true, true) => TransformKind::None,
+            (false, true) => TransformKind::Channel,
+            (true, false) => TransformKind::Drift,
+            (false, false) => TransformKind::Both,
+        }
+    }
+
+    /// load the transforms from a yml file
+    pub fn load<P: AsRef<Path>>(path: P) -> Result<Self, Error> {
+        let mut file = std::fs::File::open(&path)?;
+        Ok(
+            if let Ok(t) = yaml_serde::from_reader::<_, Self>(&file)
+                && !(t.channel.is_empty() && t.drift.is_empty())
+            {
+                t
+            } else {
+                file.seek(std::io::SeekFrom::Start(0))?;
+                yaml_serde::from_reader::<_, HashMap<String, TransformOld>>(&file)?.try_into()?
+            },
+        )
+    }
+
+    /// save the transforms to a yml file
+    pub fn save<P: AsRef<Path>>(&self, path: P) -> Result<(), Error> {
+        let file = std::fs::OpenOptions::new()
+            .write(true)
+            .create(true)
+            .truncate(true)
+            .open(path)?;
+        Ok(yaml_serde::to_writer(file, self)?)
+    }
+
+    /// helper function to eiter load transforms from file or calculate and save them
+    pub fn calculate_transforms_from<P, Q, R>(
+        transform_file: Option<P>,
+        bead_files: Option<&[Q]>,
+        drift_file: Option<R>,
+        main_channel: Option<usize>,
+        default_transform: Option<Vec<f64>>,
+    ) -> Result<Self, Error>
+    where
+        P: AsRef<Path>,
+        Q: AsRef<Path>,
+        R: AsRef<Path>,
+    {
+        if let Some(transform_file) = transform_file.as_ref()
+            && let Ok(transforms) = Self::load(transform_file.as_ref().with_extension("yml"))
+        {
+            Ok(transforms)
+        } else {
+            let default_transform = default_transform.map(Transform::new);
+            let mut transforms = Self::default();
+            if let Some(bead_files) = bead_files {
+                let main_channel = main_channel.unwrap_or(0);
+                transforms.channel = Self::calculate_channel_transforms_2d(
+                    bead_files,
+                    main_channel,
+                    default_transform,
+                )?
+                .into_iter()
+                .map(|(k, t)| (k, t.into_dyn()))
+                .collect();
+            };
+            if let Some(drift_file) = drift_file {
+                transforms.drift = Self::calculate_drift_transforms_2d(drift_file)?
+                    .into_iter()
+                    .map(|t| t.into_dyn())
+                    .collect();
+            };
+            if let Some(transform_file) = transform_file {
+                let transform_file = transform_file.as_ref();
+                transforms.save(transform_file.with_extension("yml"))?;
+                if let Some(bead_files) = bead_files {
+                    transforms.save_channel_transform_tiff(
+                        transform_file.with_extension("tif"),
+                        bead_files,
+                    )?;
+                }
+            }
+            Ok(transforms)
+        }
+    }
+
+    /// save a tiff file with the uncorrected and corrected frames of the bead files side by side
+    pub fn save_channel_transform_tiff<P, Q>(&self, file: P, bead_files: &[Q]) -> Result<(), Error>
+    where
+        P: AsRef<Path>,
+        Q: AsRef<Path>,
+    {
+        let mut transforms = self.clone();
+        transforms.drift = Vec::new();
+        let mut tif = IJTiffFile::new(file)?;
+        for (t, bead_file) in bead_files.iter().enumerate() {
+            let mut view = View::<_, DynReader>::from_path(bead_file)?.into_dyn();
+            let axes = view.axes();
+            if axes.contains(&Axis::Z) {
+                view = view.max_proj(Axis::Z)?;
+            }
+            if axes.contains(&Axis::T) {
+                view = view.max_proj(Axis::T)?;
+            }
+            let view_u = view.into_dimensionality::<Ix3>()?;
+            let view_t = view_u.with_transforms(transforms.clone());
+            for c in 0..view_u.shape().c {
+                let a = view_u.slice_cztyx(s![c, 0, 0, .., ..])?.as_array::<u16>()?;
+                let b = view_t.slice_cztyx(s![c, 0, 0, .., ..])?.as_array::<u16>()?;
+                let frame = concatenate(ndarray::Axis(1), &[a.view(), b.view()])?;
+                tif.save(frame, c, 0, t)?;
+            }
+        }
+        Ok(())
+    }
+
+    /// calculate drift transforms in yx using image file
+    pub fn calculate_drift_transforms_2d<P: AsRef<Path>>(
+        file: P,
+    ) -> Result<Vec<Transform<Ix2>>, Error> {
+        View::<_, DynReader>::from_path(file)?.calculate_drift_transforms_2d()
+    }
+
+    /// calculate drift transforms in zyx using image file
+    pub fn calculate_drift_transforms_3d<P: AsRef<Path>>(
+        file: P,
+    ) -> Result<Vec<Transform<Ix3>>, Error> {
+        View::<_, DynReader>::from_path(file)?.calculate_drift_transforms_3d()
+    }
+
+    /// calculate channel transforms in yx using image file
+    pub fn calculate_channel_transforms_2d<P: AsRef<Path>>(
+        bead_files: &[P],
+        main_channel: usize,
+        default_transform: Option<Transform<Ix2>>,
+    ) -> Result<HashMap<String, Transform<Ix2>>, Error> {
+        // TODO: autodetect main channel
+        let mut transforms = Vec::new();
+        let default_transform = default_transform.unwrap_or_default();
+        for file in bead_files {
+            let view = View::<_, DynReader>::from_path(file)?;
+            transforms.push(view.calculate_channel_transforms_2d(main_channel)?)
+        }
+        let transforms = transforms.into_iter().flat_map(|map| map.into_iter()).fold(
+            HashMap::new(),
+            |mut acc, (channel, transform)| {
+                acc.entry(channel)
+                    .or_insert_with(Vec::new)
+                    .push(transform * &default_transform);
+                acc
+            },
+        );
+        Ok(transforms
+            .iter()
+            .map(|(channel, ts)| {
+                let n = ts.len();
+                let parameters = ts
+                    .iter()
+                    .fold(vec![0.0; 6], |mut acc, t| {
+                        for (a, p) in acc.iter_mut().zip(&t.parameters) {
+                            *a += p;
+                        }
+                        acc
+                    })
+                    .iter()
+                    .map(|i| i / n as f64)
+                    .collect::<Vec<_>>();
+                let dparameters = if n == 1 {
+                    ts[0].dparameters.clone()
+                } else {
+                    ts.iter()
+                        .fold(vec![0.0; 6], |mut acc, t| {
+                            for ((a, p), dp) in acc.iter_mut().zip(&parameters).zip(&t.dparameters)
+                            {
+                                *a += (dp - p).powi(2);
+                            }
+                            acc
+                        })
+                        .iter()
+                        .map(|i| (i / n as f64).sqrt())
+                        .collect::<Vec<_>>()
+                };
+                (
+                    channel.to_string(),
+                    Transform::new(parameters).with_dparameters(dparameters),
+                )
+            })
+            .collect())
+    }
+
+    /// calculate channel transforms in zyx using image file
+    pub fn calculate_channel_transforms_3d<P: AsRef<Path>>(
+        bead_files: &[P],
+        main_channel: usize,
+        default_transform: Option<Transform<Ix3>>,
+    ) -> Result<HashMap<String, Transform<Ix3>>, Error> {
+        let mut transforms = Vec::new();
+        let default_transform = default_transform.unwrap_or_default();
+        for file in bead_files {
+            let view = View::<_, DynReader>::from_path(file)?;
+            transforms.push(view.calculate_channel_transforms_3d(main_channel)?)
+        }
+        let channels = transforms
+            .iter()
+            .flat_map(|t| t.keys())
+            .collect::<HashSet<_>>();
+        let mut average_transforms = HashMap::new();
+        for channel in channels {
+            let (s, n) = transforms
+                .iter()
+                .filter_map(|t| {
+                    if t.contains_key(channel) {
+                        Some((&t[channel] * &default_transform).matrix())
+                    } else {
+                        None
+                    }
+                })
+                .fold((Array2::zeros((4, 4)), 0), |(s, n), t| (s + t, n + 1));
+            let matrix = s / (n as f64);
+            let (s, n) = transforms
+                .iter()
+                .filter_map(|t| {
+                    if t.contains_key(channel) {
+                        Some(((&t[channel] * &default_transform).matrix() - &matrix).powi(2))
+                    } else {
+                        None
+                    }
+                })
+                .fold((Array2::zeros((4, 4)), 0), |(s, n), t| (s + t, n + 1));
+            let dmatrix = (s / n as f64).sqrt();
+            average_transforms.insert(
+                channel.clone(),
+                Transform::default()
+                    .with_matrix(matrix.view())
+                    .with_dmatrix(dmatrix.view()),
+            );
+        }
+        Ok(average_transforms)
+    }
+
+    /// get the combination of channel and drift transforms for a certain channel and time point
+    pub fn get_transform(
+        &self,
+        channel: &str,
+        t: usize,
+    ) -> Result<Option<Cow<'_, Transform<IxDyn>>>, Error> {
+        match (self.channel.get(channel), self.drift.get(t)) {
+            (Some(channel), Some(drift)) => match (channel.ndim(), drift.ndim()) {
+                (2, 2) | (3, 3) => Ok(Some(Cow::Owned(drift * channel))),
+                (2, 3) => Ok(Some(Cow::Owned(drift.add_dimension(0) * channel))),
+                (3, 2) => Ok(Some(Cow::Owned(drift * channel.add_dimension(0)))),
+                (2 | 3, n) | (n, _) => Err(Error::TransformationNotImplemented(n)),
+            },
+            (Some(channel), None) => Ok(Some(Cow::Borrowed(channel))),
+            (None, Some(drift)) => Ok(Some(Cow::Borrowed(drift))),
+            (None, None) => Ok(None),
+        }
+    }
+}
+
+impl<D: Dimension, R: Reader> View<D, R> {
+    pub fn get_transforms(&self) -> &Transforms {
+        &self.transforms
+    }
+
+    pub fn set_transforms(&mut self, transforms: Transforms) {
+        self.transforms = transforms;
+    }
+
+    /// apply transforms to a new view
+    pub fn with_transforms(&self, transforms: Transforms) -> Self {
+        let mut new = self.clone();
+        new.transforms = transforms;
+        new
+    }
+
+    /// load the transforms from a yml file and apply to a new view
+    pub fn with_transforms_from_yaml(mut self, path: &Path) -> Result<Self, Error> {
+        self.set_transforms(Transforms::load(path)?);
+        Ok(self)
+    }
+
+    /// load the transforms from a yml file and apply to the view
+    pub fn load_transforms_from_yaml(&mut self, path: &Path) -> Result<(), Error> {
+        self.set_transforms(Transforms::load(path)?);
+        Ok(())
+    }
+
+    /// calculate channel transforms in yx using this image view
+    pub fn calculate_channel_transforms_2d(
+        &self,
+        main_channel: usize,
+    ) -> Result<HashMap<String, Transform<Ix2>>, Error> {
+        let main_max = self
+            .slice_cztyx(s![main_channel, .., 0, .., ..])?
+            .max_proj(Axis::Z)?
+            .as_array::<f64>()?;
+        let mut transforms = HashMap::new();
+        for (channel, metadata) in self.metadata()?.image[0].pixels.channel.iter().enumerate() {
+            let name = metadata
+                .name
+                .clone()
+                .unwrap_or_else(|| format!("channel{}", channel));
+            if channel == main_channel {
+                transforms.insert(name, Transform::default());
+            } else {
+                let max = self
+                    .slice_cztyx(s![channel, .., 0, .., ..])?
+                    .max_proj(Axis::Z)?
+                    .as_array::<f64>()?;
+                transforms.insert(
+                    name,
+                    Transform::register_affine(main_max.view(), max.view())?,
+                );
+            }
+        }
+        Ok(transforms)
+    }
+
+    /// calculate channel transforms in zyx using this image view
+    pub fn calculate_channel_transforms_3d(
+        &self,
+        main_channel: usize,
+    ) -> Result<HashMap<String, Transform<Ix3>>, Error> {
+        let main_max = self
+            .slice_cztyx(s![main_channel, .., 0, .., ..])?
+            .as_array::<f64>()?;
+        let mut transforms = HashMap::new();
+        for (channel, metadata) in self.metadata()?.image[0].pixels.channel.iter().enumerate() {
+            let name = metadata
+                .name
+                .clone()
+                .unwrap_or_else(|| format!("channel{}", channel));
+            if channel == main_channel {
+                transforms.insert(name, Transform::default());
+            } else {
+                let max = self
+                    .slice_cztyx(s![channel, .., 0, .., ..])?
+                    .as_array::<f64>()?;
+                transforms.insert(
+                    name,
+                    Transform::register_affine(main_max.view(), max.view())?,
+                );
+            }
+        }
+        Ok(transforms)
+    }
+
+    /// calculate drift transforms in yx using this image view
+    pub fn calculate_drift_transforms_2d(&self) -> Result<Vec<Transform<Ix2>>, Error> {
+        let mut view = self.clone().into_dyn();
+        let axes = view.axes();
+        if axes.contains(&Axis::Z) {
+            view = view.max_proj(Axis::Z)?;
+        }
+        if axes.contains(&Axis::C) {
+            view = view.max_proj(Axis::C)?;
+        }
+        view = view.permute_axes(&[Axis::T, Axis::Y, Axis::X])?;
+        let n = view.shape().t;
+        let m = (n as f64).sqrt().round() as usize;
+        let t_groups = (0..=m)
+            .map(|i| (((i * n) as f64) / (m as f64)).round() as usize)
+            .collect::<Vec<_>>();
+        let t_keys = t_groups
+            .array_windows()
+            .map(|[start, end]| (start + end) / 2)
+            .collect::<Vec<_>>();
+        let t_pairs = t_keys
+            .iter()
+            .zip(t_groups.array_windows())
+            .flat_map(|(t_key, [start, end])| (*start..*end).map(|i| (*t_key, i)))
+            .chain(t_keys.array_windows().map(|[start, end]| (*start, *end)))
+            .collect::<Vec<_>>();
+        let fmaxz_keys = t_keys
+            .iter()
+            .map(|t| {
+                Ok((
+                    t,
+                    filter::gaussian_smooth(
+                        &view.slice(s![*t, .., ..])?.as_array::<f64>()?,
+                        &[5.0; 2],
+                    )?,
+                ))
+            })
+            .collect::<Result<HashMap<_, _>, Error>>()?;
+        let shifts = t_pairs
+            .into_par_iter()
+            .map(|(t_key, t)| {
+                if t_key == t {
+                    Ok(Transform::<Ix2>::from_translation(&vec![0.0; view.ndim()]))
+                } else {
+                    if let Some(fmaxz) = fmaxz_keys.get(&t_key) {
+                        Ok(Transform::<Ix2>::register_translation(
+                            fmaxz_keys[&t_key].view(),
+                            fmaxz,
+                        )?)
+                    } else {
+                        let fmaxz = filter::gaussian_smooth(
+                            &view.slice(s![t, .., ..])?.as_array::<f64>()?,
+                            &[5.0; 2],
+                        )?;
+                        Ok(Transform::<Ix2>::register_translation(
+                            fmaxz_keys[&t_key].view(),
+                            fmaxz.view(),
+                        )?)
+                    }
+                }
+            })
+            .collect::<Result<Vec<_>, Error>>()?;
+        let mut t = Transform::<Ix2>::from_translation(&[0.0; 2]);
+        let mut key_shifts = vec![t.clone()];
+        key_shifts.extend(shifts.iter().skip(n).map(|i| {
+            t = i * &t;
+            t.clone()
+        }));
+        Ok(key_shifts
+            .iter()
+            .zip(t_groups.array_windows())
+            .flat_map(|(s, [start, end])| {
+                let mut t = s.clone();
+                (*start..*end)
+                    .map(|i| {
+                        t = &shifts[i] * &t;
+                        t.clone()
+                    })
+                    .collect::<Vec<_>>()
+            })
+            .collect())
+    }
+
+    /// calculate drift transforms in zyx using this image view
+    pub fn calculate_drift_transforms_3d(&self) -> Result<Vec<Transform<Ix3>>, Error> {
+        todo!()
+    }
+
+    /// retrieve a frame and warp it according to the transform
+    pub(crate) fn get_transformed_frame(
+        &self,
+        c: usize,
+        z: usize,
+        t: usize,
+        transform: &Transform<IxDyn>,
+        order: usize,
+        default_value: &DefaultValue,
+    ) -> Result<Arc<Frame>, Error> {
+        match transform.ndim() {
+            2 => {
+                let tr = transform.clone().into_dimensionality::<Ix2>()?;
+
+                macro_rules! interpolate {
+                    ($array:expr ,$var:ident, bool) => {
+                        Ok(Arc::from(Frame::$var(
+                            tr.interpolate_par_order(
+                                &$array.mapv(|i| if i { 1.0 } else { 0.0 }),
+                                order,
+                                default_value,
+                            )?
+                            .mapv(|i| i > 0.5),
+                        )))
+                    };
+                    ($array:expr ,$var:ident, f64) => {
+                        Ok(Arc::from(Frame::$var(tr.interpolate_par_order(
+                            $array,
+                            order,
+                            default_value,
+                        )?)))
+                    };
+                    ($array:expr ,$var:ident, $t:ty) => {
+                        Ok(Arc::from(Frame::$var(
+                            tr.interpolate_par_order($array, order, default_value)?
+                                .mapv(|i| i as $t),
+                        )))
+                    };
+                }
+
+                match self.get_cached_frame(c, z, t, None)?.as_ref() {
+                    Frame::Bool(frame) => interpolate!(frame, Bool, bool),
+                    Frame::I8(frame) => interpolate!(frame, I8, i8),
+                    Frame::U8(frame) => interpolate!(frame, U8, u8),
+                    Frame::I16(frame) => interpolate!(frame, I16, i16),
+                    Frame::U16(frame) => interpolate!(frame, U16, u16),
+                    Frame::I32(frame) => interpolate!(frame, I32, i32),
+                    Frame::U32(frame) => interpolate!(frame, U32, u32),
+                    Frame::U64(frame) => interpolate!(frame, U64, u64),
+                    Frame::I64(frame) => interpolate!(frame, I64, i64),
+                    Frame::I128(frame) => interpolate!(frame, I128, i128),
+                    Frame::U128(frame) => interpolate!(frame, U128, u128),
+                    Frame::F32(frame) => interpolate!(frame, F32, f32),
+                    Frame::F64(frame) => interpolate!(frame, F64, f64),
+                }
+            }
+            3 => {
+                let tr = transform.clone().into_dimensionality::<Ix3>()?;
+                let shape = self.shape();
+                let center = shape.iter().map(|&i| (i as f64) / 2.0).collect::<Vec<_>>();
+                let zf = z as f64;
+                let yf = shape.y as f64 - 1.0;
+                let xf = shape.x as f64 - 1.0;
+                let corners = array![[zf, 0.0, 0.0], [zf, 0.0, xf], [zf, yf, 0.0], [zf, yf, xf]];
+                let extent = tr.transform_points(&corners, &center)?;
+                let z_extent = extent.slice(s![.., 2]);
+                let z_min = (z_extent
+                    .iter()
+                    .min_by(|a, b| a.total_cmp(b))
+                    .unwrap()
+                    .floor() as isize
+                    - 2)
+                .max(0) as usize;
+                let z_max = (z_extent
+                    .iter()
+                    .max_by(|a, b| a.total_cmp(b))
+                    .unwrap()
+                    .ceil() as isize
+                    + 2)
+                .min(shape.z as isize - 1)
+                .max(0) as usize;
+                let frames = (z_min..=z_max)
+                    .map(|zi| self.get_cached_frame(c, zi, t, None))
+                    .collect::<Result<Vec<_>, Error>>()?;
+
+                // TODO; check these translations
+                let u0 = tr.as_translation(&[z_min as f64, 0.0, 0.0]);
+                let u1 = tr.as_translation(&[-(z_min as f64), 0.0, 0.0]);
+                let tr = u0 * tr * u1;
+
+                match frames.first().map(|f| f.as_ref()) {
+                    Some(first) => {
+                        macro_rules! interpolate {
+                            ($var:ident, bool) => {{
+                                let views = frames
+                                    .iter()
+                                    .map(|f| match f.as_ref() {
+                                        Frame::$var(a) => Ok(a.view()),
+                                        _ => unreachable!(),
+                                    })
+                                    .collect::<Result<Vec<_>, Error>>()?;
+                                let stack = ndarray::stack(ndarray::Axis(0), &views)?;
+                                let transformed = tr.interpolate_par_order(
+                                    &stack.mapv(|i| if i { 1.0 } else { 0.0 }),
+                                    order,
+                                    default_value,
+                                )?;
+                                Ok(Arc::from(Frame::$var(
+                                    transformed.slice(s![0, .., ..]).mapv(|i| i > 0.5),
+                                )))
+                            }};
+                            ($var:ident, f64) => {{
+                                let views = frames
+                                    .iter()
+                                    .map(|f| match f.as_ref() {
+                                        Frame::$var(a) => Ok(a.view()),
+                                        _ => unreachable!(),
+                                    })
+                                    .collect::<Result<Vec<_>, Error>>()?;
+                                let stack = ndarray::stack(ndarray::Axis(0), &views)?;
+                                let transformed =
+                                    tr.interpolate_par_order(&stack, order, default_value)?;
+                                Ok(Arc::from(Frame::$var(
+                                    transformed.slice(s![0, .., ..]).to_owned(),
+                                )))
+                            }};
+                            ($var:ident, $t:ty) => {{
+                                let views = frames
+                                    .iter()
+                                    .map(|f| match f.as_ref() {
+                                        Frame::$var(a) => Ok(a.view()),
+                                        _ => unreachable!(),
+                                    })
+                                    .collect::<Result<Vec<_>, Error>>()?;
+                                let stack = ndarray::stack(ndarray::Axis(0), &views)?;
+                                let transformed =
+                                    tr.interpolate_par_order(&stack, order, default_value)?;
+                                Ok(Arc::from(Frame::$var(
+                                    transformed.slice(s![0, .., ..]).mapv(|i| i as $t),
+                                )))
+                            }};
+                        }
+                        match first {
+                            Frame::Bool(_) => interpolate!(Bool, bool),
+                            Frame::I8(_) => interpolate!(I8, i8),
+                            Frame::U8(_) => interpolate!(U8, u8),
+                            Frame::I16(_) => interpolate!(I16, i16),
+                            Frame::U16(_) => interpolate!(U16, u16),
+                            Frame::I32(_) => interpolate!(I32, i32),
+                            Frame::U32(_) => interpolate!(U32, u32),
+                            Frame::F32(_) => interpolate!(F32, f32),
+                            Frame::F64(_) => interpolate!(F64, f64),
+                            Frame::I64(_) => interpolate!(I64, i64),
+                            Frame::U64(_) => interpolate!(U64, u64),
+                            Frame::I128(_) => interpolate!(I128, i128),
+                            Frame::U128(_) => interpolate!(U128, u128),
+                        }
+                    }
+                    None => {
+                        let shape = self.shape();
+                        Ok(Arc::new(match self.pixel_type() {
+                            PixelType::Bool => Frame::Bool(Array2::default((shape.y, shape.x))),
+                            PixelType::I8 => Frame::I8(Array2::zeros((shape.y, shape.x))),
+                            PixelType::U8 => Frame::U8(Array2::zeros((shape.y, shape.x))),
+                            PixelType::I16 => Frame::I16(Array2::zeros((shape.y, shape.x))),
+                            PixelType::U16 => Frame::U16(Array2::zeros((shape.y, shape.x))),
+                            PixelType::I32 => Frame::I32(Array2::zeros((shape.y, shape.x))),
+                            PixelType::U32 => Frame::U32(Array2::zeros((shape.y, shape.x))),
+                            PixelType::F32 => Frame::F32(Array2::zeros((shape.y, shape.x))),
+                            PixelType::F64 => Frame::F64(Array2::zeros((shape.y, shape.x))),
+                            PixelType::I64 => Frame::I64(Array2::zeros((shape.y, shape.x))),
+                            PixelType::U64 => Frame::U64(Array2::zeros((shape.y, shape.x))),
+                            PixelType::I128 => Frame::I128(Array2::zeros((shape.y, shape.x))),
+                            PixelType::U128 => Frame::U128(Array2::zeros((shape.y, shape.x))),
+                            PixelType::F128 => Err(Error::F128NotSupported)?,
+                        }))
+                    }
+                }
+            }
+            d => Err(Error::TransformationNotImplemented(d)),
+        }
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use crate::axes::Axis;
+    use crate::error::Error;
+    use crate::readers::{DynReader, Reader};
+    use crate::transforms::Transforms;
+    use image_registration::transform::Transform;
+    use ndarray::s;
+    use std::collections::HashMap;
+
+    fn open(file: &str) -> Result<DynReader, Error> {
+        let path = std::env::current_dir()?
+            .join("tests")
+            .join("files")
+            .join(file);
+        DynReader::new(&path, 0, 0)
+    }
+
+    #[test]
+    fn test_transforms() -> Result<(), Box<dyn std::error::Error>> {
+        let _t = Transforms::default();
+
+        Ok(())
+    }
+
+    #[test]
+    fn load_transforms() -> Result<(), Box<dyn std::error::Error>> {
+        let t = Transforms::load("tests/files/transforms.yml")?;
+        println!("t: {:?}", t);
+        assert_eq!(t.channel["TV1"].ndim(), 2);
+        assert_eq!(t.channel["TV2"].ndim(), 2);
+        Ok(())
+    }
+
+    #[test]
+    fn load_old_transforms() -> Result<(), Box<dyn std::error::Error>> {
+        let t = Transforms::load("tests/files/transform_red.yml")?;
+        println!("t: {:?}", t);
+        assert_eq!(t.channel["TV1"].ndim(), 2);
+        assert_eq!(t.channel["TV2"].ndim(), 2);
+        Ok(())
+    }
+
+    #[test]
+    fn transforms_as_array() -> Result<(), Box<dyn std::error::Error>> {
+        let file = "czi/beads_2020_08_05__11_43_07.czi";
+        let reader = open(file)?;
+        let view = reader.view();
+        let mut channel = HashMap::new();
+        channel.insert(
+            "TV1".to_string(),
+            Transform::new(vec![1.0, 0.0, 0.0, 1.0, 0.0, 0.0]),
+        );
+        channel.insert(
+            "TV2".to_string(),
+            Transform::new(vec![
+                0.9995868322142533,
+                0.012140955373934502,
+                0.012592501668218145,
+                0.8571690967587099,
+                -5.270859902118885,
+                -5.81118547158953,
+            ]),
+        );
+        let view_t = view.with_transforms(Transforms {
+            default_value: Default::default(),
+            interpolation_order: None,
+            channel,
+            drift: Vec::new(),
+        });
+
+        let s_a = view.max_proj(Axis::Z)?.slice_cztyx(s![1, 0, .., .., ..])?;
+
+        println!(
+            "s_a: {:?}, {:?}, {:?}",
+            s_a.get_slice(),
+            s_a.get_axes(),
+            s_a.get_operations()
+        );
+
+        let s_b = view_t
+            .max_proj(Axis::Z)?
+            .slice_cztyx(s![1, 0, .., .., ..])?;
+
+        println!(
+            "s_b: {:?}, {:?}, {:?}",
+            s_b.get_slice(),
+            s_b.get_axes(),
+            s_b.get_operations()
+        );
+
+        let a = s_a.as_array::<u16>()?;
+        let b = s_b.as_array::<u16>()?;
+        assert_ne!(a, b);
+
+        // let a = view.slice_cztyx(s![1, 0, 0, .., ..])?.as_array::<u16>()?;
+        // let b = view_t.slice_cztyx(s![1, 0, 0, .., ..])?.as_array::<u16>()?;
+        // assert_ne!(a, b);
+        //
+        // let a = view
+        //     .max_proj(Axis::Z)?
+        //     .slice(s![1, 0, .., ..])?
+        //     .as_array::<u16>()?;
+        // let b = view_t
+        //     .max_proj(Axis::Z)?
+        //     .slice(s![1, 0, .., ..])?
+        //     .as_array::<u16>()?;
+        // assert_ne!(a, b);
+
+        Ok(())
+    }
+}

ndbioimage/tests/files/tiffseq/YTL1841B2-2-1_1hr_DMSO_galinduction_1/display_and_comments.txt

@@ -0,0 +1,12 @@
+{
+  "Channels": [{
+    "Name": "mScarlet_GFP-mSc-filter",
+    "DisplayMode": 1,
+    "HistogramMax": -1,
+    "Max": 523,
+    "Gamma": 1,
+    "Min": 352,
+    "Color": -1
+  }],
+  "Comments": {"Summary": "5-GFP-mScarlet/GFP/OD2.0/100x/1x/2x2binning/15s_interval/{-2,2,0/5zstack}/autofocus/mScarlet_GFP-mSc-filter/500ms/100% OR: GFP 20%/exp150ms"}
+}

ndbioimage/tests/test_slicing.py

@@ -26,8 +26,10 @@ def image(array):
 
 @pytest.mark.parametrize(
     "s",
-    combinations_with_replacement(
-        (0, -1, 1, slice(None), slice(0, 1), slice(-1, 0), slice(1, 1)), 5
+    list(
+        combinations_with_replacement(
+            (0, -1, 1, slice(None), slice(0, 1), slice(-1, 0), slice(1, 1), slice(1, 2)), 5
+        )
     ),
 )
 def test_slicing(s, image, array):

ndbioimage/tests/test_ufuncs.py

@@ -26,24 +26,26 @@ def image(array):
 
 @pytest.mark.parametrize(
     "fun_and_axis",
-    product(
-        (
-            np.sum,
-            np.nansum,
-            np.min,
-            np.nanmin,
-            np.max,
-            np.nanmax,
-            np.argmin,
-            np.argmax,
-            np.mean,
-            np.nanmean,
-            np.var,
-            np.nanvar,
-            np.std,
-            np.nanstd,
-        ),
-        (None, 0, 1, 2, 3, 4),
+    list(
+        product(
+            (
+                np.sum,
+                np.nansum,
+                np.min,
+                np.nanmin,
+                np.max,
+                np.nanmax,
+                np.argmin,
+                np.argmax,
+                np.mean,
+                np.nanmean,
+                np.var,
+                np.nanvar,
+                np.std,
+                np.nanstd,
+            ),
+            (None, 0, 1, 2, 3, 4),
+        )
     ),
 )
 def test_ufuncs(fun_and_axis, image, array):

py/ndbioimage/__init__.pyi

@@ -1,3 +0,0 @@
-from .ndbioimage_rs import Imread, Shape
-
-__all__ = ["Imread", "Shape"]

src/readers/tiff.rs

@@ -1,497 +0,0 @@
-use crate::error::Error;
-use crate::readers::{ArrayT, DynReader, Frame, PixelType, Reader, Shape};
-use ndarray::Array2;
-use ome_metadata::{Ome, ome};
-use serde::{Deserialize, Serialize};
-use std::cell::{RefCell, RefMut};
-use std::collections::{HashMap, HashSet};
-use std::fs::File;
-use std::hash::{Hash, Hasher};
-use std::path::{Path, PathBuf};
-use thread_local::ThreadLocal;
-use tiff::decoder::ifd::Value;
-use tiff::decoder::{Decoder, DecodingResult};
-use tiff::tags::Tag;
-
-/// reader for single tiff image files
-#[derive(Debug, Serialize, Deserialize)]
-pub struct TiffReader {
-    #[serde(skip)]
-    reader: ThreadLocal<RefCell<Decoder<File>>>,
-    path: PathBuf,
-    series: usize,
-    position: usize,
-    shape: Shape,
-    pixel_type: PixelType,
-    n_samples: usize,
-    p_ndim: u8,
-    planar: bool,
-}
-
-impl From<TiffReader> for DynReader {
-    fn from(value: TiffReader) -> Self {
-        DynReader::Tiff(value)
-    }
-}
-
-impl Hash for TiffReader {
-    fn hash<H: Hasher>(&self, state: &mut H) {
-        self.path.hash(state);
-        self.series.hash(state);
-        self.position.hash(state);
-    }
-}
-
-impl Clone for TiffReader {
-    fn clone(&self) -> Self {
-        TiffReader {
-            reader: ThreadLocal::default(),
-            path: self.path.clone(),
-            series: self.series,
-            position: self.position,
-            shape: self.shape.clone(),
-            pixel_type: self.pixel_type,
-            n_samples: self.n_samples,
-            p_ndim: self.p_ndim,
-            planar: self.planar,
-        }
-    }
-}
-
-impl PartialEq for TiffReader {
-    fn eq(&self, other: &Self) -> bool {
-        self.path == other.path && self.series == other.series && self.position == other.position
-    }
-}
-
-impl Eq for TiffReader {}
-
-impl TiffReader {
-    fn get_reader(&self) -> Result<RefMut<'_, Decoder<File>>, Error> {
-        self.reader
-            .get_or_try(|| {
-                let file = File::open(&self.path)?;
-                Ok(RefCell::new(Decoder::new(file)?))
-            })
-            .map(|i| i.borrow_mut())
-    }
-}
-
-fn pixel_type_from_bits_sample(bits: u16, fmt: u16) -> Result<PixelType, Error> {
-    match (fmt, bits) {
-        (1, 8) => Ok(PixelType::U8),
-        (1, 16) => Ok(PixelType::U16),
-        (1, 32) => Ok(PixelType::U32),
-        (2, 8) => Ok(PixelType::I8),
-        (2, 16) => Ok(PixelType::I16),
-        (2, 32) => Ok(PixelType::I32),
-        (3, 32) => Ok(PixelType::F32),
-        (3, 64) => Ok(PixelType::F64),
-        _ => Err(Error::Parse(format!(
-            "unsupported pixel type: SampleFormat={}, BitsPerSample={}",
-            fmt, bits
-        ))),
-    }
-}
-
-fn parse_imagej_metadata(s: &str) -> HashMap<String, String> {
-    let mut map = HashMap::new();
-    for line in s.lines() {
-        let line = line.trim();
-        if line.is_empty() || line.starts_with('#') {
-            continue;
-        }
-        if let Some(pos) = line.find('=') {
-            let key = line[..pos].trim().to_string();
-            let val = line[pos + 1..].trim().to_string();
-            map.insert(key, val);
-        }
-    }
-    map
-}
-
-impl Reader for TiffReader {
-    fn new<P>(path: P, series: usize, position: usize) -> Result<Self, Error>
-    where
-        P: AsRef<Path>,
-    {
-        let mut new = Self {
-            reader: ThreadLocal::default(),
-            path: path.as_ref().to_path_buf(),
-            series,
-            position,
-            shape: Shape::default(),
-            pixel_type: PixelType::U16,
-            n_samples: 0,
-            p_ndim: 0,
-            planar: true,
-        };
-
-        let mut reader = new.get_reader()?;
-        let width = reader.get_tag(Tag::ImageWidth)?.into_u32()? as usize;
-        let height = reader.get_tag(Tag::ImageLength)?.into_u32()? as usize;
-
-        let n_samples = reader
-            .get_tag(Tag::SamplesPerPixel)
-            .ok()
-            .and_then(|t| t.into_u16().ok())
-            .unwrap_or(1) as usize;
-
-        let bits = tag_first_u16(&mut reader, Tag::BitsPerSample).unwrap_or(8);
-        let sfmt = tag_first_u16(&mut reader, Tag::SampleFormat).unwrap_or(1);
-        let pixtype = pixel_type_from_bits_sample(bits, sfmt)?;
-
-        let planar = reader
-            .get_tag(Tag::PlanarConfiguration)
-            .ok()
-            .and_then(|t| t.into_u16().ok())
-            .unwrap_or(1)
-            == 2;
-
-        let imagej_bytes = reader
-            .get_tag(Tag::Unknown(50839))
-            .or_else(|_| reader.get_tag(Tag::ImageDescription))?;
-        let imagej_bytes = match imagej_bytes {
-            Value::Ascii(s) => s.into_bytes(),
-            other => other.into_u8_vec()?,
-        };
-        let metadata_str = String::from_utf8_lossy(&imagej_bytes);
-        let metadata_map = parse_imagej_metadata(&metadata_str);
-        let p_ndim = if n_samples > 1 { 3u8 } else { 2u8 };
-        let size_c = if p_ndim == 3 {
-            n_samples
-        } else {
-            metadata_map
-                .get("channels")
-                .and_then(|v| v.parse::<usize>().ok())
-                .unwrap_or(1)
-        };
-        let size_z = metadata_map
-            .get("slices")
-            .and_then(|v| v.parse::<usize>().ok())
-            .unwrap_or(1);
-        let size_t = metadata_map
-            .get("frames")
-            .and_then(|v| v.parse::<usize>().ok())
-            .unwrap_or(1);
-
-        let mult = if p_ndim == 3 { n_samples } else { 1 };
-        let expected = size_c * size_z * size_t / mult;
-
-        let mut n_ifds = 1;
-        while reader.more_images() {
-            reader.next_image()?;
-            n_ifds += 1;
-        }
-
-        let size_z = if n_ifds != expected && n_ifds > 0 {
-            n_ifds / (size_c * size_t / mult).max(1)
-        } else {
-            size_z
-        };
-        drop(reader);
-
-        new.shape = Shape {
-            c: size_c,
-            z: size_z,
-            t: size_t,
-            y: height,
-            x: width,
-            ..Default::default()
-        };
-        new.pixel_type = pixtype;
-        new.n_samples = n_samples;
-        new.p_ndim = p_ndim;
-        new.planar = planar;
-
-        // test reading first frame, if error another reader (bioformats) can take over
-        new.get_frame(0, 0, 0)?;
-        Ok(new)
-    }
-
-    fn metadata(&self) -> Result<Ome, Error> {
-        let mut reader = self.get_reader()?;
-        let imagej_bytes = reader
-            .get_tag(Tag::Unknown(50839))
-            .or_else(|_| reader.get_tag(Tag::ImageDescription))?;
-        let imagej_bytes = match imagej_bytes {
-            Value::Ascii(s) => s.into_bytes(),
-            other => other.into_u8_vec()?,
-        };
-        let metadata_str = String::from_utf8_lossy(&imagej_bytes);
-        let metadata_map = parse_imagej_metadata(&metadata_str);
-
-        let interval_t = metadata_map
-            .get("interval")
-            .and_then(|v| v.parse::<f64>().ok());
-        let ru = reader
-            .get_tag(Tag::ResolutionUnit)
-            .ok()
-            .and_then(|i| i.into_u16().ok())
-            .unwrap_or(2);
-        let pixel_size = reader
-            .get_tag(Tag::XResolution)
-            .ok()
-            .and_then(|i| i.into_u64_vec().ok())
-            .and_then(|i| {
-                i.first()
-                    .and_then(|&n| i.get(1).map(|&d| (n as f64) / (d as f64)))
-            })
-            .and_then(|i| {
-                if i <= 0.0 {
-                    None
-                } else if ru == 3 {
-                    Some(1e4 / i)
-                } else if ru == 2 {
-                    Some(25400.0 / i)
-                } else {
-                    Some(1.0 / i)
-                }
-            });
-        let interval_z = metadata_map
-            .get("spacing")
-            .and_then(|v| v.parse::<f64>().ok());
-
-        let ome_pixel_type = match self.pixel_type {
-            PixelType::I8 => ome::PixelType::Int8,
-            PixelType::U8 => ome::PixelType::Uint8,
-            PixelType::I16 => ome::PixelType::Int16,
-            PixelType::U16 => ome::PixelType::Uint16,
-            PixelType::I32 => ome::PixelType::Int32,
-            PixelType::U32 => ome::PixelType::Uint32,
-            PixelType::F32 => ome::PixelType::Float,
-            PixelType::F64 => ome::PixelType::Double,
-            _ => ome::PixelType::Uint16,
-        };
-
-        let mut pixels = ome::Pixels {
-            id: "Pixels:0".to_string(),
-            size_x: self.shape.x as i32,
-            size_y: self.shape.y as i32,
-            size_z: self.shape.z as i32,
-            size_c: self.shape.c as i32,
-            size_t: self.shape.t as i32,
-            dimension_order: ome::PixelsDimensionOrderType::Xyczt,
-            r#type: ome_pixel_type,
-            physical_size_x: pixel_size.map(|v| v as f32),
-            physical_size_x_unit: ome::UnitsLength::um,
-            physical_size_y: pixel_size.map(|v| v as f32),
-            physical_size_y_unit: ome::UnitsLength::um,
-            physical_size_z: interval_z.map(|v| v as f32),
-            physical_size_z_unit: ome::UnitsLength::um,
-            time_increment: interval_t.map(|v| v as f32),
-            time_increment_unit: ome::UnitsTime::s,
-            significant_bits: None,
-            interleaved: None,
-            big_endian: None,
-            channel: Vec::new(),
-            bin_data: Vec::new(),
-            tiff_data: Vec::new(),
-            metadata_only: None,
-            plane: Vec::new(),
-        };
-
-        for c in 0..self.shape.c {
-            pixels.channel.push(ome::Channel {
-                id: format!("Channel:{}", c),
-                name: Some(format!("Channel {}", c)),
-                ..Default::default()
-            });
-        }
-
-        for c in 0..self.shape.c {
-            for z in 0..self.shape.z {
-                for t in 0..self.shape.t {
-                    pixels.plane.push(ome::Plane {
-                        the_c: Some(c as i32),
-                        the_z: Some(z as i32),
-                        the_t: Some(t as i32),
-                        delta_t: interval_t.and_then(|v| {
-                            if v > 0.0 {
-                                Some((t as f64 * v) as f32)
-                            } else {
-                                None
-                            }
-                        }),
-                        ..Default::default()
-                    });
-                }
-            }
-        }
-
-        let mut ome = Ome::default();
-        ome.image.push(ome::Image {
-            id: "Image:0".to_string(),
-            name: self
-                .path
-                .file_name()
-                .map(|n| n.to_string_lossy().to_string()),
-            pixels,
-            instrument_ref: None,
-            objective_settings: None,
-            acquisition_date: None,
-            description: None,
-            experimenter_ref: None,
-            experiment_ref: None,
-            experimenter_group_ref: None,
-            imaging_environment: None,
-            stage_label: None,
-            roi_ref: Vec::new(),
-            microbeam_manipulation_ref: Vec::new(),
-            annotation_ref: Vec::new(),
-        });
-        Ok(ome)
-    }
-
-    fn get_frame(&self, c: usize, z: usize, t: usize) -> Result<Frame, Error> {
-        let (page_idx, offset, stride) = if self.p_ndim == 3 {
-            (t * self.shape.z + z, c, self.n_samples)
-        } else {
-            (c + z * self.shape.c + t * self.shape.c * self.shape.z, 0, 1)
-        };
-
-        let mut reader = self.get_reader()?;
-        reader.seek_to_image(page_idx)?;
-
-        let data: DecodingResult = if self.planar {
-            let mut result = DecodingResult::U16(Vec::new());
-            reader.read_image_to_buffer(&mut result)?;
-            result
-        } else {
-            reader.read_image()?
-        };
-
-        macro_rules! to_frame {
-            ($var:ident, $ty:ty) => {{
-                let v: Vec<$ty> = match data {
-                    DecodingResult::$var(v) => v,
-                    _ => {
-                        return Err(Error::Parse(format!(
-                            "type mismatch: expected {}",
-                            stringify!($var)
-                        )));
-                    }
-                };
-                let yx = self.shape.y * self.shape.x;
-                if self.planar {
-                    let chan: Vec<$ty> = v.into_iter().skip(c * yx).take(yx).collect();
-                    Ok(ArrayT::$var(Array2::from_shape_vec(
-                        (self.shape.y, self.shape.x),
-                        chan,
-                    )?))
-                } else if stride > 1 {
-                    let chan: Vec<$ty> = v
-                        .into_iter()
-                        .skip(offset)
-                        .step_by(stride)
-                        .take(yx)
-                        .collect();
-                    Ok(ArrayT::$var(Array2::from_shape_vec(
-                        (self.shape.y, self.shape.x),
-                        chan,
-                    )?))
-                } else {
-                    Ok(ArrayT::$var(Array2::from_shape_vec(
-                        (self.shape.y, self.shape.x),
-                        v,
-                    )?))
-                }
-            }};
-        }
-
-        match self.pixel_type {
-            PixelType::U8 => to_frame!(U8, u8),
-            PixelType::U16 => to_frame!(U16, u16),
-            PixelType::U32 => to_frame!(U32, u32),
-            PixelType::U64 => to_frame!(U64, u64),
-            PixelType::I8 => to_frame!(I8, i8),
-            PixelType::I16 => to_frame!(I16, i16),
-            PixelType::I32 => to_frame!(I32, i32),
-            PixelType::I64 => to_frame!(I64, i64),
-            PixelType::F32 => to_frame!(F32, f32),
-            PixelType::F64 => to_frame!(F64, f64),
-            _ => Err(Error::NotImplemented(format!(
-                "unsupported TIFF pixel type: {:?}",
-                self.pixel_type
-            ))),
-        }
-    }
-
-    fn path(&self) -> &Path {
-        &self.path
-    }
-    fn series(&self) -> usize {
-        self.series
-    }
-    fn position(&self) -> usize {
-        self.position
-    }
-    fn shape(&self) -> &Shape {
-        &self.shape
-    }
-    fn pixel_type(&self) -> &PixelType {
-        &self.pixel_type
-    }
-
-    fn get_available_positions<P>(_path: P, _series: usize) -> Result<HashSet<usize>, Error>
-    where
-        P: AsRef<Path>,
-    {
-        Ok(HashSet::from([0]))
-    }
-
-    fn get_available_series<P>(_path: P) -> Result<HashSet<usize>, Error>
-    where
-        P: AsRef<Path>,
-    {
-        Ok(HashSet::from([0]))
-    }
-}
-
-fn tag_first_u16(reader: &mut RefMut<Decoder<File>>, tag: Tag) -> Option<u16> {
-    let val = reader.get_tag(tag).ok()?;
-    use tiff::decoder::ifd::Value;
-    match val {
-        Value::Short(v) => Some(v),
-        Value::List(list) => list.first().and_then(|v| match v {
-            Value::Short(v) => Some(*v),
-            _ => None,
-        }),
-        _ => None,
-    }
-}
-
-#[cfg(test)]
-mod tests {
-    use super::*;
-
-    fn open(file: &str) -> Result<TiffReader, Error> {
-        let path = std::env::current_dir()?
-            .join("tests")
-            .join("files")
-            .join(file);
-        TiffReader::new(&path, 0, 0)
-    }
-
-    macro_rules! test_metadata {
-        ($($name:ident: $file:expr $(,)?)*) => {
-            $(
-                #[test]
-                fn $name() -> Result<(), Error> {
-                    let ts = open($file)?;
-                    println!("{}", ts.view().squeeze()?.summary()?);
-                    Ok(())
-                }
-            )*
-        };
-    }
-
-    test_metadata! {
-        // metadata_a: "tiff/20251014_20-Pos_000_000_mask.tif",
-        metadata_b: "tiff/20251014_20-Pos_000_000_max.tif",
-        metadata_c: "tiff/20251014_20-Pos_000_000_loc_results_Cy3.tif",
-        metadata_e: "tiff/1xp53-01-AP1.tiff",
-        metadata_f: "tiff/test.tif",
-        // metadata_g: "tiff/YTL1849A111_2023_05_04__14_46_19_cellnr_1_track.tif",
-    }
-}

src/transforms.rs

@@ -1,316 +0,0 @@
-use crate::axes::Axis;
-use crate::error::Error;
-use crate::readers::{DynReader, Reader};
-use crate::view::View;
-pub use image_registration::transform::Transform;
-use ndarray::{Array2, Dimension, Ix2, Ix3, IxDyn, s};
-use serde::{Deserialize, Serialize};
-use std::collections::{HashMap, HashSet};
-use std::io::Seek;
-use std::path::Path;
-
-#[derive(Clone, Debug, Default, PartialEq, Serialize, Deserialize)]
-struct TransformOld {
-    #[serde(alias = "CenterOfRotationPoint")]
-    center_of_rotation_point: Vec<f64>,
-    #[serde(alias = "Size")]
-    size: Vec<f64>,
-    #[serde(alias = "TransformParameters")]
-    transform_parameters: Vec<f64>,
-    #[serde(alias = "dTransformParameters")]
-    dtransform_parameters: Vec<f64>,
-}
-
-impl TryFrom<HashMap<String, TransformOld>> for Transforms {
-    type Error = Error;
-
-    fn try_from(t: HashMap<String, TransformOld>) -> Result<Self, Error> {
-        Ok(Self {
-            channel: t
-                .into_iter()
-                .map(|(n, i)| {
-                    Ok((
-                        n,
-                        Transform::new_with_center(
-                            i.transform_parameters,
-                            i.center_of_rotation_point,
-                            i.size.into_iter().map(|j| j as usize).collect(),
-                        )
-                        .with_dparameters(i.dtransform_parameters)
-                        .inverse()?,
-                    ))
-                })
-                .collect::<Result<_, Error>>()?,
-            drift: Vec::new(),
-        })
-    }
-}
-
-#[derive(Clone, Debug, Default, Eq, PartialEq, Serialize, Deserialize)]
-pub struct Transforms {
-    #[serde(skip_serializing_if = "HashMap::is_empty", default = "HashMap::new")]
-    pub channel: HashMap<String, Transform<IxDyn>>,
-    #[serde(skip_serializing_if = "Vec::is_empty", default = "Vec::new")]
-    pub drift: Vec<Transform<IxDyn>>,
-}
-
-impl Transforms {
-    pub fn load<P: AsRef<Path>>(path: P) -> Result<Self, Error> {
-        let mut file = std::fs::File::open(&path)?;
-        Ok(
-            if let Ok(t) = yaml_serde::from_reader::<_, Self>(&file)
-                && !(t.channel.is_empty() && t.drift.is_empty())
-            {
-                t
-            } else {
-                file.seek(std::io::SeekFrom::Start(0))?;
-                yaml_serde::from_reader::<_, HashMap<String, TransformOld>>(&file)?.try_into()?
-            },
-        )
-    }
-
-    pub fn save<P: AsRef<Path>>(&self, path: P) -> Result<(), Error> {
-        let file = std::fs::OpenOptions::new()
-            .write(true)
-            .create(true)
-            .truncate(true)
-            .open(path)?;
-        Ok(yaml_serde::to_writer(file, self)?)
-    }
-
-    pub fn calculate_channel_transforms_2d(
-        bead_files: &[&Path],
-        main_channel: usize,
-        default_transform: Option<Transform<Ix2>>,
-    ) -> Result<HashMap<String, Transform<Ix2>>, Error> {
-        let mut transforms = Vec::new();
-        let default_transform = default_transform.unwrap_or_default();
-        for file in bead_files {
-            let view = View::<_, DynReader>::from_path(file)?;
-            transforms.push(view.calculate_channel_transforms_2d(main_channel)?)
-        }
-        let channels = transforms
-            .iter()
-            .flat_map(|t| t.keys())
-            .collect::<HashSet<_>>();
-        let mut average_transforms = HashMap::new();
-        for channel in channels {
-            let (s, n) = transforms
-                .iter()
-                .filter_map(|t| {
-                    if t.contains_key(channel) {
-                        Some((&t[channel] * &default_transform).matrix())
-                    } else {
-                        None
-                    }
-                })
-                .fold((Array2::zeros((3, 3)), 0), |(s, n), t| (s + t, n + 1));
-            let matrix = s / (n as f64);
-            let (s, n) = transforms
-                .iter()
-                .filter_map(|t| {
-                    if t.contains_key(channel) {
-                        Some(((&t[channel] * &default_transform).matrix() - &matrix).powi(2))
-                    } else {
-                        None
-                    }
-                })
-                .fold((Array2::zeros((3, 3)), 0), |(s, n), t| (s + t, n + 1));
-            let dmatrix = (s / n as f64).sqrt();
-            average_transforms.insert(
-                channel.clone(),
-                Transform::default()
-                    .with_matrix(matrix.view())
-                    .with_dmatrix(dmatrix.view()),
-            );
-        }
-        Ok(average_transforms)
-    }
-
-    pub fn calculate_channel_transforms_3d(
-        bead_files: &[&Path],
-        main_channel: usize,
-        default_transform: Option<Transform<Ix3>>,
-    ) -> Result<HashMap<String, Transform<Ix3>>, Error> {
-        let mut transforms = Vec::new();
-        let default_transform = default_transform.unwrap_or_default();
-        for file in bead_files {
-            let view = View::<_, DynReader>::from_path(file)?;
-            transforms.push(view.calculate_channel_transforms_3d(main_channel)?)
-        }
-        let channels = transforms
-            .iter()
-            .flat_map(|t| t.keys())
-            .collect::<HashSet<_>>();
-        let mut average_transforms = HashMap::new();
-        for channel in channels {
-            let (s, n) = transforms
-                .iter()
-                .filter_map(|t| {
-                    if t.contains_key(channel) {
-                        Some((&t[channel] * &default_transform).matrix())
-                    } else {
-                        None
-                    }
-                })
-                .fold((Array2::zeros((4, 4)), 0), |(s, n), t| (s + t, n + 1));
-            let matrix = s / (n as f64);
-            let (s, n) = transforms
-                .iter()
-                .filter_map(|t| {
-                    if t.contains_key(channel) {
-                        Some(((&t[channel] * &default_transform).matrix() - &matrix).powi(2))
-                    } else {
-                        None
-                    }
-                })
-                .fold((Array2::zeros((4, 4)), 0), |(s, n), t| (s + t, n + 1));
-            let dmatrix = (s / n as f64).sqrt();
-            average_transforms.insert(
-                channel.clone(),
-                Transform::default()
-                    .with_matrix(matrix.view())
-                    .with_dmatrix(dmatrix.view()),
-            );
-        }
-        Ok(average_transforms)
-    }
-}
-
-impl<D: Dimension, R: Reader> View<D, R> {
-    #[cfg(feature = "transforms")]
-    pub(crate) fn get_transforms(&self) -> &Transforms {
-        &self.transforms
-    }
-
-    #[cfg(feature = "transforms")]
-    pub(crate) fn set_transforms(&mut self, transforms: Transforms) {
-        self.transforms = transforms;
-    }
-
-    pub fn with_transforms_from_yaml(mut self, path: &Path) -> Result<Self, Error> {
-        self.set_transforms(Transforms::load(path)?);
-        Ok(self)
-    }
-
-    pub fn load_transforms_from_yaml(&mut self, path: &Path) -> Result<(), Error> {
-        self.set_transforms(Transforms::load(path)?);
-        Ok(())
-    }
-
-    pub fn calculate_channel_transforms_2d(
-        &self,
-        main_channel: usize,
-    ) -> Result<HashMap<String, Transform<Ix2>>, Error> {
-        let main_max = self
-            .slice_cztyx(s![main_channel, .., 0, .., ..])?
-            .max_proj(Axis::Z)?
-            .as_array::<f64>()?;
-        let mut transforms = HashMap::new();
-        let shape = vec![self.shape().y, self.shape().x];
-        let center = shape
-            .iter()
-            .map(|&i| (i as f64 - 1.0) / 2.0)
-            .collect::<Vec<_>>();
-        for (channel, metadata) in self.metadata()?.image[0].pixels.channel.iter().enumerate() {
-            let name = metadata
-                .name
-                .clone()
-                .unwrap_or_else(|| format!("channel{}", channel));
-            if channel == main_channel {
-                transforms.insert(
-                    name,
-                    Transform::default()
-                        .with_shape(shape.clone())
-                        .with_center(center.clone()),
-                );
-            } else {
-                let max = self
-                    .slice_cztyx(s![channel, .., 0, .., ..])?
-                    .max_proj(Axis::Z)?
-                    .as_array::<f64>()?;
-                transforms.insert(
-                    name,
-                    Transform::register_affine(main_max.view(), max.view())?,
-                );
-            }
-        }
-        Ok(transforms)
-    }
-
-    pub fn calculate_channel_transforms_3d(
-        &self,
-        main_channel: usize,
-    ) -> Result<HashMap<String, Transform<Ix3>>, Error> {
-        let main_max = self
-            .slice_cztyx(s![main_channel, .., 0, .., ..])?
-            .as_array::<f64>()?;
-        let mut transforms = HashMap::new();
-        let shape = vec![self.shape().z, self.shape().y, self.shape().x];
-        let center = shape
-            .iter()
-            .map(|&i| (i as f64 - 1.0) / 2.0)
-            .collect::<Vec<_>>();
-        for (channel, metadata) in self.metadata()?.image[0].pixels.channel.iter().enumerate() {
-            let name = metadata
-                .name
-                .clone()
-                .unwrap_or_else(|| format!("channel{}", channel));
-            if channel == main_channel {
-                transforms.insert(
-                    name,
-                    Transform::default()
-                        .with_shape(shape.clone())
-                        .with_center(center.clone()),
-                );
-            } else {
-                let max = self
-                    .slice_cztyx(s![channel, .., 0, .., ..])?
-                    .as_array::<f64>()?;
-                transforms.insert(
-                    name,
-                    Transform::register_affine(main_max.view(), max.view())?,
-                );
-            }
-        }
-        Ok(transforms)
-    }
-
-    pub fn calculate_drift_transforms_2d(&self) -> Result<Vec<Transform<Ix2>>, Error> {
-        todo!()
-    }
-
-    pub fn calculate_drift_transforms_3d(&self) -> Result<Vec<Transform<Ix3>>, Error> {
-        todo!()
-    }
-}
-
-#[cfg(test)]
-mod tests {
-    use crate::transforms::Transforms;
-
-    #[test]
-    fn test_transforms() -> Result<(), Box<dyn std::error::Error>> {
-        let _t = Transforms::default();
-
-        Ok(())
-    }
-
-    #[test]
-    fn load_transforms() -> Result<(), Box<dyn std::error::Error>> {
-        let t = Transforms::load("tests/files/transforms.yml")?;
-        println!("t: {:?}", t);
-        assert_eq!(t.channel["TV1"].ndim(), 2);
-        assert_eq!(t.channel["TV2"].ndim(), 2);
-        Ok(())
-    }
-
-    #[test]
-    fn load_old_transforms() -> Result<(), Box<dyn std::error::Error>> {
-        let t = Transforms::load("tests/files/transform_red.yml")?;
-        println!("t: {:?}", t);
-        assert_eq!(t.channel["TV1"].ndim(), 2);
-        assert_eq!(t.channel["TV2"].ndim(), 2);
-        Ok(())
-    }
-}