@@ -243,10 +243,10 @@ def save_buffer(f: Path, buffer_name: str) -> list[Path]:
                 print("[ ] Dose not exist:", (folder), f"{' ':20}") if verbose else None
                 continue
             if (folder / buffer_name).exists():
-                import datetime
+                import datetime as dt
 
-                file_mod_time = datetime.datetime.fromtimestamp(os.path.getmtime(folder / buffer_name))
-                today = datetime.datetime.today()
+                file_mod_time = dt.datetime.fromtimestamp(os.path.getmtime(folder / buffer_name))
+                today = dt.datetime.today()
 
                 age = today - file_mod_time
                 if age.days >= int(max_age_days):
 
@@ -58,6 +58,7 @@
     np_unique,
     np_unique_withoutzero,
     np_volume,
+    np_voronoi_labels,
 )
 from TPTBox.core.vert_constants import (
     AFFINE,
@@ -1808,6 +1809,59 @@ def dilate_msk_(self, n_pixel: int = 5, labels: LABEL_REFERENCE = None, connecti
         """In-place variant of `dilate_msk`."""
         return self.dilate_msk(n_pixel=n_pixel, labels=labels, connectivity=connectivity, mask=mask, inplace=True, verbose=verbose,use_crop=use_crop,ignore_direction=ignore_direction)
 
+    def voronoi_labels(self, labels: LABEL_REFERENCE = None, max_distance: float | None = None, use_mm: bool = True, signed_background: bool = True, tie_to_zero: bool = True, inplace=False, verbose: logging = True) -> Self:
+        """Assigns every voxel the label of the nearest labelled region (a Voronoi partition).
+
+        Args:
+            labels (LABEL_REFERENCE, optional): Labels to partition by. Everything else is treated as
+                background. Defaults to None (all labels present in the segmentation).
+            max_distance (float | None, optional): Leave background voxels further than this from any
+                region at 0. In mm when `use_mm` is True, else in voxels. Defaults to None (no limit).
+            use_mm (bool, optional): If True, distances are measured in mm using this image's voxel
+                spacing, so the partition is correct on anisotropic images. If False, voxels are
+                treated as isotropic. Defaults to True.
+            signed_background (bool, optional): If True, voxels filled in from the background get the
+                negated label of their region, so the input mask stays recoverable as `out > 0`. The
+                returned NII then has a signed dtype. Defaults to True.
+            tie_to_zero (bool, optional): If True, voxels equidistant to two different regions are set
+                to 0 instead of being assigned to one of them. Defaults to True.
+            inplace (bool, optional): Whether to modify this NII in place or return a new object.
+                Defaults to False.
+            verbose (bool, optional): Whether to print a message. Defaults to True.
+
+        Returns:
+            NII: The partition, signed when `signed_background` is set.
+
+        Notes:
+            Runtime is independent of the number of labels (one Euclidean feature transform for the
+            whole volume), unlike `dilate_msk`.
+            Set `signed_background=False, tie_to_zero=False` and a `max_distance` for a plain,
+            unsigned expansion of the labels (`np_expand_labels`).
+
+        Examples:
+            >>> nii.voronoi_labels(max_distance=10)  # doctest: +SKIP
+            # every voxel within 1 cm of a labelled region takes its label
+        """
+        assert self.seg
+        log.print("voronoi labels", end="\r", verbose=verbose)
+        out = np_voronoi_labels(
+            self.get_seg_array(),
+            label_ref=labels,
+            max_distance=max_distance,
+            zoom=self.zoom if use_mm else None,
+            signed_background=signed_background,
+            tie_to_zero=tie_to_zero,
+        )
+        if not signed_background:
+            # keep the image dtype; the negated labels of a signed partition would not fit an unsigned one
+            out = out.astype(self.dtype)
+        log.print("Voronoi partition of", len(np_unique_withoutzero(np.abs(out))), "labels", verbose=verbose)
+        return self.set_array(out, inplace=inplace)
+
+    def voronoi_labels_(self, labels: LABEL_REFERENCE = None, max_distance: float | None = None, use_mm: bool = True, signed_background: bool = True, tie_to_zero: bool = True, verbose: logging = True) -> Self:
+        """In-place variant of `voronoi_labels`."""
+        return self.voronoi_labels(labels=labels, max_distance=max_distance, use_mm=use_mm, signed_background=signed_background, tie_to_zero=tie_to_zero, inplace=True, verbose=verbose)
+
 
     def fill_holes(self, labels: LABEL_REFERENCE = None, slice_wise_dim: int|str | None = None, verbose:logging=False, inplace=False,use_crop=True) -> Self:  # noqa: ARG002
         """Fills holes in segmentation.
@@ -2996,6 +3050,126 @@ def relabel_by_position_(self, axis: DIRECTIONS = "I", offset: int = 0, verbose:
         """In-place version of :meth:`relabel_by_position`."""
         return self.relabel_by_position(axis=axis, offset=offset, inplace=True, verbose=verbose)
 
+    def _split_axis(self, axis: DIRECTIONS | int) -> tuple[int, bool]:
+        """Resolves ``axis`` to an array axis and whether that direction runs with increasing index."""
+        if isinstance(axis, int):
+            return axis, True
+        # get_axis falls back to the opposite letter, so recover which way the axis actually runs.
+        return self.get_axis(axis), axis in self.orientation
+
+    def _split_label_profile(self, labels: LABEL_REFERENCE, ax: int) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
+        """Returns the seg array, the mask of ``labels`` and its per-slice voxel count along ``ax``."""
+        arr = self.get_seg_array()
+        selected = arr > 0 if labels is None else np_isin(arr, labels)
+        counts = selected.sum(axis=tuple(i for i in range(arr.ndim) if i != ax))
+        return arr, selected, counts
+
+    def _split_label_at_index(self, arr: np.ndarray, selected: np.ndarray, ax: int, region: slice, offset: int, inplace: bool) -> Self:
+        """Adds ``offset`` to the selected voxels that lie inside ``region`` along axis ``ax``."""
+        part = np.zeros_like(selected)
+        part[(slice(None),) * ax + (region,)] = True
+        part &= selected
+        if offset != 0 and part.any():
+            arr = arr.astype(np.promote_types(arr.dtype, np.min_scalar_type(int(arr.max()) + offset)), copy=False)
+            arr[part] += offset
+        return self.set_array(arr, inplace=inplace)
+
+    def partition_segmentation_by_volume(
+        self,
+        labels: LABEL_REFERENCE = None,
+        axis: DIRECTIONS | int = "P",
+        fraction: float = 0.5,
+        offset: int = 100,
+        inplace: bool = False,
+    ) -> Self:
+        """Splits labels along an axis at the plane that cuts off a given fraction of their volume.
+
+        Slices are walked inwards from the ``axis``-most end of the mask, summing voxels until
+        ``fraction`` of the total volume is covered; that leading slab gets ``offset`` added to its
+        label values, the rest is left untouched. With the default ``axis="P"`` and
+        ``fraction=1/3`` the posterior third (by volume, not by extent) is split off.
+
+        The direction is resolved against :attr:`orientation`, so no reorientation round-trip is
+        needed. Passing an ``int`` instead selects that array axis, pointing at high indices.
+
+        Args:
+            labels (LABEL_REFERENCE, optional): Label(s) to split. Defaults to None (every label).
+            axis (DIRECTIONS | int, optional): Anatomical direction the split-off part lies in,
+                or a raw array axis. Defaults to ``"P"``.
+            fraction (float, optional): Volume fraction that ends up in the split-off part, in
+                ``[0, 1]``. Defaults to 0.5.
+            offset (int, optional): Added to the labels of the split-off part. Defaults to 100.
+            inplace (bool, optional): If True, modifies this NII in place. Defaults to False.
+
+        Returns:
+            NII: The segmentation with the split-off part relabeled to ``label + offset``.
+
+        Examples:
+            >>> ethmoid.partition_segmentation_by_volume([7, 8], "P", 1 / 3)  # doctest: +SKIP
+            # posterior third of labels 7 and 8 becomes 107 and 108
+        """
+        assert 0 <= fraction <= 1, f"fraction has to be in [0,1], got {fraction}"
+        ax, forward = self._split_axis(axis)
+        arr, selected, counts = self._split_label_profile(labels, ax)
+        total = int(counts.sum())
+        if total == 0:
+            return self if inplace else self.copy()
+        n = len(counts)
+        # walk inwards from the axis-most end until the slab holds fraction of the volume
+        order = counts[::-1] if forward else counts
+        target = fraction * total
+        taken = 0 if target <= 0 else min(int(np.searchsorted(np.cumsum(order), target)) + 1, n)
+        region = slice(n - taken, None) if forward else slice(None, taken)
+        return self._split_label_at_index(arr, selected, ax, region, offset, inplace)
+
+    def partition_segmentation_by_center_of_mass(
+        self,
+        labels: LABEL_REFERENCE = None,
+        axis: DIRECTIONS | int = "L",
+        offset: int = 100,
+        inplace: bool = False,
+    ) -> Self:
+        """Splits labels along an axis at their center of mass.
+
+        The center of mass of ``labels`` is projected onto ``axis`` and rounded to the nearest
+        slice; everything from there towards ``axis`` gets ``offset`` added to its label values.
+        With the default ``axis="L"`` this separates the left half of a structure from its right
+        half, and flipping the letter to ``"R"`` splits off the other side. The slice the plane
+        lands on always goes to the ``axis`` side, so the two halves overlap in that one slice
+        when the center of mass sits exactly on a voxel center.
+
+        The direction is resolved against :attr:`orientation` and the rounding is done in the
+        anatomical frame, so the split is the same no matter which orientation the image is
+        stored in and no reorientation round-trip is needed. Passing an ``int`` instead selects
+        that array axis, pointing at high indices.
+
+        Args:
+            labels (LABEL_REFERENCE, optional): Label(s) to split. Defaults to None (every label).
+            axis (DIRECTIONS | int, optional): Anatomical direction the split-off part lies in,
+                or a raw array axis. Defaults to ``"L"``.
+            offset (int, optional): Added to the labels of the split-off part. Defaults to 100.
+            inplace (bool, optional): If True, modifies this NII in place. Defaults to False.
+
+        Returns:
+            NII: The segmentation with the split-off part relabeled to ``label + offset``.
+
+        Examples:
+            >>> ethmoid.partition_segmentation_by_center_of_mass([7, 15], "L")  # doctest: +SKIP
+            # left half of labels 7 and 15 becomes 107 and 115
+        """
+        ax, forward = self._split_axis(axis)
+        arr, selected, counts = self._split_label_profile(labels, ax)
+        total = int(counts.sum())
+        if total == 0:
+            return self if inplace else self.copy()
+        n = len(counts)
+        com = float((counts * np.arange(n)).sum()) / total
+        # round in the anatomical frame, not the array frame, so the result does not depend on
+        # whether the axis happens to be stored forwards or backwards
+        c = int(np.clip(np.rint(com if forward else n - 1 - com), 0, n))
+        region = slice(c, None) if forward else slice(None, n - c)
+        return self._split_label_at_index(arr, selected, ax, region, offset, inplace)
+
     def filter_connected_components_by_bbox_chain(
         self,
         margin_mm: float = 0.0,
 
@@ -24,6 +24,8 @@
     distance_transform_edt,
     gaussian_filter,
     generate_binary_structure,
+    maximum_filter,
+    minimum_filter,
     sobel,
 )
 from scipy.spatial.distance import cdist
@@ -598,6 +600,145 @@ def np_dilate_msk_euclid(arr: np.ndarray, n_pixel: int = 3, use_crop=True, label
     return arr
 
 
+def _np_voronoi_ties(arr: UINTARRAY, indices: np.ndarray, nearest: np.ndarray, zoom: Sequence[float] | None) -> np.ndarray:
+    """Finds the voxels that are equidistant to two differently labelled regions.
+
+    ``distance_transform_edt`` picks one nearest source per voxel and never reports that the choice
+    was arbitrary, so the ties have to be recovered afterwards. Only voxels next to a change in the
+    nearest-label assignment can sit on a boundary between two regions, so the search is restricted
+    to those and then each of their neighbours' nearest sources is tested for being at exactly the
+    same distance. Any neighbour's source is a real foreground voxel, so a hit is always a genuine
+    tie -- the test can miss one, but never invent one.
+
+    Args:
+        arr (UINTARRAY): The labelled input array.
+        indices (np.ndarray): Feature transform of ``arr``, shape ``(ndim, *arr.shape)``.
+        nearest (np.ndarray): ``arr`` gathered at ``indices``.
+        zoom (Sequence[float] | None): Voxel spacing used for the distances.
+
+    Returns:
+        np.ndarray: Boolean mask of the tied voxels.
+    """
+    ties = np.zeros(arr.shape, dtype=bool)
+    boundary = maximum_filter(nearest, size=3) != minimum_filter(nearest, size=3)
+    pts = np.argwhere(boundary)
+    if len(pts) == 0:
+        return ties
+    sampling = np.ones(arr.ndim) if zoom is None else np.asarray(zoom, dtype=float)
+    limit = np.asarray(arr.shape) - 1
+    key = (slice(None), *pts.T)
+    dist2 = (((indices[key].T - pts) * sampling) ** 2).sum(1)
+    label = nearest[tuple(pts.T)]
+    tied = np.zeros(len(pts), dtype=bool)
+    for offset in itertools.product((-1, 0, 1), repeat=arr.ndim):
+        if not any(offset):
+            continue
+        neighbour = np.clip(pts + offset, 0, limit)
+        source = indices[(slice(None), *neighbour.T)]
+        cand2 = (((source.T - pts) * sampling) ** 2).sum(1)
+        tied |= (arr[tuple(source)] != label) & (np.abs(cand2 - dist2) <= 1e-9 * np.maximum(dist2, 1.0))
+    ties[tuple(pts[tied].T)] = True
+    return ties
+
+
+def np_voronoi_labels(
+    arr: UINTARRAY,
+    label_ref: LABEL_REFERENCE = None,
+    max_distance: float | None = None,
+    zoom: Sequence[float] | None = None,
+    signed_background: bool = True,
+    tie_to_zero: bool = True,
+) -> INTARRAY:
+    """Assigns every voxel the label of the nearest labelled region (a Voronoi partition).
+
+    Computed with a single Euclidean feature transform of the background: the nearest background
+    voxel of ``arr == 0`` is exactly the nearest labelled voxel of ``arr``, so one pass yields the
+    assignment for the whole volume, independent of how many labels there are.
+
+    By default the partition is signed: voxels that were already labelled keep their positive
+    label, voxels filled in from the background carry the *negative* label of the region they were
+    assigned to, so the original mask can still be recovered from the result. Voxels that are
+    exactly equidistant to two different regions are set to 0 rather than being broken arbitrarily.
+
+    For a bounded expansion that overwrites the input in place and crops to the mask first, see
+    :func:`np_dilate_msk_euclid`; this function is the unbounded, spacing-aware variant.
+
+    Args:
+        arr (UINTARRAY): Labelled input array with background 0.
+        label_ref (int | list[int] | None, optional): Labels to partition by. Everything else is
+            treated as background. Defaults to None (all labels found in ``arr``).
+        max_distance (float | None, optional): Leave background voxels further than this from any
+            region at 0. In the units of ``zoom``, i.e. voxels when ``zoom`` is None.
+            Defaults to None (no limit).
+        zoom (Sequence[float] | None, optional): Voxel spacing, so distances are measured in mm on
+            anisotropic images. Defaults to None (isotropic voxels).
+        signed_background (bool, optional): If True, background voxels get the negated label of
+            their region. Defaults to True.
+        tie_to_zero (bool, optional): If True, voxels equidistant to two different regions are set
+            to 0 instead of being assigned to one of them. Defaults to True.
+
+    Returns:
+        INTARRAY: The partition, signed when ``signed_background`` is set.
+
+    Examples:
+        >>> arr = np.zeros((5, 1, 1), dtype=np.uint8)
+        >>> arr[0] = 1
+        >>> arr[4] = 2
+        >>> np_voronoi_labels(arr).ravel().tolist()
+        [1, -1, 0, -2, 2]
+    """
+    assert 2 <= arr.ndim <= 3, f"expected 2D or 3D, but got {arr.ndim}"
+    assert np.min(arr) >= 0, f"expected non-negative labels, got {np.min(arr)}"
+    assert zoom is None or len(zoom) == arr.ndim, f"zoom has to have {arr.ndim} entries, got {zoom}"
+    if label_ref is not None:
+        arr = np.where(np_isin(arr, label_ref), arr, 0)
+    background = arr == 0
+    if background.all():
+        return np.zeros(arr.shape, dtype=np.int8 if signed_background else arr.dtype)
+
+    want_distance = max_distance is not None
+    result = distance_transform_edt(background, sampling=zoom, return_distances=want_distance, return_indices=True)
+    distances, indices = result if want_distance else (None, result)
+    nearest = arr[tuple(indices)]
+
+    out = nearest.astype(np.min_scalar_type(-int(arr.max()) - 1) if signed_background else arr.dtype, copy=True)
+    if signed_background:
+        out[background] = -out[background]
+    if want_distance:
+        out[background & (distances > max_distance)] = 0
+    if tie_to_zero:
+        out[_np_voronoi_ties(arr, indices, nearest, zoom)] = 0
+    return out
+
+
+def np_expand_labels(
+    arr: UINTARRAY,
+    distance: float | None = None,
+    label_ref: LABEL_REFERENCE = None,
+    zoom: Sequence[float] | None = None,
+) -> UINTARRAY:
+    """Grows every label into the surrounding background, up to ``distance``.
+
+    Each background voxel within reach takes the label of the nearest region; labelled voxels are
+    left alone, so labels never eat into each other. This is :func:`np_voronoi_labels` with the
+    plain, unsigned semantics -- see there for the signed partition and the handling of voxels that
+    are equidistant to two regions.
+
+    Args:
+        arr (UINTARRAY): Labelled input array with background 0.
+        distance (float | None, optional): How far to grow, in the units of ``zoom``, i.e. voxels
+            when ``zoom`` is None. Defaults to None (fill the whole volume).
+        label_ref (int | list[int] | None, optional): Labels to grow. Everything else is treated as
+            background. Defaults to None (all labels found in ``arr``).
+        zoom (Sequence[float] | None, optional): Voxel spacing, so ``distance`` is in mm on
+            anisotropic images. Defaults to None (isotropic voxels).
+
+    Returns:
+        UINTARRAY: A copy of ``arr`` with the labels grown.
+    """
+    return np_voronoi_labels(arr, label_ref=label_ref, max_distance=distance, zoom=zoom, signed_background=False, tie_to_zero=False)
+
+
 def np_dilate_msk(
     arr: np.ndarray,
     label_ref: LABEL_REFERENCE = None,