From c7e96fc69f5cb909188854a24f1c9ffaca9ca762 Mon Sep 17 00:00:00 2001
From: uermel <utz.ermel@czii.org>
Date: Mon, 5 Oct 2026 18:03:46 -0700
Subject: [PATCH 1/3] chore: build against the copick filament stack (drop at
 release)

copick-utils' filament, instance and pick-identity work uses APIs from the open
copick PRs #479, #480, #492 and #494 (filament objects, pick identity, the
Filaments entity with curves, instance and panoptic segmentations). Until
those are released, depend on their stack's tip. Replace this with a version
floor once copick publishes them.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 pyproject.toml | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/pyproject.toml b/pyproject.toml
index 1683ae9..c2c13c8 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -18,7 +18,9 @@ classifiers = [
 ]
 dynamic = ["version"]
 dependencies = [
-    "copick>=1.26.0",
+    # TODO:remove before release -- build against the unreleased copick filament stack (copick#479/#480/#492/#494);
+    # restore "copick>=<release with those PRs>" when it is published.
+    "copick @ git+https://github.com/copick/copick@uermel/panoptic-segmentations",
     "click",
     "click-option-group",
     "numpy",

From c0fe02e0f3723c1249b835d96d0e3d2e630c6bf1 Mon Sep 17 00:00:00 2001
From: uermel <utz.ermel@czii.org>
Date: Mon, 5 Oct 2026 18:09:02 -0700
Subject: [PATCH 2/3] fix: pick tools keep identity and place particles at
 their centre

picksin, picksout and clippicks rebuilt their output from picks.numpy()/
from_numpy(), which dropped every pick's instance_id and score. For a filament
object that loses which filament a pick belongs to. They also tested the
pick's location, not its particle centre: copick defines the centre as the
location plus the transform's translation (a shift, e.g. from a RELION
import). Now they test the centre and copy the surviving points unchanged and
in order, so orientation, shift, instance_id and score survive.

Also:
- An empty result writes an empty pick set. Before, the run was reported as an
  error and an earlier output in the same place was left on disk.
- An instance segmentation works as a reference (any non-zero ID is inside),
  and a panoptic one is tested on its label channel.
- clippicks --ref-seg built its distance field on the segmentation's (z, y, x)
  grid but indexed it with (x, y, z) indices, which is wrong for any non-cubic
  volume. It now indexes (z, y, x) and rounds to the nearest voxel, matching
  copick's voxel-to-Angstrom mapping (voxel i at i * spacing).
- picks2seg, picks2mesh, picks2sphere, picks2ellipsoid, picks2plane,
  picks2surface and io.readers.coordinates use particle centres.

The shared helpers are in copick_utils.util.picks. Covered by
tests/test_pick_tools.py; every case fails on the previous code.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/copick_utils/cli/clippicks.py             |   5 +
 src/copick_utils/cli/picksin.py               |   5 +
 src/copick_utils/cli/picksout.py              |   5 +
 .../converters/ellipsoid_from_picks.py        |   8 +-
 .../converters/mesh_from_picks.py             |   8 +-
 .../converters/plane_from_picks.py            |   8 +-
 .../converters/segmentation_from_picks.py     |  13 +-
 .../converters/sphere_from_picks.py           |   8 +-
 .../converters/surface_from_picks.py          |   8 +-
 src/copick_utils/io/readers.py                |  24 +-
 .../logical/distance_operations.py            |  55 ++--
 src/copick_utils/logical/point_operations.py  | 237 ++++++++----------
 src/copick_utils/util/picks.py                | 106 ++++++++
 tests/conftest.py                             |  43 ++++
 tests/test_pick_tools.py                      | 162 ++++++++++++
 15 files changed, 483 insertions(+), 212 deletions(-)
 create mode 100644 src/copick_utils/util/picks.py
 create mode 100644 tests/conftest.py
 create mode 100644 tests/test_pick_tools.py

diff --git a/src/copick_utils/cli/clippicks.py b/src/copick_utils/cli/clippicks.py
index f05defd..29df8ae 100644
--- a/src/copick_utils/cli/clippicks.py
+++ b/src/copick_utils/cli/clippicks.py
@@ -65,6 +65,11 @@ def clippicks(
     Exactly one reference must be supplied via `--ref-mesh`, `--ref-seg`, or
     `--ref-tomogram`. Distances are measured in angstroms.
 
+    Each pick is measured at its particle centre (its location plus the shift stored in its
+    transform). Kept picks are written unchanged and in their original order, so their
+    orientations, instance IDs (e.g. filament IDs) and scores survive. If no pick is kept, an
+    empty pick set is written.
+
     URI Format:
 
         \b
diff --git a/src/copick_utils/cli/picksin.py b/src/copick_utils/cli/picksin.py
index 4287a41..0239d18 100644
--- a/src/copick_utils/cli/picksin.py
+++ b/src/copick_utils/cli/picksin.py
@@ -50,6 +50,11 @@ def picksin(
     whose coordinates fall inside it are written to the output. This is the inverse of
     `copick logical picksout`, which keeps the picks that fall outside the volume instead.
 
+    Each pick is tested at its particle centre (its location plus the shift stored in its
+    transform). Kept picks are written unchanged and in their original order, so their
+    orientations, instance IDs (e.g. filament IDs) and scores survive. If no pick is kept, an
+    empty pick set is written.
+
     URI Format:
 
         \b
diff --git a/src/copick_utils/cli/picksout.py b/src/copick_utils/cli/picksout.py
index ae720d5..ed42bcf 100644
--- a/src/copick_utils/cli/picksout.py
+++ b/src/copick_utils/cli/picksout.py
@@ -52,6 +52,11 @@ def picksout(
     whose coordinates fall outside it are written to the output. This is the inverse of
     `copick logical picksin`, which keeps the picks that fall inside the volume instead.
 
+    Each pick is tested at its particle centre (its location plus the shift stored in its
+    transform). Kept picks are written unchanged and in their original order, so their
+    orientations, instance IDs (e.g. filament IDs) and scores survive. If no pick is kept, an
+    empty pick set is written.
+
     URI Format:
 
         \b
diff --git a/src/copick_utils/converters/ellipsoid_from_picks.py b/src/copick_utils/converters/ellipsoid_from_picks.py
index e073537..661d0d2 100644
--- a/src/copick_utils/converters/ellipsoid_from_picks.py
+++ b/src/copick_utils/converters/ellipsoid_from_picks.py
@@ -11,6 +11,7 @@
     validate_points,
 )
 from copick_utils.converters.lazy_converter import create_lazy_batch_converter
+from copick_utils.util.picks import pick_centres
 
 if TYPE_CHECKING:
     from copick.models import CopickMesh, CopickRun
@@ -317,9 +318,10 @@ def create_ellipsoid_from_points(cluster_points):
 
 
 def ellipsoid_from_picks_standard(picks, run, object_name, session_id, user_id, voxel_spacing=None, **kwargs):
-    """Adapt ellipsoid_from_picks to the lazy-batch contract: take a CopickPicks and extract its points."""
-    pos, _ = picks.numpy()
-    if pos is None or len(pos) == 0:
+    """Adapt ellipsoid_from_picks to the lazy-batch contract: take a CopickPicks and extract its particle centres
+    (``location`` plus the transform's shift)."""
+    pos = pick_centres(picks)
+    if len(pos) == 0:
         logger.warning("Could not load pick data")
         return None
     return ellipsoid_from_picks(pos, run, object_name, session_id, user_id, **kwargs)
diff --git a/src/copick_utils/converters/mesh_from_picks.py b/src/copick_utils/converters/mesh_from_picks.py
index 8677c1b..8f8a5c0 100644
--- a/src/copick_utils/converters/mesh_from_picks.py
+++ b/src/copick_utils/converters/mesh_from_picks.py
@@ -10,6 +10,7 @@
     validate_points,
 )
 from copick_utils.converters.lazy_converter import create_lazy_batch_converter
+from copick_utils.util.picks import pick_centres
 
 if TYPE_CHECKING:
     from copick.models import CopickMesh, CopickRun
@@ -189,9 +190,10 @@ def create_mesh_from_points(cluster_points):
 
 # Lazy batch converter for new architecture
 def mesh_from_picks_standard(picks, run, object_name, session_id, user_id, voxel_spacing=None, **kwargs):
-    """Adapt mesh_from_picks to the lazy-batch contract: take a CopickPicks and extract its points."""
-    pos, _ = picks.numpy()
-    if pos is None or len(pos) == 0:
+    """Adapt mesh_from_picks to the lazy-batch contract: take a CopickPicks and extract its particle centres
+    (``location`` plus the transform's shift)."""
+    pos = pick_centres(picks)
+    if len(pos) == 0:
         logger.warning("Could not load pick data")
         return None
     return mesh_from_picks(pos, run, object_name, session_id, user_id, **kwargs)
diff --git a/src/copick_utils/converters/plane_from_picks.py b/src/copick_utils/converters/plane_from_picks.py
index d8f6271..cc363d9 100644
--- a/src/copick_utils/converters/plane_from_picks.py
+++ b/src/copick_utils/converters/plane_from_picks.py
@@ -11,6 +11,7 @@
     validate_points,
 )
 from copick_utils.converters.lazy_converter import create_lazy_batch_converter
+from copick_utils.util.picks import pick_centres
 
 if TYPE_CHECKING:
     from copick.models import CopickMesh, CopickRun
@@ -232,9 +233,10 @@ def create_plane_from_points(cluster_points):
 
 # Lazy batch converter for new architecture
 def plane_from_picks_standard(picks, run, object_name, session_id, user_id, voxel_spacing=None, **kwargs):
-    """Adapt plane_from_picks to the lazy-batch contract: take a CopickPicks and extract its points."""
-    pos, _ = picks.numpy()
-    if pos is None or len(pos) == 0:
+    """Adapt plane_from_picks to the lazy-batch contract: take a CopickPicks and extract its particle centres
+    (``location`` plus the transform's shift)."""
+    pos = pick_centres(picks)
+    if len(pos) == 0:
         logger.warning("Could not load pick data")
         return None
     return plane_from_picks(pos, run, object_name, session_id, user_id, **kwargs)
diff --git a/src/copick_utils/converters/segmentation_from_picks.py b/src/copick_utils/converters/segmentation_from_picks.py
index ce75397..095ffc7 100644
--- a/src/copick_utils/converters/segmentation_from_picks.py
+++ b/src/copick_utils/converters/segmentation_from_picks.py
@@ -5,6 +5,7 @@
 from copick.util.log import get_logger
 
 from copick_utils.converters.lazy_converter import create_lazy_batch_converter
+from copick_utils.util.picks import pick_centres
 
 if TYPE_CHECKING:
     from copick.models import CopickObject, CopickPicks, CopickRun, CopickSegmentation
@@ -54,14 +55,10 @@ def get_relative_target_coordinates(center, delta, shape):
     radius_voxel = max(radius / voxel_spacing, 1)
     delta = int(np.ceil(radius_voxel))
 
-    # Paint each pick as a sphere
-    for point in pick.points:
-        # Convert the pick's location from angstroms to voxel units
-        cx, cy, cz = (
-            point.location.x / voxel_spacing,
-            point.location.y / voxel_spacing,
-            point.location.z / voxel_spacing,
-        )
+    # Paint each pick as a sphere at its particle centre (location plus the transform's shift)
+    for centre in pick_centres(pick):
+        # Convert the particle centre from angstroms to voxel units
+        cx, cy, cz = centre / voxel_spacing
 
         # Calculate subarray bounds
         xLow, xHigh = get_relative_target_coordinates(cx, delta, seg_volume.shape[2])
diff --git a/src/copick_utils/converters/sphere_from_picks.py b/src/copick_utils/converters/sphere_from_picks.py
index ca81158..b4bc879 100644
--- a/src/copick_utils/converters/sphere_from_picks.py
+++ b/src/copick_utils/converters/sphere_from_picks.py
@@ -11,6 +11,7 @@
     validate_points,
 )
 from copick_utils.converters.lazy_converter import create_lazy_batch_converter
+from copick_utils.util.picks import pick_centres
 
 if TYPE_CHECKING:
     from copick.models import CopickMesh, CopickRun
@@ -287,9 +288,10 @@ def sphere_from_picks(
 
 # Lazy batch converter for new architecture
 def sphere_from_picks_standard(picks, run, object_name, session_id, user_id, voxel_spacing=None, **kwargs):
-    """Adapt sphere_from_picks to the lazy-batch contract: take a CopickPicks and extract its points."""
-    pos, _ = picks.numpy()
-    if pos is None or len(pos) == 0:
+    """Adapt sphere_from_picks to the lazy-batch contract: take a CopickPicks and extract its particle centres
+    (``location`` plus the transform's shift)."""
+    pos = pick_centres(picks)
+    if len(pos) == 0:
         logger.warning("Could not load pick data")
         return None
     return sphere_from_picks(pos, run, object_name, session_id, user_id, **kwargs)
diff --git a/src/copick_utils/converters/surface_from_picks.py b/src/copick_utils/converters/surface_from_picks.py
index 2c33399..3f17026 100644
--- a/src/copick_utils/converters/surface_from_picks.py
+++ b/src/copick_utils/converters/surface_from_picks.py
@@ -13,6 +13,7 @@
     validate_points,
 )
 from copick_utils.converters.lazy_converter import create_lazy_batch_converter
+from copick_utils.util.picks import pick_centres
 
 if TYPE_CHECKING:
     from copick.models import CopickMesh, CopickRun
@@ -318,9 +319,10 @@ def create_surface_from_points(cluster_points):
 
 # Lazy batch converter for new architecture
 def surface_from_picks_standard(picks, run, object_name, session_id, user_id, voxel_spacing=None, **kwargs):
-    """Adapt surface_from_picks to the lazy-batch contract: take a CopickPicks and extract its points."""
-    pos, _ = picks.numpy()
-    if pos is None or len(pos) == 0:
+    """Adapt surface_from_picks to the lazy-batch contract: take a CopickPicks and extract its particle centres
+    (``location`` plus the transform's shift)."""
+    pos = pick_centres(picks)
+    if len(pos) == 0:
         logger.warning("Could not load pick data")
         return None
     return surface_from_picks(pos, run, object_name, session_id, user_id, **kwargs)
diff --git a/src/copick_utils/io/readers.py b/src/copick_utils/io/readers.py
index 118f299..27df162 100644
--- a/src/copick_utils/io/readers.py
+++ b/src/copick_utils/io/readers.py
@@ -1,6 +1,7 @@
-import numpy as np
 from copick.util.uri import resolve_copick_objects
 
+from copick_utils.util.picks import pick_centres
+
 
 def tomogram(run, voxel_size: float = 10, algorithm: str = "wbp", raise_error: bool = False, verbose=True):
     """Read a tomogram from a copick run as a NumPy array.
@@ -143,7 +144,8 @@ def coordinates(
     """Read pick coordinates from a copick run as an array of voxel indices.
 
     Looks up the picks for `name` (optionally filtered by `user_id` and
-    `session_id`) and returns their locations divided by `voxel_size`. If several
+    `session_id`) and returns their particle centres (location plus the transform's
+    shift) divided by `voxel_size`. If several
     pick sets match, the first is used; if none match, reports the available picks.
 
     Args:
@@ -194,19 +196,5 @@ def coordinates(
                 f"Defaulting to loading:\n {picks[0]}\n",
             )
 
-    points = picks[0].points
-
-    # Initialize an array to store the coordinates
-    nPoints = len(picks[0].points)  # Number of points retrieved
-    coordinates = np.zeros([len(picks[0].points), 3])  # Create an empty array to hold the (z, y, x) coordinates
-
-    # Iterate over all points and convert their locations to coordinates in voxel space
-    for ii in range(nPoints):
-        coordinates[ii,] = [
-            points[ii].location.z / voxel_size,  # Scale z-coordinate by voxel size
-            points[ii].location.y / voxel_size,  # Scale y-coordinate by voxel size
-            points[ii].location.x / voxel_size,
-        ]  # Scale x-coordinate by voxel size
-
-    # Return the array of coordinates
-    return coordinates
+    # Particle centres (location plus the transform's shift), converted to (z, y, x) voxel coordinates
+    return pick_centres(picks[0])[:, ::-1] / voxel_size
diff --git a/src/copick_utils/logical/distance_operations.py b/src/copick_utils/logical/distance_operations.py
index f2dd207..2425ce0 100644
--- a/src/copick_utils/logical/distance_operations.py
+++ b/src/copick_utils/logical/distance_operations.py
@@ -10,6 +10,7 @@
     store_mesh_with_stats,
 )
 from copick_utils.converters.lazy_converter import create_lazy_batch_converter
+from copick_utils.util.picks import point_centres, segmentation_label_volume, store_pick_subset, voxel_indices
 
 if TYPE_CHECKING:
     from copick.models import CopickMesh, CopickPicks, CopickRun, CopickSegmentation
@@ -586,6 +587,10 @@ def limit_picks_by_distance(
     """
     Limit picks to those within a certain distance of a reference surface.
 
+    A pick is measured at its particle centre (``location`` plus the transform's shift). Kept picks are written
+    unchanged and in their original order: transform, ``instance_id`` and ``score`` survive. If no pick passes, an
+    empty pick set is written.
+
     Args:
         picks: CopickPicks to limit
         run: CopickRun object
@@ -610,13 +615,14 @@ def limit_picks_by_distance(
                 "One of reference_mesh, reference_segmentation, or reference_tomogram_info must be provided",
             )
 
-        # Load pick data
-        points, transforms = picks.numpy()
-        if points is None or len(points) == 0:
+        # Load pick data; picks are tested at their particle centres (location plus the transform's shift)
+        points = picks.points or []
+        if len(points) == 0:
             logger.error("Could not load pick data")
             return None
 
-        pick_positions = points[:, :3]  # Use only x, y, z coordinates
+        pick_positions = point_centres(points)
+        within_or_beyond = "beyond" if invert else "within"
 
         # Handle tomogram boundary reference (direct distance calculation, no distance field)
         if reference_tomogram_info is not None:
@@ -627,27 +633,19 @@ def limit_picks_by_distance(
             # Apply invert logic: default keeps within distance, invert keeps beyond
             final_valid = pick_distances > max_distance if invert else pick_distances <= max_distance
 
-            if not np.any(final_valid):
-                within_or_beyond = "beyond" if invert else "within"
-                logger.warning(f"No picks {within_or_beyond} {max_distance} units of tomogram boundary")
-                return None
-
         # Handle segmentation or mesh reference (uses distance field)
         else:
             # Create distance field from reference
             if reference_segmentation is not None:
-                ref_seg_array = reference_segmentation.numpy()
+                ref_seg_array = segmentation_label_volume(reference_segmentation)
                 if ref_seg_array is None or ref_seg_array.size == 0:
                     logger.error("Could not load reference segmentation data")
                     return None
 
-                # Use reference segmentation's coordinate space
+                # The field is on the reference segmentation's own (z, y, x) grid; voxel i sits at i * spacing.
                 field_voxel_spacing = reference_segmentation.voxel_size
                 distance_field = _create_distance_field_from_segmentation(ref_seg_array, field_voxel_spacing)
-
-                # Convert pick coordinates to voxel indices in reference segmentation space
-                pick_voxel_coords = pick_positions / field_voxel_spacing
-                pick_voxel_indices = np.floor(pick_voxel_coords).astype(int)
+                pick_voxel_indices = voxel_indices(pick_positions, field_voxel_spacing)
 
             else:  # reference_mesh is not None
                 ref_mesh = reference_mesh.mesh
@@ -674,7 +672,7 @@ def limit_picks_by_distance(
                 field_size = coord_bounds[1] - coord_bounds[0]
                 field_shape = np.ceil(field_size / field_voxel_spacing).astype(int)
 
-                # Create distance field from mesh in this coordinate space.
+                # Create distance field from mesh in this coordinate space (axis order x, y, z).
                 # Pass origin=coord_bounds[0] so the field aligns with the pick-index origin below.
                 distance_field = _create_distance_field_from_mesh(
                     ref_mesh,
@@ -693,10 +691,6 @@ def limit_picks_by_distance(
                 axis=1,
             )
 
-            if not np.any(valid_picks):
-                logger.warning("No picks within distance field bounds")
-                return None
-
             # Get distances for valid picks
             valid_indices = pick_voxel_indices[valid_picks]
             pick_distances = distance_field[valid_indices[:, 0], valid_indices[:, 1], valid_indices[:, 2]]
@@ -708,22 +702,15 @@ def limit_picks_by_distance(
             final_valid = np.zeros(len(points), dtype=bool)
             final_valid[valid_picks] = distance_valid
 
-            if not np.any(final_valid):
-                within_or_beyond = "beyond" if invert else "within"
-                logger.warning(f"No picks {within_or_beyond} {max_distance} units of reference surface")
-                return None
-
-        # Filter picks
-        valid_points = points[final_valid]
-        valid_transforms = transforms[final_valid] if transforms is not None else None
+        if not np.any(final_valid):
+            logger.warning(
+                f"No picks {within_or_beyond} {max_distance} units of the reference; writing an empty pick set",
+            )
 
-        # Create output picks
-        output_picks = run.new_picks(object_name, session_id, user_id, exist_ok=True)
-        output_picks.from_numpy(positions=valid_points, transforms=valid_transforms)
-        output_picks.store()
+        # Kept picks are written unchanged and in order (transform, instance_id and score survive)
+        output_picks = store_pick_subset(points, final_valid, run, object_name, session_id, user_id)
 
-        stats = {"points_created": len(valid_points)}
-        within_or_beyond = "beyond" if invert else "within"
+        stats = {"points_created": int(np.count_nonzero(final_valid))}
         logger.info(f"Limited picks to {stats['points_created']} points {within_or_beyond} {max_distance} units")
         return output_picks, stats
 
diff --git a/src/copick_utils/logical/point_operations.py b/src/copick_utils/logical/point_operations.py
index 7d5a2bb..1c65831 100644
--- a/src/copick_utils/logical/point_operations.py
+++ b/src/copick_utils/logical/point_operations.py
@@ -7,6 +7,7 @@
 from copick.util.log import get_logger
 
 from copick_utils.converters.lazy_converter import create_lazy_batch_converter
+from copick_utils.util.picks import point_centres, segmentation_label_volume, store_pick_subset, voxel_indices
 
 if TYPE_CHECKING:
     from copick.models import CopickMesh, CopickPicks, CopickRun, CopickSegmentation
@@ -55,36 +56,83 @@ def _check_points_in_segmentation(
     Check which points are inside a segmentation volume.
 
     Args:
-        points: Array of points to check (N, 3) in physical coordinates
-        segmentation_array: Binary segmentation array
+        points: Array of points to check (N, 3) in physical (x, y, z) coordinates
+        segmentation_array: Segmentation array (z, y, x); any non-zero voxel is inside (binary, multilabel, or an
+            instance segmentation's IDs)
         voxel_spacing: Spacing between voxels
 
     Returns:
         Boolean array indicating which points are inside the segmentation
     """
-    # Convert points to voxel coordinates
-    voxel_coords = np.round(points / voxel_spacing).astype(int)
-
-    # Check bounds
-    valid_bounds = (
-        (voxel_coords[:, 0] >= 0)
-        & (voxel_coords[:, 0] < segmentation_array.shape[2])
-        & (voxel_coords[:, 1] >= 0)
-        & (voxel_coords[:, 1] < segmentation_array.shape[1])
-        & (voxel_coords[:, 2] >= 0)
-        & (voxel_coords[:, 2] < segmentation_array.shape[0])
-    )
+    zyx = voxel_indices(points, voxel_spacing)
+    valid_bounds = np.all((zyx >= 0) & (zyx < np.array(segmentation_array.shape)), axis=1)
 
     inside = np.zeros(len(points), dtype=bool)
+    valid = zyx[valid_bounds]
+    if len(valid) > 0:
+        inside[valid_bounds] = segmentation_array[valid[:, 0], valid[:, 1], valid[:, 2]] > 0
+    return inside
 
-    # Check only points within bounds
-    valid_coords = voxel_coords[valid_bounds]
-    if len(valid_coords) > 0:
-        # Check if voxels are non-zero (inside segmentation)
-        voxel_values = segmentation_array[valid_coords[:, 2], valid_coords[:, 1], valid_coords[:, 0]]
-        inside[valid_bounds] = voxel_values > 0
 
-    return inside
+def _inside_reference(
+    centres: np.ndarray,
+    reference_mesh: Optional["CopickMesh"],
+    reference_segmentation: Optional["CopickSegmentation"],
+) -> Optional[np.ndarray]:
+    """Which particle centres lie inside the reference mesh or segmentation (None if the reference cannot be read)."""
+    if reference_mesh is not None:
+        ref_mesh = reference_mesh.mesh
+        if ref_mesh is None:
+            logger.error("Could not load reference mesh data")
+            return None
+
+        if isinstance(ref_mesh, tm.Scene):
+            if len(ref_mesh.geometry) == 0:
+                logger.error("Reference mesh is empty")
+                return None
+            ref_mesh = tm.util.concatenate(list(ref_mesh.geometry.values()))
+
+        return _check_points_in_mesh(centres, ref_mesh)
+
+    ref_seg_array = segmentation_label_volume(reference_segmentation)
+    if ref_seg_array is None or ref_seg_array.size == 0:
+        logger.error("Could not load reference segmentation data")
+        return None
+    return _check_points_in_segmentation(centres, ref_seg_array, reference_segmentation.voxel_size)
+
+
+def _filter_picks_by_reference(
+    picks: "CopickPicks",
+    run: "CopickRun",
+    object_name: str,
+    session_id: str,
+    user_id: str,
+    reference_mesh: Optional["CopickMesh"],
+    reference_segmentation: Optional["CopickSegmentation"],
+    keep_inside: bool,
+) -> Optional[Tuple["CopickPicks", Dict[str, int]]]:
+    """Keep the picks whose particle centres are inside (or outside) a reference, unchanged and in order."""
+    if reference_mesh is None and reference_segmentation is None:
+        raise ValueError("Either reference_mesh or reference_segmentation must be provided")
+
+    points = picks.points or []
+    if len(points) == 0:
+        logger.error("Could not load pick data")
+        return None
+
+    inside_mask = _inside_reference(point_centres(points), reference_mesh, reference_segmentation)
+    if inside_mask is None:
+        return None
+
+    keep = inside_mask if keep_inside else ~inside_mask
+    where = "inside" if keep_inside else "outside"
+    if not np.any(keep):
+        logger.warning(f"No picks found {where} reference volume; writing an empty pick set")
+
+    output_picks = store_pick_subset(points, keep, run, object_name, session_id, user_id)
+    stats = {"points_created": int(np.count_nonzero(keep))}
+    logger.info(f"Kept {stats['points_created']} of {len(points)} picks {where} reference volume")
+    return output_picks, stats
 
 
 def picks_inclusion_by_mesh(
@@ -100,10 +148,15 @@ def picks_inclusion_by_mesh(
     """
     Filter picks to include only those inside a reference mesh or segmentation.
 
+    A pick is tested at its particle centre (``location`` plus the transform's shift). Kept picks are written
+    unchanged and in their original order: transform, ``instance_id`` and ``score`` survive. If no pick is inside,
+    an empty pick set is written.
+
     Args:
         picks: CopickPicks to filter
         reference_mesh: Reference CopickMesh (either this or reference_segmentation must be provided)
-        reference_segmentation: Reference CopickSegmentation
+        reference_segmentation: Reference CopickSegmentation (binary, multilabel, instance, or the label channel
+            of a panoptic segmentation; any non-zero voxel is inside)
         run: CopickRun object
         object_name: Name for the output picks
         session_id: Session ID for the output picks
@@ -115,61 +168,16 @@ def picks_inclusion_by_mesh(
         Stats dict contains 'points_created'.
     """
     try:
-        if reference_mesh is None and reference_segmentation is None:
-            raise ValueError("Either reference_mesh or reference_segmentation must be provided")
-
-        # Load pick data
-        points, transforms = picks.numpy()
-        if points is None or len(points) == 0:
-            logger.error("Could not load pick data")
-            return None
-
-        pick_positions = points[:, :3]  # Use only x, y, z coordinates
-
-        # Check which points are inside the reference
-        if reference_mesh is not None:
-            ref_mesh = reference_mesh.mesh
-            if ref_mesh is None:
-                logger.error("Could not load reference mesh data")
-                return None
-
-            if isinstance(ref_mesh, tm.Scene):
-                if len(ref_mesh.geometry) == 0:
-                    logger.error("Reference mesh is empty")
-                    return None
-                ref_mesh = tm.util.concatenate(list(ref_mesh.geometry.values()))
-
-            inside_mask = _check_points_in_mesh(pick_positions, ref_mesh)
-
-        else:  # reference_segmentation is not None
-            ref_seg_array = reference_segmentation.numpy()
-            if ref_seg_array is None or ref_seg_array.size == 0:
-                logger.error("Could not load reference segmentation data")
-                return None
-
-            inside_mask = _check_points_in_segmentation(
-                pick_positions,
-                ref_seg_array,
-                reference_segmentation.voxel_size,
-            )
-
-        if not np.any(inside_mask):
-            logger.warning("No picks found inside reference volume")
-            return None
-
-        # Filter picks to include only those inside
-        included_points = points[inside_mask]
-        included_transforms = transforms[inside_mask] if transforms is not None else None
-
-        # Create output picks
-        output_picks = run.new_picks(object_name, session_id, user_id, exist_ok=True)
-        output_picks.from_numpy(positions=included_points, transforms=included_transforms)
-        output_picks.store()
-
-        stats = {"points_created": len(included_points)}
-        logger.info(f"Included {stats['points_created']} picks inside reference volume")
-        return output_picks, stats
-
+        return _filter_picks_by_reference(
+            picks,
+            run,
+            object_name,
+            session_id,
+            user_id,
+            reference_mesh,
+            reference_segmentation,
+            keep_inside=True,
+        )
     except Exception as e:
         logger.error(f"Error filtering picks by inclusion: {e}")
         return None
@@ -188,10 +196,15 @@ def picks_exclusion_by_mesh(
     """
     Filter picks to exclude those inside a reference mesh or segmentation.
 
+    A pick is tested at its particle centre (``location`` plus the transform's shift). Kept picks are written
+    unchanged and in their original order: transform, ``instance_id`` and ``score`` survive. If every pick is
+    inside, an empty pick set is written.
+
     Args:
         picks: CopickPicks to filter
         reference_mesh: Reference CopickMesh (either this or reference_segmentation must be provided)
-        reference_segmentation: Reference CopickSegmentation
+        reference_segmentation: Reference CopickSegmentation (binary, multilabel, instance, or the label channel
+            of a panoptic segmentation; any non-zero voxel is inside)
         run: CopickRun object
         object_name: Name for the output picks
         session_id: Session ID for the output picks
@@ -203,66 +216,16 @@ def picks_exclusion_by_mesh(
         Stats dict contains 'points_created'.
     """
     try:
-        if reference_mesh is None and reference_segmentation is None:
-            raise ValueError("Either reference_mesh or reference_segmentation must be provided")
-
-        # Load pick data
-        points, transforms = picks.numpy()
-        if points is None or len(points) == 0:
-            logger.error("Could not load pick data")
-            return None
-
-        pick_positions = points[:, :3]  # Use only x, y, z coordinates
-
-        # Check which points are inside the reference
-        if reference_mesh is not None:
-            ref_mesh = reference_mesh.mesh
-            if ref_mesh is None:
-                logger.error("Could not load reference mesh data")
-                return None
-
-            if isinstance(ref_mesh, tm.Scene):
-                if len(ref_mesh.geometry) == 0:
-                    logger.error("Reference mesh is empty")
-                    return None
-                ref_mesh = tm.util.concatenate(list(ref_mesh.geometry.values()))
-
-            inside_mask = _check_points_in_mesh(pick_positions, ref_mesh)
-
-        else:  # reference_segmentation is not None
-            ref_seg_array = reference_segmentation.numpy()
-            if ref_seg_array is None or ref_seg_array.size == 0:
-                logger.error("Could not load reference segmentation data")
-                return None
-
-            inside_mask = _check_points_in_segmentation(
-                pick_positions,
-                ref_seg_array,
-                reference_segmentation.voxel_size,
-            )
-
-        # Invert mask to exclude points inside
-        outside_mask = ~inside_mask
-
-        if not np.any(outside_mask):
-            logger.warning("No picks found outside reference volume")
-            return None
-
-        # Filter picks to exclude those inside
-        excluded_points = points[outside_mask]
-        excluded_transforms = transforms[outside_mask] if transforms is not None else None
-
-        # Create output picks
-        output_picks = run.new_picks(object_name, session_id, user_id, exist_ok=True)
-        output_picks.from_numpy(positions=excluded_points, transforms=excluded_transforms)
-        output_picks.store()
-
-        stats = {"points_created": len(excluded_points)}
-        logger.info(
-            f"Excluded {len(points) - stats['points_created']} picks inside reference volume, kept {stats['points_created']} picks",
+        return _filter_picks_by_reference(
+            picks,
+            run,
+            object_name,
+            session_id,
+            user_id,
+            reference_mesh,
+            reference_segmentation,
+            keep_inside=False,
         )
-        return output_picks, stats
-
     except Exception as e:
         logger.error(f"Error filtering picks by exclusion: {e}")
         return None
diff --git a/src/copick_utils/util/picks.py b/src/copick_utils/util/picks.py
new file mode 100644
index 0000000..825022f
--- /dev/null
+++ b/src/copick_utils/util/picks.py
@@ -0,0 +1,106 @@
+"""Helpers that keep a pick's identity intact while tools filter, paint or measure picks.
+
+A pick's particle centre is its ``location`` plus the translation of its ``transformation`` (copick's geometry
+documentation, "Shifts"). Tools that place or test particles use that centre, not ``location`` alone.
+
+Tools that keep a subset of picks copy the surviving points themselves instead of rebuilding them from arrays, so each
+point keeps its transform (rotation and shift), ``instance_id`` and ``score``, and the points keep their order. For
+filament objects the order is the order along the filament, and ``instance_id`` is the filament.
+"""
+
+from typing import TYPE_CHECKING, Sequence
+
+import numpy as np
+
+if TYPE_CHECKING:
+    from copick.models import CopickPicks, CopickPoint, CopickRun, CopickSegmentation
+
+
+def point_centres(points: Sequence["CopickPoint"]) -> np.ndarray:
+    """The particle centres of picks, in Angstrom.
+
+    Args:
+        points: Copick points.
+
+    Returns:
+        (N, 3) array of ``location + transformation[:3, 3]`` in (x, y, z) order.
+    """
+    if len(points) == 0:
+        return np.zeros((0, 3), dtype=float)
+    locations = np.array([[p.location.x, p.location.y, p.location.z] for p in points], dtype=float)
+    shifts = np.array([np.asarray(p.transformation, dtype=float)[:3, 3] for p in points], dtype=float)
+    return locations + shifts
+
+
+def pick_centres(picks: "CopickPicks") -> np.ndarray:
+    """The particle centres of a pick set, in Angstrom.
+
+    Args:
+        picks: Copick pick set.
+
+    Returns:
+        (N, 3) array of ``location + transformation[:3, 3]`` in (x, y, z) order.
+    """
+    return point_centres(picks.points or [])
+
+
+def store_pick_subset(
+    points: Sequence["CopickPoint"],
+    keep: np.ndarray,
+    run: "CopickRun",
+    object_name: str,
+    session_id: str,
+    user_id: str,
+) -> "CopickPicks":
+    """Store the kept points, unchanged and in their original order, as a pick set.
+
+    An empty selection stores an empty pick set, so an earlier result in the same place does not survive a run that
+    kept nothing.
+
+    Args:
+        points: The input points.
+        keep: Boolean mask over ``points``.
+        run: Run to write into.
+        object_name: Object name of the output pick set.
+        session_id: Session ID of the output pick set.
+        user_id: User ID of the output pick set.
+
+    Returns:
+        The stored pick set.
+    """
+    keep = np.asarray(keep, dtype=bool)
+    if keep.shape != (len(points),):
+        raise ValueError(f"The selection has shape {keep.shape}, but there are {len(points)} points.")
+    kept = [p.model_copy(deep=True) for p, k in zip(points, keep) if k]
+    output = run.new_picks(object_name, session_id, user_id, exist_ok=True)
+    output.points = kept
+    output.store()
+    return output
+
+
+def segmentation_label_volume(segmentation: "CopickSegmentation") -> np.ndarray:
+    """The volume a tool tests positions against: the label channel of a panoptic segmentation, otherwise the
+    segmentation itself (binary, multilabel or instance; non-zero means inside).
+
+    Args:
+        segmentation: Copick segmentation.
+
+    Returns:
+        (Z, Y, X) array.
+    """
+    if getattr(segmentation, "is_panoptic", False):
+        return segmentation.numpy(channel="label")
+    return segmentation.numpy()
+
+
+def voxel_indices(positions: np.ndarray, voxel_spacing: float) -> np.ndarray:
+    """The (z, y, x) indices of the voxels holding positions, under copick's mapping of voxel ``i`` to ``i * s``.
+
+    Args:
+        positions: (N, 3) positions in Angstrom, (x, y, z) order.
+        voxel_spacing: Voxel spacing in Angstrom.
+
+    Returns:
+        (N, 3) integer array in (z, y, x) order.
+    """
+    return np.round(np.asarray(positions, dtype=float) / float(voxel_spacing)).astype(int)[:, ::-1]
diff --git a/tests/conftest.py b/tests/conftest.py
new file mode 100644
index 0000000..8faadb1
--- /dev/null
+++ b/tests/conftest.py
@@ -0,0 +1,43 @@
+"""Shared fixtures: a small filesystem copick project in a temporary directory."""
+
+import json
+
+import copick
+import pytest
+
+OBJECTS = [
+    {
+        "name": "microtubule",
+        "is_particle": True,
+        "label": 1,
+        "radius": 120.0,
+        "metadata": {"copick": {"filament": {"polar": True}}},
+    },
+    {"name": "ribosome", "is_particle": True, "label": 2, "radius": 150.0},
+    {"name": "membrane", "is_particle": False, "label": 3},
+]
+
+
+@pytest.fixture
+def copick_root(tmp_path):
+    """A copick root with one run (``run1``), a filament object, a particle object and a segmentation object."""
+    config = {
+        "config_type": "filesystem",
+        "name": "test",
+        "description": "copick-utils test project",
+        "version": "1.0.0",
+        "pickable_objects": OBJECTS,
+        "overlay_root": f"local://{tmp_path / 'overlay'}",
+        "overlay_fs_args": {"auto_mkdir": True},
+    }
+    path = tmp_path / "copick_config.json"
+    path.write_text(json.dumps(config))
+    root = copick.from_file(str(path))
+    root.new_run("run1")
+    return root
+
+
+@pytest.fixture
+def run(copick_root):
+    """The fixture project's single run."""
+    return copick_root.get_run("run1")
diff --git a/tests/test_pick_tools.py b/tests/test_pick_tools.py
new file mode 100644
index 0000000..7a545db
--- /dev/null
+++ b/tests/test_pick_tools.py
@@ -0,0 +1,162 @@
+"""Pick tools keep each pick's identity and place or test it at its particle centre (location plus shift)."""
+
+import numpy as np
+import pytest
+from copick_utils.converters.segmentation_from_picks import from_picks
+from copick_utils.io.readers import coordinates
+from copick_utils.logical.distance_operations import limit_picks_by_distance
+from copick_utils.logical.point_operations import picks_exclusion_by_mesh, picks_inclusion_by_mesh
+from scipy.spatial.transform import Rotation
+
+VS = 10.0
+
+
+def _transforms(n, shifts=None, seed=0):
+    rng = np.random.default_rng(seed)
+    t = np.tile(np.eye(4), (n, 1, 1))
+    t[:, :3, :3] = Rotation.random(n, random_state=rng).as_matrix()
+    if shifts is not None:
+        t[:, :3, 3] = shifts
+    return t
+
+
+def _picks(run, positions, transforms, ids, scores, session="in"):
+    picks = run.new_picks("microtubule", session, "test", exist_ok=True)
+    picks.from_numpy(np.asarray(positions, float), transforms, instance_ids=np.asarray(ids), scores=np.asarray(scores))
+    return picks
+
+
+def _segmentation(run, array, name="membrane", session="ref", **kwargs):
+    run.new_voxel_spacing(VS, exist_ok=True)
+    seg = run.new_segmentation(VS, name, session, user_id="test", exist_ok=True, **kwargs)
+    seg.from_numpy(array)
+    return seg
+
+
+def _state(points):
+    return [
+        (
+            (p.location.x, p.location.y, p.location.z),
+            np.round(np.asarray(p.transformation), 9).tolist(),
+            p.instance_id,
+            p.score,
+        )
+        for p in points
+    ]
+
+
+@pytest.fixture
+def box_reference(run):
+    """A segmentation whose voxels x, y, z in [0, 10) are inside (shape z=12, y=14, x=16, not cubic)."""
+    array = np.zeros((12, 14, 16), dtype=np.uint8)
+    array[:10, :10, :10] = 1
+    return _segmentation(run, array)
+
+
+def test_picksin_keeps_identity_order_and_tests_centres(run, box_reference):
+    # Five picks along a filament; the third has its location outside the box but its centre inside, the fourth
+    # the opposite. Positions are in Angstrom.
+    positions = np.array([[10, 10, 10], [30, 30, 30], [120, 50, 50], [50, 50, 50], [80, 80, 80]], float)
+    shifts = np.zeros((5, 3))
+    shifts[2] = [-80, 0, 0]  # centre (40, 50, 50): inside
+    shifts[3] = [100, 0, 0]  # centre (150, 50, 50): outside
+    transforms = _transforms(5, shifts)
+    picks = _picks(run, positions, transforms, ids=[1, 1, 1, 1, 2], scores=[0.9, 0.8, 0.7, 0.6, 0.5])
+
+    out, stats = picks_inclusion_by_mesh(
+        picks,
+        run,
+        "microtubule",
+        "kept",
+        "test",
+        reference_segmentation=box_reference,
+    )
+
+    expected = [p for i, p in enumerate(picks.points) if i in (0, 1, 2, 4)]
+    assert stats["points_created"] == 4
+    assert _state(out.points) == _state(expected)
+
+    out, stats = picks_exclusion_by_mesh(
+        picks,
+        run,
+        "microtubule",
+        "dropped",
+        "test",
+        reference_segmentation=box_reference,
+    )
+    assert _state(out.points) == _state([picks.points[3]])
+
+
+def test_empty_result_replaces_an_earlier_output(run, box_reference):
+    picks = _picks(run, [[10, 10, 10]], _transforms(1), ids=[3], scores=[1.0])
+    picks_inclusion_by_mesh(picks, run, "microtubule", "result", "test", reference_segmentation=box_reference)
+    assert len(run.get_picks("microtubule", "test", "result")[0].points) == 1
+
+    # Every pick is now inside, so nothing survives exclusion: the earlier result must not be left behind.
+    out, stats = picks_exclusion_by_mesh(
+        picks,
+        run,
+        "microtubule",
+        "result",
+        "test",
+        reference_segmentation=box_reference,
+    )
+    assert stats["points_created"] == 0
+    run.refresh_picks()
+    stored = run.get_picks("microtubule", "test", "result")
+    assert len(stored) == 1 and len(stored[0].points) == 0
+
+
+def test_instance_reference_counts_every_instance_as_inside(run):
+    array = np.zeros((12, 14, 16), dtype=np.uint16)
+    array[2:4, 2:4, 2:4] = 1
+    array[6:8, 6:8, 6:8] = 7
+    instances = _segmentation(run, array, name="microtubule", session="inst", is_instance=True)
+    picks = _picks(run, [[30, 30, 30], [70, 70, 70], [100, 100, 100]], _transforms(3), ids=[1, 2, 3], scores=[1, 1, 1])
+    out, stats = picks_inclusion_by_mesh(picks, run, "microtubule", "kept", "test", reference_segmentation=instances)
+    assert [p.instance_id for p in out.points] == [1, 2]
+
+
+def test_clippicks_reference_segmentation_is_indexed_zyx(run):
+    # Non-cubic reference (z=4, y=20, x=60) with one voxel set at x=50, y=10, z=2. With x/z swapped, the pick near
+    # it would index out of bounds and be dropped.
+    array = np.zeros((4, 20, 60), dtype=np.uint8)
+    array[2, 10, 50] = 1
+    reference = _segmentation(run, array, session="clip")
+    positions = [[500, 100, 20], [100, 100, 20]]  # 0 and 400 Angstrom from the voxel
+    shifts = np.array([[0, 0, 0], [380, 0, 0]], float)  # the second's centre is 20 Angstrom away
+    picks = _picks(run, positions, _transforms(2, shifts), ids=[4, 5], scores=[0.3, 0.4])
+
+    out, stats = limit_picks_by_distance(
+        picks,
+        run,
+        "microtubule",
+        "clipped",
+        "test",
+        reference_segmentation=reference,
+        max_distance=50.0,
+    )
+    assert stats["points_created"] == 2
+    assert _state(out.points) == _state(picks.points)
+
+    out, stats = limit_picks_by_distance(
+        picks,
+        run,
+        "microtubule",
+        "far",
+        "test",
+        reference_segmentation=reference,
+        max_distance=50.0,
+        invert=True,
+    )
+    assert stats["points_created"] == 0 and len(out.points) == 0
+
+
+def test_painting_and_coordinates_use_centres(run):
+    picks = _picks(run, [[30, 40, 50]], _transforms(1, np.array([[20, 0, -10]])), ids=[1], scores=[1.0])
+    volume = from_picks(picks, np.zeros((10, 10, 10), np.uint8), radius=5.0, label_value=4, voxel_spacing=VS)
+    assert volume[4, 4, 5] == 4  # centre (50, 40, 40) Angstrom -> voxel z=4, y=4, x=5
+    assert volume[5, 4, 3] == 0  # location (30, 40, 50) is not painted
+
+    zyx = coordinates(run, "microtubule", "test", "in", voxel_size=VS)
+    np.testing.assert_allclose(zyx, [[4.0, 4.0, 5.0]])

From ef1da6fe3d756c27fb9126a1653a815107e3e2f6 Mon Sep 17 00:00:00 2001
From: uermel <utz.ermel@czii.org>
Date: Tue, 6 Oct 2026 00:30:57 -0700
Subject: [PATCH 3/3] chore: require copick 1.28.0

copick 1.28.0 ships the filament stack this builds on (copick#479, #480,
#492, #494), so the build-against pin to the unreleased branch goes.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 pyproject.toml | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/pyproject.toml b/pyproject.toml
index c2c13c8..a7c868c 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -18,9 +18,7 @@ classifiers = [
 ]
 dynamic = ["version"]
 dependencies = [
-    # TODO:remove before release -- build against the unreleased copick filament stack (copick#479/#480/#492/#494);
-    # restore "copick>=<release with those PRs>" when it is published.
-    "copick @ git+https://github.com/copick/copick@uermel/panoptic-segmentations",
+    "copick>=1.28.0",
     "click",
     "click-option-group",
     "numpy",