diff --git a/docs/api.rst b/docs/api.rst
index 1e4d3ef..b756957 100644
--- a/docs/api.rst
+++ b/docs/api.rst
@@ -159,37 +159,42 @@ Region-level streaming
 ~~~~~~~~~~~~~~~~~~~~~~~
 
 The encoder-level API above operates on tiles you have already read and
-normalized. To extract dense grids over many slide regions at once, slide2vec
-provides a higher-level streaming primitive,
-``iter_regions_dense`` (from ``slide2vec.runtime.dense_regions``), that wraps
-spacing-aware region reads, padding, and encoding into a single generator:
+normalized. To extract dense grids over the regions of an hs2p ``TilingResult``,
+slide2vec provides a higher-level streaming primitive,
+``iter_regions_dense`` (from ``slide2vec.runtime.dense_regions``), that wraps the
+region reads, padding, and encoding into a single generator:
 
 .. code-block:: python
 
    from slide2vec.runtime.dense_regions import iter_regions_dense
 
+   # ``tiling_result`` is an hs2p TilingResult (from annotation/tissue sampling);
+   # it already resolved spacing -> level (read_level / read_tile_size_px /
+   # requested_tile_size_px) and carries the region coordinates and slide path.
    for grid in iter_regions_dense(
        model=model,
        device=model.device,
-       wsi=wsi,
-       coordinates=[(0, 0), (4096, 0)],
-       requested_spacing_um=0.5,
-       target_size=2048,
+       tiling_result=tiling_result,
+       num_workers=4,
    ):
        print(grid.shape)  # (d, grid_h, grid_w), float32
 
-For each ``(x, y)`` coordinate (a level-0 top-left location), the region is read
-**spacing-aware** at ``requested_spacing_um`` to ``target_size``, run through the
-encoder's normalization-only ``get_dense_transform``, padded on the bottom/right
-up to the encoder's patch multiple, and encoded into a ``(d, grid_h, grid_w)``
-token grid. The ``wsi`` argument is any object exposing
-``read_region_at_spacing(location, requested_spacing_um, size, *, tolerance,
-interpolation)``, so the loop runs offline in tests with a fake reader.
+Reads go through the **same shared batched reader the pooled path uses**
+(``WSIRegionReader``, cuCIM ``read_regions(num_workers=…)``). For each region the
+reader reads ``read_tile_size_px`` at ``read_level`` and, when that differs from
+``requested_tile_size_px``, **area-resizes** to the supervision size reusing
+hs2p's own ``resize_array(..., "area")`` — the identical operation the tiling
+planner assumed — so the read pixels are unchanged. Each region is then run
+through the encoder's normalization-only ``get_dense_transform``, padded on the
+bottom/right up to the encoder's patch multiple, and encoded into a
+``(d, grid_h, grid_w)`` token grid. The low-level backend is opened lazily via
+``slide2vec.data.tile_reader._open_wsi_backend``, so the loop runs offline in
+tests by monkeypatching that seam with a fake backend.
 
 Streaming contract:
 
 - Grids are yielded **one per coordinate, in coordinate order**; an empty
-  ``coordinates`` sequence yields nothing.
+  ``tiling_result`` yields nothing.
 - Regions are read and encoded one ``batch_size`` chunk at a time, so resident
   host memory is bounded by ``batch_size`` rather than by the slide's coordinate
   count — there is no per-slide accumulation.
diff --git a/slide2vec/data/tile_reader.py b/slide2vec/data/tile_reader.py
index 07a8836..9532a4d 100644
--- a/slide2vec/data/tile_reader.py
+++ b/slide2vec/data/tile_reader.py
@@ -9,6 +9,7 @@
 from hs2p import TilingResult
 from hs2p.utils.stderr import run_with_filtered_stderr
 from hs2p.wsi.streaming.plans import build_supertile_index
+from hs2p.wsi.wsi import resize_array
 from slide2vec.utils.log_utils import suppress_c_stderr
 
 
@@ -278,7 +279,17 @@ def _run_batch():
 
 
 class WSIRegionReader:
-    """Random-access region reader for hierarchical extraction."""
+    """Random-access region reader for hierarchical and dense extraction.
+
+    Reads ``region_size_px`` at ``read_level`` (batched via cucim ``read_regions`` when the
+    backend is ``"cucim"``, else serial ``read_region``). When ``resize_to_px`` is set and
+    differs from ``region_size_px``, each region is resized to ``(resize_to_px, resize_to_px)``
+    with ``interpolation`` — reusing hs2p's own :func:`resize_array`. This is the geometry
+    resize the dense path needs (``"area"`` downscale of a coarser-level read to the
+    supervision tile size), applied identically to what ``read_region_at_spacing`` does, and
+    left off (``resize_to_px=None``) for the hierarchical unfold path, which needs the full
+    region.
+    """
 
     def __init__(
         self,
@@ -289,6 +300,8 @@ def __init__(
         backend: str = "cucim",
         num_cucim_workers: int = 4,
         gpu_decode: bool = False,
+        resize_to_px: int | None = None,
+        interpolation: str = "area",
     ):
         self._image_path = str(image_path)
         self._backend = backend
@@ -296,6 +309,11 @@ def __init__(
         self._gpu_decode = gpu_decode
         self._read_level = int(read_level)
         self._region_size_px = int(region_size_px)
+        self._resize_to_px = int(resize_to_px) if resize_to_px is not None else None
+        self._interpolation = interpolation
+        self._out_size_px = (
+            self._resize_to_px if self._resize_to_px is not None else self._region_size_px
+        )
         self._reader = None
 
     def _ensure_open(self) -> None:
@@ -321,13 +339,24 @@ def _read_regions_batch(self, locations: list[tuple[int, int]]) -> list[np.ndarr
             for loc in locations
         ]
 
+    def _prepare_region(self, region: np.ndarray) -> np.ndarray:
+        """RGB-slice a read region and area-resize it to the target size when configured."""
+        arr = np.asarray(region)[:, :, :3]
+        if self._resize_to_px is not None and self._resize_to_px != self._region_size_px:
+            arr = resize_array(
+                arr,
+                (self._resize_to_px, self._resize_to_px),
+                interpolation=self._interpolation,
+            )
+        return arr
+
     def read_batch_with_timing(
         self,
         locations: list[tuple[int, int]],
     ) -> tuple[torch.Tensor, dict[str, float]]:
         if not locations:
             return (
-                torch.empty((0, 3, self._region_size_px, self._region_size_px), dtype=torch.uint8),
+                torch.empty((0, 3, self._out_size_px, self._out_size_px), dtype=torch.uint8),
                 {"reader_open_ms": 0.0, "reader_read_ms": 0.0},
             )
         stderr_context = suppress_c_stderr() if self._backend == "cucim" else nullcontext()
@@ -339,7 +368,7 @@ def read_batch_with_timing(
             read_start = time.perf_counter()
             regions = self._read_regions_batch(locations)
         reader_read_ms = (time.perf_counter() - read_start) * 1000.0
-        batch = np.stack([np.asarray(region)[:, :, :3] for region in regions], axis=0)
+        batch = np.stack([self._prepare_region(region) for region in regions], axis=0)
         tensor = torch.from_numpy(batch).permute(0, 3, 1, 2).contiguous()
         return tensor, {"reader_open_ms": reader_open_ms, "reader_read_ms": reader_read_ms}
 
diff --git a/slide2vec/inference.py b/slide2vec/inference.py
index 36af5b9..9a14b5b 100644
--- a/slide2vec/inference.py
+++ b/slide2vec/inference.py
@@ -600,7 +600,7 @@ def embed_tiles(
                     tiling_result=tiling_result,
                     image_path=slide.image_path,
                     mask_path=slide.mask_path,
-                    backend=tiling.resolve_slide_backend(resolved_preprocessing, tiling_result),
+                    backend=tiling.resolve_slide_backend(resolved_preprocessing.backend, tiling_result),
                     preprocessing=resolved_preprocessing,
                 ),
                 annotation=embedding.tiling_result_annotation(tiling_result),
@@ -620,7 +620,7 @@ def embed_tiles(
                 image_path=slide.image_path,
                 mask_path=slide.mask_path,
                 tile_size_lv0=int(tiling_result.tile_size_lv0),
-                backend=tiling.resolve_slide_backend(resolved_preprocessing, tiling_result),
+                backend=tiling.resolve_slide_backend(resolved_preprocessing.backend, tiling_result),
             )
             artifact = embedding.write_tile_embedding_artifact(
                 slide.sample_id,
diff --git a/slide2vec/runtime/cpu_budget.py b/slide2vec/runtime/cpu_budget.py
index 5688f2f..5407303 100644
--- a/slide2vec/runtime/cpu_budget.py
+++ b/slide2vec/runtime/cpu_budget.py
@@ -50,7 +50,7 @@ def log_on_the_fly_worker_override_once(
     if not preprocessing.on_the_fly or preprocessing.read_tiles_from is not None:
         return
     if not any(
-        resolve_slide_backend(preprocessing, tiling_result) == "cucim"
+        resolve_slide_backend(preprocessing.backend, tiling_result) == "cucim"
         for tiling_result in tiling_results
     ):
         return
diff --git a/slide2vec/runtime/dense_regions.py b/slide2vec/runtime/dense_regions.py
index f413297..3cecdcb 100644
--- a/slide2vec/runtime/dense_regions.py
+++ b/slide2vec/runtime/dense_regions.py
@@ -1,25 +1,28 @@
-"""Dense ``(d, h, w)`` grid extraction over **slide regions at coordinates**.
+"""Dense ``(d, h, w)`` grid extraction over the **regions of an hs2p ``TilingResult``**.
 
 The dense counterpart of the pooled coordinate path (``compute_tile_embeddings_for_slide``
 → ``run_forward_pass`` → ``encode_tiles``): instead of pooling each region to one vector,
-each sampled ROI is read **spacing-aware** from the slide, run through the encoder's
-normalization-only dense transform (``get_dense_transform`` — NOT the pooled transform,
-which crops), padded up to the encoder's patch multiple, and encoded via
-``encode_tiles_dense`` into a ``(d, grid_h, grid_w)`` token grid. ``iter_regions_dense``
-**streams** these grids — yielding one per coordinate, in coordinate order, holding at most
-one ``batch_size`` chunk resident — so host memory is bounded by ``batch_size`` rather than
-by a slide's ROI count.
+each sampled ROI is read from the slide, run through the encoder's normalization-only dense
+transform (``get_dense_transform`` — NOT the pooled transform, which crops), padded up to the
+encoder's patch multiple, and encoded via ``encode_tiles_dense`` into a ``(d, grid_h, grid_w)``
+token grid. ``iter_regions_dense`` **streams** these grids — yielding one per coordinate, in
+coordinate order, holding at most one ``batch_size`` chunk resident — so host memory is bounded
+by ``batch_size`` rather than by a slide's ROI count.
 
 This is the extraction half of soma's slide-manifest segmentation path: slide2vec reads
-regions + encodes (it already owns the region reader and the dense encode); soma sources
-the ROI coordinates (hs2p annotation sampling) and persists/caches the grids. It mirrors
+regions + encodes (it already owns the region reader and the dense encode); soma sources the
+ROIs (hs2p annotation sampling → a ``TilingResult``) and persists/caches the grids. It mirrors
 the pooled split exactly — extraction here, caching in soma.
 
-Region reads are spacing-aware via hs2p (:meth:`hs2p.wsi.wsi.WSI.read_region_at_spacing`):
-the finest pyramid level ``<=`` the requested µm/px is read and downscaled to the exact
-``target_size`` (``area`` for images), so the token grid registers against a mask read at
-the same spacing. The ``wsi`` is injected (any object exposing ``read_region_at_spacing``),
-so the loop is unit-testable offline with a fake reader + a random-weight encoder.
+Reads go through the **same shared batched reader the pooled path uses**
+(:class:`slide2vec.data.tile_reader.WSIRegionReader`, cuCIM ``read_regions(num_workers=…)``),
+driven by the ``TilingResult`` that already resolved the spacing→level plan. The reader reads
+``read_tile_size_px`` at ``read_level`` and, when ``read_tile_size_px != requested_tile_size_px``,
+area-resizes to ``requested_tile_size_px`` reusing hs2p's own ``resize_array(..., "area")`` —
+the identical operation the legacy ``read_region_at_spacing`` performed, so the read pixels are
+unchanged. The token grid registers against a mask read at the same spacing. The offline seam
+is :func:`slide2vec.data.tile_reader._open_wsi_backend` (monkeypatch it with a fake backend to
+run the loop with a random-weight encoder and canned regions).
 
 Both dense modes run through one primitive (:func:`~slide2vec.runtime.dense_sliding.encode_dense_sliding`):
 ``window_size=None`` is a single whole-tile forward (byte-identical to the legacy
@@ -32,16 +35,18 @@
 from __future__ import annotations
 
 from dataclasses import dataclass
-from typing import Callable, Iterator, Sequence
+from typing import Callable, Iterator
 
 import numpy as np
 import torch
 import torch.nn.functional as F
 from PIL import Image
 
+from slide2vec.data.tile_reader import WSIRegionReader
 from slide2vec.runtime.dense_sliding import encode_dense_sliding
 from slide2vec.runtime.model_settings import output_torch_dtype, resolve_output_precision
 from slide2vec.runtime.slide_encode import slide_encode_autocast_ctx
+from slide2vec.runtime.tiling import resolve_slide_backend
 
 
 def _resolve_output_dtype(output_dtype: "torch.dtype | None", precision: str) -> "torch.dtype":
@@ -169,11 +174,9 @@ def iter_regions_dense(
     *,
     model,
     device: torch.device | str,
-    wsi,
-    coordinates: Sequence[tuple[int, int]],
-    requested_spacing_um: float,
-    target_size: int | tuple[int, int],
-    tolerance: float = 0.05,
+    tiling_result,
+    backend: str | None = None,
+    num_workers: int = 4,
     pad_mode: str = "reflect",
     image_pad_value: float | None = None,
     window_size: int | None = None,
@@ -186,7 +189,18 @@ def iter_regions_dense(
     output_dtype: "torch.dtype | None" = None,
     dense_transform: Callable | None = None,
 ) -> Iterator[np.ndarray]:
-    """Stream slide regions at ``coordinates`` into dense grids, one per coordinate.
+    """Stream the regions of ``tiling_result`` into dense grids, one per coordinate.
+
+    Driven by an hs2p ``TilingResult`` — which already resolved the spacing→level plan
+    (``read_level`` / ``read_tile_size_px`` / ``requested_tile_size_px`` via
+    ``plan_spacing_read`` at tiling time) — so this loop does **not** re-plan per read. It
+    reads through the shared batched :class:`~slide2vec.data.tile_reader.WSIRegionReader`
+    (cuCIM ``read_regions(num_workers=…)``), which reads ``read_tile_size_px`` at
+    ``read_level`` and area-resizes to ``requested_tile_size_px`` when they differ (hs2p
+    ``resize_array(..., "area")`` — the identical op the legacy ``read_region_at_spacing``
+    performed, so the read pixels are unchanged). Each region is then run through the
+    normalization-only dense transform, padded up to the encoder's patch multiple, and
+    encoded.
 
     Yields one ``(d, grid_h, grid_w)`` grid per coordinate, in coordinate order, in the
     model's compute ``precision`` by default (fp16 runs yield fp16 grids; see
@@ -194,22 +208,25 @@ def iter_regions_dense(
     resident host memory is bounded by ``batch_size`` rather than by a slide's ROI count
     (the loop holds at most one batch of grids resident — no per-slide accumulation).
 
-    Injectable core: takes a constructed dense-capable ``model`` (with
-    ``encode_tiles_dense`` / ``encode_tiles_attention`` / ``patch_size`` /
-    ``get_dense_transform``) and a ``wsi`` exposing
-    ``read_region_at_spacing(location, requested_spacing_um, size, *, tolerance,
-    interpolation)``, so it runs offline in tests with random weights + a fake reader.
+    Offline seam: the low-level backend is opened lazily via
+    :func:`slide2vec.data.tile_reader._open_wsi_backend`; monkeypatch it with a fake backend
+    (serving canned region arrays through ``read_regions`` / ``read_region``) to run this loop
+    with a random-weight ``model`` and no real slide.
 
     Arguments are validated and geometry is resolved **eagerly** (before any region is
     read): an invalid ``pad_mode`` or ``feature_kind`` raises at the call site, not on the
     first ``next()``. Iteration itself is lazy — reads advance one batch at a time.
 
     Args:
-        coordinates: ``(x, y)`` top-left locations in **level-0** pixel space (the hs2p
-            tiling convention; passed straight to ``read_region_at_spacing``).
-        requested_spacing_um: µm/px to read each region at.
-        target_size: supervision tile size (int or ``(h, w)``); the region is read at this
-            size at ``requested_spacing_um`` and the token grid registers to it.
+        tiling_result: hs2p ``TilingResult`` carrying ``tiles.x`` / ``tiles.y`` (level-0
+            top-left coordinates), ``read_level``, ``read_tile_size_px``,
+            ``requested_tile_size_px``, ``image_path`` and ``backend``. The supervision tile
+            size is ``requested_tile_size_px``; the token grid registers to it.
+        backend: reader backend override (``"cucim"`` / ``"openslide"`` / ``"vips"`` /
+            ``"asap"``). ``None`` (default) resolves from ``tiling_result.backend`` the same
+            way the pooled path does.
+        num_workers: cuCIM read parallelism (``read_regions(num_workers=…)``); size it via
+            :func:`slide2vec.runtime.cpu_budget.resolve_on_the_fly_num_workers`.
         window_size: encoder field-of-view chunk fed through the backbone per forward.
             ``None`` (default) is one whole-tile forward, byte-identical to the
             whole-region encode; a value smaller than the encoded tile slides the encoder
@@ -223,8 +240,8 @@ def iter_regions_dense(
             bfloat16). Pass e.g. ``torch.float32`` to force a lossless cache regardless of
             precision; an explicit ``torch.bfloat16`` is rejected (cannot cross ``.numpy()``).
 
-    Yields grids in coordinate order in ``output_dtype``; empty ``coordinates`` yields nothing.
-    ``feature_kind`` selects ``encode_tiles_dense`` (patch grid) vs
+    Yields grids in coordinate order in ``output_dtype``; an empty ``tiling_result`` yields
+    nothing. ``feature_kind`` selects ``encode_tiles_dense`` (patch grid) vs
     ``encode_tiles_attention`` (CLS-attention grid); both produce a ``(C, gh, gw)`` grid and
     share this path. Each yielded grid is a standalone contiguous copy, so it does not pin
     the rest of its batch's memory alive.
@@ -232,7 +249,12 @@ def iter_regions_dense(
     if pad_mode not in _PAD_MODES:
         raise ValueError(f"unsupported pad_mode {pad_mode!r}; expected one of {sorted(_PAD_MODES)}")
     resolved_output_dtype = _resolve_output_dtype(output_dtype, precision)
-    geometry = compute_dense_geometry(target_size=target_size, patch_size=model.patch_size)
+    read_level = int(tiling_result.read_level)
+    read_tile_size_px = int(tiling_result.read_tile_size_px)
+    requested_tile_size_px = int(tiling_result.requested_tile_size_px)
+    geometry = compute_dense_geometry(
+        target_size=requested_tile_size_px, patch_size=model.patch_size
+    )
     if dense_transform is None:
         dense_transform = model.get_dense_transform()
     encode_fn = _resolve_encode_fn(
@@ -242,18 +264,27 @@ def iter_regions_dense(
         attention_include_registers=attention_include_registers,
     )
     target_h, target_w = geometry.target_size
-    coords = [(int(x), int(y)) for x, y in coordinates]
+    x = np.asarray(tiling_result.x)
+    y = np.asarray(tiling_result.y)
+    coords = [(int(x[i]), int(y[i])) for i in range(x.shape[0])]
     step = max(1, int(batch_size))
 
-    def _read_padded(location: tuple[int, int]) -> torch.Tensor:
-        region = wsi.read_region_at_spacing(
-            location,
-            float(requested_spacing_um),
-            (target_w, target_h),  # hs2p size is (width, height)
-            tolerance=float(tolerance),
-            interpolation="area",
-        )
-        region = np.ascontiguousarray(np.asarray(region)[..., :3])
+    # Resolve the backend the way the pooled path does: an explicit ``backend`` wins,
+    # otherwise fall back to the TilingResult's own backend ("auto" resolution).
+    requested_backend = "auto" if backend is None else backend
+    resolved_backend = resolve_slide_backend(requested_backend, tiling_result)
+    reader = WSIRegionReader(
+        tiling_result.image_path,
+        read_level=read_level,
+        region_size_px=read_tile_size_px,
+        backend=resolved_backend,
+        num_cucim_workers=int(num_workers),
+        resize_to_px=requested_tile_size_px,  # area-resizes iff read != requested
+        interpolation="area",
+    )
+
+    def _transform_and_pad(region_hwc: np.ndarray, location: tuple[int, int]) -> torch.Tensor:
+        region = np.ascontiguousarray(np.asarray(region_hwc)[..., :3])
         tensor = torch.as_tensor(dense_transform(Image.fromarray(region))).as_subclass(torch.Tensor)
         if tensor.ndim != 3:
             raise ValueError(
@@ -262,8 +293,8 @@ def _read_padded(location: tuple[int, int]) -> torch.Tensor:
         if tuple(int(s) for s in tensor.shape[-2:]) != (target_h, target_w):
             raise ValueError(
                 f"region at {location} is {tuple(int(s) for s in tensor.shape[-2:])} after the dense "
-                f"transform, but target_size is {(target_h, target_w)}. The dense transform must be "
-                "normalization-only (no resize/crop)."
+                f"transform, but requested_tile_size_px is {(target_h, target_w)}. The dense transform "
+                "must be normalization-only (no resize/crop)."
             )
         return pad_image_to_encoded(
             tensor, geometry, pad_mode=pad_mode, image_pad_value=image_pad_value
@@ -273,9 +304,14 @@ def _stream() -> Iterator[np.ndarray]:
         with torch.inference_mode(), slide_encode_autocast_ctx(device, precision):
             for start in range(0, len(coords), step):
                 chunk = coords[start : start + step]
-                batch = torch.stack([_read_padded(loc) for loc in chunk]).to(
-                    device, non_blocking=True
-                )
+                # One batched read per chunk (cuCIM read_regions when cucim; serial read
+                # otherwise) through the shared reader, already RGB-sliced and area-resized
+                # to requested_tile_size_px. (B, 3, target_h, target_w) uint8.
+                region_batch, _timing = reader.read_batch_with_timing(chunk)
+                region_np = region_batch.permute(0, 2, 3, 1).numpy()  # (B, H, W, 3) uint8
+                batch = torch.stack(
+                    [_transform_and_pad(region_np[i], chunk[i]) for i in range(region_np.shape[0])]
+                ).to(device, non_blocking=True)
                 # Every batch goes through the one windowed primitive: window_size=None
                 # short-circuits to a single whole-tile forward (byte-identical to the
                 # whole-region encode), so there is no separate whole-region branch.
diff --git a/slide2vec/runtime/embedding_persist.py b/slide2vec/runtime/embedding_persist.py
index f1b8305..f04c042 100644
--- a/slide2vec/runtime/embedding_persist.py
+++ b/slide2vec/runtime/embedding_persist.py
@@ -89,7 +89,7 @@ def persist_embedded_slide(
                 image_path=embedded_slide.image_path,
                 mask_path=embedded_slide.mask_path,
                 tile_size_lv0=embedded_slide.tile_size_lv0,
-                backend=resolve_slide_backend(preprocessing, tiling_result),
+                backend=resolve_slide_backend(preprocessing.backend, tiling_result),
             ),
         )
         return None, None
@@ -103,7 +103,7 @@ def persist_embedded_slide(
                 tiling_result=tiling_result,
                 image_path=embedded_slide.image_path,
                 mask_path=embedded_slide.mask_path,
-                backend=resolve_slide_backend(preprocessing, tiling_result),
+                backend=resolve_slide_backend(preprocessing.backend, tiling_result),
                 preprocessing=preprocessing,
             ),
             annotation=annotation,
@@ -122,7 +122,7 @@ def persist_embedded_slide(
                 image_path=embedded_slide.image_path,
                 mask_path=embedded_slide.mask_path,
                 tile_size_lv0=embedded_slide.tile_size_lv0,
-                backend=resolve_slide_backend(preprocessing, tiling_result),
+                backend=resolve_slide_backend(preprocessing.backend, tiling_result),
             ),
         )
     slide_artifact = None
diff --git a/slide2vec/runtime/embedding_pipeline.py b/slide2vec/runtime/embedding_pipeline.py
index e47d35e..1861e0d 100644
--- a/slide2vec/runtime/embedding_pipeline.py
+++ b/slide2vec/runtime/embedding_pipeline.py
@@ -95,7 +95,7 @@ def compute_tile_embeddings_for_slide(
             return torch.empty((0, int(feature_dim)), dtype=torch.float32)
     supertile_reorder = None
     if preprocessing.on_the_fly and preprocessing.read_tiles_from is None:
-        resolved_backend = resolve_slide_backend(preprocessing, tiling_result)
+        resolved_backend = resolve_slide_backend(preprocessing.backend, tiling_result)
         collate_fn = OnTheFlyBatchTileCollator(
             image_path=slide.image_path,
             tiling_result=tiling_result,
@@ -139,7 +139,7 @@ def compute_tile_embeddings_for_slide(
     dataset = TileIndexDataset(resolved_indices)
     batch_preprocessor = build_batch_preprocessor(loaded, tiling_result)
     loader_kwargs = embedding_dataloader_kwargs(loaded, execution)
-    resolved_backend = resolve_slide_backend(preprocessing, tiling_result)
+    resolved_backend = resolve_slide_backend(preprocessing.backend, tiling_result)
     if preprocessing.on_the_fly and preprocessing.read_tiles_from is None and resolved_backend == "cucim":
         effective_num_workers, _ = resolve_on_the_fly_num_workers(
             preprocessing.num_cucim_workers,
@@ -211,7 +211,7 @@ def compute_hierarchical_embeddings_for_slide(
         subtile_index_within_region=index.subtile_index_within_region,
         read_region_size_px=int(geometry["read_region_size_px"]),
         read_tile_size_px=int(geometry["read_tile_size_px"]),
-        backend=resolve_slide_backend(preprocessing, tiling_result),
+        backend=resolve_slide_backend(preprocessing.backend, tiling_result),
         num_cucim_workers=preprocessing.num_cucim_workers,
         gpu_decode=preprocessing.gpu_decode,
     )
@@ -221,7 +221,7 @@ def compute_hierarchical_embeddings_for_slide(
         requested_tile_size_px=int(geometry["requested_tile_size_px"]),
     )
     loader_kwargs = embedding_dataloader_kwargs(loaded, execution)
-    resolved_backend = resolve_slide_backend(preprocessing, tiling_result)
+    resolved_backend = resolve_slide_backend(preprocessing.backend, tiling_result)
     if resolved_backend == "cucim":
         effective_num_workers, _ = resolve_on_the_fly_num_workers(
             preprocessing.num_cucim_workers,
@@ -292,7 +292,7 @@ def compute_hierarchical_embedding_shard_for_slide(
         subtile_index_within_region=index.subtile_index_within_region,
         read_region_size_px=int(geometry["read_region_size_px"]),
         read_tile_size_px=int(geometry["read_tile_size_px"]),
-        backend=resolve_slide_backend(preprocessing, tiling_result),
+        backend=resolve_slide_backend(preprocessing.backend, tiling_result),
         num_cucim_workers=preprocessing.num_cucim_workers,
         gpu_decode=preprocessing.gpu_decode,
     )
@@ -302,7 +302,7 @@ def compute_hierarchical_embedding_shard_for_slide(
         requested_tile_size_px=int(geometry["requested_tile_size_px"]),
     )
     loader_kwargs = embedding_dataloader_kwargs(loaded, execution)
-    resolved_backend = resolve_slide_backend(preprocessing, tiling_result)
+    resolved_backend = resolve_slide_backend(preprocessing.backend, tiling_result)
     if resolved_backend == "cucim":
         effective_num_workers, _ = resolve_on_the_fly_num_workers(
             preprocessing.num_cucim_workers,
diff --git a/slide2vec/runtime/patient_pipeline.py b/slide2vec/runtime/patient_pipeline.py
index 46a8823..63eb252 100644
--- a/slide2vec/runtime/patient_pipeline.py
+++ b/slide2vec/runtime/patient_pipeline.py
@@ -77,7 +77,7 @@ def run_patient_pipeline(
                     image_path=slide.image_path,
                     mask_path=slide.mask_path,
                     tile_size_lv0=int(tiling_result.tile_size_lv0),
-                    backend=resolve_slide_backend(preprocessing, tiling_result),
+                    backend=resolve_slide_backend(preprocessing.backend, tiling_result),
                 ),
             )
             tile_artifacts.append(tile_artifact)
diff --git a/slide2vec/runtime/process_list.py b/slide2vec/runtime/process_list.py
index d3f5d4c..44986b5 100644
--- a/slide2vec/runtime/process_list.py
+++ b/slide2vec/runtime/process_list.py
@@ -71,7 +71,7 @@ def write_zero_tile_embedding_sidecars(
                     tiling_result=tiling_result,
                     image_path=slide.image_path,
                     mask_path=slide.mask_path,
-                    backend=resolve_slide_backend(preprocessing, tiling_result),
+                    backend=resolve_slide_backend(preprocessing.backend, tiling_result),
                     preprocessing=preprocessing,
                 ),
                 annotation=tiling_result_annotation(tiling_result),
@@ -90,7 +90,7 @@ def write_zero_tile_embedding_sidecars(
                 image_path=slide.image_path,
                 mask_path=slide.mask_path,
                 tile_size_lv0=int(tiling_result.tile_size_lv0),
-                backend=resolve_slide_backend(preprocessing, tiling_result),
+                backend=resolve_slide_backend(preprocessing.backend, tiling_result),
             ),
         )
 
diff --git a/slide2vec/runtime/tiling.py b/slide2vec/runtime/tiling.py
index 6d93020..cf316be 100644
--- a/slide2vec/runtime/tiling.py
+++ b/slide2vec/runtime/tiling.py
@@ -18,10 +18,9 @@ def resolve_tiling_backend(preprocessing: PreprocessingConfig | None) -> str:
     return preprocessing.backend
 
 
-def resolve_slide_backend(preprocessing: PreprocessingConfig | None, tiling_result) -> str:
-    backend = resolve_tiling_backend(preprocessing)
-    if backend != "auto":
-        return backend
+def resolve_slide_backend(requested_backend: str, tiling_result) -> str:
+    if requested_backend != "auto":
+        return requested_backend
     resolved_backend = tiling_result.backend if hasattr(tiling_result, "backend") else None
     if isinstance(resolved_backend, str) and resolved_backend and resolved_backend != "auto":
         return resolved_backend
diff --git a/tests/test_dense_regions.py b/tests/test_dense_regions.py
index a5438b6..84a7d58 100644
--- a/tests/test_dense_regions.py
+++ b/tests/test_dense_regions.py
@@ -1,11 +1,19 @@
 """Tests for dense grid extraction over slide regions: ``iter_regions_dense``.
 
-Fully offline (``pretrained=False`` random weights) + an injected fake reader, so no
-weights, no real WSI. ``iter_regions_dense`` is a streaming generator: it yields one
-``(d, grid_h, grid_w)`` grid per coordinate in coordinate order, holding at most one batch
-resident. Checks (1) grid shapes over a batch of coordinates, (2) that each yielded grid is
-byte-identical to a direct ``transform → pad → encode`` of the same region (both feature
-kinds), (3) streaming/laziness via a call-counting reader, and (4) eager validation.
+Fully offline (``pretrained=False`` random weights) with the low-level WSI backend faked,
+so no weights and no real slide. ``iter_regions_dense`` is driven by an hs2p ``TilingResult``
+(the tiling planner already resolved spacing→level) and reads through the shared batched
+``WSIRegionReader``. The offline seam is ``slide2vec.data.tile_reader._open_wsi_backend`` —
+monkeypatched to return a fake backend serving canned region arrays (``read_regions`` for the
+cucim batched path, ``read_region`` for the serial path).
+
+``iter_regions_dense`` is a streaming generator: it yields one ``(d, grid_h, grid_w)`` grid
+per coordinate in coordinate order, holding at most one ``batch_size`` chunk resident. Checks
+(1) grid shapes / coordinate order, (2) byte-identity to a direct ``transform → pad → encode``
+of the same region (both feature kinds, whole + sliding window), (3) streaming/laziness via a
+read-counting fake, (4) eager validation before any read, (5) the area-resize path when
+``read_tile_size_px != requested_tile_size_px``, and (6) batch-invariance (composition
+irrelevant; only ``B`` matters — see docs/adr/0002).
 """
 
 from __future__ import annotations
@@ -16,6 +24,10 @@
 torch = pytest.importorskip("torch")
 timm = pytest.importorskip("timm")
 
+from hs2p.tiling.result import TileGeometry, TilingResult  # noqa: E402
+from hs2p.wsi.wsi import resize_array  # noqa: E402
+
+from slide2vec.data import tile_reader  # noqa: E402
 from slide2vec.encoders.base import TimmTileEncoder  # noqa: E402
 from slide2vec.runtime.dense_regions import (  # noqa: E402
     _resolve_output_dtype,
@@ -31,38 +43,137 @@ def _encoder(**kwargs) -> TimmTileEncoder:
                            dynamic_img_size=True, **kwargs)
 
 
-class _FakeWSI:
-    """Returns a deterministic RGB region per location (so reads are reproducible)."""
+def _canned_region(location, size) -> np.ndarray:
+    """Deterministic RGB region for a location, of the requested ``(width, height)``."""
+    width, height = int(size[0]), int(size[1])
+    x, y = int(location[0]), int(location[1])
+    rng = np.random.default_rng(abs(hash((x, y))) % (2**32))
+    return rng.integers(0, 256, size=(height, width, 3), dtype=np.uint8)
+
+
+class _FakeBackend:
+    """Serves deterministic region arrays and records every read.
+
+    Implements both the cucim batched ``read_regions`` and the serial ``read_region``
+    the shared reader dispatches to, so a single fake covers both backends.
+    """
+
+    def __init__(self) -> None:
+        self.read_regions_calls: list[list[tuple[int, int]]] = []
+        self.read_region_calls: list[tuple[int, int]] = []
+
+    @property
+    def locations_read(self) -> list[tuple[int, int]]:
+        flat = [loc for batch in self.read_regions_calls for loc in batch]
+        return flat + list(self.read_region_calls)
+
+    def read_regions(self, locations, level, size, num_workers):
+        locs = [(int(x), int(y)) for x, y in locations]
+        self.read_regions_calls.append(locs)
+        return [_canned_region(loc, size) for loc in locs]
 
-    def __init__(self, *, target_h: int, target_w: int):
-        self._target_h = target_h
-        self._target_w = target_w
-        self.calls: list[tuple] = []
+    def read_region(self, location, level, size):
+        loc = (int(location[0]), int(location[1]))
+        self.read_region_calls.append(loc)
+        return _canned_region(loc, size)
 
-    def read_region_at_spacing(self, location, requested_spacing_um, size, *, tolerance, interpolation):
-        self.calls.append((tuple(location), requested_spacing_um, tuple(size), tolerance, interpolation))
-        width, height = size
-        x, y = location
-        rng = np.random.default_rng(abs(hash((int(x), int(y)))) % (2**32))
-        return rng.integers(0, 256, size=(height, width, 3), dtype=np.uint8)
+
+class _FakeBackendFactory:
+    """Stands in for ``_open_wsi_backend``: hands out one shared fake, counts opens."""
+
+    def __init__(self) -> None:
+        self.open_count = 0
+        self.backend = _FakeBackend()
+
+    def __call__(self, image_path, backend, gpu_decode):
+        self.open_count += 1
+        return self.backend
+
+
+@pytest.fixture
+def fake_backend(monkeypatch) -> _FakeBackendFactory:
+    factory = _FakeBackendFactory()
+    monkeypatch.setattr(tile_reader, "_open_wsi_backend", factory)
+    return factory
+
+
+def _make_tiling_result(
+    coords,
+    *,
+    requested_tile_size_px,
+    read_tile_size_px=None,
+    read_level=0,
+    requested_spacing_um=0.5,
+    tolerance=0.05,
+    backend="cucim",
+    image_path="fake.tif",
+) -> TilingResult:
+    """Build a minimal ``TilingResult`` carrying just the fields the dense read needs.
+
+    ``read_tile_size_px`` defaults to ``requested_tile_size_px`` (the no-resize case).
+    """
+    if read_tile_size_px is None:
+        read_tile_size_px = requested_tile_size_px
+    x = np.asarray([c[0] for c in coords], dtype=np.int64)
+    y = np.asarray([c[1] for c in coords], dtype=np.int64)
+    tiles = TileGeometry(
+        x=x,
+        y=y,
+        tissue_fractions=np.ones(len(coords), dtype=np.float32),
+        requested_tile_size_px=int(requested_tile_size_px),
+        requested_spacing_um=float(requested_spacing_um),
+        read_level=int(read_level),
+        read_tile_size_px=int(read_tile_size_px),
+        read_spacing_um=float(requested_spacing_um),
+        tile_size_lv0=int(read_tile_size_px),
+        is_within_tolerance=True,
+        base_spacing_um=float(requested_spacing_um),
+        slide_dimensions=[100000, 100000],
+        level_downsamples=[1.0],
+        overlap=0.0,
+        min_tissue_fraction=0.0,
+    )
+    return TilingResult(
+        tiles=tiles,
+        sample_id="fake",
+        image_path=image_path,
+        backend=backend,
+        requested_backend=backend,
+        tolerance=float(tolerance),
+        step_px_lv0=int(read_tile_size_px),
+        tissue_method="none",
+        requested_seg_downsample=1,
+        seg_downsample=1,
+        seg_level=0,
+        seg_spacing_um=float(requested_spacing_um),
+        seg_sthresh=0,
+        seg_sthresh_up=255,
+        seg_mthresh=0,
+        seg_close=0,
+        ref_tile_size_px=int(requested_tile_size_px),
+        a_t=0.0,
+        a_h=0.0,
+        filter_white=False,
+        filter_black=False,
+        white_threshold=255,
+        black_threshold=0,
+        fraction_threshold=0.0,
+    )
 
 
 @pytest.mark.parametrize("feature_kind", ["patch_features", "cls_attention"])
 @pytest.mark.parametrize("window_size", [None, 32], ids=["whole", "window32"])
-def test_iter_regions_dense_yields_grid_per_coordinate_in_order(window_size, feature_kind):
+def test_iter_regions_dense_yields_grid_per_coordinate_in_order(fake_backend, window_size, feature_kind):
     enc = _encoder()
     target_size = 64  # patch 16 -> grid 4x4, no padding
-    wsi = _FakeWSI(target_h=target_size, target_w=target_size)
     coords = [(0, 0), (64, 0), (0, 64)]
+    result = _make_tiling_result(coords, requested_tile_size_px=target_size)
 
     grids = list(
         iter_regions_dense(
             model=enc,
             device="cpu",
-            wsi=wsi,
-            coordinates=coords,
-            requested_spacing_um=0.5,
-            target_size=target_size,
+            tiling_result=result,
             window_size=window_size,
             feature_kind=feature_kind,
             batch_size=2,
@@ -78,38 +189,47 @@ def test_iter_regions_dense_yields_grid_per_coordinate_in_order(window_size, fea
         assert grid.dtype == np.float32
         assert grid.flags["C_CONTIGUOUS"]
         assert grid.base is None  # standalone copy, not a view pinning a batch
-    # Reads went through read_region_at_spacing at (target_w, target_h), area interp, level-0 coords.
-    assert [c[0] for c in wsi.calls] == [(0, 0), (64, 0), (0, 64)]
-    assert all(c[2] == (target_size, target_size) and c[4] == "area" for c in wsi.calls)
+    # Reads went through the shared batched reader, in coordinate order.
+    assert fake_backend.backend.locations_read == [(0, 0), (64, 0), (0, 64)]
 
 
-def test_iter_regions_dense_pads_non_multiple_target():
+def test_iter_regions_dense_pads_non_multiple_target(fake_backend):
     enc = _encoder()
     target_size = 60  # padded up to 64 -> grid 4x4
-    wsi = _FakeWSI(target_h=target_size, target_w=target_size)
-    grids = list(iter_regions_dense(
-        model=enc, device="cpu", wsi=wsi, coordinates=[(0, 0)],
-        requested_spacing_um=0.5, target_size=target_size,
-    ))
+    result = _make_tiling_result([(0, 0)], requested_tile_size_px=target_size)
+    grids = list(iter_regions_dense(model=enc, device="cpu", tiling_result=result))
     assert len(grids) == 1
     assert grids[0].shape == (enc.encode_dim, 4, 4)
 
 
-def _reference_grid(enc, loc, *, target_size, feature_kind, window_size=None, overlap=0.0):
-    """Hand-rolled transform → pad → encode of one region, for parity checks.
-
-    ``window_size=None`` is the direct whole-tile forward (the byte-identity anchor for
-    the whole-region path); a ``window_size`` routes the padded tile through the same
-    windowed primitive ``iter_regions_dense`` uses, so the seam stays exactly identical.
+def _reference_grid(
+    enc,
+    loc,
+    *,
+    requested_tile_size_px,
+    read_tile_size_px=None,
+    feature_kind,
+    window_size=None,
+    overlap=0.0,
+):
+    """Hand-rolled read → (area-resize) → transform → pad → encode of one region.
+
+    Mirrors ``iter_regions_dense`` exactly: reads the same canned region at
+    ``read_tile_size_px`` and area-resizes to ``requested_tile_size_px`` when they differ
+    (reusing hs2p ``resize_array``), so the pixels are identical.
     """
     from PIL import Image
 
-    geometry = compute_dense_geometry(target_size=target_size, patch_size=enc.patch_size)
+    if read_tile_size_px is None:
+        read_tile_size_px = requested_tile_size_px
+    geometry = compute_dense_geometry(target_size=requested_tile_size_px, patch_size=enc.patch_size)
     transform = enc.get_dense_transform()
-    ref_wsi = _FakeWSI(target_h=target_size, target_w=target_size)
-    region = ref_wsi.read_region_at_spacing(
-        loc, 0.5, (target_size, target_size), tolerance=0.05, interpolation="area"
-    )
+    region = _canned_region(loc, (read_tile_size_px, read_tile_size_px))[:, :, :3]
+    if read_tile_size_px != requested_tile_size_px:
+        region = resize_array(
+            region, (requested_tile_size_px, requested_tile_size_px), interpolation="area"
+        )
+    region = np.ascontiguousarray(region)
     tensor = torch.as_tensor(transform(Image.fromarray(region))).as_subclass(torch.Tensor)
     padded = pad_image_to_encoded(tensor, geometry, pad_mode="reflect", image_pad_value=None)
     batch = padded.unsqueeze(0)
@@ -130,7 +250,7 @@ def _reference_grid(enc, loc, *, target_size, feature_kind, window_size=None, ov
 
 @pytest.mark.parametrize("feature_kind", ["patch_features", "cls_attention"])
 @pytest.mark.parametrize("window_size", [None, 32], ids=["whole", "window32"])
-def test_iter_regions_dense_matches_direct_encode(window_size, feature_kind):
+def test_iter_regions_dense_matches_direct_encode(fake_backend, window_size, feature_kind):
     """Each yielded grid is byte-identical to a hand-rolled transform+pad+encode.
 
     ``window_size=None`` pins the whole-region path against a direct encode; a smaller
@@ -138,88 +258,157 @@ def test_iter_regions_dense_matches_direct_encode(window_size, feature_kind):
     """
     enc = _encoder()
     target_size = 64
-    wsi = _FakeWSI(target_h=target_size, target_w=target_size)
     coords = [(0, 0), (128, 256)]
+    result = _make_tiling_result(coords, requested_tile_size_px=target_size)
 
     grids = list(iter_regions_dense(
-        model=enc, device="cpu", wsi=wsi, coordinates=coords,
-        requested_spacing_um=0.5, target_size=target_size,
+        model=enc, device="cpu", tiling_result=result,
         window_size=window_size, feature_kind=feature_kind,
     ))
 
     assert len(grids) == len(coords)
     for grid, loc in zip(grids, coords):
         ref = _reference_grid(
-            enc, loc, target_size=target_size, feature_kind=feature_kind,
+            enc, loc, requested_tile_size_px=target_size, feature_kind=feature_kind,
             window_size=window_size,
         )
         assert grid.shape == ref.shape
         np.testing.assert_array_equal(grid, ref)
 
 
-def test_iter_regions_dense_empty_coordinates_yields_nothing():
+@pytest.mark.parametrize("feature_kind", ["patch_features", "cls_attention"])
+def test_iter_regions_dense_area_resizes_when_read_differs_from_requested(fake_backend, feature_kind):
+    """When ``read_tile_size_px != requested_tile_size_px`` the region is area-resized.
+
+    The shared reader reads at ``read_tile_size_px`` and area-resizes to
+    ``requested_tile_size_px`` (hs2p ``resize_array``); the grid must match a reference
+    that does the same, and the reader must have requested the *read* size from the slide.
+    """
     enc = _encoder()
-    wsi = _FakeWSI(target_h=64, target_w=64)
+    requested = 64
+    read = 96  # coarser level read, downscaled to the supervision size
+    coords = [(0, 0), (512, 512)]
+    result = _make_tiling_result(
+        coords, requested_tile_size_px=requested, read_tile_size_px=read
+    )
+
     grids = list(iter_regions_dense(
-        model=enc, device="cpu", wsi=wsi, coordinates=[],
-        requested_spacing_um=0.5, target_size=64,
+        model=enc, device="cpu", tiling_result=result, feature_kind=feature_kind,
     ))
+
+    assert len(grids) == len(coords)
+    for grid, loc in zip(grids, coords):
+        ref = _reference_grid(
+            enc, loc, requested_tile_size_px=requested, read_tile_size_px=read,
+            feature_kind=feature_kind,
+        )
+        assert grid.shape == (enc.encode_dim if feature_kind == "patch_features" else grid.shape[0], 4, 4)
+        np.testing.assert_array_equal(grid, ref)
+
+
+def test_iter_regions_dense_serial_backend_reads_region(fake_backend):
+    """A non-cucim backend reads through the serial ``read_region`` path, same output."""
+    enc = _encoder()
+    target_size = 64
+    coords = [(0, 0), (64, 0)]
+    result = _make_tiling_result(coords, requested_tile_size_px=target_size, backend="openslide")
+
+    grids = list(iter_regions_dense(model=enc, device="cpu", tiling_result=result))
+
+    assert len(grids) == 2
+    # The serial path was taken (read_region, not the batched read_regions).
+    assert fake_backend.backend.read_region_calls == [(0, 0), (64, 0)]
+    assert fake_backend.backend.read_regions_calls == []
+    for grid, loc in zip(grids, coords):
+        ref = _reference_grid(enc, loc, requested_tile_size_px=target_size, feature_kind="patch_features")
+        np.testing.assert_array_equal(grid, ref)
+
+
+def test_iter_regions_dense_empty_coordinates_yields_nothing(fake_backend):
+    enc = _encoder()
+    result = _make_tiling_result([], requested_tile_size_px=64)
+    grids = list(iter_regions_dense(model=enc, device="cpu", tiling_result=result))
     assert grids == []
-    assert wsi.calls == []
+    assert fake_backend.backend.locations_read == []
+    assert fake_backend.open_count == 0  # nothing read -> slide never opened
 
 
 @pytest.mark.parametrize("feature_kind", ["patch_features", "cls_attention"])
 @pytest.mark.parametrize("window_size", [None, 32], ids=["whole", "window32"])
-def test_iter_regions_dense_streams_one_batch_at_a_time(window_size, feature_kind):
+def test_iter_regions_dense_streams_one_batch_at_a_time(fake_backend, window_size, feature_kind):
     """Reads advance one batch at a time; first grids land before all coords are read.
 
     The streaming/laziness contract is independent of the dense mode, so it holds for
-    both the whole-tile and sliding-window paths and both feature kinds.
+    both the whole-tile and sliding-window paths and both feature kinds. Reads are counted
+    by total locations pulled from the slide (one batched read per chunk).
     """
     enc = _encoder()
     target_size = 64
-    wsi = _FakeWSI(target_h=target_size, target_w=target_size)
     coords = [(0, 0), (64, 0), (0, 64), (64, 64), (128, 0)]  # 5 coords, batches of [2, 2, 1]
+    result = _make_tiling_result(coords, requested_tile_size_px=target_size)
 
     gen = iter_regions_dense(
-        model=enc, device="cpu", wsi=wsi, coordinates=coords,
-        requested_spacing_um=0.5, target_size=target_size,
+        model=enc, device="cpu", tiling_result=result,
         window_size=window_size, feature_kind=feature_kind, batch_size=2,
     )
 
-    assert wsi.calls == []  # iteration is lazy: building the generator reads nothing
+    backend = fake_backend.backend
+    assert backend.locations_read == []  # iteration is lazy: building the generator reads nothing
 
     first = next(gen)
     assert first.shape[1:] == (4, 4)
     # First grid is yielded after only the first batch (2 of 5) has been read.
-    assert len(wsi.calls) == 2
+    assert len(backend.locations_read) == 2
     next(gen)
-    assert len(wsi.calls) == 2  # second grid comes from the already-read first batch
+    assert len(backend.locations_read) == 2  # second grid comes from the already-read first batch
     next(gen)
-    assert len(wsi.calls) == 4  # third grid forces the next batch to be read
+    assert len(backend.locations_read) == 4  # third grid forces the next batch to be read
 
     rest = list(gen)
     assert len(rest) == 2
-    assert len(wsi.calls) == len(coords)  # total reads never exceed the coordinate count
+    assert len(backend.locations_read) == len(coords)  # total reads never exceed the coordinate count
+
+
+@pytest.mark.parametrize("feature_kind", ["patch_features", "cls_attention"])
+def test_iter_regions_dense_is_batch_invariant(fake_backend, feature_kind):
+    """Composition is irrelevant: only ``B`` matters, not how coords are grouped (adr/0002).
+
+    The same coordinate yields the same grid whether streamed one-per-batch or all at once.
+    """
+    enc = _encoder()
+    target_size = 64
+    coords = [(0, 0), (64, 0), (0, 64), (64, 64), (128, 0)]
+    result = _make_tiling_result(coords, requested_tile_size_px=target_size)
+
+    def _run(batch_size):
+        return list(iter_regions_dense(
+            model=enc, device="cpu", tiling_result=result,
+            feature_kind=feature_kind, batch_size=batch_size,
+        ))
+
+    per_one = _run(1)
+    all_at_once = _run(len(coords))
+    assert len(per_one) == len(all_at_once) == len(coords)
+    for g1, g_all, loc in zip(per_one, all_at_once, coords):
+        # Cosine >= 1 - 1e-4 per grid position (docs/adr/0002 tolerance).
+        cos = np.sum(g1 * g_all) / (np.linalg.norm(g1) * np.linalg.norm(g_all) + 1e-12)
+        assert cos >= 1 - 1e-4, f"batch composition changed the grid at {loc}: cos={cos}"
 
 
 @pytest.mark.parametrize(
     "kwargs", [{"pad_mode": "bogus"}, {"feature_kind": "bogus"}], ids=["pad_mode", "feature_kind"]
 )
-def test_iter_regions_dense_validates_eagerly_before_any_read(kwargs):
+def test_iter_regions_dense_validates_eagerly_before_any_read(fake_backend, kwargs):
     """Invalid pad mode / feature kind raise at the call site, before any region is read."""
     enc = _encoder()
-    target_size = 64
-    wsi = _FakeWSI(target_h=target_size, target_w=target_size)
+    result = _make_tiling_result([(0, 0)], requested_tile_size_px=64)
 
     with pytest.raises(ValueError):
         # The raise must come from the call itself, not from iterating the result — a
         # single ``def … yield`` would wrongly defer validation to the first ``next()``.
-        iter_regions_dense(
-            model=enc, device="cpu", wsi=wsi, coordinates=[(0, 0)],
-            requested_spacing_um=0.5, target_size=target_size, **kwargs,
-        )
-    assert wsi.calls == []
+        iter_regions_dense(model=enc, device="cpu", tiling_result=result, **kwargs)
+    assert fake_backend.backend.locations_read == []
+    assert fake_backend.open_count == 0
 
 
 @pytest.mark.parametrize(
@@ -238,27 +427,24 @@ def test_resolve_output_dtype_defaults_follow_precision(precision, expected):
 
 
 @pytest.mark.parametrize("dtype,np_dtype", [(torch.float16, np.float16), (torch.float32, np.float32)])
-def test_iter_regions_dense_honours_output_dtype(dtype, np_dtype):
+def test_iter_regions_dense_honours_output_dtype(fake_backend, dtype, np_dtype):
     """An explicit output_dtype materializes the grids in that dtype, deterministically."""
     enc = _encoder()
-    target_size = 64
-    wsi = _FakeWSI(target_h=target_size, target_w=target_size)
+    result = _make_tiling_result([(0, 0)], requested_tile_size_px=64)
     grids = list(iter_regions_dense(
-        model=enc, device="cpu", wsi=wsi, coordinates=[(0, 0)],
-        requested_spacing_um=0.5, target_size=target_size, output_dtype=dtype,
+        model=enc, device="cpu", tiling_result=result, output_dtype=dtype,
     ))
     assert len(grids) == 1
     assert grids[0].dtype == np_dtype
 
 
-def test_iter_regions_dense_rejects_bfloat16_output_eagerly():
+def test_iter_regions_dense_rejects_bfloat16_output_eagerly(fake_backend):
     """output_dtype=bfloat16 (uncrossable by .numpy()) raises at the call site, no read."""
     enc = _encoder()
-    target_size = 64
-    wsi = _FakeWSI(target_h=target_size, target_w=target_size)
+    result = _make_tiling_result([(0, 0)], requested_tile_size_px=64)
     with pytest.raises(ValueError):
         iter_regions_dense(
-            model=enc, device="cpu", wsi=wsi, coordinates=[(0, 0)],
-            requested_spacing_um=0.5, target_size=target_size, output_dtype=torch.bfloat16,
+            model=enc, device="cpu", tiling_result=result, output_dtype=torch.bfloat16,
         )
-    assert wsi.calls == []
+    assert fake_backend.backend.locations_read == []
+    assert fake_backend.open_count == 0
diff --git a/tests/test_regression_inference.py b/tests/test_regression_inference.py
index d2dfda9..78a28ff 100644
--- a/tests/test_regression_inference.py
+++ b/tests/test_regression_inference.py
@@ -2719,13 +2719,19 @@ def test_record_slide_metadata_in_process_list_adds_backend_columns(monkeypatch,
     assert recorded.loc[0, "backend"] == "asap"
 
 
-def test_resolve_slide_backend_uses_tiling_result_backend_for_auto():
+@pytest.mark.parametrize(
+    ("requested_backend", "tiling_result", "expected"),
+    [
+        pytest.param("auto", SimpleNamespace(backend="cucim"), "cucim", id="auto-uses-result"),
+        pytest.param("auto", SimpleNamespace(backend="asap"), "asap", id="auto-keeps-asap"),
+        pytest.param("auto", SimpleNamespace(), "asap", id="auto-falls-back-to-asap"),
+        pytest.param("cucim", SimpleNamespace(backend="asap"), "cucim", id="explicit-wins"),
+    ],
+)
+def test_resolve_slide_backend(requested_backend, tiling_result, expected):
     import slide2vec.runtime.tiling as runtime_tiling
 
-    assert runtime_tiling.resolve_slide_backend(replace(DEFAULT_PREPROCESSING, backend="auto"), SimpleNamespace(backend="cucim")) == "cucim"
-    assert runtime_tiling.resolve_slide_backend(replace(DEFAULT_PREPROCESSING, backend="auto"), SimpleNamespace(backend="asap")) == "asap"
-    assert runtime_tiling.resolve_slide_backend(replace(DEFAULT_PREPROCESSING, backend="auto"), SimpleNamespace()) == "asap"
-    assert runtime_tiling.resolve_slide_backend(replace(DEFAULT_PREPROCESSING, backend="cucim"), SimpleNamespace(backend="asap")) == "cucim"
+    assert runtime_tiling.resolve_slide_backend(requested_backend, tiling_result) == expected
 
 
 def test_preload_asap_wholeslidedata_suppresses_noisy_import(monkeypatch, capfd):