Tile and Stitching

docs/examples/run_tile_hcs_well.py

@@ -0,0 +1,157 @@
+"""
+Tile + Blend an HCS Well (Multi-FOV)
+=====================================
+
+Create a synthetic HCS plate with 1 well and 4 FOVs arranged in a 2x2
+grid with physical overlap, then composite the FOVs into a single mosaic
+and tile+blend with ``map_tiles``.
+
+This demonstrates the full pipeline:
+FOV compositing (``Well.to_xarray``) → tiling (``Slicer``) → blending (``map_tiles``).
+"""
+
+# %%
+import os
+import warnings
+from tempfile import TemporaryDirectory
+
+import numpy as np
+
+from iohub.ngff import open_ome_zarr
+from iohub.ngff.models import TransformationMeta
+from iohub.tile import Slicer, map_tiles
+
+warnings.filterwarnings("ignore")
+
+# %%
+# Create a synthetic HCS plate
+# ------------------------------
+# 1 well ("A/1") with 4 FOVs in a 2x2 grid.
+# Each FOV is 1t x 1c x 1z x 32y x 32x.
+# FOVs overlap by 8 pixels (~25%) in both Y and X.
+#
+# Layout (pixel coordinates):
+#
+# .. code-block:: text
+#
+#     FOV 0: y=[0,32), x=[0,32)       FOV 1: y=[0,32), x=[24,56)
+#     FOV 2: y=[24,56), x=[0,32)      FOV 3: y=[24,56), x=[24,56)
+#
+#     Mosaic: 56 x 56 pixels (with 8px overlap strips between FOVs)
+
+tmp_dir = TemporaryDirectory()
+plate_path = os.path.join(tmp_dir.name, "plate.zarr")
+
+rng = np.random.default_rng(123)
+
+# Pixel size (um/px) — use 1.0 for clean coordinate alignment
+pixel_size = 1.0
+
+# Grid step: 24 px → 8 px overlap per FOV pair (32 - 24 = 8)
+grid_step = 24
+
+# FOV grid positions: (row_idx, col_idx) → pixel origin (y, x)
+fov_grid = {
+    "000": (0, 0),
+    "001": (0, 1),
+    "010": (1, 0),
+    "011": (1, 1),
+}
+
+with open_ome_zarr(plate_path, layout="hcs", mode="w-", channel_names=["GFP"]) as plate:
+    for fov_name, (row_idx, col_idx) in fov_grid.items():
+        pos = plate.create_position("A", "1", fov_name)
+        data = rng.random((1, 1, 1, 32, 32), dtype=np.float32)
+        pos.create_image("0", data, chunks=(1, 1, 1, 32, 32))
+
+        # Set physical scale and translation so FOVs are placed on a grid
+        y_offset = row_idx * grid_step * pixel_size
+        x_offset = col_idx * grid_step * pixel_size
+        pos.set_transform(
+            "0",
+            [
+                TransformationMeta(
+                    type="scale",
+                    scale=[1.0, 1.0, 1.0, pixel_size, pixel_size],
+                ),
+                TransformationMeta(
+                    type="translation",
+                    translation=[0.0, 0.0, 0.0, y_offset, x_offset],
+                ),
+            ],
+        )
+
+print(f"Created plate at {plate_path}")
+
+# %%
+# Open and composite the well
+# -----------------------------
+# ``Well.to_xarray()`` composites all 4 FOVs into one mosaic.
+
+plate = open_ome_zarr(plate_path, mode="r")
+_, well = next(plate.wells())
+mosaic = well.to_xarray(compositor="mean")
+
+print(f"Mosaic shape: {mosaic.shape}")
+print(f"Mosaic Y range: [{float(mosaic.y[0]):.2f}, {float(mosaic.y[-1]):.2f}] um")
+print(f"Mosaic X range: [{float(mosaic.x[0]):.2f}, {float(mosaic.x[-1]):.2f}] um")
+
+# %%
+# Inspect the tiling
+# --------------------
+
+slicer = Slicer(mosaic, tile_size={"y": 24, "x": 24}, overlap={"y": 4, "x": 4})
+print(f"\n{slicer}")
+print(f"Tiles: {len(slicer)}")
+print(f"Overlap edges: {slicer.graph.number_of_edges()}")
+
+# %%
+# Tile, process, and blend
+# --------------------------
+# Apply a function to each tile of the mosaic and blend back.
+
+
+def process(tile):
+    """Example: double the intensity."""
+    return tile * 2
+
+
+result = map_tiles(
+    mosaic,
+    fn=process,
+    tile_size={"y": 24, "x": 24},
+    overlap={"y": 4, "x": 4},
+    weights="gaussian",
+)
+print(f"\nResult shape: {result.shape}")
+print(f"Lazy: {hasattr(result.data, 'dask')}")
+
+# %%
+# Verify the result
+# -------------------
+
+values = result.values
+expected = mosaic.values * 2
+np.testing.assert_allclose(values, expected, atol=1e-4)
+print("Round-trip check: PASSED")
+
+# %%
+# With overlap caching
+# ----------------------
+
+result_cached = map_tiles(
+    mosaic,
+    fn=process,
+    tile_size={"y": 24, "x": 24},
+    overlap={"y": 4, "x": 4},
+    weights="gaussian",
+    cache="persist",
+)
+np.testing.assert_allclose(result_cached.values, expected, atol=1e-4)
+print("Cached round-trip: PASSED")
+
+# %%
+# Clean up
+
+plate.close()
+tmp_dir.cleanup()

docs/examples/run_tile_single_fov.py

@@ -0,0 +1,142 @@
+"""
+Tile + Blend a Single FOV
+==========================
+
+Create a synthetic OME-Zarr FOV, then tile it with overlap,
+apply a processing function to each tile, and blend the results
+back into a single mosaic using ``map_tiles`` (xarray-native)
+and ``tile_and_assemble`` (zarr output).
+"""
+
+# %%
+import os
+import warnings
+from tempfile import TemporaryDirectory
+
+import numpy as np
+
+from iohub.ngff import open_ome_zarr
+from iohub.tile import Slicer, map_tiles, tile_and_assemble
+
+warnings.filterwarnings("ignore")
+
+# %%
+# Create a synthetic single-FOV OME-Zarr
+# ----------------------------------------
+# 1 timepoint, 2 channels, 4 Z-slices, 64x128 YX.
+
+tmp_dir = TemporaryDirectory()
+fov_path = os.path.join(tmp_dir.name, "fov.zarr")
+
+rng = np.random.default_rng(42)
+raw = rng.random((1, 2, 4, 64, 128), dtype=np.float32)
+
+with open_ome_zarr(fov_path, layout="fov", mode="w-", channel_names=["GFP", "DAPI"]) as dataset:
+    dataset.create_image("0", raw, chunks=(1, 1, 4, 64, 128))
+    dataset.set_scale("0", "y", 0.325)
+    dataset.set_scale("0", "x", 0.325)
+
+print(f"Created FOV at {fov_path}")
+
+# %%
+# Open and inspect the data
+# --------------------------
+
+pos = open_ome_zarr(fov_path, mode="r")
+data = pos.to_xarray()
+print(f"Shape: {data.shape}  dims: {data.dims}")
+print(f"Y range: [{float(data.y[0]):.2f}, {float(data.y[-1]):.2f}] um")
+print(f"X range: [{float(data.x[0]):.2f}, {float(data.x[-1]):.2f}] um")
+
+# %%
+# Inspect the Slicer
+# --------------------
+# See how tiles are laid out with overlap.
+
+slicer = Slicer(data, tile_size={"y": 32, "x": 64}, overlap={"y": 8, "x": 16})
+print(slicer)
+print(f"Neighborhood graph: {slicer.graph.number_of_edges()} overlap edges")
+
+# %%
+# map_tiles: xarray-native (no zarr output)
+# -------------------------------------------
+# Tile, apply a function, blend back. Result stays lazy until ``.values``.
+
+
+def my_algorithm(tile):
+    """Example: scale by 2 and add 1."""
+    return tile * 2 + 1
+
+
+result = map_tiles(
+    data,
+    fn=my_algorithm,
+    tile_size={"y": 32, "x": 64},
+    overlap={"y": 8, "x": 16},
+    weights="gaussian",
+)
+print(f"Result shape: {result.shape}, lazy: {hasattr(result.data, 'dask')}")
+print(f"Coords preserved: c={list(result.c.values)}")
+
+# Trigger computation and verify
+values = result.values
+expected = raw * 2 + 1
+np.testing.assert_allclose(values, expected, atol=1e-4)
+print("Round-trip check: PASSED")
+
+# %%
+# map_tiles with overlap caching
+# --------------------------------
+# ``cache="persist"`` pre-loads overlap strips so they aren't read twice.
+# ``cache="bfs"`` reorders tile processing for cache locality.
+
+result_cached = map_tiles(
+    data,
+    fn=my_algorithm,
+    tile_size={"y": 32, "x": 64},
+    overlap={"y": 8, "x": 16},
+    weights="gaussian",
+    cache="persist",
+)
+np.testing.assert_allclose(result_cached.values, expected, atol=1e-4)
+print("Cached round-trip: PASSED")
+
+# %%
+# tile_and_assemble: zarr output
+# --------------------------------
+# Same pipeline, but writes to zarr on disk.
+
+out_path = os.path.join(tmp_dir.name, "result.zarr")
+result_zarr = tile_and_assemble(
+    data,
+    fn=my_algorithm,
+    tile_size={"y": 32, "x": 64},
+    output=out_path,
+    overlap={"y": 8, "x": 16},
+    weights="gaussian",
+)
+print(f"Output zarr: {out_path}")
+np.testing.assert_allclose(result_zarr.values, expected, atol=1e-4)
+print("Zarr round-trip: PASSED")
+
+# %%
+# Identity round-trip with different blenders
+# -----------------------------------------------
+# Verify that blending is correct: ``fn=identity`` recovers the original.
+
+for blender in ["uniform", "gaussian", "distance"]:
+    r = map_tiles(
+        data,
+        fn=lambda t: t,
+        tile_size={"y": 32, "x": 64},
+        overlap={"y": 8, "x": 16},
+        weights=blender,
+    )
+    maxerr = float(np.max(np.abs(r.values - raw)))
+    print(f"  {blender:10s} identity max error: {maxerr:.2e}")
+
+# %%
+# Clean up
+
+pos.close()
+tmp_dir.cleanup()