diff --git a/bdpy/dataform/__init__.py b/bdpy/dataform/__init__.py
index 30d66a60..fe9ea43a 100644
--- a/bdpy/dataform/__init__.py
+++ b/bdpy/dataform/__init__.py
@@ -5,6 +5,7 @@
 """
 
 from .datastore import *
+from .feature_hdf5 import FeatureWriter, convert_features_to_hdf5, save_features
 from .features import *
 from .kvs import SQLite3KeyValueStore
 from .pd import *
diff --git a/bdpy/dataform/_feature_chunking.py b/bdpy/dataform/_feature_chunking.py
new file mode 100644
index 00000000..24eaf742
--- /dev/null
+++ b/bdpy/dataform/_feature_chunking.py
@@ -0,0 +1,146 @@
+"""Chunk-shape policy for chunked feature storage.
+
+HDF5 reads whole chunks even when only a few elements of a chunk are selected,
+so the chunk shape - not the file format alone - is what decides the cost of a
+partial read. This module turns an array shape and dtype into a chunk shape
+under a byte budget, rather than hard-coding an extent along any one axis.
+
+This file is a part of BdPy.
+"""
+
+import math
+from typing import Sequence, Tuple
+
+import numpy as np
+
+__all__ = ["DEFAULT_TARGET_CHUNK_BYTES", "choose_chunk_shape"]
+
+#: Default per-chunk byte budget. h5py's own guidance puts a useful chunk
+#: somewhere between 10 KiB and 1 MiB; features are large and read in slabs, so
+#: we sit at the top of that range.
+DEFAULT_TARGET_CHUNK_BYTES = 1 << 20
+
+
+def _prod(values: Sequence[int]) -> int:
+    out = 1
+    for v in values:
+        out *= int(v)
+    return out
+
+
+def choose_chunk_shape(
+    shape: Sequence[int],
+    dtype: np.dtype,
+    target_bytes: int = DEFAULT_TARGET_CHUNK_BYTES,
+    n_samples_known: bool = True,
+) -> Tuple[int, ...]:
+    """Choose an HDF5 chunk shape for a feature array under a byte budget.
+
+    The two axes callers actually slice are axis 0 (samples) and axis 1 (the
+    outermost feature axis, typically channels or units), so both are chunked.
+    Trailing axes - the spatial axes of a convolutional feature map, say
+    ``(13, 13)`` - are kept whole: they are small and are never sliced in
+    practice, and splitting them would only shrink the chunks.
+
+    Within the budget the two chunked axes are split roughly evenly in element
+    count, so that neither a sample-wise read nor a channel-wise read degenerates
+    into reading the entire dataset. Each extent is then shrunk to the smallest
+    one needing the same number of chunks, which keeps HDF5 from padding the
+    edge chunks and inflating the file.
+
+    Parameters
+    ----------
+    shape : sequence of int
+        Full array shape, ``(n_samples, *feature_shape)``.
+    dtype : numpy.dtype
+        Array dtype; only its itemsize matters.
+    target_bytes : int, optional
+        Size budget for one chunk, in bytes
+        (default: ``DEFAULT_TARGET_CHUNK_BYTES``). It is an upper bound only
+        when the budget can be met at all: the trailing axes are kept whole, so
+        a shape whose trailing axes alone exceed `target_bytes` necessarily
+        produces a larger chunk.
+    n_samples_known : bool, optional
+        Whether ``shape[0]`` is the final number of samples. Pass ``False`` when
+        writing to a resizable dataset whose sample count is not yet known, so
+        that the sample-axis extent is not capped by the current size.
+
+    Returns
+    -------
+    tuple of int
+        Chunk shape, same length as ``shape``. Every element is at least 1, and
+        no larger than the corresponding entry of ``shape`` for axes of
+        non-zero size (a layer with zero samples still needs a chunk extent of
+        at least 1 on the sample axis, since HDF5 forbids a zero chunk dim).
+
+    Raises
+    ------
+    ValueError
+        If `shape` has fewer than two axes, any feature axis has size 0, the
+        sample axis is negative, or `target_bytes` is not positive.
+
+    Examples
+    --------
+    >>> choose_chunk_shape((1200, 1000), np.dtype(np.float32))
+    (400, 500)
+    >>> choose_chunk_shape((1200, 256, 13, 13), np.dtype(np.float32))
+    (39, 37, 13, 13)
+    >>> choose_chunk_shape((50, 1000), np.dtype(np.float32))
+    (50, 1000)
+    """
+    shape = tuple(int(s) for s in shape)
+    if len(shape) < 2:
+        raise ValueError(
+            "Feature arrays need at least a sample axis and a feature axis; "
+            "got shape {}".format(shape)
+        )
+    if target_bytes < 1:
+        raise ValueError("target_bytes must be positive, got {}".format(target_bytes))
+    if shape[0] < 0:
+        raise ValueError("sample axis cannot be negative, got {}".format(shape))
+    if any(s < 1 for s in shape[1:]):
+        # A zero-width feature axis carries nothing and cannot be chunked.
+        # A zero-sample layer, by contrast, is meaningful and is allowed below.
+        raise ValueError(
+            "feature axes must all be non-empty, got shape {}".format(shape)
+        )
+
+    itemsize = np.dtype(dtype).itemsize
+
+    # The whole array fits in one chunk: nothing to gain from splitting it.
+    # Clamp to at least 1 per axis -- a zero-sample layer is legal but HDF5
+    # rejects a chunk with a zero dimension.
+    if n_samples_known and _prod(shape) * itemsize <= target_bytes:
+        return tuple(max(1, s) for s in shape)
+
+    # Bytes taken by one (sample, feature-0) cell, i.e. one element of the two
+    # chunked axes with the trailing axes kept whole.
+    inner = max(1, _prod(shape[2:])) * itemsize
+    budget_cells = max(1, target_bytes // inner)
+
+    # Split the budget evenly between the two chunked axes, then let each axis
+    # absorb the headroom the other did not need: a short sample axis should not
+    # force a narrow feature chunk, and vice versa.
+    c1 = min(shape[1], max(1, math.isqrt(budget_cells)))
+    c0 = max(1, budget_cells // c1)
+    if n_samples_known and shape[0] > 0:
+        c0 = min(shape[0], c0)
+    c1 = min(shape[1], max(c1, budget_cells // c0))
+
+    # HDF5 allocates whole chunks, so a chunk extent that divides its axis
+    # unevenly pads the edge chunks and inflates the file. Shrinking the extent
+    # to the smallest one that needs the same number of chunks removes most of
+    # that padding, and can only lower the chunk size, so the budget still holds.
+    c1 = _snap(shape[1], c1)
+    if n_samples_known and shape[0] > 0:
+        c0 = _snap(shape[0], c0)
+
+    return (c0, c1, *shape[2:])
+
+
+def _snap(length: int, chunk: int) -> int:
+    """Smallest extent covering `length` in the same number of chunks as `chunk`."""
+    if length < 1 or chunk >= length:
+        return chunk
+    n_chunks = -(-length // chunk)  # ceil
+    return -(-length // n_chunks)   # ceil
diff --git a/bdpy/dataform/_feature_store.py b/bdpy/dataform/_feature_store.py
new file mode 100644
index 00000000..d57c7ff5
--- /dev/null
+++ b/bdpy/dataform/_feature_store.py
@@ -0,0 +1,909 @@
+"""Storage backends for DNN feature collections.
+
+``Features`` used to be welded to one on-disk layout: a directory per layer
+holding one ``.mat`` file per stimulus. That layout has a file boundary only on
+the sample axis, so reading a few channels still costs a full read of every
+stimulus. This module puts a backend interface between ``Features`` and the
+bytes, so the legacy layout and the chunked HDF5 layout are interchangeable and
+a future backend needs no change to ``Features`` itself.
+
+Every backend answers the same question -- "give me these labels, sliced this
+way along the feature axes" -- and differs only in how much it has to read to do
+it. The legacy backend loads whole stimulus files and slices in memory; the HDF5
+backend pushes the slice down to h5py and reads only the chunks it covers.
+
+This file is a part of BdPy.
+"""
+
+import glob
+import os
+from abc import ABC, abstractmethod
+from functools import partial
+from multiprocessing import Pool
+from typing import TYPE_CHECKING, Dict, Iterator, List, Optional, Sequence, Tuple, Union
+
+import h5py
+import numpy as np
+
+from . import _mat_v73
+from ._feature_chunking import DEFAULT_TARGET_CHUNK_BYTES
+
+__all__ = [
+    "FeatureStore",
+    "HDF5FeatureStore",
+    "MatFeatureStore",
+    "detect_format",
+]
+
+#: Root attribute marking a file as bdpy chunked feature storage.
+FORMAT_ATTR = "bdpy_format"
+#: Root attribute holding the schema version.
+FORMAT_VERSION_ATTR = "bdpy_format_version"
+#: Value of ``FORMAT_ATTR`` for feature files.
+FORMAT_NAME = "features"
+#: Highest schema version this module can read.
+SUPPORTED_FORMAT_VERSION = 1
+#: Dataset holding the feature array, shape ``(n_samples, *feature_shape)``.
+FEATURES_DATASET = "features"
+#: Dataset holding the stimulus labels, shape ``(n_samples,)``.
+LABELS_DATASET = "labels"
+#: Extension of a chunked feature file.
+HDF5_EXT = "h5"
+
+# A slice spec for the feature axes. Deliberately narrow: basic indexing only,
+# i.e. what ``numpy.s_[128:256]`` or ``numpy.s_[8:16, 1:4]`` produces. See
+# _validate_feature_slice for what is rejected and why.
+# `types.EllipsisType` is 3.10+, and a bare "ellipsis" forward reference does not
+# resolve at runtime (it is not a Python name), which makes
+# typing.get_type_hints() raise. typeshed does define builtins.ellipsis, so the
+# checker gets the real type and the interpreter gets an equivalent object.
+if TYPE_CHECKING:
+    from builtins import ellipsis
+else:
+    ellipsis = type(Ellipsis)
+BasicIndexer = Union[slice, int, np.integer, ellipsis]
+FeatureSlice = Union[BasicIndexer, Tuple[BasicIndexer, ...], None]
+
+
+def _load_array_with_key(key: str, path: str) -> np.ndarray:
+    """Load one array from a ``.mat`` file. Shared with :mod:`bdpy.dataform.features`.
+
+    v5 ``.mat`` via scipy, v7.3 (HDF5) via h5py; avoids hdf5storage on the load
+    path, which breaks under NumPy 2.0 (see :mod:`bdpy.dataform._mat_v73`).
+    `key` comes first so ``partial(..., 'feat')`` is a ``Pool.map`` callable.
+    """
+    return _mat_v73.loadmat_key(path, key)
+
+
+def _determine_num_parallel(num_files: int) -> int:
+    # NOTE: optimal number of parallel processes is not clear. It could depend
+    # on several factors such as the number of files, the size of files, the
+    # number of cores, etc. For now, we use a simple heuristic based on the
+    # number of files.
+    num_parallel: int
+    if num_files < 16:
+        num_parallel = 1
+    elif num_files < 64:
+        num_parallel = 16
+    else:
+        num_parallel = 64
+    return num_parallel
+
+
+def _normalize_feature_slice(feature_slice: FeatureSlice) -> Tuple:
+    """Normalize and validate a feature-axis slice spec.
+
+    Returns a tuple of per-axis indexers; ``None`` means "no slice" and gives
+    the empty tuple. Note this makes ``numpy.newaxis`` (which *is* ``None``)
+    inexpressible, so a bare newaxis reads as "no slice"; inside a tuple it is
+    rejected explicitly.
+    """
+    if feature_slice is None:
+        return ()
+    indexers = feature_slice if isinstance(feature_slice, tuple) else (feature_slice,)
+    _validate_feature_slice(indexers)
+    return indexers
+
+
+def _validate_feature_slice(indexers: Tuple) -> None:
+    """Reject anything outside basic, forward indexing.
+
+    `feature_slice` exists to read part of a layer, not to reimplement NumPy
+    indexing. Keeping it to basic forward indexing is what lets every backend
+    mean the same thing by it: h5py handles slices with a positive step
+    (negative start/stop included), integers and one ``Ellipsis`` exactly as
+    NumPy does, but rejects a negative step outright. Allowing more would make
+    the legacy and HDF5 backends diverge, or require reimplementing NumPy's
+    semantics to keep them together.
+
+    Parameters
+    ----------
+    indexers : tuple
+        Per-axis indexers, as returned by :func:`_normalize_feature_slice`.
+
+    Raises
+    ------
+    ValueError
+        If any indexer is outside the supported set.
+    """
+    advice = (
+        " feature_slice supports basic forward indexing only (slices with a "
+        "positive step, integers and a single Ellipsis). Read without "
+        "feature_slice and index the resulting array with NumPy instead."
+    )
+    if sum(ix is Ellipsis for ix in indexers) > 1:
+        raise ValueError("feature_slice may contain at most one Ellipsis." + advice)
+    for ix in indexers:
+        # bool is a subclass of int in Python, so it must be rejected *before*
+        # the integer check or True would silently be read as index 1.
+        if isinstance(ix, (bool, np.bool_)):
+            raise ValueError(
+                "feature_slice does not support boolean indexing." + advice
+            )
+        if ix is Ellipsis or isinstance(ix, (int, np.integer)):
+            continue
+        if isinstance(ix, slice):
+            for part, value in (
+                ("start", ix.start), ("stop", ix.stop), ("step", ix.step),
+            ):
+                if value is None:
+                    continue
+                if isinstance(value, (bool, np.bool_)) or not isinstance(
+                    value, (int, np.integer)
+                ):
+                    raise ValueError(
+                        "feature_slice slice {} must be an integer or None, "
+                        "got {!r}.".format(part, value) + advice
+                    )
+            if ix.step is not None and ix.step < 1:
+                raise ValueError(
+                    "feature_slice does not support a step below 1 "
+                    "(got {}).".format(ix.step) + advice
+                )
+            continue
+        if ix is None:
+            raise ValueError(
+                "feature_slice does not support numpy.newaxis." + advice
+            )
+        raise ValueError(
+            "feature_slice does not support {} indexing.".format(
+                type(ix).__name__
+            )
+            + advice
+        )
+
+
+def _normalize_axis(axis: int, shape: Tuple[int, ...]) -> int:
+    """Return `axis` as a non-negative axis of an array of `shape`.
+
+    Raises
+    ------
+    ValueError
+        If `axis` is out of range for `shape`.
+    """
+    normalized = axis + len(shape) if axis < 0 else axis
+    if not 0 <= normalized < len(shape):
+        raise ValueError(
+            "axis {} is out of range for selected shape {}".format(axis, shape)
+        )
+    return normalized
+
+
+def _check_size(size: int) -> None:
+    """Reject a slab size below 1."""
+    if size < 1:
+        raise ValueError("size must be positive, got {}".format(size))
+
+
+def _partition(length: int, size: int) -> Iterator[slice]:
+    """Split ``range(length)`` into consecutive slices of at most `size`.
+
+    The slices cover the range exactly once, in order; only the last may be
+    shorter. ``length == 0`` yields nothing. `size` must already be positive.
+
+    A generator rather than a list: a materialized one would cost an object per
+    slab, which for a large axis at size 1 is millions of them.
+    """
+    for start in range(0, length, size):
+        yield slice(start, min(start + size, length))
+
+
+def _axis_index(axis: int, sl: slice) -> Tuple[slice, ...]:
+    """Index tuple that applies `sl` to `axis` and leaves earlier axes whole."""
+    return (slice(None),) * axis + (sl,)
+
+
+def _default_iter_size(
+    shape: Tuple[int, ...], axis: int, itemsize: int, extent: Optional[int]
+) -> int:
+    """Default slab size along `axis` of an array of `shape`.
+
+    The on-disk chunk `extent` along the axis when there is one, so iteration is
+    chunk-aligned; otherwise as many elements as fit in the byte budget. Always
+    at least 1 and at most the axis length (when that is positive).
+    """
+    if extent is not None:
+        return max(1, min(extent, shape[axis]))
+    other = 1
+    for i, n in enumerate(shape):
+        if i != axis:
+            other *= int(n)
+    per_element = max(1, other * itemsize)
+    # Read the module constant at call time so it can be patched in tests.
+    return max(1, min(shape[axis], DEFAULT_TARGET_CHUNK_BYTES // per_element))
+
+
+def _iter_array_slabs(
+    array: np.ndarray, axis: int, slices: Iterator[slice]
+) -> Iterator[Tuple[slice, np.ndarray]]:
+    """Yield ``(sl, view)`` for each of `slices` applied to `axis` of `array`.
+
+    The slabs are views, so holding one keeps `array` alive.
+    """
+    for sl in slices:
+        yield sl, array[_axis_index(axis, sl)]
+
+
+class FeatureStore(ABC):
+    """Interface to a collection of DNN features on disk.
+
+    A store exposes layers and stimulus labels, and reads features for a set of
+    labels with an optional slice along the feature axes. Implementations differ
+    only in how much I/O a given request costs.
+    """
+
+    @property
+    @abstractmethod
+    def layers(self) -> List[str]:
+        """List of DNN layers held by this store."""
+
+    @property
+    @abstractmethod
+    def labels(self) -> List[str]:
+        """List of stimulus labels, in the store's own order."""
+
+    @abstractmethod
+    def shape(self, layer: str) -> Tuple[int, ...]:
+        """Full shape ``(n_samples, *feature_shape)`` of `layer`."""
+
+    @abstractmethod
+    def dtype(self, layer: str) -> np.dtype:
+        """Dtype of the features in `layer`."""
+
+    def read(
+        self,
+        layer: str,
+        labels: Optional[Sequence[str]] = None,
+        feature_slice: FeatureSlice = None,
+    ) -> np.ndarray:
+        """Read features from `layer`.
+
+        `feature_slice` is validated here, before any I/O, so every backend
+        accepts and rejects exactly the same specs. Backends implement
+        :meth:`_read` and must not override this method.
+
+        Parameters
+        ----------
+        layer : str
+            DNN layer.
+        labels : sequence of str, optional
+            Stimulus labels to read. ``None`` reads every label in store order.
+            Otherwise rows come back **in the order given**, and repeated labels
+            yield repeated rows.
+        feature_slice : slice, int, Ellipsis or tuple of those, optional
+            Index applied to the feature axes (axes 1 and up), as produced by
+            ``numpy.s_[...]``. Basic forward indexing only: slices with a
+            positive step and integer bounds, integers, and at most one
+            ``Ellipsis``. Fancy indexing, a negative step, a non-integer slice
+            bound, booleans, and ``numpy.newaxis`` *inside a tuple* raise
+            ``ValueError``. A bare ``numpy.newaxis`` is ``None``, the "no
+            slice" default, so it reads the whole feature tensor.
+
+        Returns
+        -------
+        numpy.ndarray
+            Array of shape ``(n_labels, *sliced_feature_shape)``.
+        """
+        return self._read(layer, labels, _normalize_feature_slice(feature_slice))
+
+    @abstractmethod
+    def _read(
+        self,
+        layer: str,
+        labels: Optional[Sequence[str]],
+        indexers: Tuple,
+    ) -> np.ndarray:
+        """Backend read: `indexers` is already validated, ``()`` for no slice."""
+
+    def chunk_extent(self, layer: str, axis: int) -> Optional[int]:
+        """On-disk chunk extent of `layer` along `axis`, if the store has one.
+
+        Returns ``None`` when the backend has no chunking to align with, in
+        which case :meth:`iter_chunks` falls back to a byte budget.
+        """
+        return None
+
+    def iter_chunks(
+        self,
+        layer: str,
+        labels: Optional[Sequence[str]] = None,
+        axis: int = 1,
+        size: Optional[int] = None,
+    ) -> Iterator[Tuple[slice, np.ndarray]]:
+        """Iterate over `layer` in slabs along `axis`.
+
+        Yields ``(sl, block)`` pairs, where `sl` is the slice applied to `axis`
+        of the selected array and `block` is that slab. The slice is yielded so
+        that callers can place results back without tracking offsets themselves.
+
+        There is deliberately no `feature_slice` here. Composing an arbitrary
+        slice with the per-slab slice means reimplementing NumPy's index
+        arithmetic, which is both the source of subtle wrong results and more
+        generality than reading a layer in slabs needs. Slice the blocks as they
+        come out instead.
+
+        Parameters
+        ----------
+        layer : str
+            DNN layer.
+        labels : sequence of str, optional
+            Stimulus labels, as in :meth:`read`.
+        axis : int, optional
+            Axis to iterate over. Axis 0 is the sample axis; the default, axis
+            1, is the outermost feature axis, which is the one partial feature
+            reads are about.
+        size : int, optional
+            Number of elements per slab. Defaults to the on-disk chunk extent
+            along `axis` when the backend has one, so that iteration is
+            chunk-aligned and every element is read exactly once; otherwise a
+            size derived from the same byte budget is used.
+
+        Yields
+        ------
+        tuple of (slice, numpy.ndarray)
+            The slice applied to `axis`, and the corresponding slab.
+
+        Notes
+        -----
+        How much is resident at once depends on the backend. Along the sample
+        axis every backend reads one slab at a time. Along a feature axis only
+        a backend with real partial reads can do so: the legacy ``.mat``
+        layout has no file boundary there, so it reads the selection once and
+        yields views into it, and the blocks keep that array alive for as long
+        as they are held.
+        """
+        selected_shape = self._selected_shape(layer, labels)
+        axis = _normalize_axis(axis, selected_shape)
+
+        length = selected_shape[axis]
+        if size is None:
+            size = _default_iter_size(
+                selected_shape,
+                axis,
+                np.dtype(self.dtype(layer)).itemsize,
+                self.chunk_extent(layer, axis),
+            )
+        _check_size(size)
+
+        if length == 0:
+            # Nothing to iterate. Returning here rather than falling into the
+            # loop keeps a backend that reads the selection up front from
+            # reading anything the per-slab path would not have read either.
+            return
+
+        slices = _partition(length, size)
+        if axis == 0:
+            all_labels = list(self.labels) if labels is None else list(labels)
+            for sl in slices:
+                yield sl, self.read(layer, all_labels[sl])
+        else:
+            yield from self._iter_feature_slabs(layer, labels, axis, slices)
+
+    def _iter_feature_slabs(
+        self,
+        layer: str,
+        labels: Optional[Sequence[str]],
+        axis: int,
+        slices: Iterator[slice],
+    ) -> Iterator[Tuple[slice, np.ndarray]]:
+        """Yield the feature-axis slabs of `layer` named by `slices`.
+
+        This is the seam for a backend that cannot read a feature axis
+        partially. `axis` is already normalized and validated to a feature axis
+        (``>= 1``), and `slices` partitions that axis, so an implementation
+        decides only *how* a slab is obtained -- never which slabs there are.
+
+        The default reads each slab on its own, which is what a backend with
+        real partial reads wants.
+        """
+        for sl in slices:
+            # axis >= 1 indexes the feature axes, where axis 1 of the array is
+            # entry 0 of the feature-slice tuple. Go through read(), the single
+            # entry point that validates.
+            yield sl, self.read(layer, labels, _axis_index(axis - 1, sl))
+
+    def _selected_shape(
+        self,
+        layer: str,
+        labels: Optional[Sequence[str]] = None,
+    ) -> Tuple[int, ...]:
+        """Shape the selection would have, without reading the data."""
+        full = self.shape(layer)
+        n_samples = full[0] if labels is None else len(labels)
+        return (n_samples, *full[1:])
+
+
+
+class MatFeatureStore(FeatureStore):
+    """Legacy per-stimulus ``.mat`` feature directory.
+
+    Layout is ``<dpath>/<layer>/<label>.<ext>``, each file holding one sample
+    under `key` with a leading sample axis. There is no file boundary on the
+    feature axes, so `feature_slice` is applied after the full stimulus files
+    have been loaded and concatenated -- the read cost is the same as before,
+    and the slice only saves the caller from doing it themselves. For the same
+    reason :meth:`iter_chunks` along a feature axis reads the selection once
+    and cuts the slabs out of it in memory; only sample-axis iteration is
+    bounded by the slab.
+
+    Parameters
+    ----------
+    dpath : str
+        Feature directory.
+    ext : str, optional
+        Feature file extension (default: ``'mat'``).
+    key : str, optional
+        Variable name inside each file (default: ``'feat'``).
+    """
+
+    def __init__(self, dpath: str, ext: str = "mat", key: str = "feat"):
+        self._dpath = dpath
+        self._ext = ext
+        self._key = key
+        self._layers = self._collect_layers()
+        self._labels = self._collect_labels()
+        self._file_table: Dict[str, Dict[str, str]] = {
+            layer: {
+                label: os.path.join(dpath, layer, label + "." + ext)
+                for label in self._labels
+            }
+            for layer in self._layers
+        }
+
+    @property
+    def layers(self) -> List[str]:
+        return self._layers
+
+    @property
+    def labels(self) -> List[str]:
+        return self._labels
+
+    def path(self, layer: str, label: str) -> str:
+        """Path of the file holding `label` in `layer`."""
+        return self._file_table[layer][label]
+
+    def shape(self, layer: str) -> Tuple[int, ...]:
+        if not self._labels:
+            raise RuntimeError("No features found in {}".format(self._dpath))
+        sample = _load_array_with_key(self._key, self.path(layer, self._labels[0]))
+        return (len(self._labels), *sample.shape[1:])
+
+    def dtype(self, layer: str) -> np.dtype:
+        if not self._labels:
+            raise RuntimeError("No features found in {}".format(self._dpath))
+        sample = _load_array_with_key(self._key, self.path(layer, self._labels[0]))
+        return sample.dtype
+
+    def _read(
+        self,
+        layer: str,
+        labels: Optional[Sequence[str]],
+        indexers: Tuple,
+    ) -> np.ndarray:
+        if labels is None:
+            labels = self._labels
+        paths = [self.path(layer, label) for label in labels]
+
+        num_parallel = _determine_num_parallel(len(paths))
+        load = partial(_load_array_with_key, self._key)
+        if num_parallel == 1:
+            arrays = list(map(load, paths))
+        else:
+            with Pool(processes=num_parallel) as pool:
+                arrays = pool.map(load, paths)
+        features = np.concatenate(arrays, axis=0)
+
+        if indexers:
+            features = features[(slice(None), *indexers)]
+        return features
+
+    def _iter_feature_slabs(
+        self,
+        layer: str,
+        labels: Optional[Sequence[str]],
+        axis: int,
+        slices: Iterator[slice],
+    ) -> Iterator[Tuple[slice, np.ndarray]]:
+        """Read the selection once, then cut every slab out of it in memory.
+
+        This layout has no file boundary on the feature axes, so one slab
+        already costs a read of every selected stimulus file. Reading per slab
+        therefore re-read the whole selection once per slab, making a k-slab
+        iteration k times the cost of a single read.
+
+        Nothing is kept on the store, so this is a per-call strategy rather
+        than a cache. The selection is a plain array held for the duration of
+        the iteration, and the yielded slabs are views into it -- exactly as
+        they were when each slab was read separately -- so holding a slab keeps
+        the whole selection alive. Sample-axis iteration is deliberately left
+        to the base class, where each slab reads only its own files and a layer
+        larger than memory can still be processed end to end.
+        """
+        yield from _iter_array_slabs(self.read(layer, labels), axis, slices)
+
+    def _collect_layers(self) -> List[str]:
+        return sorted(
+            d
+            for d in os.listdir(self._dpath)
+            if os.path.isdir(os.path.join(self._dpath, d))
+        )
+
+    def _collect_labels(self) -> List[str]:
+        labels: List[str] = []
+        for layer in self._layers:
+            layer_dir = os.path.join(self._dpath, layer)
+            layer_dir = layer_dir.replace("[", "[[]")  # Use glob.escape for Python 3.4 or later
+            files = glob.glob(os.path.join(layer_dir, "*." + self._ext))
+            labels_t = sorted(
+                os.path.splitext(os.path.basename(f))[0] for f in files
+            )
+            if not labels:
+                labels = labels_t
+            elif labels != labels_t:
+                raise RuntimeError("Invalid feature file in %s " % self._dpath)
+        return labels
+
+
+class HDF5FeatureStore(FeatureStore):
+    """Chunked HDF5 feature directory, one ``<layer>.h5`` file per layer.
+
+    Each file holds ``/features`` with shape ``(n_samples, *feature_shape)`` and
+    ``/labels`` with the matching stimulus labels. ``/features`` is explicitly
+    chunked, so a slice along the feature axes reads only the chunks it covers
+    instead of the whole layer.
+
+    Parameters
+    ----------
+    dpath : str
+        Directory holding ``<layer>.h5`` files.
+    """
+
+    def __init__(self, dpath: str):
+        self._dpath = dpath
+        self._layers = sorted(
+            os.path.splitext(os.path.basename(p))[0]
+            for p in glob.glob(os.path.join(glob.escape(dpath), "*." + HDF5_EXT))
+        )
+        if not self._layers:
+            raise RuntimeError("No .{} feature file found in {}".format(HDF5_EXT, dpath))
+        self._labels: List[str] = []
+        self._label_index: Dict[str, int] = {}
+        self._validate_all()
+
+    @property
+    def layers(self) -> List[str]:
+        return self._layers
+
+    @property
+    def labels(self) -> List[str]:
+        return self._labels
+
+    def path(self, layer: str) -> str:
+        """Path of the file holding `layer`."""
+        return os.path.join(self._dpath, layer + "." + HDF5_EXT)
+
+    def shape(self, layer: str) -> Tuple[int, ...]:
+        with self._open(layer) as f:
+            return tuple(int(n) for n in f[FEATURES_DATASET].shape)
+
+    def dtype(self, layer: str) -> np.dtype:
+        with self._open(layer) as f:
+            dtype: np.dtype = np.dtype(f[FEATURES_DATASET].dtype)
+        return dtype
+
+    def chunk_extent(self, layer: str, axis: int) -> Optional[int]:
+        with self._open(layer) as f:
+            chunks = f[FEATURES_DATASET].chunks
+        if chunks is None or axis >= len(chunks):
+            return None
+        return int(chunks[axis])
+
+    def _read(
+        self,
+        layer: str,
+        labels: Optional[Sequence[str]],
+        indexers: Tuple,
+    ) -> np.ndarray:
+        with self._open(layer) as f:
+            dset = f[FEATURES_DATASET]
+
+            if labels is None:
+                return self._read_rows(dset, slice(None), indexers)
+
+            rows = self._row_indices(labels)
+            if rows.size == dset.shape[0] and np.array_equal(
+                rows, np.arange(dset.shape[0])
+            ):
+                # Asking for every row in order: a plain slice reads
+                # contiguously, where a full index list would not.
+                return self._read_rows(dset, slice(None), indexers)
+            # h5py wants a strictly increasing index list and allows at most one
+            # fancy index per selection. Read each distinct row once in order,
+            # then restore the caller's order (and any repeats) with `inverse`.
+            uniq, inverse = np.unique(rows, return_inverse=True)
+            block = self._read_rows(dset, [int(i) for i in uniq], indexers)
+            if uniq.size == rows.size and np.array_equal(uniq, rows):
+                return block
+            return block[inverse]
+
+    @staticmethod
+    def _read_rows(
+        dset: h5py.Dataset, rows: Union[slice, List[int]], indexers: Tuple
+    ) -> np.ndarray:
+        """Read `rows` from `dset`, pushing `indexers` down into the h5py selection."""
+        if not indexers:
+            return np.asarray(dset[rows])
+        # Basic forward indexers are not fancy, so they ride along in the same
+        # selection as the row list and h5py reads only the chunks they cover.
+        # _validate_feature_slice has already rejected anything else, which is
+        # what keeps this to a single line.
+        return np.asarray(dset[(rows, *indexers)])
+
+    def _row_indices(self, labels: Sequence[str]) -> np.ndarray:
+        try:
+            return np.array([self._label_index[label] for label in labels], dtype=int)
+        except KeyError as exc:
+            raise KeyError(
+                "Label {} not found in {}".format(exc.args[0], self._dpath)
+            ) from None
+
+    def _open(self, layer: str) -> h5py.File:
+        return h5py.File(self.path(layer), "r")
+
+    def _validate_all(self) -> None:
+        """Validate every layer file and pin down the shared label sequence.
+
+        A row index built from one layer is used against every layer's
+        ``/features``, so the layers must agree on the labels *and on their
+        order*; otherwise a label lookup would silently return another
+        stimulus' row. The legacy ``.mat`` backend enforces the same invariant
+        in :meth:`MatFeatureStore._collect_labels`.
+        """
+        reference: Optional[List[str]] = None
+        reference_layer = ""
+        for layer in self._layers:
+            with self._open(layer) as f:
+                labels = _validate_format(f, self.path(layer))
+            if reference is None:
+                reference, reference_layer = labels, layer
+            elif labels != reference:
+                raise RuntimeError(
+                    _label_mismatch_message(
+                        self._dpath, reference_layer, reference, layer, labels
+                    )
+                )
+
+        assert reference is not None  # self._layers is non-empty
+        self._labels = reference
+        self._label_index = {label: i for i, label in enumerate(reference)}
+
+
+def _validate_format(f: h5py.File, path: str) -> List[str]:
+    """Check that `f` is well-formed bdpy feature storage, and return its labels.
+
+    Everything a reader relies on is checked in one place: the format marker,
+    a version in the range this build understands, both datasets present with
+    the expected rank, one label per feature row, and no duplicate labels.
+
+    Duplicate labels matter more than they look. Labels are mapped to row
+    indices once, so a repeated label would silently resolve every occurrence
+    to the last row and drop the earlier one. The legacy ``.mat`` layout cannot
+    express duplicates at all -- the file name *is* the label -- so rejecting
+    them keeps the two layouts equivalent.
+
+    Parameters
+    ----------
+    f : h5py.File
+        Open feature file.
+    path : str
+        Its path, for error messages.
+
+    Returns
+    -------
+    list of str
+        The file's labels, already decoded.
+
+    Raises
+    ------
+    RuntimeError
+        If any of the above does not hold.
+    """
+    fmt = f.attrs.get(FORMAT_ATTR)
+    if isinstance(fmt, bytes):
+        fmt = fmt.decode("utf-8")
+    if fmt != FORMAT_NAME:
+        raise RuntimeError(
+            "{} is not bdpy feature storage ({} = {!r}); expected {!r}".format(
+                path, FORMAT_ATTR, fmt, FORMAT_NAME
+            )
+        )
+
+    if FORMAT_VERSION_ATTR not in f.attrs:
+        raise RuntimeError(
+            "{} has no {} attribute, so it is not valid bdpy feature "
+            "storage.".format(path, FORMAT_VERSION_ATTR)
+        )
+    raw_version = f.attrs[FORMAT_VERSION_ATTR]
+    # Require an actual integer scalar rather than coercing: int(1.5) would
+    # silently round a malformed version down to a supported one. bool is
+    # checked first because it is a subclass of int.
+    if isinstance(raw_version, (bool, np.bool_)) or not isinstance(
+        raw_version, (int, np.integer)
+    ):
+        raise RuntimeError(
+            "{} has a malformed {} ({!r}); expected an integer.".format(
+                path, FORMAT_VERSION_ATTR, raw_version
+            )
+        )
+    version = int(raw_version)
+    if version < 1:
+        raise RuntimeError(
+            "{} declares feature storage format version {}; versions start at "
+            "1.".format(path, version)
+        )
+    if version > SUPPORTED_FORMAT_VERSION:
+        raise RuntimeError(
+            "{} uses feature storage format version {}, but this version of "
+            "bdpy supports up to version {}. Please upgrade bdpy.".format(
+                path, version, SUPPORTED_FORMAT_VERSION
+            )
+        )
+
+    for name in (FEATURES_DATASET, LABELS_DATASET):
+        if name not in f:
+            raise RuntimeError("{} has no /{} dataset".format(path, name))
+
+    features, labels_dset = f[FEATURES_DATASET], f[LABELS_DATASET]
+    if features.ndim < 2:
+        raise RuntimeError(
+            "{}: /{} must have a sample axis and at least one feature axis, "
+            "got shape {}".format(path, FEATURES_DATASET, features.shape)
+        )
+    if labels_dset.ndim != 1:
+        raise RuntimeError(
+            "{}: /{} must be one-dimensional, got shape {}".format(
+                path, LABELS_DATASET, labels_dset.shape
+            )
+        )
+
+    labels = _decode_labels(labels_dset[()])
+    if len(labels) != features.shape[0]:
+        raise RuntimeError(
+            "{}: /{} has {} rows but /{} has {} entries".format(
+                path, FEATURES_DATASET, features.shape[0],
+                LABELS_DATASET, len(labels),
+            )
+        )
+    duplicates = _duplicates(labels)
+    if duplicates:
+        raise RuntimeError(
+            "{}: /{} contains duplicate labels ({}). Labels identify rows, so "
+            "they must be unique.".format(
+                path, LABELS_DATASET, ", ".join(repr(d) for d in duplicates[:3])
+            )
+        )
+    return labels
+
+
+def _duplicates(labels: Sequence[str]) -> List[str]:
+    """Labels appearing more than once, in first-seen order."""
+    seen: Dict[str, int] = {}
+    for label in labels:
+        seen[label] = seen.get(label, 0) + 1
+    return [label for label, n in seen.items() if n > 1]
+
+
+def _label_mismatch_message(
+    dpath: str,
+    reference_layer: str,
+    reference: List[str],
+    layer: str,
+    labels: List[str],
+) -> str:
+    """Describe *how* two layers' label sequences differ, not just that they do."""
+    if len(labels) != len(reference):
+        detail = "{} has {} labels, {} has {}".format(
+            reference_layer, len(reference), layer, len(labels)
+        )
+    elif sorted(labels) == sorted(reference):
+        first = next(
+            i for i, (a, b) in enumerate(zip(reference, labels)) if a != b
+        )
+        detail = (
+            "same labels in a different order; first difference at index {}: "
+            "{} has {!r}, {} has {!r}".format(
+                first, reference_layer, reference[first], layer, labels[first]
+            )
+        )
+    else:
+        only_ref = sorted(set(reference) - set(labels))[:3]
+        only_this = sorted(set(labels) - set(reference))[:3]
+        detail = "different labels; only in {}: {}, only in {}: {}".format(
+            reference_layer, only_ref or "-", layer, only_this or "-"
+        )
+    return (
+        "Inconsistent labels across layers in {}: {}. Every layer must hold the "
+        "same stimulus labels in the same order.".format(dpath, detail)
+    )
+
+
+def _decode_labels(raw: np.ndarray) -> List[str]:
+    """Decode an HDF5 string dataset into a list of str (h5py 3.x gives bytes)."""
+    return [
+        item.decode("utf-8") if isinstance(item, bytes) else str(item)
+        for item in np.asarray(raw).ravel().tolist()
+    ]
+
+
+def detect_format(dpath: str, ext: str = "mat") -> str:
+    """Work out which storage layout `dpath` holds.
+
+    The two layouts are told apart by the files that are actually there, not by
+    the mere presence of a subdirectory: a legacy tree is one whose
+    subdirectories really contain ``*.<ext>`` feature files, so an unrelated
+    subdirectory next to ``<layer>.h5`` files does not turn a chunked directory
+    into a legacy one.
+
+    A directory holding both is ambiguous and raises rather than picking one
+    silently -- guessing wrong means reading different data than the caller
+    meant. Pass an explicit `format` to resolve it.
+
+    Parameters
+    ----------
+    dpath : str
+        Feature directory.
+    ext : str, optional
+        Extension of the legacy per-stimulus files (default: ``'mat'``).
+
+    Returns
+    -------
+    str
+        ``'hdf5'`` or ``'mat'``.
+
+    Raises
+    ------
+    RuntimeError
+        If the directory holds both layouts, or neither.
+    """
+    escaped = glob.escape(dpath)
+    has_h5 = bool(glob.glob(os.path.join(escaped, "*." + HDF5_EXT)))
+    has_legacy = any(
+        glob.glob(os.path.join(glob.escape(os.path.join(dpath, d)), "*." + ext))
+        for d in os.listdir(dpath)
+        if os.path.isdir(os.path.join(dpath, d))
+    )
+
+    if has_h5 and has_legacy:
+        raise RuntimeError(
+            "{} holds both chunked .{} files and a legacy .{} tree, so the "
+            "storage format is ambiguous. Pass format='hdf5' or format='mat' "
+            "to choose one.".format(dpath, HDF5_EXT, ext)
+        )
+    if has_h5:
+        return "hdf5"
+    if has_legacy:
+        return "mat"
+    raise RuntimeError(
+        "No features found in {}: expected either <layer>.{} files or "
+        "<layer>/<label>.{} subdirectories.".format(dpath, HDF5_EXT, ext)
+    )
diff --git a/bdpy/dataform/feature_hdf5.py b/bdpy/dataform/feature_hdf5.py
new file mode 100644
index 00000000..02fb8af2
--- /dev/null
+++ b/bdpy/dataform/feature_hdf5.py
@@ -0,0 +1,624 @@
+"""Chunked HDF5 storage for DNN features.
+
+The legacy feature layout writes one ``.mat`` file per stimulus, so the only
+file boundary is the sample axis and a partial read along the feature axes is
+impossible. This module writes the chunked layout instead: one ``<layer>.h5``
+per layer, holding the whole layer in a dataset that is explicitly chunked on
+both the sample axis and the outermost feature axis.
+
+Schema version 1::
+
+    <layer>.h5
+        /features   (n_samples, *feature_shape)
+        /labels     (n_samples,)   variable-length UTF-8
+        attrs:
+            bdpy_format         = "features"
+            bdpy_format_version = 1
+
+One file per layer rather than one file for everything: layers differ in shape
+and in optimal chunk shape, and separate files keep regeneration, copying and
+parallel writing per-layer.
+
+Compression is off by default. The point of this format is partial-read latency,
+and every compressed chunk costs a decompression on the way out; pass
+``compression=`` when size matters more than speed.
+
+This file is a part of BdPy.
+
+Examples
+--------
+Write a whole layer at once::
+
+    save_features('features/conv5.h5', array, labels)
+
+Write incrementally, as feature extraction produces one stimulus at a time::
+
+    with FeatureWriter('features/conv5.h5', feature_shape=(256, 13, 13),
+                       dtype=np.float32) as writer:
+        for label, feature in extract():
+            writer.append(feature, label)
+
+Migrate an existing ``.mat`` tree::
+
+    convert_features_to_hdf5('features_mat', 'features_h5')
+"""
+
+import errno
+import os
+import uuid
+import weakref
+from types import TracebackType
+from typing import Iterable, Optional, Sequence, Set, Tuple, Type
+
+import h5py
+import numpy as np
+
+from ._feature_chunking import DEFAULT_TARGET_CHUNK_BYTES, choose_chunk_shape
+from ._feature_store import (
+    FEATURES_DATASET,
+    FORMAT_ATTR,
+    FORMAT_NAME,
+    FORMAT_VERSION_ATTR,
+    HDF5_EXT,
+    LABELS_DATASET,
+    SUPPORTED_FORMAT_VERSION,
+    MatFeatureStore,
+    _duplicates,
+)
+
+__all__ = [
+    "FeatureWriter",
+    "convert_features_to_hdf5",
+    "save_features",
+]
+
+#: Chunk length of the resizable labels dataset.
+_LABELS_CHUNK = 64
+
+
+def _string_dtype() -> np.dtype:
+    dtype: np.dtype = h5py.string_dtype(encoding="utf-8")
+    return dtype
+
+
+def _write_header(f: h5py.File) -> None:
+    """Stamp the format marker and version. The layer name is the file name."""
+    f.attrs[FORMAT_ATTR] = FORMAT_NAME
+    f.attrs[FORMAT_VERSION_ATTR] = SUPPORTED_FORMAT_VERSION
+
+
+#: Suffix of the scratch file a write is staged in.
+_PARTIAL_SUFFIX = ".partial"
+#: Cap on the copied basename so the staged name stays inside NAME_MAX.
+_MAX_STEM = 64
+
+
+def _staging_path(path: str) -> str:
+    """Name a scratch file next to `path`, invisible to feature-directory scans.
+
+    The name is dotted and does not end in ``.h5``, so neither
+    ``glob('*.h5')`` in :class:`~bdpy.dataform._feature_store.HDF5FeatureStore`
+    nor :func:`~bdpy.dataform._feature_store.detect_format` picks it up as a
+    layer while a write is in flight. The random component lets two processes
+    stage the same layer without colliding.
+    """
+    directory, name = os.path.split(os.path.abspath(path))
+    return os.path.join(
+        directory,
+        ".{}.{}{}".format(name[:_MAX_STEM], uuid.uuid4().hex, _PARTIAL_SUFFIX),
+    )
+
+
+def _prepare_target(path: str, overwrite: bool) -> str:
+    """Refuse to clobber, then name a scratch file next to `path`.
+
+    Writing goes to a scratch file in the *same directory*, so publishing it
+    (see :func:`_promote`) is a single atomic step on the same filesystem:
+    either the finished file appears at `path`, or nothing does. A half-written
+    file must never be left where a reader -- or the converter's skip-if-exists
+    check -- would take it for a complete one.
+
+    This check fails early, before any work is done; :func:`_promote` repeats
+    it when publishing, since another writer may publish `path` meanwhile.
+
+    The caller opens the returned path with mode ``"x"``, which both creates it
+    exclusively and lets HDF5 apply the process umask, so published files keep
+    the permissions they had before staging existed.
+    """
+    if os.path.exists(path) and not overwrite:
+        raise _exists_error(path)
+    _makedirs_for(path)
+    return _staging_path(path)
+
+
+def _exists_error(path: str) -> FileExistsError:
+    return FileExistsError(
+        "{} already exists. Pass overwrite=True to replace it.".format(path)
+    )
+
+
+#: errno values with which os.link reports that the filesystem cannot hard-link.
+_NO_HARDLINK_ERRNOS = frozenset(
+    {errno.EPERM, errno.ENOTSUP, errno.EOPNOTSUPP, errno.ENOSYS}
+)
+
+
+def _promote(tmp_path: str, path: str, overwrite: bool) -> None:
+    """Publish the finished scratch file at `path`.
+
+    With `overwrite`, an atomic rename replaces whatever is there. Without it,
+    the file is published with :func:`os.link`, which fails if `path` exists:
+    checking for a file and publishing ours are one atomic step, so a file that
+    another writer published while this one was staging is never clobbered.
+
+    On a filesystem that cannot hard-link, this falls back to checking for
+    `path` and then renaming, which leaves a short window between the two.
+
+    Raises
+    ------
+    FileExistsError
+        If `path` exists and `overwrite` is False.
+    """
+    if overwrite:
+        os.replace(tmp_path, path)
+        return
+    try:
+        os.link(tmp_path, path)
+    except FileExistsError:
+        raise _exists_error(path) from None
+    except OSError as exc:
+        if exc.errno not in _NO_HARDLINK_ERRNOS:
+            raise
+        if os.path.exists(path):
+            raise _exists_error(path) from None
+        os.replace(tmp_path, path)
+    else:
+        _discard(tmp_path)
+
+
+def _discard(tmp_path: str) -> None:
+    """Remove the scratch file, ignoring an already-vanished one."""
+    try:
+        os.remove(tmp_path)
+    except OSError:
+        pass
+
+
+def save_features(
+    path: str,
+    features: np.ndarray,
+    labels: Sequence[str],
+    chunks: Optional[Tuple[int, ...]] = None,
+    target_chunk_bytes: int = DEFAULT_TARGET_CHUNK_BYTES,
+    compression: Optional[str] = None,
+    dtype: Optional[np.dtype] = None,
+    overwrite: bool = False,
+) -> None:
+    """Write a whole layer to a chunked HDF5 feature file.
+
+    Parameters
+    ----------
+    path : str
+        Output file. An existing file -- including one another writer
+        publishes meanwhile -- is kept unless `overwrite` is True. The write is
+        atomic: on failure nothing is left at `path`.
+    features : numpy.ndarray
+        Feature array of shape ``(n_samples, *feature_shape)``.
+    labels : sequence of str
+        One stimulus label per sample, in the same order as `features`.
+    chunks : tuple of int, optional
+        Explicit chunk shape. Defaults to :func:`choose_chunk_shape`.
+    target_chunk_bytes : int, optional
+        Per-chunk byte budget used when `chunks` is not given.
+    compression : str, optional
+        h5py compression filter, e.g. ``'lzf'`` or ``'gzip'``. ``None``
+        (default) stores the features uncompressed, which keeps partial reads
+        fast.
+    dtype : numpy.dtype, optional
+        Dtype to store. Defaults to the dtype of `features`.
+    overwrite : bool, optional
+        Replace an existing file at `path` (default: False).
+
+    Raises
+    ------
+    ValueError
+        If `labels` does not have one entry per sample, or `features` has fewer
+        than two axes.
+    FileExistsError
+        If `path` exists or is published meanwhile, and `overwrite` is False.
+    """
+    features = np.asarray(features)
+    if features.ndim < 2:
+        raise ValueError(
+            "features need a sample axis and at least one feature axis; "
+            "got shape {}".format(features.shape)
+        )
+    labels = list(labels)
+    if len(labels) != features.shape[0]:
+        raise ValueError(
+            "got {} labels for {} samples".format(len(labels), features.shape[0])
+        )
+
+    # One write path, not two: writing a whole layer is the degenerate case of
+    # writing it incrementally, and keeping them separate is what let a bug
+    # (a zero-sample layer) land in only one of them.
+    with FeatureWriter(
+        path,
+        feature_shape=features.shape[1:],
+        dtype=features.dtype if dtype is None else dtype,
+        n_samples=features.shape[0],
+        chunks=chunks,
+        target_chunk_bytes=target_chunk_bytes,
+        compression=compression,
+        overwrite=overwrite,
+    ) as writer:
+        writer.extend(features, labels)
+
+
+class FeatureWriter:
+    """Incremental writer for a chunked HDF5 feature file.
+
+    Feature extraction produces one stimulus at a time, so the datasets are
+    created resizable (``maxshape=(None, *feature_shape)``) and grown in blocks
+    as samples arrive.
+
+    Writing goes to a temporary file next to `path`; :meth:`close` moves it into
+    place and :meth:`abort` throws it away, so a failed write leaves nothing at
+    `path`. **Prefer the context manager**, which aborts when the block raises --
+    a bare ``try/finally: writer.close()`` would publish a half-written file.
+
+    Parameters
+    ----------
+    path : str
+        Output file. An existing file -- including one another writer
+        publishes meanwhile -- is kept unless `overwrite` is True, and nothing
+        is written there until :meth:`close` succeeds. On a filesystem without
+        hard links, a short window remains in which a file published at the
+        same moment can still be replaced.
+    feature_shape : sequence of int
+        Shape of a single sample's features, without the sample axis.
+    dtype : numpy.dtype
+        Dtype to store.
+    n_samples : int, optional
+        Expected number of samples, if known. Used only to pick a chunk shape.
+    chunks : tuple of int, optional
+        Explicit chunk shape, overriding the byte budget.
+    target_chunk_bytes : int, optional
+        Per-chunk byte budget used when `chunks` is not given.
+    compression : str, optional
+        h5py compression filter. ``None`` (default) stores uncompressed.
+    overwrite : bool, optional
+        Replace an existing file at `path` (default: False).
+
+    Raises
+    ------
+    FileExistsError
+        If `path` exists, or (from :meth:`close`) was published meanwhile,
+        and `overwrite` is False.
+
+    Examples
+    --------
+    >>> with FeatureWriter('conv5.h5', (256, 13, 13), np.float32) as w:  # doctest: +SKIP
+    ...     w.append(feature, 'img0001')
+    """
+
+    def __init__(
+        self,
+        path: str,
+        feature_shape: Sequence[int],
+        dtype: np.dtype,
+        n_samples: Optional[int] = None,
+        chunks: Optional[Tuple[int, ...]] = None,
+        target_chunk_bytes: int = DEFAULT_TARGET_CHUNK_BYTES,
+        compression: Optional[str] = None,
+        overwrite: bool = False,
+    ):
+        self._feature_shape = tuple(int(s) for s in feature_shape)
+        if not self._feature_shape:
+            raise ValueError("feature_shape must have at least one axis")
+        self._dtype = np.dtype(dtype)
+        self._n = 0
+        # Labels identify rows, so they must be unique across the whole file --
+        # not just within one call. A repeated label would make every lookup
+        # resolve to the last row and silently drop the earlier one.
+        self._seen: Set[str] = set()
+
+        if chunks is None:
+            chunks = choose_chunk_shape(
+                (n_samples if n_samples is not None else 0, *self._feature_shape),
+                self._dtype,
+                target_bytes=target_chunk_bytes,
+                n_samples_known=n_samples is not None,
+            )
+        # Fail before doing any work if the target is occupied, then write to
+        # a scratch file so `path` stays untouched until close() succeeds.
+        self._path = path
+        self._overwrite = overwrite
+        self._aborted = False
+        self._tmp_path: Optional[str] = _prepare_target(path, overwrite)
+        self._file: Optional[h5py.File] = h5py.File(self._tmp_path, "x")
+        # If the writer is dropped without close() or abort(), the scratch file
+        # would otherwise linger next to the output. weakref.finalize (rather
+        # than __del__) runs exactly once, is detachable at publish time, and
+        # closes over the path instead of the writer, so it can never resurrect
+        # the object or outlive a name it no longer owns.
+        self._finalizer = weakref.finalize(self, _discard, self._tmp_path)
+        _write_header(self._file)
+        self._features = self._file.create_dataset(
+            FEATURES_DATASET,
+            shape=(0, *self._feature_shape),
+            maxshape=(None, *self._feature_shape),
+            dtype=self._dtype,
+            chunks=chunks,
+            compression=compression,
+        )
+        self._labels = self._file.create_dataset(
+            LABELS_DATASET,
+            shape=(0,),
+            maxshape=(None,),
+            dtype=_string_dtype(),
+            chunks=(_LABELS_CHUNK,),
+        )
+
+    @property
+    def n_samples(self) -> int:
+        """Number of samples written so far."""
+        return self._n
+
+    def append(self, feature: np.ndarray, label: str) -> None:
+        """Append a single sample.
+
+        Parameters
+        ----------
+        feature : numpy.ndarray
+            One sample, either ``feature_shape`` or ``(1, *feature_shape)``.
+        label : str
+            Its stimulus label.
+        """
+        feature = np.asarray(feature)
+        if feature.shape == self._feature_shape:
+            feature = feature[np.newaxis]
+        self.extend(feature, [label])
+
+    def extend(self, features: np.ndarray, labels: Sequence[str]) -> None:
+        """Append several samples at once.
+
+        Parameters
+        ----------
+        features : numpy.ndarray
+            Samples of shape ``(n, *feature_shape)``.
+        labels : sequence of str
+            One label per sample.
+
+        Raises
+        ------
+        ValueError
+            If the feature shape does not match the writer's, or the label count
+            does not match the sample count.
+        RuntimeError
+            If the writer is already closed.
+        """
+        if self._file is None:
+            raise RuntimeError("FeatureWriter is closed")
+        features = np.asarray(features)
+        labels = list(labels)
+        if features.shape[1:] != self._feature_shape:
+            raise ValueError(
+                "expected features of shape (n, {}), got {}".format(
+                    ", ".join(str(s) for s in self._feature_shape), features.shape
+                )
+            )
+        if len(labels) != features.shape[0]:
+            raise ValueError(
+                "got {} labels for {} samples".format(len(labels), features.shape[0])
+            )
+        # Check before writing anything, so a rejected batch leaves the file
+        # exactly as it was.
+        self._reject_duplicates(labels)
+        if not labels:
+            return
+
+        new_n = self._n + len(labels)
+        self._features.resize(new_n, axis=0)
+        self._labels.resize(new_n, axis=0)
+        self._features[self._n:new_n] = features.astype(self._dtype, copy=False)
+        self._labels[self._n:new_n] = labels
+        self._seen.update(labels)
+        self._n = new_n
+
+    def _reject_duplicates(self, labels: Sequence[str]) -> None:
+        """Refuse labels repeated within this batch or already written."""
+        within = _duplicates(labels)
+        if within:
+            raise ValueError(
+                "duplicate labels in this batch: {}. Labels identify rows, so "
+                "they must be unique.".format(
+                    ", ".join(repr(d) for d in within[:3])
+                )
+            )
+        already = [label for label in labels if label in self._seen]
+        if already:
+            raise ValueError(
+                "labels already written to {}: {}. Labels identify rows, so "
+                "they must be unique.".format(
+                    self._path, ", ".join(repr(d) for d in already[:3])
+                )
+            )
+
+    def close(self) -> None:
+        """Finish the file and move it into place. Idempotent.
+
+        This is the commit: call it only once the data is complete, because it
+        publishes whatever has been written so far. To throw a partial file
+        away, call :meth:`abort` instead. ``try: ... finally: writer.close()``
+        is therefore wrong -- it would publish a truncated file -- so use the
+        writer as a context manager, which routes to :meth:`abort` when the
+        body raises.
+
+        Raises
+        ------
+        RuntimeError
+            If the writer was already aborted.
+        FileExistsError
+            If the target was published meanwhile; the writer is then aborted.
+        """
+        if self._aborted:
+            raise RuntimeError(
+                "FeatureWriter was aborted; nothing was written to {}".format(
+                    self._path
+                )
+            )
+        if self._file is not None:
+            self._file.close()
+            self._file = None
+        if self._tmp_path is not None:
+            try:
+                _promote(self._tmp_path, self._path, self._overwrite)
+            except FileExistsError:
+                self.abort()
+                raise
+            self._tmp_path = None
+            self._finalizer.detach()  # the scratch file is now the output
+
+    def abort(self) -> None:
+        """Discard the partial file without touching the target path.
+
+        Idempotent, and a no-op once :meth:`close` has published the file --
+        it must never remove a published output.
+        """
+        if self._file is not None:
+            self._file.close()
+            self._file = None
+        if self._tmp_path is not None:
+            self._finalizer()  # runs _discard exactly once
+            self._tmp_path = None
+            self._aborted = True
+
+    def __enter__(self) -> "FeatureWriter":
+        """Enter the context manager."""
+        return self
+
+    def __exit__(
+        self,
+        exc_type: Optional[Type[BaseException]],
+        exc: Optional[BaseException],
+        tb: Optional[TracebackType],
+    ) -> None:
+        """Publish the file on success, discard it if the block raised."""
+        if self._tmp_path is None:
+            return  # the body already called close() or abort() itself
+        if exc_type is None:
+            self.close()
+        else:
+            self.abort()
+
+
+def convert_features_to_hdf5(
+    src_dir: str,
+    dst_dir: str,
+    layers: Optional[Iterable[str]] = None,
+    ext: str = "mat",
+    key: str = "feat",
+    overwrite: bool = False,
+    batch_size: int = 64,
+    target_chunk_bytes: int = DEFAULT_TARGET_CHUNK_BYTES,
+    compression: Optional[str] = None,
+    verbose: bool = False,
+) -> None:
+    """Convert a legacy per-stimulus feature directory to chunked HDF5.
+
+    Reads `src_dir` through the same backend :class:`~bdpy.dataform.Features`
+    uses for the legacy layout, so the output is by construction what the legacy
+    reader sees. Samples are streamed in batches, so a layer is never held in
+    memory in full.
+
+    Parameters
+    ----------
+    src_dir : str
+        Legacy directory, ``<src_dir>/<layer>/<label>.<ext>``.
+    dst_dir : str
+        Output directory; ``<dst_dir>/<layer>.h5`` is written per layer. Created
+        if missing. Each layer is written atomically, so a failed conversion
+        leaves no file behind and can simply be re-run.
+    layers : iterable of str, optional
+        Layers to convert. Defaults to every layer found in `src_dir`.
+    ext : str, optional
+        Extension of the legacy files (default: ``'mat'``).
+    key : str, optional
+        Variable name inside the legacy files (default: ``'feat'``).
+    overwrite : bool, optional
+        Overwrite an existing output file instead of skipping it
+        (default: False). Without it, a layer that another process publishes
+        while this one is converting it is skipped as well.
+    batch_size : int, optional
+        Number of stimulus files read per batch. Must be positive.
+    target_chunk_bytes : int, optional
+        Per-chunk byte budget.
+    compression : str, optional
+        h5py compression filter. ``None`` (default) stores uncompressed.
+    verbose : bool, optional
+        Print progress.
+
+    Raises
+    ------
+    ValueError
+        If `batch_size` is not positive.
+    KeyError
+        If a requested layer is not present in `src_dir`.
+    """
+    if batch_size < 1:
+        # range(0, n, -1) is empty, so a non-positive batch size would write no
+        # samples at all and then publish that empty file as a finished layer.
+        raise ValueError("batch_size must be positive, got {}".format(batch_size))
+
+    store = MatFeatureStore(src_dir, ext=ext, key=key)
+    selected = list(store.layers) if layers is None else list(layers)
+    missing = [layer for layer in selected if layer not in store.layers]
+    if missing:
+        raise KeyError(
+            "Layer(s) {} not found in {}".format(", ".join(missing), src_dir)
+        )
+
+    os.makedirs(dst_dir, exist_ok=True)
+    all_labels = store.labels
+
+    for layer in selected:
+        out_path = os.path.join(dst_dir, layer + "." + HDF5_EXT)
+        if os.path.exists(out_path) and not overwrite:
+            if verbose:
+                print("{} already exists. Skipped.".format(out_path))
+            continue
+
+        full_shape = store.shape(layer)
+        # The context manager is load-bearing: if a batch fails, it discards the
+        # partial file instead of publishing it. A published partial file would
+        # be taken for a finished one by the skip check above, and the layer
+        # would stay silently truncated across re-runs.
+        try:
+            with FeatureWriter(
+                out_path,
+                feature_shape=full_shape[1:],
+                dtype=store.dtype(layer),
+                n_samples=full_shape[0],
+                target_chunk_bytes=target_chunk_bytes,
+                compression=compression,
+                overwrite=overwrite,
+            ) as writer:
+                for start in range(0, len(all_labels), batch_size):
+                    batch = all_labels[start:start + batch_size]
+                    writer.extend(store.read(layer, batch), batch)
+        except FileExistsError:
+            # Only reachable without overwrite: another process published this
+            # layer while we were converting it. Treat it like the check above;
+            # the writer has already discarded our scratch file.
+            if verbose:
+                print("{} already exists. Skipped.".format(out_path))
+            continue
+
+        if verbose:
+            print("Saved {} ({} samples).".format(out_path, full_shape[0]))
+
+
+def _makedirs_for(path: str) -> None:
+    parent = os.path.dirname(os.path.abspath(path))
+    os.makedirs(parent, exist_ok=True)
diff --git a/bdpy/dataform/features.py b/bdpy/dataform/features.py
index 1a181f5b..610dafc5 100644
--- a/bdpy/dataform/features.py
+++ b/bdpy/dataform/features.py
@@ -15,12 +15,26 @@
 import warnings
 from functools import partial
 from multiprocessing import Pool
-from typing import Any, Dict, List, Optional, Union
+from typing import Any, Dict, Iterator, List, Optional, Sequence, Tuple, Union
 
 import hdf5storage
 import numpy as np
 
 from . import _mat_v73
+from ._feature_store import (
+    FeatureSlice,
+    FeatureStore,
+    HDF5FeatureStore,
+    MatFeatureStore,
+    _check_size,
+    _default_iter_size,
+    _determine_num_parallel,
+    _iter_array_slabs,
+    _load_array_with_key,
+    _normalize_axis,
+    _partition,
+    detect_format,
+)
 
 # Deprecation notice emitted by the MATLAB-compatible write path. The actual
 # switch to bdpy-native plain HDF5 (and the drop of the hdf5storage write
@@ -34,36 +48,32 @@
 )
 
 
-def _load_array_with_key(key: str, path: str) -> np.ndarray:
-    # v5 .mat via scipy, v7.3 (HDF5) via h5py; avoids hdf5storage on the load
-    # path, which breaks under NumPy 2.0 (see bdpy/dataform/_mat_v73.py).
-    return _mat_v73.loadmat_key(path, key)
-
-
-def _determine_num_parallel(num_files: int) -> int:
-    # NOTE: optimal number of parallel processes is not clear. It could depend
-    # on several factors such as the number of files, the size of files, the
-    # number of cores, etc. For now, we use a simple heuristic based on the
-    # number of files.
-    num_parallel: int
-    if num_files < 16:
-        num_parallel = 1
-    elif num_files < 64:
-        num_parallel = 16
-    else:
-        num_parallel = 64
-    return num_parallel
-
-
 class Features(object):
     """DNN features class.
 
+    Reads a collection of DNN features from disk. Two on-disk layouts are
+    supported and behave identically through this API:
+
+    - the legacy per-stimulus layout, ``<dpath>/<layer>/<label>.mat``;
+    - chunked HDF5 storage, ``<dpath>/<layer>.h5`` (see
+      :mod:`bdpy.dataform.feature_hdf5`), which additionally supports partial
+      reads along the feature axes without loading a whole layer.
+
     Parameters
     ----------
     dpath: str or list
         (List of) DNN feature directory(ies)
     ext: str
-        DNN feature file extension (default: mat)
+        DNN feature file extension of the legacy per-stimulus layout
+        (default: mat). Ignored by chunked HDF5 storage.
+    feature_index: str, optional
+        Path to a ``.mat`` file holding a per-layer unit index under the key
+        ``index``. When given, features are flattened and restricted to those
+        units.
+    format: str
+        Storage format: ``'auto'`` (default), ``'mat'`` or ``'hdf5'``. With
+        ``'auto'`` each directory is inspected independently, so a mix of
+        layouts across `dpath` entries is fine.
 
     Attributes
     ----------
@@ -77,17 +87,20 @@ class Features(object):
 
     def __init__(
             self, dpath: Union[str, List[str]] = [],
-            ext: str = 'mat', feature_index: Optional[str] = None
+            ext: str = 'mat', feature_index: Optional[str] = None,
+            format: str = 'auto'
         ):
         if not isinstance(dpath, list):
             dpath = [dpath]
         self.__dpath = dpath
+        self.__format = format
 
-        self.__feature_file_table: Dict[str, Dict[str, str]] = {}  # Stimulus feature file tables
+        self.__stores: List[FeatureStore] = []
+        self.__label_store: Dict[str, FeatureStore] = {}  # label -> owning store
         self.__labels: List[str] = []  # Stimulus labels
         self.__index: List[int] = []  # Stimulus index (one-based)
         self.__layers: List[str] = []  # DNN layers
-        self.__collect_feature_files(ext=ext)
+        self.__collect_features(ext=ext)
 
         self.__c_feature_name: Optional[str] = None  # Loaded layer
         self.__features: Optional[np.ndarray] = None  # Loaded features
@@ -126,7 +139,26 @@ def layers(self):
     def feature_index(self):
         return self.__feature_index
 
-    def get(self, layer: str, label: Union[str, List[str], None] = None) -> np.ndarray:
+    def shape(self, layer: str) -> Tuple[int, ...]:
+        """Return the full shape of `layer` without reading the features.
+
+        Parameters
+        ----------
+        layer: str
+            DNN layer
+
+        Returns
+        -------
+        tuple of int
+            ``(n_samples, *feature_shape)``
+        """
+        feature_shape = self.__stores[0].shape(layer)[1:]
+        return (len(self.__labels), *feature_shape)
+
+    def get(
+            self, layer: str, label: Union[str, List[str], None] = None,
+            feature_slice: FeatureSlice = None
+        ) -> np.ndarray:
         """Return features in `layer`.
 
         Parameters
@@ -134,40 +166,158 @@ def get(self, layer: str, label: Union[str, List[str], None] = None) -> np.ndarr
         layer: str
             DNN layer
         label: str or list
-            Sample label(s)
+            Sample label(s). Rows come back in the order given.
+        feature_slice: slice, int, Ellipsis or tuple of those, optional
+            Index applied to the feature axes (axes 1 and up), as produced by
+            ``numpy.s_[...]``. Basic forward indexing only: slices with a
+            positive step and integer bounds, integers, and at most one
+            ``Ellipsis``. Fancy indexing, a negative step, a non-integer slice
+            bound, booleans, and ``numpy.newaxis`` *inside a tuple* raise
+            ``ValueError``. A bare ``numpy.newaxis`` is ``None``, the "no
+            slice" default, so it reads the whole feature tensor. With chunked
+            HDF5 storage this is a genuine partial read; with the legacy layout
+            the files are loaded in full and then sliced.
 
         Returns
         -------
         numpy.ndarray, shape=(n_samples, shape_layers)
             DNN features
+
+        Examples
+        --------
+        >>> features.get('conv5')                                  # doctest: +SKIP
+        >>> features.get('conv5', label=['img0001', 'img0002'])     # doctest: +SKIP
+        >>> features.get('conv5', feature_slice=np.s_[128:256])     # doctest: +SKIP
+
+        Raises
+        ------
+        ValueError
+            If `feature_slice` is combined with a unit index (`feature_index`).
+            The index addresses the flattened *full* feature space, so it is
+            meaningless against an already-sliced array.
         """
-        if label is None:
+        if feature_slice is not None and self.__feat_index_table is not None:
+            # The unit index addresses the flattened full feature space, while
+            # feature_slice has already narrowed it. Applying one to the other
+            # would silently select the wrong units (or raise IndexError),
+            # so refuse the combination instead of guessing.
+            raise ValueError(
+                'feature_slice cannot be combined with feature_index: the unit '
+                'index addresses the full feature space, not a slice of it. '
+                'Call get() without feature_slice and slice the result, or '
+                'construct Features without feature_index.'
+            )
+
+        if label is None and feature_slice is None:
             return self.get_features(layer)
 
-        if isinstance(label, str):
+        if label is None:
+            labels = self.__labels
+        elif isinstance(label, str):
             labels = [label]
         else:
-            labels = label
+            labels = list(label)
 
-        features: np.ndarray
-        num_labels = len(labels)
-        num_parallel = _determine_num_parallel(num_labels)
-        path_iterator = map(lambda label: self.__feature_file_table[layer][label], labels)
-        if num_parallel == 1:
-            features = np.concatenate(list(map(partial(_load_array_with_key, 'feat'), path_iterator)), axis=0)
-        else:
-            with Pool(processes=num_parallel) as pool:
-                features = np.concatenate(pool.map(partial(_load_array_with_key, 'feat'), path_iterator), axis=0)
+        features = self.__read(layer, labels, feature_slice)
+        return self.__apply_feature_index(features, layer)
 
-        if self.__feat_index_table is not None:
-            # Select features by index
-            self.__feature_index = self.__feat_index_table[layer]
-            n_sample = features.shape[0]
-            n_feat = np.array(features.shape[1:]).prod()
+    def iter_chunks(
+            self, layer: str, label: Union[str, List[str], None] = None,
+            axis: int = 1, size: Optional[int] = None
+        ) -> Iterator[Tuple[slice, np.ndarray]]:
+        """Iterate over `layer` in slabs along `axis`.
 
-            features = features.reshape([n_sample, n_feat], order='C')[:, self.__feature_index]
+        This is the streaming counterpart of :meth:`get`: with chunked HDF5
+        storage it holds only one slab at a time, so a layer far larger than RAM
+        can be processed end to end.
 
-        return features
+        There is no `feature_slice` here on purpose -- slice the blocks as they
+        come out instead. Composing an arbitrary slice with the per-slab slice
+        means reimplementing NumPy's index arithmetic, which is more generality
+        than reading a layer in slabs needs.
+
+        Parameters
+        ----------
+        layer: str
+            DNN layer
+        label: str or list, optional
+            Sample label(s). ``None`` iterates over every label.
+        axis: int
+            Axis of the selected array to iterate over. Axis 0 is the sample
+            axis; the default, axis 1, is the outermost feature axis.
+        size: int, optional
+            Elements per slab. Defaults to the on-disk chunk extent along
+            `axis`, so that each element is read exactly once.
+
+        Notes
+        -----
+        The one-slab-at-a-time guarantee is a property of chunked HDF5 storage
+        and of `axis=0` on any layout. The legacy `.mat` layout has no file
+        boundary on the feature axes, so iterating one reads the selection once
+        -- the same peak as a single :meth:`get` -- and yields views into it,
+        which keep that array alive for as long as a block is held. It is still
+        read once rather than once per slab.
+
+        The guarantee also assumes the requested labels all live in a single
+        feature directory, which is the usual case. When they are spread across
+        several `dpath` entries, no single store can stream them, so the
+        selection is read in full and then split; without `size`, the slabs
+        then follow the byte budget, as on a store without on-disk chunking.
+
+        Yields
+        ------
+        tuple of (slice, numpy.ndarray)
+            The slice applied to `axis`, and the corresponding slab.
+
+        Raises
+        ------
+        ValueError
+            If a unit index (`feature_index`) is in use, which flattens the
+            feature axes and so has no meaningful per-axis iteration; if `axis`
+            is out of range; or if `size` is given and is not positive.
+
+        Examples
+        --------
+        >>> for sl, block in features.iter_chunks('conv5'):  # doctest: +SKIP
+        ...     out[:, sl] = transform(block)
+        """
+        if self.__feat_index_table is not None:
+            raise ValueError(
+                'iter_chunks is not supported together with feature_index, '
+                'which flattens the feature axes'
+            )
+
+        if label is None:
+            labels = self.__labels
+        elif isinstance(label, str):
+            labels = [label]
+        else:
+            labels = list(label)
+
+        # Validate up front with the same helpers the stores use, so that both
+        # dispatch paths below reject the same inputs the same way, and the
+        # fallback rejects them before reading the selection.
+        selected_shape = (len(labels), *self.shape(layer)[1:])
+        axis = _normalize_axis(axis, selected_shape)
+        if size is not None:
+            _check_size(size)
+
+        store = self.__store_for(labels)
+        if store is not None:
+            yield from store.iter_chunks(layer, labels, axis=axis, size=size)
+            return
+
+        # Labels span several directories, so no single store can stream them.
+        # Fall back to slicing a full read; the default size follows the same
+        # byte budget as a store without on-disk chunking.
+        features = self.__read(layer, labels)
+        if size is None:
+            size = _default_iter_size(
+                features.shape, axis, features.dtype.itemsize, extent=None
+            )
+        yield from _iter_array_slabs(
+            features, axis, _partition(features.shape[axis], size)
+        )
 
     def statistic(self, statistic: str = 'mean', layer: Optional[str] = None):
         # NOTE: return type is ambiguous. currently, it is inferred as Unkown | Any
@@ -220,79 +370,84 @@ def get_features(self, layer: str) -> np.ndarray:
             assert isinstance(self.__features, np.ndarray)
             return self.__features  # self.__features could be None
 
-        features: np.ndarray
-        num_labels = len(self.__labels)
-        num_parallel = _determine_num_parallel(num_labels)
-        path_iterator = map(lambda label: self.__feature_file_table[layer][label], self.__labels)
-        if num_parallel == 1:
-            features = np.concatenate(list(map(partial(_load_array_with_key, 'feat'), path_iterator)), axis=0)
-        else:
-            with Pool(processes=num_parallel) as pool:
-                features = np.concatenate(pool.map(partial(_load_array_with_key, 'feat'), path_iterator), axis=0)
-        self.__features = features
-
+        self.__features = self.__read(layer, self.__labels)
         self.__c_feature_name = layer
-
-        if self.__feat_index_table is not None:
-            # Select features by index
-            self.__feature_index = self.__feat_index_table[layer]
-            n_sample = self.__features.shape[0]
-            n_feat = np.array(self.__features.shape[1:]).prod()
-
-            self.__features = self.__features.reshape([n_sample, n_feat], order='C')[:, self.__feature_index]
+        self.__features = self.__apply_feature_index(self.__features, layer)
 
         return self.__features
 
-    def __collect_feature_files(self, ext='mat'):
-        dpath_lst = self.__dpath
-
-        # List-up layers and stimulus labels
-        label_dir = {}
-        for dpath in dpath_lst:
-            # List-up layers
-            self.__layers = self.__get_layers(dpath)
+    def __read(
+            self, layer: str, labels: Sequence[str],
+            feature_slice: FeatureSlice = None
+        ) -> np.ndarray:
+        """Read `labels` from `layer`, dispatching to the owning store(s)."""
+        store = self.__store_for(labels)
+        if store is not None:
+            return store.read(layer, labels, feature_slice)
+        # Labels are spread over several directories: read each run of
+        # consecutive labels from its own store, preserving the caller's order.
+        blocks: List[np.ndarray] = []
+        run: List[str] = []
+        run_store: Optional[FeatureStore] = None
+        for label in labels:
+            owner = self.__label_store[label]
+            if owner is not run_store and run:
+                assert run_store is not None
+                blocks.append(run_store.read(layer, run, feature_slice))
+                run = []
+            run_store = owner
+            run.append(label)
+        if run:
+            assert run_store is not None
+            blocks.append(run_store.read(layer, run, feature_slice))
+        return np.concatenate(blocks, axis=0)
+
+    def __store_for(self, labels: Sequence[str]) -> Optional[FeatureStore]:
+        """The single store owning every label, or None if they are spread out."""
+        if len(self.__stores) == 1:
+            return self.__stores[0]
+        owners = {id(self.__label_store[label]) for label in labels}
+        if len(owners) == 1:
+            return self.__label_store[labels[0]]
+        return None
 
-            # List-up stimulus labels
-            labels_in_dir = self.__get_labels(dpath, self.__layers, ext=ext)
-            label_dir.update({label: dpath for label in labels_in_dir})
-            self.__labels += labels_in_dir
+    def __apply_feature_index(self, features: np.ndarray, layer: str) -> np.ndarray:
+        if self.__feat_index_table is None:
+            return features
+        # Select features by index
+        self.__feature_index = self.__feat_index_table[layer]
+        n_sample = features.shape[0]
+        n_feat = np.array(features.shape[1:]).prod()
+        return features.reshape([n_sample, n_feat], order='C')[:, self.__feature_index]
+
+    def __collect_features(self, ext='mat'):
+        for dpath in self.__dpath:
+            store = self.__make_store(dpath, ext)
+            if self.__layers and store.layers != self.__layers:
+                raise RuntimeError('Invalid layers in %s' % dpath)
+            self.__layers = store.layers
+            self.__stores.append(store)
+            for label in store.labels:
+                self.__label_store[label] = store
+            self.__labels += store.labels
 
         # NOTE: type incompatibility here. Is it OK to cast to list?
         self.__index = np.arange(len(self.__labels)) + 1
 
-        # List-up feature files
-        for lay in self.__layers:
-            self.__feature_file_table.update(
-                {
-                    lay:
-                    {
-                        label:
-                         os.path.join(label_dir[label], lay, label + '.' + ext)
-                         for label in self.__labels
-                    }
-                })
-
         return None
 
-    def __get_layers(self, dpath: str):
-        layers = sorted([d for d in os.listdir(dpath) if os.path.isdir(os.path.join(dpath, d))])
-        if self.__layers and (layers != self.__layers):
-            raise RuntimeError('Invalid layers in %s' % dpath)
-        return layers
-
-    def __get_labels(self, dpath: str, layers: List[str], ext: str = 'mat'):
-        labels: List[str] = []
-        for lay in layers:
-            lay_dir = os.path.join(dpath, lay)
-            lay_dir = lay_dir.replace('[', '[[]') # Use glob.escape for Python 3.4 or later
-            files = glob.glob(os.path.join(lay_dir, '*.' + ext))
-            labels_t = sorted([os.path.splitext(os.path.basename(f))[0] for f in files])
-            if not labels:
-                labels = labels_t
-            else:
-                if labels != labels_t:
-                    raise RuntimeError('Invalid feature file in %s ' % dpath)
-        return labels
+    def __make_store(self, dpath: str, ext: str) -> FeatureStore:
+        fmt = self.__format
+        if fmt == 'auto':
+            fmt = detect_format(dpath, ext=ext)
+        if fmt == 'hdf5':
+            return HDF5FeatureStore(dpath)
+        if fmt == 'mat':
+            return MatFeatureStore(dpath, ext=ext)
+        raise ValueError(
+            "Unknown feature storage format {!r}; expected "
+            "'auto', 'mat' or 'hdf5'".format(fmt)
+        )
 
 
 class DecodedFeatures(object):
diff --git a/docs/dataform_features.md b/docs/dataform_features.md
index cd07c19c..cccb82a1 100644
--- a/docs/dataform_features.md
+++ b/docs/dataform_features.md
@@ -51,4 +51,181 @@ decoded_features.statistic('mean', layer='fc8', subject='sub-01', roi='VC')
 
 ```
 
+## Chunked HDF5 feature storage
 
+The layout above stores one file per stimulus, so the only file boundary is the
+sample axis: reading a few channels of a layer still costs a full read of every
+stimulus. For large DNN features, `Features` can instead read **chunked HDF5**
+storage, one file per layer:
+
+```
+features/
+  conv1_1.h5    # /features (n_stimuli, *feature_shape) + /labels (n_stimuli,)
+  conv2_1.h5
+  fc8.h5
+```
+
+`/features` is explicitly chunked along both the sample axis and the outermost
+feature axis, so a slice reads only the chunks it covers.
+
+### Reading
+
+Nothing changes for existing code -- `Features` detects the layout per
+directory, and `.mat` trees keep working exactly as before:
+
+``` python
+features = Features('/path/to/features')   # either layout
+feat = features.get(layer='conv5')
+```
+
+Detection keys off the files that are actually present, so an unrelated
+subdirectory next to `<layer>.h5` files does not make a directory look legacy.
+A directory holding *both* layouts is ambiguous and raises rather than picking
+one silently; pass `format='mat'` or `format='hdf5'` to resolve it.
+
+### Partial reads
+
+`feature_slice` selects along the feature axes (axis 1 and up). On chunked HDF5
+this is a real partial read; on the legacy layout the files are loaded in full
+and then sliced, so the same code works on both.
+
+``` python
+import numpy as np
+
+# Channels 128-255 only, without loading the rest of the layer
+feat = features.get(layer='conv5', feature_slice=np.s_[128:256])
+
+# Combine with label selection; rows come back in the order given
+feat = features.get(
+    layer='conv5',
+    label=['stimulus-0003', 'stimulus-0001'],
+    feature_slice=np.s_[128:256],
+)
+
+# NOTE: feature_slice cannot be combined with a unit index (`feature_index`).
+# The index addresses the flattened full feature space, so applying it to an
+# already-sliced array would select the wrong units; the combination raises
+# ValueError instead.
+```
+
+`feature_slice` is **basic forward indexing**: slices with a positive step
+(negative `start`/`stop` are fine), integers, a single `Ellipsis`, and tuples of
+those. Anything else -- a negative step, a non-integer slice bound, fancy
+indexing with a list or array, booleans, or `np.newaxis` *inside a tuple* --
+raises `ValueError` on *both* backends. (A bare `np.newaxis` is simply `None`,
+which is the "no slice" default, so it reads the whole feature tensor rather
+than raising.) The restriction is what lets the two layouts mean the same thing
+by the same index; read without `feature_slice` and index the result with NumPy
+when you need more.
+
+``` python
+features.get('conv5', feature_slice=np.s_[128:256])     # ok
+features.get('conv5', feature_slice=np.s_[8:16, 1:4])   # ok
+features.get('conv5', feature_slice=np.s_[::-1])        # ValueError
+features.get('conv5', feature_slice=np.s_[[3, 1, 7]])   # ValueError
+features.get('conv5', feature_slice=slice(1.5, 3))      # ValueError
+
+# Full shape without reading anything
+n_stimuli, *feature_shape = features.shape('conv5')
+```
+
+To stream a layer that does not fit in memory, iterate. `iter_chunks` yields
+`(slice, block)` so results can be placed back without tracking offsets, and by
+default uses the on-disk chunk extent, so every element is read exactly once:
+
+``` python
+out = np.empty(features.shape('conv5'))
+for sl, block in features.iter_chunks('conv5', axis=1):
+    out[:, sl] = transform(block)
+```
+
+`axis=0` iterates over stimuli instead of features. `iter_chunks` takes no
+`feature_slice` -- slice each `block` as it comes out instead.
+
+Only one slab is resident at a time for chunked HDF5, and for `axis=0` on either
+layout. On the legacy `.mat` layout there is no file boundary on the feature
+axes, so iterating one reads the selection once -- the same peak as a single
+`get()` -- and yields views into it, which keep that array alive for as long as a
+block is held. It is still read once rather than once per slab, which is what
+makes `iter_chunks` usable on that layout at all. The same full-read fallback
+applies when the requested labels are spread across several `dpath` entries,
+since no single store can stream them.
+
+### Writing
+
+Write a whole layer at once:
+
+``` python
+from bdpy.dataform import save_features
+
+save_features('features/conv5.h5', array, labels)
+```
+
+Writing never clobbers: an existing file raises `FileExistsError` unless you pass
+`overwrite=True`. Writes are also atomic -- the file is built in a hidden
+`.<name>.<id>.partial` sibling and moved into place only once it is complete --
+so a failed write leaves nothing behind and can simply be re-run, and an
+overwrite keeps the old file readable until the moment it is replaced. (A
+process killed outright cannot clean up after itself, so a stray `.partial` file
+may be left for you to delete.)
+
+Or incrementally, which is what feature extraction needs since it produces one
+stimulus at a time:
+
+``` python
+from bdpy.dataform import FeatureWriter
+
+with FeatureWriter('features/conv5.h5', feature_shape=(256, 13, 13),
+                   dtype=np.float32) as writer:
+    for label, feature in extract():
+        writer.append(feature, label)
+```
+
+### Migrating existing features
+
+``` python
+from bdpy.dataform import convert_features_to_hdf5
+
+convert_features_to_hdf5('/path/to/features_mat', '/path/to/features_h5')
+```
+
+The converter reads through the same backend `Features` uses for the legacy
+layout, and streams in batches, so a layer is never held in memory in full.
+Existing `<layer>.h5` files are skipped unless `overwrite=True`, and a layer that
+fails to convert leaves no file behind, so re-running picks up where it stopped.
+
+All layers in a directory must hold the same stimulus labels in the same order,
+and labels must be unique within a layer; both are rejected when the directory
+is opened, as they are for the legacy layout (where the file name *is* the
+label, so duplicates cannot arise). Writing a duplicate label raises too.
+
+### Chunk shape
+
+Chunk shape, not the file format alone, is what decides the cost of a partial
+read: HDF5 reads whole chunks even when only a few of their elements are
+selected. By default bdpy derives it from a 1 MiB budget, splitting it between
+the sample axis and the outermost feature axis and keeping trailing spatial axes
+whole -- e.g. `(1200, 256, 13, 13)` float32 becomes chunks of
+`(39, 37, 13, 13)`. Extents are chosen to divide each axis as evenly as
+possible, so HDF5 does not pad the edge chunks.
+
+That means reading a *single* stimulus costs one chunk row (~39 stimuli of I/O),
+which matters for `FeaturesDataset`-style per-sample access. Tune it if your
+access pattern is skewed:
+
+``` python
+save_features(path, array, labels, target_chunk_bytes=256 * 1024)  # smaller chunks
+save_features(path, array, labels, chunks=(1, 256, 13, 13))        # per-stimulus
+```
+
+Compression is off by default, since every compressed chunk costs a
+decompression on the way out. Note that this makes the new format *larger* on
+disk than a legacy `.mat` tree, whose files `hdf5storage` compresses by default.
+Pass `compression='gzip'` or `compression='lzf'` when size matters more than
+read speed.
+
+### Format
+
+Files are marked with root attributes `bdpy_format="features"` and
+`bdpy_format_version=1`, and are rejected on read if those are missing or
+newer than the running bdpy understands.
diff --git a/tests/dataform/test_feature_chunking.py b/tests/dataform/test_feature_chunking.py
new file mode 100644
index 00000000..3dfb1b20
--- /dev/null
+++ b/tests/dataform/test_feature_chunking.py
@@ -0,0 +1,124 @@
+import unittest
+
+import numpy as np
+
+from bdpy.dataform._feature_chunking import (
+    DEFAULT_TARGET_CHUNK_BYTES,
+    choose_chunk_shape,
+)
+
+
+def _nbytes(chunk, dtype):
+    return int(np.prod(chunk)) * np.dtype(dtype).itemsize
+
+
+class TestChooseChunkShape(unittest.TestCase):
+    """Chunk shape decides partial-read cost, so pin the policy down."""
+
+    shapes = [
+        (1200, 1000),
+        (1200, 256, 13, 13),
+        (50, 1000),
+        (1, 256, 13, 13),
+        (3, 4096, 7, 7),
+        (1000000, 1),
+        (2, 1000000),
+        (17, 63, 5),
+    ]
+
+    def test_invariants(self):
+        for dtype in (np.float32, np.float64, np.int16):
+            for shape in self.shapes:
+                chunk = choose_chunk_shape(shape, np.dtype(dtype))
+                with self.subTest(shape=shape, dtype=dtype):
+                    self.assertEqual(len(chunk), len(shape))
+                    # Never zero; never larger than the dataset on axes that
+                    # actually have a size. A zero-sample layer still needs a
+                    # chunk extent of at least 1 -- HDF5 forbids a zero dim.
+                    for c, s in zip(chunk, shape):
+                        self.assertGreaterEqual(c, 1)
+                        if s > 0:
+                            self.assertLessEqual(c, s)
+
+    def test_chunk_stays_within_budget_when_the_budget_is_reachable(self):
+        # The budget is an upper bound only when some chunk can meet it. The
+        # current policy keeps the trailing axes whole, so a shape whose
+        # trailing axes alone exceed the budget cannot: assert only where the
+        # smallest chunk that policy can produce does fit.
+        dtype = np.dtype(np.float32)
+        for shape in self.shapes:
+            chunk = choose_chunk_shape(shape, dtype)
+            smallest = _nbytes((1, 1, *shape[2:]), dtype)
+            with self.subTest(shape=shape):
+                if smallest <= DEFAULT_TARGET_CHUNK_BYTES:
+                    self.assertLessEqual(
+                        _nbytes(chunk, dtype), DEFAULT_TARGET_CHUNK_BYTES
+                    )
+
+    def test_single_sample(self):
+        chunk = choose_chunk_shape((1, 4096, 7, 7), np.dtype(np.float32))
+        self.assertEqual(chunk[0], 1)
+
+    def test_target_bytes_is_respected(self):
+        dtype = np.dtype(np.float32)
+        shape = (1200, 256, 13, 13)
+        small = choose_chunk_shape(shape, dtype, target_bytes=64 * 1024)
+        large = choose_chunk_shape(shape, dtype, target_bytes=4 * 1024 * 1024)
+        self.assertLessEqual(_nbytes(small, dtype), 64 * 1024)
+        self.assertLessEqual(_nbytes(large, dtype), 4 * 1024 * 1024)
+        self.assertLess(_nbytes(small, dtype), _nbytes(large, dtype))
+
+    def test_unknown_sample_count_ignores_the_current_size(self):
+        # The contract of n_samples_known=False is exactly that the current
+        # sample-axis size is not consulted, so two calls that differ only in
+        # that size must agree. Stated this way the test does not depend on
+        # which extent the policy happens to pick.
+        dtype = np.dtype(np.float32)
+        empty = choose_chunk_shape((0, 256, 13, 13), dtype, n_samples_known=False)
+        huge = choose_chunk_shape(
+            (10 ** 6, 256, 13, 13), dtype, n_samples_known=False
+        )
+        self.assertEqual(empty, huge)
+
+    def test_both_sliceable_axes_are_chunked(self):
+        # The promise of issue #144: reading a slice of channels must not have
+        # to read every channel, which requires the feature axis to be split as
+        # well as the sample axis. Asserted on a representative layer, since a
+        # small one legitimately fits in a single chunk.
+        shape = (1200, 256, 13, 13)
+        chunk = choose_chunk_shape(shape, np.dtype(np.float32))
+        self.assertLess(chunk[0], shape[0])
+        self.assertLess(chunk[1], shape[1])
+
+    def test_deterministic(self):
+        shape = (1200, 256, 13, 13)
+        dtype = np.dtype(np.float32)
+        self.assertEqual(
+            choose_chunk_shape(shape, dtype), choose_chunk_shape(shape, dtype)
+        )
+
+    def test_zero_sample_layer_still_gets_a_usable_chunk(self):
+        # HDF5 forbids a zero chunk dimension, so an empty layer needs >= 1.
+        chunk = choose_chunk_shape((0, 8), np.dtype(np.float32))
+        self.assertEqual(chunk[0], 1)
+        self.assertGreaterEqual(chunk[1], 1)
+
+    def test_rejects_empty_feature_axis(self):
+        for bad in ((5, 0), (5, 3, 0)):
+            with self.subTest(shape=bad):
+                with self.assertRaises(ValueError):
+                    choose_chunk_shape(bad, np.dtype(np.float32))
+
+    def test_rejects_negative_sample_axis(self):
+        with self.assertRaises(ValueError):
+            choose_chunk_shape((-1, 4), np.dtype(np.float32))
+
+    def test_rejects_bad_input(self):
+        with self.assertRaises(ValueError):
+            choose_chunk_shape((100,), np.dtype(np.float32))
+        with self.assertRaises(ValueError):
+            choose_chunk_shape((100, 10), np.dtype(np.float32), target_bytes=0)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/tests/dataform/test_feature_hdf5.py b/tests/dataform/test_feature_hdf5.py
new file mode 100644
index 00000000..8a30b1c1
--- /dev/null
+++ b/tests/dataform/test_feature_hdf5.py
@@ -0,0 +1,695 @@
+import errno
+import gc
+import glob
+import os
+import stat
+import tempfile
+import unittest
+import warnings
+from unittest import mock
+
+import h5py
+import numpy as np
+from numpy.testing import assert_array_equal
+
+from bdpy.dataform import Features
+from bdpy.dataform._feature_store import (
+    FORMAT_ATTR,
+    FORMAT_NAME,
+    FORMAT_VERSION_ATTR,
+    SUPPORTED_FORMAT_VERSION,
+    HDF5FeatureStore,
+    MatFeatureStore,
+)
+from bdpy.dataform.feature_hdf5 import (
+    FeatureWriter,
+    convert_features_to_hdf5,
+    save_features,
+)
+
+from .test_features import prepare_mat_features
+
+
+def _scratch_files(directory):
+    return [n for n in os.listdir(directory) if n.endswith('.partial')]
+
+
+def _current_umask():
+    mask = os.umask(0)
+    os.umask(mask)
+    return mask
+
+
+class TestSaveFeatures(unittest.TestCase):
+    def setUp(self):
+        self.tmpdir = tempfile.TemporaryDirectory()
+        self.labels = ['img%04d' % i for i in range(12)]
+        self.data = np.random.rand(12, 32, 5, 5).astype(np.float32)
+        self.path = os.path.join(self.tmpdir.name, 'conv5.h5')
+
+    def tearDown(self):
+        self.tmpdir.cleanup()
+
+    def test_writes_schema_v1(self):
+        save_features(self.path, self.data, self.labels)
+
+        with h5py.File(self.path, 'r') as f:
+            self.assertEqual(f.attrs[FORMAT_ATTR], FORMAT_NAME)
+            self.assertEqual(
+                int(f.attrs[FORMAT_VERSION_ATTR]), SUPPORTED_FORMAT_VERSION
+            )
+            assert_array_equal(f['features'][()], self.data)
+            self.assertEqual(
+                [s.decode('utf-8') for s in f['labels'][()]], self.labels
+            )
+
+    def test_features_dataset_is_chunked_and_uncompressed(self):
+        # Chunking is the whole point; compression would tax every partial read.
+        save_features(self.path, self.data, self.labels)
+        with h5py.File(self.path, 'r') as f:
+            self.assertIsNotNone(f['features'].chunks)
+            self.assertIsNone(f['features'].compression)
+
+    def test_explicit_chunks_and_compression(self):
+        save_features(
+            self.path, self.data, self.labels,
+            chunks=(4, 8, 5, 5), compression='gzip',
+        )
+        with h5py.File(self.path, 'r') as f:
+            self.assertEqual(f['features'].chunks, (4, 8, 5, 5))
+            self.assertEqual(f['features'].compression, 'gzip')
+        # Compression must not change what comes back out.
+        assert_array_equal(
+            HDF5FeatureStore(self.tmpdir.name).read('conv5'), self.data
+        )
+
+    def test_dtype_is_preserved(self):
+        for dtype in (np.float32, np.float64, np.int32):
+            with self.subTest(dtype=dtype):
+                path = os.path.join(self.tmpdir.name, 'l_%s.h5' % np.dtype(dtype).name)
+                save_features(path, self.data.astype(dtype), self.labels)
+                with h5py.File(path, 'r') as f:
+                    self.assertEqual(f['features'].dtype, np.dtype(dtype))
+
+    def test_rejects_mismatched_labels(self):
+        with self.assertRaises(ValueError):
+            save_features(self.path, self.data, self.labels[:-1])
+
+    def test_rejects_1d_features(self):
+        with self.assertRaises(ValueError):
+            save_features(self.path, np.arange(5), ['a'] * 5)
+
+
+class TestFeatureWriter(unittest.TestCase):
+    def setUp(self):
+        self.tmpdir = tempfile.TemporaryDirectory()
+        self.path = os.path.join(self.tmpdir.name, 'conv5.h5')
+        self.labels = ['img%04d' % i for i in range(150)]
+        self.data = np.random.rand(150, 16, 3, 3).astype(np.float32)
+
+    def tearDown(self):
+        self.tmpdir.cleanup()
+
+    def test_append_round_trip(self):
+        # 150 samples crosses several resize blocks, which is where an
+        # incremental writer would lose data if the resize were wrong.
+        with FeatureWriter(self.path, (16, 3, 3), np.float32) as writer:
+            for label, feature in zip(self.labels, self.data):
+                writer.append(feature, label)
+            self.assertEqual(writer.n_samples, len(self.labels))
+
+        store = HDF5FeatureStore(self.tmpdir.name)
+        self.assertEqual(store.labels, self.labels)
+        assert_array_equal(store.read('conv5'), self.data)
+
+    def test_append_accepts_leading_sample_axis(self):
+        with FeatureWriter(self.path, (16, 3, 3), np.float32) as writer:
+            writer.append(self.data[0][np.newaxis], self.labels[0])
+            writer.append(self.data[1], self.labels[1])
+        assert_array_equal(
+            HDF5FeatureStore(self.tmpdir.name).read('conv5'), self.data[:2]
+        )
+
+    def test_extend(self):
+        with FeatureWriter(self.path, (16, 3, 3), np.float32) as writer:
+            writer.extend(self.data[:40], self.labels[:40])
+            writer.extend(self.data[40:], self.labels[40:])
+        assert_array_equal(
+            HDF5FeatureStore(self.tmpdir.name).read('conv5'), self.data
+        )
+
+    def test_writes_valid_header(self):
+        with FeatureWriter(self.path, (16, 3, 3), np.float32) as writer:
+            writer.append(self.data[0], self.labels[0])
+        with h5py.File(self.path, 'r') as f:
+            self.assertEqual(f.attrs[FORMAT_ATTR], FORMAT_NAME)
+            self.assertIsNotNone(f['features'].chunks)
+
+    def test_rejects_wrong_shape(self):
+        with FeatureWriter(self.path, (16, 3, 3), np.float32) as writer:
+            with self.assertRaises(ValueError):
+                writer.extend(np.zeros((2, 8, 3, 3)), ['a', 'b'])
+
+    def test_rejects_label_count_mismatch(self):
+        with FeatureWriter(self.path, (16, 3, 3), np.float32) as writer:
+            with self.assertRaises(ValueError):
+                writer.extend(np.zeros((2, 16, 3, 3)), ['a'])
+
+    def test_write_after_close_raises(self):
+        writer = FeatureWriter(self.path, (16, 3, 3), np.float32)
+        writer.close()
+        writer.close()  # idempotent
+        with self.assertRaises(RuntimeError):
+            writer.append(self.data[0], self.labels[0])
+
+
+class TestConvertFeaturesToHDF5(unittest.TestCase):
+    """The converter must reproduce exactly what the legacy reader sees."""
+
+    def setUp(self):
+        warnings.simplefilter('ignore', FutureWarning)
+        self.tmpdir = tempfile.TemporaryDirectory()
+        self.matdir = os.path.join(self.tmpdir.name, 'mat')
+        self.h5dir = os.path.join(self.tmpdir.name, 'h5')
+        os.makedirs(self.matdir)
+        self.labels = ['img%04d' % i for i in range(10)]
+        self.stacked = prepare_mat_features(
+            self.matdir, ['conv5', 'fc8'], self.labels,
+            [(1, 16, 3, 3), (1, 50)],
+        )
+
+    def tearDown(self):
+        self.tmpdir.cleanup()
+
+    def test_converted_features_match_source(self):
+        convert_features_to_hdf5(self.matdir, self.h5dir)
+
+        from_mat = Features(self.matdir)
+        from_h5 = Features(self.h5dir)
+        self.assertEqual(from_h5.layers, from_mat.layers)
+        self.assertEqual(from_h5.labels, from_mat.labels)
+        for layer in from_mat.layers:
+            with self.subTest(layer=layer):
+                assert_array_equal(from_h5.get(layer), self.stacked[layer])
+                assert_array_equal(from_h5.get(layer), from_mat.get(layer))
+
+    def test_selected_layers_only(self):
+        convert_features_to_hdf5(self.matdir, self.h5dir, layers=['fc8'])
+        self.assertEqual(Features(self.h5dir).layers, ['fc8'])
+
+    def test_unknown_layer_raises(self):
+        with self.assertRaises(KeyError):
+            convert_features_to_hdf5(self.matdir, self.h5dir, layers=['nope'])
+
+    def test_existing_file_is_skipped_unless_overwrite(self):
+        convert_features_to_hdf5(self.matdir, self.h5dir)
+        marker = os.path.join(self.h5dir, 'fc8.h5')
+        mtime = os.path.getmtime(marker)
+        os.utime(marker, (mtime - 100, mtime - 100))
+
+        convert_features_to_hdf5(self.matdir, self.h5dir)
+        self.assertEqual(os.path.getmtime(marker), mtime - 100)
+
+        convert_features_to_hdf5(self.matdir, self.h5dir, overwrite=True)
+        self.assertGreater(os.path.getmtime(marker), mtime - 100)
+
+    def test_layer_published_by_another_process_meanwhile_is_skipped(self):
+        # Two converters on the same destination: the other one publishes
+        # conv5 while this one is still reading it. This one must keep the
+        # other's file, discard its own, and go on to the next layer.
+        other = self.stacked['conv5'][::-1].copy()
+        out_path = os.path.join(self.h5dir, 'conv5.h5')
+        original_read = MatFeatureStore.read
+
+        def read_then_race(store, layer, labels=None, feature_slice=None):
+            if layer == 'conv5' and not os.path.exists(out_path):
+                save_features(out_path, other, self.labels)
+            return original_read(store, layer, labels, feature_slice)
+
+        with mock.patch.object(MatFeatureStore, 'read', read_then_race):
+            convert_features_to_hdf5(self.matdir, self.h5dir)
+
+        converted = Features(self.h5dir)
+        assert_array_equal(converted.get('conv5'), other)
+        assert_array_equal(converted.get('fc8'), self.stacked['fc8'])
+        self.assertEqual(_scratch_files(self.h5dir), [])
+
+    def test_small_batches_produce_the_same_file(self):
+        # Batching is an implementation detail; it must not affect the result.
+        convert_features_to_hdf5(self.matdir, self.h5dir, batch_size=3)
+        assert_array_equal(Features(self.h5dir).get('conv5'), self.stacked['conv5'])
+
+
+class TestSchemaValidation(unittest.TestCase):
+    """Everything a reader depends on is checked when the file is opened."""
+
+    def setUp(self):
+        self.tmpdir = tempfile.TemporaryDirectory()
+        self.path = os.path.join(self.tmpdir.name, 'conv5.h5')
+        save_features(self.path, np.zeros((3, 4)), ['a', 'b', 'c'])
+
+    def tearDown(self):
+        self.tmpdir.cleanup()
+
+    def _mutate(self, fn):
+        with h5py.File(self.path, 'a') as f:
+            fn(f)
+
+    def test_missing_version_is_rejected(self):
+        # Previously read as version 0 and accepted, contradicting the docs.
+        self._mutate(lambda f: f.attrs.__delitem__(FORMAT_VERSION_ATTR))
+        with self.assertRaises(RuntimeError):
+            HDF5FeatureStore(self.tmpdir.name)
+
+    def test_malformed_version_is_rejected(self):
+        # Anything that is not an integer scalar, including a float that would
+        # round down to a supported version under int().
+        for bad in ('v1', 1.5, 1.0, True, np.True_):
+            with self.subTest(version=bad):
+                self._mutate(
+                    lambda f, b=bad: f.attrs.__setitem__(FORMAT_VERSION_ATTR, b)
+                )
+                with self.assertRaises(RuntimeError):
+                    HDF5FeatureStore(self.tmpdir.name)
+
+    def test_version_below_one_is_rejected(self):
+        for bad in (0, -1):
+            with self.subTest(version=bad):
+                self._mutate(lambda f, b=bad: f.attrs.__setitem__(FORMAT_VERSION_ATTR, b))
+                with self.assertRaises(RuntimeError):
+                    HDF5FeatureStore(self.tmpdir.name)
+
+    def test_future_version_is_rejected(self):
+        self._mutate(
+            lambda f: f.attrs.__setitem__(
+                FORMAT_VERSION_ATTR, SUPPORTED_FORMAT_VERSION + 1
+            )
+        )
+        with self.assertRaises(RuntimeError) as ctx:
+            HDF5FeatureStore(self.tmpdir.name)
+        self.assertIn('upgrade bdpy', str(ctx.exception))
+
+    def test_row_count_mismatch_is_rejected(self):
+        def shrink(f):
+            del f['labels']
+            f.create_dataset('labels', data=['a', 'b'],
+                             dtype=h5py.string_dtype(encoding='utf-8'))
+        self._mutate(shrink)
+        with self.assertRaises(RuntimeError) as ctx:
+            HDF5FeatureStore(self.tmpdir.name)
+        self.assertIn('rows', str(ctx.exception))
+
+    def test_wrong_ndim_is_rejected(self):
+        def flatten_features(f):
+            del f['features']
+            f.create_dataset('features', data=np.zeros(3))
+        self._mutate(flatten_features)
+        with self.assertRaises(RuntimeError):
+            HDF5FeatureStore(self.tmpdir.name)
+
+        save_features(self.path, np.zeros((2, 4)), ['a', 'b'], overwrite=True)
+
+        def widen_labels(f):
+            del f['labels']
+            f.create_dataset('labels', data=[['a'], ['b']],
+                             dtype=h5py.string_dtype(encoding='utf-8'))
+        self._mutate(widen_labels)
+        with self.assertRaises(RuntimeError):
+            HDF5FeatureStore(self.tmpdir.name)
+
+    def test_duplicate_labels_are_rejected_on_read(self):
+        # Labels map to row indices, so a repeat would resolve both rows to the
+        # last one and silently drop the first. The .mat layout cannot express
+        # duplicates at all, so neither should this one.
+        def duplicate(f):
+            del f['labels']
+            f.create_dataset('labels', data=['a', 'a', 'c'],
+                             dtype=h5py.string_dtype(encoding='utf-8'))
+        self._mutate(duplicate)
+        with self.assertRaises(RuntimeError) as ctx:
+            HDF5FeatureStore(self.tmpdir.name)
+        self.assertIn('duplicate', str(ctx.exception))
+
+
+class TestDuplicateLabelsOnWrite(unittest.TestCase):
+    def setUp(self):
+        self.tmpdir = tempfile.TemporaryDirectory()
+        self.path = os.path.join(self.tmpdir.name, 'conv5.h5')
+
+    def tearDown(self):
+        self.tmpdir.cleanup()
+
+    def test_save_features_rejects_duplicates(self):
+        with self.assertRaises(ValueError):
+            save_features(self.path, np.zeros((2, 4)), ['a', 'a'])
+        self.assertFalse(os.path.exists(self.path))
+
+    def test_writer_rejects_duplicates_within_one_call(self):
+        with FeatureWriter(self.path, (4,), np.float32) as writer:
+            with self.assertRaises(ValueError) as ctx:
+                writer.extend(np.zeros((2, 4)), ['a', 'a'])
+            self.assertIn('this batch', str(ctx.exception))
+            writer.abort()
+
+    def test_writer_rejects_duplicates_across_calls(self):
+        with FeatureWriter(self.path, (4,), np.float32) as writer:
+            writer.append(np.zeros(4), 'a')
+            with self.assertRaises(ValueError) as ctx:
+                writer.append(np.ones(4), 'a')
+            self.assertIn('already written', str(ctx.exception))
+            writer.abort()
+
+    def test_rejected_batch_leaves_earlier_rows_intact(self):
+        # The check runs before writing, so a refused batch is a no-op.
+        with FeatureWriter(self.path, (4,), np.float32) as writer:
+            writer.extend(np.ones((2, 4)), ['a', 'b'])
+            with self.assertRaises(ValueError):
+                writer.extend(np.zeros((2, 4)), ['c', 'a'])
+            self.assertEqual(writer.n_samples, 2)
+            writer.extend(np.full((1, 4), 7.0), ['c'])
+        store = HDF5FeatureStore(self.tmpdir.name)
+        self.assertEqual(store.labels, ['a', 'b', 'c'])
+
+
+class TestZeroSampleLayer(unittest.TestCase):
+    """Both write paths agree on an empty layer."""
+
+    def setUp(self):
+        self.tmpdir = tempfile.TemporaryDirectory()
+
+    def tearDown(self):
+        self.tmpdir.cleanup()
+
+    def test_save_features_accepts_zero_samples(self):
+        path = os.path.join(self.tmpdir.name, 'conv5.h5')
+        save_features(path, np.empty((0, 8), dtype=np.float32), [])
+        store = HDF5FeatureStore(self.tmpdir.name)
+        self.assertEqual(store.labels, [])
+        self.assertEqual(store.read('conv5').shape, (0, 8))
+
+    def test_writer_accepts_zero_samples(self):
+        path = os.path.join(self.tmpdir.name, 'fc8.h5')
+        with FeatureWriter(path, (8,), np.float32):
+            pass
+        self.assertEqual(HDF5FeatureStore(self.tmpdir.name).read('fc8').shape, (0, 8))
+
+
+class TestAtomicWriteAndOverwrite(unittest.TestCase):
+    """A failed write must leave nothing behind, and must never clobber.
+
+    The converter treats an existing <layer>.h5 as finished and skips it, so a
+    half-written file published at the target path would make a truncated layer
+    permanent across re-runs.
+    """
+
+    def setUp(self):
+        warnings.simplefilter('ignore', FutureWarning)
+        self.tmpdir = tempfile.TemporaryDirectory()
+        self.matdir = os.path.join(self.tmpdir.name, 'mat')
+        self.h5dir = os.path.join(self.tmpdir.name, 'h5')
+        os.makedirs(self.matdir)
+        self.labels = ['img%04d' % i for i in range(10)]
+        self.stacked = prepare_mat_features(
+            self.matdir, ['conv5'], self.labels, [(1, 16, 3, 3)]
+        )
+        self.path = os.path.join(self.tmpdir.name, 'conv5.h5')
+        self.data = np.random.rand(10, 8).astype(np.float32)
+
+    def tearDown(self):
+        self.tmpdir.cleanup()
+
+    def _leftovers(self, directory):
+        return sorted(os.listdir(directory))
+
+    # --- overwrite guard -------------------------------------------------
+
+    def test_save_features_refuses_to_clobber(self):
+        save_features(self.path, self.data, self.labels)
+        with self.assertRaises(FileExistsError):
+            save_features(self.path, self.data * 2, self.labels)
+        # The original must be untouched.
+        assert_array_equal(
+            HDF5FeatureStore(self.tmpdir.name).read('conv5'), self.data
+        )
+
+    def test_save_features_overwrite_replaces(self):
+        save_features(self.path, self.data, self.labels)
+        save_features(self.path, self.data * 2, self.labels, overwrite=True)
+        assert_array_equal(
+            HDF5FeatureStore(self.tmpdir.name).read('conv5'), self.data * 2
+        )
+
+    def test_writer_refuses_to_clobber(self):
+        save_features(self.path, self.data, self.labels)
+        with self.assertRaises(FileExistsError):
+            FeatureWriter(self.path, (8,), np.float32)
+        assert_array_equal(
+            HDF5FeatureStore(self.tmpdir.name).read('conv5'), self.data
+        )
+
+    def test_writer_overwrite_replaces(self):
+        save_features(self.path, self.data, self.labels)
+        with FeatureWriter(self.path, (8,), np.float32, overwrite=True) as w:
+            w.extend(self.data * 3, self.labels)
+        assert_array_equal(
+            HDF5FeatureStore(self.tmpdir.name).read('conv5'), self.data * 3
+        )
+
+    # --- atomicity -------------------------------------------------------
+
+    def test_target_is_absent_until_close(self):
+        writer = FeatureWriter(self.path, (8,), np.float32)
+        writer.append(self.data[0], self.labels[0])
+        self.assertFalse(os.path.exists(self.path))
+        writer.close()
+        self.assertTrue(os.path.exists(self.path))
+
+    def test_abort_leaves_nothing(self):
+        writer = FeatureWriter(self.path, (8,), np.float32)
+        writer.append(self.data[0], self.labels[0])
+        writer.abort()
+        self.assertFalse(os.path.exists(self.path))
+        self.assertEqual(self._leftovers(self.tmpdir.name), ['mat'])
+
+    def test_context_manager_aborts_on_exception(self):
+        with self.assertRaises(ZeroDivisionError):
+            with FeatureWriter(self.path, (8,), np.float32) as writer:
+                writer.append(self.data[0], self.labels[0])
+                raise ZeroDivisionError
+        self.assertFalse(os.path.exists(self.path))
+        self.assertEqual(self._leftovers(self.tmpdir.name), ['mat'])
+
+    def test_close_and_abort_are_idempotent(self):
+        writer = FeatureWriter(self.path, (8,), np.float32)
+        writer.extend(self.data, self.labels)
+        writer.close()
+        writer.close()
+        writer.abort()  # must not delete the published file
+        self.assertTrue(os.path.exists(self.path))
+        assert_array_equal(
+            HDF5FeatureStore(self.tmpdir.name).read('conv5'), self.data
+        )
+
+    def test_failed_save_features_leaves_nothing(self):
+        # Mismatched labels raise after the target has been prepared.
+        with self.assertRaises(ValueError):
+            save_features(self.path, self.data, self.labels[:-1])
+        self.assertFalse(os.path.exists(self.path))
+        self.assertEqual(self._leftovers(self.tmpdir.name), ['mat'])
+
+    def test_published_file_is_group_readable(self):
+        # Staging must not tighten permissions: these files live on shared lab
+        # storage, and 0600 would lock collaborators out.
+        save_features(self.path, self.data, self.labels)
+        mode = stat.S_IMODE(os.stat(self.path).st_mode)
+        expected = 0o666 & ~_current_umask()
+        self.assertEqual(mode, expected)
+
+    def test_staging_file_is_not_mistaken_for_a_layer(self):
+        # A directory being written into is also a directory someone may read.
+        # The scratch file must be invisible to the *.h5 scan, or a concurrent
+        # reader would either invent a bogus layer or fail validation.
+        writer = FeatureWriter(os.path.join(self.h5dir, 'conv5.h5'), (8,), np.float32)
+        try:
+            writer.append(self.data[0], self.labels[0])
+            self.assertEqual(glob.glob(os.path.join(self.h5dir, '*.h5')), [])
+            with self.assertRaises(RuntimeError):
+                HDF5FeatureStore(self.h5dir)  # "no .h5 feature file found"
+        finally:
+            writer.abort()
+
+    def test_failed_overwrite_preserves_the_previous_file(self):
+        # The old file stays readable for the whole write and is swapped only at
+        # the end, so a failure mid-overwrite loses neither old nor new.
+        save_features(self.path, self.data, self.labels)
+        with self.assertRaises(ZeroDivisionError):
+            with FeatureWriter(self.path, (8,), np.float32, overwrite=True) as w:
+                w.extend(self.data * 9, self.labels)
+                raise ZeroDivisionError
+        assert_array_equal(
+            HDF5FeatureStore(self.tmpdir.name).read('conv5'), self.data
+        )
+
+    def test_dropped_writer_does_not_leak_a_scratch_file(self):
+        writer = FeatureWriter(self.path, (8,), np.float32)
+        writer.append(self.data[0], self.labels[0])
+        del writer
+        gc.collect()
+        self.assertEqual(self._leftovers(self.tmpdir.name), ['mat'])
+
+    # --- another writer publishes first ---------------------------------
+
+    def test_close_does_not_clobber_a_file_published_meanwhile(self):
+        # The overwrite check in __init__ is not enough on its own: another
+        # writer can publish the same path while this one is still staging.
+        writer = FeatureWriter(self.path, (8,), np.float32)
+        writer.extend(self.data, self.labels)
+        save_features(self.path, self.data * 2, self.labels)
+        with self.assertRaises(FileExistsError):
+            writer.close()
+        assert_array_equal(
+            HDF5FeatureStore(self.tmpdir.name).read('conv5'), self.data * 2
+        )
+        self.assertEqual(_scratch_files(self.tmpdir.name), [])
+        # The failed publish leaves the writer aborted, not half-closed.
+        with self.assertRaises(RuntimeError):
+            writer.close()
+
+    def test_publishes_without_hard_links(self):
+        # A filesystem that cannot hard-link falls back to check-then-rename.
+        with mock.patch('os.link', side_effect=OSError(errno.EPERM, 'no links')):
+            save_features(self.path, self.data, self.labels)
+        assert_array_equal(
+            HDF5FeatureStore(self.tmpdir.name).read('conv5'), self.data
+        )
+        self.assertEqual(_scratch_files(self.tmpdir.name), [])
+
+    def test_without_hard_links_a_file_published_meanwhile_is_kept(self):
+        writer = FeatureWriter(self.path, (8,), np.float32)
+        writer.extend(self.data, self.labels)
+        save_features(self.path, self.data * 2, self.labels)
+        with mock.patch('os.link', side_effect=OSError(errno.EPERM, 'no links')):
+            with self.assertRaises(FileExistsError):
+                writer.close()
+        assert_array_equal(
+            HDF5FeatureStore(self.tmpdir.name).read('conv5'), self.data * 2
+        )
+        self.assertEqual(_scratch_files(self.tmpdir.name), [])
+
+    def test_close_after_abort_is_refused(self):
+        # Silently doing nothing would let a caller believe it published.
+        writer = FeatureWriter(self.path, (8,), np.float32)
+        writer.abort()
+        with self.assertRaises(RuntimeError):
+            writer.close()
+        self.assertFalse(os.path.exists(self.path))
+
+    def test_zero_samples_round_trips(self):
+        # An empty layer is representable; refusing it would add a failure mode
+        # on the success path for something the caller can check themselves.
+        with FeatureWriter(self.path, (8,), np.float32):
+            pass
+        store = HDF5FeatureStore(self.tmpdir.name)
+        self.assertEqual(store.labels, [])
+        self.assertEqual(store.read('conv5').shape, (0, 8))
+
+    # --- the converter ---------------------------------------------------
+
+    def _fail_on_second_batch(self):
+        """Patch MatFeatureStore.read to blow up partway through a layer."""
+        original = MatFeatureStore.read
+        state = {'calls': 0}
+
+        def flaky(store, layer, labels=None, feature_slice=None):
+            state['calls'] += 1
+            if state['calls'] > 1:
+                raise OSError('simulated read failure')
+            return original(store, layer, labels, feature_slice)
+
+        return original, flaky
+
+    def test_non_positive_batch_size_is_refused(self):
+        # range(0, n, -1) is empty, so this used to write nothing and then
+        # publish the empty file as a finished layer.
+        for bad in (-1, 0):
+            with self.subTest(batch_size=bad):
+                with self.assertRaises(ValueError):
+                    convert_features_to_hdf5(
+                        self.matdir, self.h5dir, batch_size=bad
+                    )
+                self.assertFalse(os.path.exists(self.h5dir))
+
+    def test_failed_conversion_leaves_no_file(self):
+        original, flaky = self._fail_on_second_batch()
+        MatFeatureStore.read = flaky
+        try:
+            with self.assertRaises(OSError):
+                convert_features_to_hdf5(
+                    self.matdir, self.h5dir, batch_size=3
+                )
+        finally:
+            MatFeatureStore.read = original
+
+        out_path = os.path.join(self.h5dir, 'conv5.h5')
+        self.assertFalse(os.path.exists(out_path))
+        # No scratch file left behind either.
+        self.assertEqual(self._leftovers(self.h5dir), [])
+
+    def test_rerun_after_failure_succeeds(self):
+        # The regression this guards: a published partial file would be taken
+        # for a finished one by the skip check and never repaired.
+        original, flaky = self._fail_on_second_batch()
+        MatFeatureStore.read = flaky
+        try:
+            with self.assertRaises(OSError):
+                convert_features_to_hdf5(
+                    self.matdir, self.h5dir, batch_size=3
+                )
+        finally:
+            MatFeatureStore.read = original
+
+        convert_features_to_hdf5(self.matdir, self.h5dir, batch_size=3)
+        assert_array_equal(
+            Features(self.h5dir).get('conv5'), self.stacked['conv5']
+        )
+
+
+class TestFormatValidation(unittest.TestCase):
+    def setUp(self):
+        self.tmpdir = tempfile.TemporaryDirectory()
+        self.path = os.path.join(self.tmpdir.name, 'conv5.h5')
+
+    def tearDown(self):
+        self.tmpdir.cleanup()
+
+    def test_foreign_hdf5_is_rejected(self):
+        # A file that merely happens to have /features is not bdpy storage.
+        with h5py.File(self.path, 'w') as f:
+            f.create_dataset('features', data=np.zeros((2, 3)))
+            f.create_dataset('labels', data=['a', 'b'],
+                             dtype=h5py.string_dtype(encoding='utf-8'))
+        with self.assertRaises(RuntimeError):
+            HDF5FeatureStore(self.tmpdir.name)
+
+    def test_future_version_is_rejected_with_a_clear_message(self):
+        save_features(self.path, np.zeros((2, 3)), ['a', 'b'])
+        with h5py.File(self.path, 'a') as f:
+            f.attrs[FORMAT_VERSION_ATTR] = SUPPORTED_FORMAT_VERSION + 1
+        with self.assertRaises(RuntimeError) as ctx:
+            HDF5FeatureStore(self.tmpdir.name)
+        self.assertIn('upgrade bdpy', str(ctx.exception))
+
+    def test_missing_dataset_is_rejected(self):
+        save_features(self.path, np.zeros((2, 3)), ['a', 'b'])
+        with h5py.File(self.path, 'a') as f:
+            del f['labels']
+        with self.assertRaises(RuntimeError):
+            HDF5FeatureStore(self.tmpdir.name)
+
+    def test_empty_directory_is_rejected(self):
+        with self.assertRaises(RuntimeError):
+            HDF5FeatureStore(self.tmpdir.name)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/tests/dataform/test_feature_store.py b/tests/dataform/test_feature_store.py
new file mode 100644
index 00000000..41bba404
--- /dev/null
+++ b/tests/dataform/test_feature_store.py
@@ -0,0 +1,689 @@
+"""Backend equivalence for feature storage.
+
+The legacy .mat layout and chunked HDF5 storage must be indistinguishable
+through the Features API: same labels, same order, same arrays, for every
+combination of label selection and feature slicing. These tests run each query
+against both backends and compare.
+"""
+
+import os
+import tempfile
+import unittest
+import warnings
+from unittest import mock
+
+import h5py
+import numpy as np
+from numpy.testing import assert_array_equal
+
+from bdpy.dataform import Features, convert_features_to_hdf5, save_features
+from bdpy.dataform import _feature_store
+from bdpy.dataform._feature_store import (
+    FeatureStore,
+    HDF5FeatureStore,
+    MatFeatureStore,
+    _axis_index,
+    _check_size,
+    _default_iter_size,
+    _iter_array_slabs,
+    _normalize_axis,
+    _normalize_feature_slice,
+    _partition,
+    detect_format,
+)
+
+from .test_features import prepare_mat_features
+
+LAYERS = ['conv5', 'fc8']
+SHAPES = [(1, 24, 5, 5), (1, 60)]
+LABELS = ['img%04d' % i for i in range(11)]
+
+# Every way a caller might ask for features. Each is run against both backends.
+QUERIES = [
+    ('all', 'conv5', None, None),
+    ('single label', 'conv5', 'img0003', None),
+    ('label list', 'conv5', ['img0005', 'img0001', 'img0009'], None),
+    ('unsorted labels', 'conv5', ['img0009', 'img0000'], None),
+    ('repeated labels', 'conv5', ['img0009', 'img0001', 'img0009', 'img0000'], None),
+    ('feature slice', 'conv5', None, np.s_[8:16]),
+    ('slice with step', 'conv5', None, np.s_[::2]),
+    ('multi-axis slice', 'conv5', None, np.s_[8:16, 1:4, :]),
+    ('integer index', 'conv5', None, np.s_[5]),
+    ('ellipsis', 'conv5', None, np.s_[..., 1:3]),
+    ('labels and slice', 'conv5', ['img0009', 'img0001', 'img0009'], np.s_[8:16]),
+    ('2d layer', 'fc8', None, np.s_[10:40]),
+    ('2d layer with labels', 'fc8', ['img0002', 'img0000'], np.s_[10:40]),
+]
+
+
+# feature_slice specs outside basic forward indexing. Validation happens once,
+# in FeatureStore.read, so the table is checked against the validator itself;
+# TestFeatureSliceValidation checks that every backend goes through it.
+REJECTED_FEATURE_SLICES = [
+    ('negative step', np.s_[::-1]),
+    ('reversed with bounds', np.s_[16:4:-1]),
+    ('zero step', slice(None, None, 0)),
+    ('list', np.s_[[3, 1, 7]]),
+    ('ndarray', np.array([1, 2])),
+    ('range', range(3)),
+    ('bool', True),
+    ('numpy bool', np.True_),
+    ('newaxis in tuple', (np.newaxis, slice(None))),
+    ('two ellipses', (Ellipsis, Ellipsis)),
+    ('float start', slice(1.5, 3)),
+    ('float stop', slice(1, 3.5)),
+    ('float step', slice(None, None, 2.0)),
+    ('bool start', slice(True, 3)),
+]
+
+
+class TestIterationHelpers(unittest.TestCase):
+    """The index arithmetic behind iter_chunks, on plain in-memory values."""
+
+    def test_normalize_axis(self):
+        shape = (11, 24, 5, 5)
+        for axis, expected in ((0, 0), (3, 3), (-1, 3), (-4, 0)):
+            with self.subTest(axis=axis):
+                self.assertEqual(_normalize_axis(axis, shape), expected)
+
+    def test_normalize_axis_rejects_out_of_range(self):
+        shape = (11, 24, 5, 5)
+        for axis, message in (
+            (4, 'axis 4 is out of range for selected shape (11, 24, 5, 5)'),
+            (-5, 'axis -5 is out of range for selected shape (11, 24, 5, 5)'),
+        ):
+            with self.subTest(axis=axis):
+                with self.assertRaises(ValueError) as ctx:
+                    _normalize_axis(axis, shape)
+                self.assertEqual(str(ctx.exception), message)
+
+    def test_check_size(self):
+        _check_size(1)
+        for bad in (0, -1):
+            with self.subTest(size=bad):
+                with self.assertRaises(ValueError):
+                    _check_size(bad)
+
+    def test_partition(self):
+        cases = [
+            ('exact', 6, 3, [slice(0, 3), slice(3, 6)]),
+            ('partial last', 7, 3, [slice(0, 3), slice(3, 6), slice(6, 7)]),
+            ('size above length', 4, 10, [slice(0, 4)]),
+            ('size one', 3, 1, [slice(0, 1), slice(1, 2), slice(2, 3)]),
+            ('empty', 0, 3, []),
+        ]
+        for name, length, size, expected in cases:
+            with self.subTest(case=name):
+                self.assertEqual(list(_partition(length, size)), expected)
+
+    def test_axis_index(self):
+        sl = slice(2, 4)
+        self.assertEqual(_axis_index(0, sl), (sl,))
+        self.assertEqual(_axis_index(2, sl), (slice(None), slice(None), sl))
+
+    def test_default_size_follows_the_chunk_extent(self):
+        shape = (3, 10, 4)
+        self.assertEqual(_default_iter_size(shape, 1, 8, extent=4), 4)
+        # Never longer than the axis, never below one.
+        self.assertEqual(_default_iter_size(shape, 1, 8, extent=50), 10)
+        self.assertEqual(_default_iter_size(shape, 1, 8, extent=0), 1)
+
+    def test_default_size_without_extent_follows_the_byte_budget(self):
+        # One element along axis 1 spans 3 * 4 items of 8 bytes = 96 bytes.
+        shape = (3, 10, 4)
+        with mock.patch.object(_feature_store, 'DEFAULT_TARGET_CHUNK_BYTES', 300):
+            self.assertEqual(_default_iter_size(shape, 1, 8, extent=None), 3)
+        with mock.patch.object(_feature_store, 'DEFAULT_TARGET_CHUNK_BYTES', 10):
+            self.assertEqual(_default_iter_size(shape, 1, 8, extent=None), 1)
+        with mock.patch.object(_feature_store, 'DEFAULT_TARGET_CHUNK_BYTES', 10 ** 6):
+            self.assertEqual(_default_iter_size(shape, 1, 8, extent=None), 10)
+
+    def test_iter_array_slabs(self):
+        array = np.arange(2 * 7 * 3).reshape(2, 7, 3)
+        slabs = list(_iter_array_slabs(array, 1, _partition(7, 3)))
+        self.assertEqual([sl for sl, _ in slabs], list(_partition(7, 3)))
+        for sl, block in slabs:
+            assert_array_equal(block, array[:, sl])
+        assert_array_equal(np.concatenate([b for _, b in slabs], axis=1), array)
+
+    def test_rejected_feature_slices(self):
+        for name, spec in REJECTED_FEATURE_SLICES:
+            with self.subTest(index=name):
+                with self.assertRaises(ValueError):
+                    _normalize_feature_slice(spec)
+
+
+class _BackendPair(unittest.TestCase):
+    """Builds the same features as a .mat tree and as chunked HDF5."""
+
+    def setUp(self):
+        warnings.simplefilter('ignore', FutureWarning)
+        self.tmpdir = tempfile.TemporaryDirectory()
+        self.matdir = os.path.join(self.tmpdir.name, 'mat')
+        self.h5dir = os.path.join(self.tmpdir.name, 'h5')
+        os.makedirs(self.matdir)
+        self.stacked = prepare_mat_features(self.matdir, LAYERS, LABELS, SHAPES)
+        convert_features_to_hdf5(self.matdir, self.h5dir)
+        self.from_mat = Features(self.matdir)
+        self.from_h5 = Features(self.h5dir)
+
+    def tearDown(self):
+        self.tmpdir.cleanup()
+
+
+class TestBackendEquivalence(_BackendPair):
+    def test_metadata_matches(self):
+        self.assertEqual(self.from_h5.layers, self.from_mat.layers)
+        self.assertEqual(self.from_h5.labels, self.from_mat.labels)
+        assert_array_equal(self.from_h5.index, self.from_mat.index)
+        for layer in LAYERS:
+            self.assertEqual(self.from_h5.shape(layer), self.from_mat.shape(layer))
+
+    def test_queries_match_across_backends(self):
+        for name, layer, label, feature_slice in QUERIES:
+            with self.subTest(query=name):
+                from_mat = self.from_mat.get(
+                    layer, label=label, feature_slice=feature_slice
+                )
+                from_h5 = self.from_h5.get(
+                    layer, label=label, feature_slice=feature_slice
+                )
+                self.assertEqual(from_h5.shape, from_mat.shape)
+                assert_array_equal(from_h5, from_mat)
+
+    def test_queries_match_ground_truth(self):
+        # Both backends agreeing is not enough if both are wrong the same way.
+        for name, layer, label, feature_slice in QUERIES:
+            with self.subTest(query=name):
+                if label is None:
+                    rows = slice(None)
+                elif isinstance(label, str):
+                    rows = [LABELS.index(label)]
+                else:
+                    rows = [LABELS.index(s) for s in label]
+                expected = self.stacked[layer][rows]
+                if feature_slice is not None:
+                    indexers = (
+                        feature_slice
+                        if isinstance(feature_slice, tuple)
+                        else (feature_slice,)
+                    )
+                    expected = expected[(slice(None), *indexers)]
+                assert_array_equal(
+                    self.from_h5.get(layer, label=label, feature_slice=feature_slice),
+                    expected,
+                )
+
+    def test_statistic_matches(self):
+        for statistic in ('mean', 'std', 'std, ddof=0'):
+            for layer in LAYERS:
+                with self.subTest(statistic=statistic, layer=layer):
+                    assert_array_equal(
+                        self.from_h5.statistic(statistic, layer=layer),
+                        self.from_mat.statistic(statistic, layer=layer),
+                    )
+
+
+class TestMatIterChunksReadsOnce(_BackendPair):
+    """The legacy backend must read a feature axis once, not once per slab.
+
+    ``MatFeatureStore.read`` loads every selected stimulus file, so reading a
+    slab at a time made a k-slab iteration cost k full reads of the layer
+    (measured: 1088 s against 5 s for one read plus in-memory slicing).
+    """
+
+    def setUp(self):
+        super().setUp()
+        self.store = MatFeatureStore(self.matdir)
+        self.reads = []
+        original = self.store.read
+
+        def spy(*args, **kwargs):
+            self.reads.append((args, kwargs))
+            return original(*args, **kwargs)
+
+        self.store.read = spy
+
+    def test_feature_axis_reads_the_selection_once(self):
+        slab_counts = {}
+        for size in (1, 3, 4, 24, 100, None):
+            with self.subTest(size=size):
+                self.reads.clear()
+                blocks = list(self.store.iter_chunks('conv5', axis=1, size=size))
+                self.assertEqual(len(self.reads), 1)
+                slab_counts[size] = len(blocks)
+        # Not vacuous: the slab count really does vary, and the worst case is
+        # the 24 reads this test exists to prevent.
+        self.assertEqual(slab_counts[1], 24)
+        self.assertEqual(slab_counts[100], 1)
+        self.assertGreater(len(set(slab_counts.values())), 1)
+
+    def test_sample_axis_still_reads_one_slab_at_a_time(self):
+        # Preloading here would destroy the one guarantee the legacy backend
+        # can make: a layer larger than memory streams along the sample axis.
+        blocks = list(self.store.iter_chunks('conv5', axis=0, size=4))
+        self.assertEqual(len(blocks), 3)
+        self.assertEqual(len(self.reads), 3)
+        for (args, _), (sl, _) in zip(self.reads, blocks):
+            self.assertEqual(list(args[1]), LABELS[sl])
+
+    def test_nothing_is_cached_between_calls(self):
+        for call in (1, 2):
+            with self.subTest(call=call):
+                self.reads.clear()
+                list(self.store.iter_chunks('conv5', axis=1, size=4))
+                self.assertEqual(len(self.reads), 1)
+
+    def test_writing_to_one_block_does_not_disturb_another(self):
+        # Slabs are views into one array now, but they partition the axis, so
+        # they cannot alias each other.
+        blocks = list(self.from_mat.iter_chunks('conv5', axis=1, size=4))
+        blocks[0][1][:] = 0
+        assert_array_equal(blocks[1][1], self.stacked['conv5'][:, 4:8])
+
+    def test_validation_precedes_selection_read(self):
+        cases = (
+            ({'axis': 9}, 'axis 9 is out of range for selected shape (11, 24, 5, 5)'),
+            ({'axis': -9}, 'axis -9 is out of range for selected shape (11, 24, 5, 5)'),
+            ({'size': 0}, 'size must be positive, got 0'),
+        )
+        for kwargs, message in cases:
+            with self.subTest(**kwargs):
+                self.reads.clear()
+                with self.assertRaises(ValueError) as ctx:
+                    list(self.store.iter_chunks('conv5', **kwargs))
+                self.assertEqual(str(ctx.exception), message)
+                # Reading the selection up front must not outrun validation.
+                self.assertEqual(len(self.reads), 0)
+                with self.assertRaises(ValueError) as from_h5:
+                    list(HDF5FeatureStore(self.h5dir).iter_chunks('conv5', **kwargs))
+                self.assertEqual(str(from_h5.exception), message)
+
+
+class TestIterChunks(_BackendPair):
+    def test_sample_axis_reassembles_to_a_full_read(self):
+        # The sample axis is iterated by label, not through a feature slice, so
+        # it is a separate path from the one test_feature_axis_blocks_match_ground_truth
+        # covers.
+        for features in (self.from_mat, self.from_h5):
+            with self.subTest(store=type(features._Features__stores[0]).__name__):
+                blocks = list(features.iter_chunks('conv5', axis=0, size=4))
+                joined = np.concatenate([b for _, b in blocks], axis=0)
+                assert_array_equal(joined, self.stacked['conv5'])
+
+    def test_blocks_can_be_sliced_by_the_caller(self):
+        # iter_chunks deliberately takes no feature_slice; slicing the blocks as
+        # they come out is the supported way to combine the two.
+        wanted = np.zeros_like(self.stacked['conv5'][:, :, 1:2])
+        for sl, block in self.from_h5.iter_chunks('conv5', axis=1, size=6):
+            wanted[:, sl] = block[:, :, 1:2]
+        assert_array_equal(wanted, self.stacked['conv5'][:, :, 1:2])
+
+    def test_default_size_is_chunk_aligned(self):
+        store = HDF5FeatureStore(self.h5dir)
+        extent = store.chunk_extent('conv5', 1)
+        sizes = {
+            sl.stop - sl.start
+            for sl, _ in self.from_h5.iter_chunks('conv5', axis=1)
+        }
+        # Every block but the last is a whole chunk along the axis.
+        self.assertLessEqual(max(sizes), max(1, extent))
+
+    def test_multi_directory_validates_like_a_single_one(self):
+        # The multi-store path does not go through FeatureStore.iter_chunks, so
+        # it used to skip these checks entirely: size=-1 yielded nothing and a
+        # bad axis raised IndexError instead of ValueError.
+        other = os.path.join(self.tmpdir.name, 'other')
+        os.makedirs(other)
+        other_labels = ['other%04d' % i for i in range(3)]
+        for layer, shape in zip(LAYERS, SHAPES):
+            save_features(
+                os.path.join(other, layer + '.h5'),
+                np.random.rand(len(other_labels), *shape[1:]),
+                other_labels,
+            )
+        spread = Features([self.h5dir, other])
+        both = [LABELS[0], other_labels[0]]
+        # Spanning two directories is what takes the fallback path.
+        self.assertEqual(len(spread.get('conv5', label=both)), 2)
+
+        for kwargs in ({'size': -1}, {'size': 0}, {'axis': 9}, {'axis': -9}):
+            with self.subTest(**kwargs):
+                with self.assertRaises(ValueError):
+                    list(spread.iter_chunks('conv5', label=both, **kwargs))
+                with self.assertRaises(ValueError):
+                    list(self.from_h5.iter_chunks('conv5', **kwargs))
+
+    def test_feature_axis_blocks_match_ground_truth(self):
+        # What each backend owns is cutting a feature-axis slab out of what it
+        # reads; the axis and size arithmetic is checked in
+        # TestIterationHelpers. Each size leaves a partial last slab, and the
+        # labels repeat and are out of order.
+        labels = ['img0009', 'img0001', 'img0009']
+        cases = [('conv5', 1, 5), ('conv5', 2, 2), ('conv5', 3, 2), ('fc8', 1, 7)]
+        for features in (self.from_mat, self.from_h5):
+            for layer, axis, size in cases:
+                expected = self.stacked[layer][[LABELS.index(s) for s in labels]]
+                with self.subTest(store=type(features._Features__stores[0]).__name__,
+                                  layer=layer, axis=axis):
+                    blocks = list(features.iter_chunks(
+                        layer, label=labels, axis=axis, size=size))
+                    self.assertEqual(
+                        [sl for sl, _ in blocks],
+                        list(_partition(expected.shape[axis], size)),
+                    )
+                    for sl, block in blocks:
+                        assert_array_equal(block, expected[_axis_index(axis, sl)])
+
+    def test_multi_directory_default_size_follows_the_byte_budget(self):
+        # Spanning two directories takes the in-memory fallback, which used to
+        # yield the whole axis as one slab when size was None. It now splits
+        # the way a store without on-disk chunking does. The budget is shrunk
+        # so that the small fixture spans several slabs.
+        other_dir = tempfile.TemporaryDirectory()
+        self.addCleanup(other_dir.cleanup)
+        other = other_dir.name
+        rng = np.random.default_rng(1)
+        for layer, shape in zip(LAYERS, SHAPES):
+            data = rng.random((1, *shape[1:]))
+            save_features(os.path.join(other, layer + '.h5'), data, ['other0000'])
+            if layer == 'conv5':
+                other_data = data
+        spread = Features([self.matdir, other])
+        both = [LABELS[0], 'other0000']
+        expected = np.concatenate([self.stacked['conv5'][:1], other_data])
+
+        with mock.patch.object(_feature_store, 'DEFAULT_TARGET_CHUNK_BYTES', 1000):
+            from_spread = list(spread.iter_chunks('conv5', label=both, axis=1))
+            # The .mat store has no chunk extent, so it uses the same budget.
+            from_store = list(self.from_mat.iter_chunks('conv5', label=LABELS[:2], axis=1))
+
+        self.assertGreater(len(from_spread), 1)
+        self.assertEqual([sl for sl, _ in from_spread], [sl for sl, _ in from_store])
+        for sl, block in from_spread:
+            assert_array_equal(block, expected[:, sl])
+
+    def test_negative_axis_is_accepted(self):
+        by_negative = list(self.from_h5.iter_chunks('conv5', axis=-3, size=4))
+        by_positive = list(self.from_h5.iter_chunks('conv5', axis=1, size=4))
+        self.assertEqual(len(by_negative), len(by_positive))
+        for (sl_a, a), (sl_b, b) in zip(by_negative, by_positive):
+            self.assertEqual(sl_a, sl_b)
+            assert_array_equal(a, b)
+
+
+class TestPartialReads(_BackendPair):
+    def test_slice_is_pushed_down_to_h5py(self):
+        # The point of the format: h5py must receive the narrowed selection, so
+        # that only the covering chunks are read, rather than a full read that
+        # NumPy then slices.
+        seen = []
+        original = h5py.Dataset.__getitem__
+
+        def spy(dataset, args):
+            seen.append(args)
+            return original(dataset, args)
+
+        h5py.Dataset.__getitem__ = spy
+        try:
+            self.from_h5.get('conv5', feature_slice=np.s_[8:16])
+        finally:
+            h5py.Dataset.__getitem__ = original
+
+        self.assertEqual(len(seen), 1)
+        rows, feature_index = seen[0]
+        self.assertEqual(rows, slice(None))
+        self.assertEqual(feature_index, slice(8, 16, None))
+
+    def test_labels_are_read_once_and_in_order(self):
+        # Duplicated and out-of-order labels must reach h5py as a single
+        # strictly increasing index list, which is all h5py accepts.
+        seen = []
+        original = h5py.Dataset.__getitem__
+
+        def spy(dataset, args):
+            seen.append(args)
+            return original(dataset, args)
+
+        h5py.Dataset.__getitem__ = spy
+        try:
+            out = self.from_h5.get('conv5', label=['img0009', 'img0001', 'img0009'])
+        finally:
+            h5py.Dataset.__getitem__ = original
+
+        self.assertEqual(len(seen), 1)
+        rows = seen[0] if not isinstance(seen[0], tuple) else seen[0][0]
+        self.assertEqual(rows, [1, 9])
+        assert_array_equal(out, self.stacked['conv5'][[9, 1, 9]])
+
+    def test_dataset_is_chunked(self):
+        # This fixture is small enough to fit in one chunk, so only the fact
+        # that chunking is enabled can be asserted here. That both sliceable
+        # axes are actually split is a property of the chunk-shape policy and
+        # is tested against a representative large shape in
+        # tests/dataform/test_feature_chunking.py.
+        with h5py.File(os.path.join(self.h5dir, 'conv5.h5'), 'r') as f:
+            self.assertIsNotNone(f['features'].chunks)
+
+
+class TestFeatureSliceValidation(_BackendPair):
+    """feature_slice is basic forward indexing, identically on both backends.
+
+    The point of narrowing it: whatever is accepted must mean the same thing on
+    both layouts, and whatever is rejected must be rejected by both. Previously
+    a negative step worked on .mat and raised inside h5py, which is exactly the
+    divergence this table guards against.
+    """
+
+    accepted = [
+        ('plain slice', np.s_[4:16]),
+        ('negative stop', np.s_[:-1]),
+        ('negative start', np.s_[-8:]),
+        ('positive step', np.s_[::2]),
+        ('integer', np.s_[5]),
+        ('ellipsis', np.s_[..., 1:3]),
+        ('multi-axis', np.s_[4:16, 1:3]),
+        ('open slice', np.s_[:]),
+    ]
+
+    def test_accepted_agree_across_backends_and_numpy(self):
+        for name, spec in self.accepted:
+            with self.subTest(index=name):
+                from_mat = self.from_mat.get('conv5', feature_slice=spec)
+                from_h5 = self.from_h5.get('conv5', feature_slice=spec)
+                indexers = spec if isinstance(spec, tuple) else (spec,)
+                expected = self.stacked['conv5'][(slice(None), *indexers)]
+                assert_array_equal(from_mat, expected)
+                assert_array_equal(from_h5, expected)
+
+    def test_backends_do_not_override_read(self):
+        # FeatureStore.read is where feature_slice is validated; a backend that
+        # overrode it could skip validation, which no spy on _read would see.
+        for store_cls in (MatFeatureStore, HDF5FeatureStore):
+            with self.subTest(store=store_cls.__name__):
+                self.assertIs(store_cls.read, FeatureStore.read)
+
+    def test_rejected_by_both_backends_before_backend_read(self):
+        # The full table of rejected specs is checked in TestIterationHelpers;
+        # here each backend only has to route through the shared validation,
+        # and reject before the backend's _read runs.
+        for store in (MatFeatureStore(self.matdir), HDF5FeatureStore(self.h5dir)):
+            with self.subTest(store=type(store).__name__):
+                with mock.patch.object(store, '_read') as backend_read:
+                    with self.assertRaises(ValueError):
+                        store.read('conv5', feature_slice=np.s_[::-1])
+                backend_read.assert_not_called()
+
+    def test_bool_is_not_read_as_an_integer(self):
+        # isinstance(True, int) is True in Python, so a naive integer check
+        # would silently read True as index 1.
+        with self.assertRaises(ValueError):
+            self.from_h5.get('conv5', feature_slice=True)
+
+    def test_rejection_message_points_at_the_alternative(self):
+        with self.assertRaises(ValueError) as ctx:
+            self.from_h5.get('conv5', feature_slice=np.s_[[1, 2]])
+        self.assertIn('NumPy', str(ctx.exception))
+
+    def test_iter_chunks_takes_no_feature_slice(self):
+        for features in (self.from_mat, self.from_h5):
+            with self.assertRaises(TypeError):
+                list(features.iter_chunks('conv5', feature_slice=np.s_[0:4]))
+
+
+class TestFeatureSliceAnnotation(unittest.TestCase):
+    def test_type_hints_resolve(self):
+        # FeatureSlice once contained an "ellipsis" forward reference, which is
+        # not a Python name, so any tool resolving annotations blew up.
+        import typing
+
+        from bdpy.dataform._feature_store import FeatureStore
+
+        typing.get_type_hints(FeatureStore.read)
+        typing.get_type_hints(MatFeatureStore.read)
+
+
+class TestCrossLayerLabelConsistency(unittest.TestCase):
+    """Every layer must hold the same labels in the same order.
+
+    A row index built from one layer is used against every layer's /features,
+    so a layer whose labels are ordered differently would make a label lookup
+    silently return another stimulus' row. The legacy .mat backend enforces the
+    same invariant, and the HDF5 backend must match it.
+    """
+
+    def setUp(self):
+        warnings.simplefilter('ignore', FutureWarning)
+        self.tmpdir = tempfile.TemporaryDirectory()
+
+    def tearDown(self):
+        self.tmpdir.cleanup()
+
+    def _write(self, layer, labels, n_features=8):
+        data = np.random.rand(len(labels), n_features)
+        save_features(os.path.join(self.tmpdir.name, layer + '.h5'), data, labels)
+        return data
+
+    def test_consistent_layers_are_accepted(self):
+        self._write('conv5', LABELS)
+        self._write('fc8', LABELS)
+        self.assertEqual(Features(self.tmpdir.name).labels, LABELS)
+
+    def test_different_order_is_rejected(self):
+        self._write('conv5', LABELS)
+        self._write('fc8', list(reversed(LABELS)))
+        with self.assertRaises(RuntimeError) as ctx:
+            HDF5FeatureStore(self.tmpdir.name)
+        message = str(ctx.exception)
+        self.assertIn('different order', message)
+        self.assertIn('fc8', message)
+
+    def test_different_content_is_rejected(self):
+        self._write('conv5', LABELS)
+        self._write('fc8', ['other%04d' % i for i in range(len(LABELS))])
+        with self.assertRaises(RuntimeError) as ctx:
+            HDF5FeatureStore(self.tmpdir.name)
+        self.assertIn('different labels', str(ctx.exception))
+
+    def test_different_length_is_rejected(self):
+        self._write('conv5', LABELS)
+        self._write('fc8', LABELS[:-2])
+        with self.assertRaises(RuntimeError) as ctx:
+            HDF5FeatureStore(self.tmpdir.name)
+        self.assertIn('labels', str(ctx.exception))
+
+    def test_features_construction_also_rejects(self):
+        # The guard must fire through the public entry point too.
+        self._write('conv5', LABELS)
+        self._write('fc8', list(reversed(LABELS)))
+        with self.assertRaises(RuntimeError):
+            Features(self.tmpdir.name)
+
+    def test_legacy_backend_rejects_the_same_way(self):
+        # Parity with the .mat layout, which raises when layers disagree.
+        matdir = os.path.join(self.tmpdir.name, 'mat')
+        os.makedirs(matdir)
+        prepare_mat_features(matdir, ['conv5'], LABELS, [(1, 8)])
+        prepare_mat_features(matdir, ['fc8'], LABELS[:-2], [(1, 8)])
+        with self.assertRaises(RuntimeError):
+            Features(matdir)
+
+
+class TestFormatDetection(unittest.TestCase):
+    def setUp(self):
+        warnings.simplefilter('ignore', FutureWarning)
+        self.tmpdir = tempfile.TemporaryDirectory()
+        self.matdir = os.path.join(self.tmpdir.name, 'mat')
+        self.h5dir = os.path.join(self.tmpdir.name, 'h5')
+        os.makedirs(self.matdir)
+        self.stacked = prepare_mat_features(self.matdir, LAYERS, LABELS, SHAPES)
+        convert_features_to_hdf5(self.matdir, self.h5dir)
+
+    def tearDown(self):
+        self.tmpdir.cleanup()
+
+    def test_detects_each_layout(self):
+        self.assertEqual(detect_format(self.matdir), 'mat')
+        self.assertEqual(detect_format(self.h5dir), 'hdf5')
+
+    def test_unrelated_subdirectory_does_not_change_detection(self):
+        # Detection keys off real feature files, not the bare existence of a
+        # subdirectory, so a stray directory cannot flip a chunked tree to .mat.
+        os.makedirs(os.path.join(self.h5dir, 'notes'))
+        self.assertEqual(detect_format(self.h5dir), 'hdf5')
+        self.assertEqual(Features(self.h5dir).layers, LAYERS)
+
+    def test_both_layouts_present_is_ambiguous(self):
+        # Guessing here means silently reading different data than intended.
+        for layer in LAYERS:
+            os.makedirs(os.path.join(self.h5dir, layer))
+            open(os.path.join(self.h5dir, layer, 'x.mat'), 'w').close()
+        with self.assertRaises(RuntimeError) as ctx:
+            detect_format(self.h5dir)
+        self.assertIn('ambiguous', str(ctx.exception))
+        with self.assertRaises(RuntimeError):
+            Features(self.h5dir)
+        # An explicit format resolves it.
+        self.assertEqual(Features(self.h5dir, format='hdf5').layers, LAYERS)
+
+    def test_empty_directory_is_rejected(self):
+        empty = os.path.join(self.tmpdir.name, 'empty')
+        os.makedirs(empty)
+        with self.assertRaises(RuntimeError):
+            detect_format(empty)
+
+    def test_explicit_format_overrides_detection(self):
+        self.assertIsInstance(
+            Features(self.h5dir, format='hdf5')._Features__stores[0],
+            HDF5FeatureStore,
+        )
+        self.assertIsInstance(
+            Features(self.matdir, format='mat')._Features__stores[0],
+            MatFeatureStore,
+        )
+
+    def test_unknown_format_raises(self):
+        with self.assertRaises(ValueError):
+            Features(self.h5dir, format='parquet')
+
+    def test_mixed_directories(self):
+        # One dpath per layout, read through a single Features.
+        h5_only = os.path.join(self.tmpdir.name, 'h5b')
+        os.makedirs(h5_only)
+        other_labels = ['other%04d' % i for i in range(4)]
+        for layer, shape in zip(LAYERS, SHAPES):
+            data = np.random.rand(len(other_labels), *shape[1:])
+            save_features(
+                os.path.join(h5_only, layer + '.h5'), data, other_labels
+            )
+
+        features = Features([self.matdir, h5_only])
+        self.assertEqual(features.labels, LABELS + other_labels)
+        # A query spanning both directories keeps the caller's order.
+        got = features.get('conv5', label=['other0002', 'img0003', 'other0000'])
+        assert_array_equal(got[1], self.stacked['conv5'][3])
+        self.assertEqual(got.shape[0], 3)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/tests/dataform/test_features.py b/tests/dataform/test_features.py
index 8b8bbd4c..2bd68d58 100644
--- a/tests/dataform/test_features.py
+++ b/tests/dataform/test_features.py
@@ -13,13 +13,13 @@
 from bdpy.dataform.features import Features, save_feature
 
 
-def _prepare_mock_data(
+def prepare_mat_features(
         tmpdir: str,
         mock_layer_names: List[str],
         mock_image_names: List[str],
         mock_shapes: List[Tuple[int, ...]]
     ) -> dict:
-    """Prepare mock data for testing.
+    """Write a legacy per-stimulus feature directory for testing.
 
     Files are written as MATLAB v7.3 (HDF5) so the loader's h5py path (used for
     NumPy 2.0 compatibility) is exercised. The stacked arrays are returned so
@@ -57,7 +57,7 @@ def setUp(self):
         ]
         self.mock_shapes = [(1, 1000), (1, 256, 13, 13)]
         self.feature_dir = tempfile.TemporaryDirectory()
-        stacked = _prepare_mock_data(
+        stacked = prepare_mat_features(
             self.feature_dir.name,
             self.mock_layer_names,
             self.mock_image_names,
@@ -122,6 +122,132 @@ def test_features_get_label(self):
         )
 
 
+class TestFeaturesPartialRead(unittest.TestCase):
+    """feature_slice / iter_chunks on the legacy .mat backend.
+
+    The legacy layout cannot read partially, so these only check that the
+    result is the same as slicing a full read -- which is what makes the API
+    backend-independent.
+    """
+
+    def setUp(self):
+        self.layers = ['fc8', 'conv5']
+        self.labels = [
+            'n01443537_22563',
+            'n01443537_22564',
+            'n01677366_18182',
+            'n04572121_3262',
+        ]
+        self.shapes = [(1, 100), (1, 32, 5, 5)]
+        self.feature_dir = tempfile.TemporaryDirectory()
+        self.stacked = prepare_mat_features(
+            self.feature_dir.name, self.layers, self.labels, self.shapes
+        )
+
+    def tearDown(self):
+        self.feature_dir.cleanup()
+
+    def test_shape_without_reading(self):
+        feat = Features(self.feature_dir.name)
+        self.assertEqual(feat.shape('conv5'), self.stacked['conv5'].shape)
+        self.assertEqual(feat.shape('fc8'), self.stacked['fc8'].shape)
+
+    def test_feature_slice(self):
+        feat = Features(self.feature_dir.name)
+        assert_array_equal(
+            feat.get('conv5', feature_slice=np.s_[8:16]),
+            self.stacked['conv5'][:, 8:16],
+        )
+        assert_array_equal(
+            feat.get('fc8', feature_slice=np.s_[10:50]),
+            self.stacked['fc8'][:, 10:50],
+        )
+
+    def test_feature_slice_with_labels(self):
+        feat = Features(self.feature_dir.name)
+        labels = [self.labels[2], self.labels[0]]
+        assert_array_equal(
+            feat.get('conv5', label=labels, feature_slice=np.s_[8:16]),
+            self.stacked['conv5'][[2, 0]][:, 8:16],
+        )
+
+    def test_iter_chunks_reassembles(self):
+        feat = Features(self.feature_dir.name)
+        blocks = list(feat.iter_chunks('conv5', axis=1, size=7))
+        assert_array_equal(
+            np.concatenate([b for _, b in blocks], axis=1), self.stacked['conv5']
+        )
+
+    def test_unsliced_get_still_uses_the_layer_cache(self):
+        # The cache holds a whole layer; a sliced read must not replace it.
+        feat = Features(self.feature_dir.name)
+        cached = feat.get_features('conv5')
+        feat.get('conv5', feature_slice=np.s_[0:2])
+        self.assertIs(feat.get_features('conv5'), cached)
+
+
+class TestFeaturesFeatureIndex(unittest.TestCase):
+    """Unit-index selection.
+
+    NOTE: the happy path is not covered here. `feature_index` files are struct
+    .mat files, which the current reader (_mat_v73.loadmat_key) cannot load --
+    it handles dense arrays only, so a struct raises TypeError. That is a
+    pre-existing regression from the hdf5storage -> h5py read-path change
+    (issue #106), independent of the storage backends, and is left untouched
+    here rather than silently changed. Only the unambiguous case is asserted.
+    """
+
+    def setUp(self):
+        self.labels = ['img0001', 'img0002', 'img0003']
+        self.feature_dir = tempfile.TemporaryDirectory()
+        self.stacked = prepare_mat_features(
+            self.feature_dir.name, ['fc8'], self.labels, [(1, 20)]
+        )
+
+    def tearDown(self):
+        self.feature_dir.cleanup()
+
+    def _features_with_index(self):
+        """Features with a unit index installed.
+
+        NOTE: the index table is injected directly rather than loaded from a
+        file. Loading is broken on `dev` (struct .mat files are unreadable by
+        _mat_v73.loadmat_key) and fixing that is out of scope here, but the
+        guards below must still be covered.
+        """
+        feat = Features(self.feature_dir.name)
+        feat._Features__feat_index_table = {'fc8': np.array([0, 5, 11, 19])}
+        return feat
+
+    def test_feature_slice_with_feature_index_is_refused(self):
+        # The unit index addresses the flattened FULL feature space, so applying
+        # it to an already-sliced array would silently pick the wrong units.
+        feat = self._features_with_index()
+        with self.assertRaises(ValueError) as ctx:
+            feat.get('fc8', feature_slice=np.s_[0:10])
+        self.assertIn('feature_index', str(ctx.exception))
+
+        # ... including when labels are also given.
+        with self.assertRaises(ValueError):
+            feat.get('fc8', label=self.labels[0], feature_slice=np.s_[0:10])
+
+    def test_unsliced_get_still_works_with_feature_index(self):
+        # The guard must not break the supported path.
+        feat = self._features_with_index()
+        assert_array_equal(
+            feat.get('fc8'), self.stacked['fc8'][:, [0, 5, 11, 19]]
+        )
+
+    def test_iter_chunks_with_feature_index_is_refused(self):
+        feat = self._features_with_index()
+        with self.assertRaises(ValueError):
+            list(feat.iter_chunks('fc8'))
+
+    def test_missing_index_file_raises(self):
+        with self.assertRaises(RuntimeError):
+            Features(self.feature_dir.name, feature_index='/no/such/file.mat')
+
+
 class TestSaveFeature(unittest.TestCase):
     def test_save_feature_warns_future(self):
         # save_feature writes a MATLAB-compatible v7.3 .mat file, a path that is