Skip to content

Volume filtering

patchworks.filter_labels_by_size(store_path: str, component: str, min_voxels: 'int | None' = None, max_voxels: 'int | None' = None, *, relabel: bool = True) -> 'tuple[int, int]'

Drop label objects outside [min_voxels, max_voxels], in place.

Two-pass streaming scan (see module docstring) -- the array never has to fit in RAM.

Parameters:

Name Type Description Default
store_path str

Path to the zarr store containing the label array.

 required
component str

Array name inside the store to filter in place.

 required
min_voxels int

Objects with fewer voxels than this are zeroed (dropped). None (default) sets no lower bound. Use :func:min_voxels_for_volume to derive this from a physical volume and calibration.

 None
max_voxels int

Objects with more voxels than this are zeroed (dropped) -- e.g. a segmentation artifact where several objects merged into one giant blob. None (default) sets no upper bound. Use :func:max_voxels_for_volume to derive this from a physical volume and calibration.

 None
relabel bool

Renumber the surviving objects to a contiguous 1..N range in the same LUT that drops the out-of-range ones (default True) -- otherwise the removed ids leave permanent gaps and survivors keep their original ids.

 True

Returns:

Type Description
tuple of int

(n_kept, n_removed).

Examples:

>>> import zarr
>>> root = zarr.open_group("labels.zarr", mode="w")
>>> root.create_array(
...     "labels", shape=(4, 4), chunks=(4, 4), dtype="int32"
... )[:] = [
...     [0, 1, 1, 0],
...     [0, 1, 1, 0],
...     [0, 0, 0, 2],
...     [0, 0, 0, 0],
... ]
>>> filter_labels_by_size("labels.zarr", "labels", min_voxels=2)
(1, 1)
Source code in src/patchworks/_volume_filter.py 
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
def filter_labels_by_size(
    store_path: str,
    component: str,
    min_voxels: "int | None" = None,
    max_voxels: "int | None" = None,
    *,
    relabel: bool = True,
) -> "tuple[int, int]":
    """Drop label objects outside ``[min_voxels, max_voxels]``, in place.

    Two-pass streaming scan (see module docstring) -- the array never has
    to fit in RAM.

    Parameters
    ----------
    store_path : str
        Path to the zarr store containing the label array.
    component : str
        Array name inside the store to filter in place.
    min_voxels : int, optional
        Objects with fewer voxels than this are zeroed (dropped). ``None``
        (default) sets no lower bound. Use :func:`min_voxels_for_volume` to
        derive this from a physical volume and calibration.
    max_voxels : int, optional
        Objects with more voxels than this are zeroed (dropped) -- e.g. a
        segmentation artifact where several objects merged into one giant
        blob. ``None`` (default) sets no upper bound. Use
        :func:`max_voxels_for_volume` to derive this from a physical volume
        and calibration.
    relabel : bool, optional
        Renumber the surviving objects to a contiguous ``1..N`` range in
        the same LUT that drops the out-of-range ones (default ``True``)
        -- otherwise the removed ids leave permanent gaps and survivors
        keep their original ids.

    Returns
    -------
    tuple of int
        ``(n_kept, n_removed)``.

    Examples
    --------
    >>> import zarr
    >>> root = zarr.open_group("labels.zarr", mode="w")  # doctest: +SKIP
    >>> root.create_array(
    ...     "labels", shape=(4, 4), chunks=(4, 4), dtype="int32"
    ... )[:] = [
    ...     [0, 1, 1, 0],
    ...     [0, 1, 1, 0],
    ...     [0, 0, 0, 2],
    ...     [0, 0, 0, 0],
    ... ]  # doctest: +SKIP
    >>> filter_labels_by_size("labels.zarr", "labels", min_voxels=2)  # doctest: +SKIP
    (1, 1)
    """
    if min_voxels is None and max_voxels is None:
        raise ValueError(
            "filter_labels_by_size needs min_voxels, max_voxels, or both"
        )

    root = zarr.open_group(store_path, mode="r+")
    z = root[component]
    slices = _chunk_slices(z.shape, z.chunks)

    counts: "dict[int, int]" = {}
    for sl in slices:
        ids, n = np.unique(np.asarray(z[sl]), return_counts=True)
        for label_id, count in zip(ids.tolist(), n.tolist()):
            if label_id == 0:
                continue
            counts[label_id] = counts.get(label_id, 0) + count

    kept = sorted(
        i
        for i, c in counts.items()
        if (min_voxels is None or c >= min_voxels)
        and (max_voxels is None or c <= max_voxels)
    )
    n_kept = len(kept)
    n_removed = len(counts) - n_kept

    # Sized to the largest id *seen*, not just the largest surviving one --
    # a removed object's id can still exceed every kept id and must stay
    # in bounds so the LUT gather below maps it to 0 rather than indexing
    # past the end.
    max_label = max(counts) if counts else 0
    if max_label > _LUT_WARN_THRESHOLD:
        logger.warning(
            "filter_labels_by_size: max_label=%d -> LUT size ~%.0f MB.",
            max_label,
            max_label * 8 / 1024**2,
        )
    lut = np.zeros(max_label + 1, dtype=np.int64)
    if kept:
        lut[kept] = np.arange(1, n_kept + 1) if relabel else np.asarray(kept)

    max_out = n_kept if relabel else max_label
    out_dtype = np.uint16 if max_out < np.iinfo(np.uint16).max else np.uint32
    for sl in slices:
        block = np.asarray(z[sl])
        z[sl] = lut[block].astype(out_dtype)

    bounds = "-".join(
        str(v) if v is not None else "" for v in (min_voxels, max_voxels)
    )
    logger.info(
        "filter_labels_by_size: dropped %d/%d object(s) outside [%s] voxels, "
        "%d remain",
        n_removed,
        len(counts),
        bounds,
        n_kept,
    )
    return n_kept, n_removed

patchworks.voxel_volume(voxel_size: 'dict[str, float]') -> float

Physical volume of one voxel, from a per-axis calibration.

Axes missing from voxel_size are treated as 1.0 -- e.g. a 2-D calibration with no z gives an area, not a bogus volume shrunk by a fake axis. Units follow whatever voxel_size is in (micrometers for :func:patchworks.plugins.ome_zarr.read_pixel_size).

Parameters:

Name Type Description Default
voxel_size dict

Per-axis physical size, e.g. {"z": .., "y": .., "x": ..}.

 required

Returns:

Type Description
float

Product of the given axis sizes.

Examples:

>>> voxel_volume({"z": 0.24, "y": 0.10833, "x": 0.10833})
0.0028164933359999997
Source code in src/patchworks/_volume_filter.py 
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
def voxel_volume(voxel_size: "dict[str, float]") -> float:
    """Physical volume of one voxel, from a per-axis calibration.

    Axes missing from *voxel_size* are treated as 1.0 -- e.g. a 2-D
    calibration with no ``z`` gives an area, not a bogus volume shrunk by a
    fake axis. Units follow whatever *voxel_size* is in (micrometers for
    :func:`patchworks.plugins.ome_zarr.read_pixel_size`).

    Parameters
    ----------
    voxel_size : dict
        Per-axis physical size, e.g. ``{"z": .., "y": .., "x": ..}``.

    Returns
    -------
    float
        Product of the given axis sizes.

    Examples
    --------
    >>> voxel_volume({"z": 0.24, "y": 0.10833, "x": 0.10833})
    0.0028164933359999997
    """
    vol = 1.0
    for size in voxel_size.values():
        vol *= size
    return vol

patchworks.min_voxels_for_volume(min_volume: float, voxel_size: 'dict[str, float]') -> int

Convert a physical volume threshold to a voxel count.

Rounds up: an object must reach min_volume to survive, so a partial voxel's worth of extra volume should not tip it over the line.

Parameters:

Name Type Description Default
min_volume float

Minimum object volume to keep, in the same physical units as voxel_size (micrometers³ for an NGFF calibration).

 required
voxel_size dict

Per-axis physical size -- see :func:voxel_volume.

 required

Returns:

Type Description
int

Minimum voxel count for an object to survive filtering.

Examples:

>>> min_voxels_for_volume(5.0, {"z": 0.24, "y": 0.10833, "x": 0.10833})
1776
Source code in src/patchworks/_volume_filter.py 
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
def min_voxels_for_volume(
    min_volume: float, voxel_size: "dict[str, float]"
) -> int:
    """Convert a physical volume threshold to a voxel count.

    Rounds up: an object must reach *min_volume* to survive, so a partial
    voxel's worth of extra volume should not tip it over the line.

    Parameters
    ----------
    min_volume : float
        Minimum object volume to keep, in the same physical units as
        *voxel_size* (micrometers³ for an NGFF calibration).
    voxel_size : dict
        Per-axis physical size -- see :func:`voxel_volume`.

    Returns
    -------
    int
        Minimum voxel count for an object to survive filtering.

    Examples
    --------
    >>> min_voxels_for_volume(5.0, {"z": 0.24, "y": 0.10833, "x": 0.10833})
    1776
    """
    return math.ceil(min_volume / voxel_volume(voxel_size))

patchworks.max_voxels_for_volume(max_volume: float, voxel_size: 'dict[str, float]') -> int

Convert a physical volume threshold to a voxel count.

Rounds down, the mirror image of :func:min_voxels_for_volume: an object must not exceed max_volume, so a voxel count whose volume would tip past it must not survive.

Parameters:

Name Type Description Default
max_volume float

Maximum object volume to keep, in the same physical units as voxel_size (micrometers³ for an NGFF calibration).

 required
voxel_size dict

Per-axis physical size -- see :func:voxel_volume.

 required

Returns:

Type Description
int

Maximum voxel count for an object to survive filtering.

Examples:

>>> max_voxels_for_volume(5.0, {"z": 0.24, "y": 0.10833, "x": 0.10833})
1775
Source code in src/patchworks/_volume_filter.py 
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
def max_voxels_for_volume(
    max_volume: float, voxel_size: "dict[str, float]"
) -> int:
    """Convert a physical volume threshold to a voxel count.

    Rounds down, the mirror image of :func:`min_voxels_for_volume`: an
    object must not *exceed* *max_volume*, so a voxel count whose volume
    would tip past it must not survive.

    Parameters
    ----------
    max_volume : float
        Maximum object volume to keep, in the same physical units as
        *voxel_size* (micrometers³ for an NGFF calibration).
    voxel_size : dict
        Per-axis physical size -- see :func:`voxel_volume`.

    Returns
    -------
    int
        Maximum voxel count for an object to survive filtering.

    Examples
    --------
    >>> max_voxels_for_volume(5.0, {"z": 0.24, "y": 0.10833, "x": 0.10833})
    1775
    """
    return math.floor(max_volume / voxel_volume(voxel_size))