diff --git a/src/copick_utils/cli/combine_labels.py b/src/copick_utils/cli/combine_labels.py index baac9f7..638d916 100644 --- a/src/copick_utils/cli/combine_labels.py +++ b/src/copick_utils/cli/combine_labels.py @@ -20,6 +20,15 @@ @add_run_names_option @optgroup.group("\nInput Options", help="Options related to the input segmentations.") @add_input_option("segmentation") +@optgroup.option( + "--instances", + "instances_uri", + type=str, + default=None, + help="Instance segmentations to add to a panoptic output (-o ...?panoptic=true), e.g. " + '"microtubule:trace/1@10.0" (patterns allowed; ?instance=true is implied). Each one adds its object\'s ' + "label and its instance IDs.", +) @add_workers_option @optgroup.group("\nOutput Options", help="Options related to the output segmentation.") @add_output_option("segmentation", default_tool="combine") @@ -28,6 +37,7 @@ def combine( config, run_names, input_uri, + instances_uri, workers, output_uri, debug, @@ -43,11 +53,17 @@ def combine( When multiple inputs overlap, the lowest label value wins. This resolution is deterministic and reproducible, and overlapping voxels are logged as warnings. + With an output URI ending in `?panoptic=true`, the result is a panoptic segmentation + instead: the inputs fill its label channel (regions without instances, such as + membranes), and the instance segmentations given with `--instances` add their objects' + labels and their instance IDs. An instance ID is kept where its object's label won. + URI Format: \b Segmentations: name:user_id/session_id@voxel_spacing Use glob/regex patterns to match multiple segmentations per run. + Panoptic output: append ?panoptic=true to the output URI Examples: @@ -63,6 +79,11 @@ def combine( # Combine for specific runs copick process combine -r run1 -r run2 -i "*:user1/session@10.0" -o "labels:combine/0" + \b + # Membranes as regions plus microtubule instances, as one panoptic segmentation + copick process combine -i "membrane:data-portal/*@10.0" --instances "microtubule:trace/1@10.0" \\ + -o "cell:combine/0@10.0?panoptic=true" + See Also: \b @@ -100,6 +121,7 @@ def combine( config=task_config, run_names=run_names_list, workers=workers, + instances_uri=instances_uri, ) successful = sum(1 for result in results.values() if result and result.get("processed", 0) > 0) diff --git a/src/copick_utils/cli/filter_components.py b/src/copick_utils/cli/filter_components.py index 9e74a05..5135988 100644 --- a/src/copick_utils/cli/filter_components.py +++ b/src/copick_utils/cli/filter_components.py @@ -53,6 +53,19 @@ help="Keep only the N largest connected components by voxel count (e.g. 1 = keep only the " "single largest). Applied in addition to any --min-size/--max-size filter.", ) +@optgroup.option( + "--min-skeleton-length", + type=float, + default=None, + help="Minimum skeleton length of a component (or instance) to keep, in the unit set by --length-unit " + "(optional). Drops specks and short blobs of filamentous structures; see seg-stats --skeleton.", +) +@optgroup.option( + "--length-unit", + type=click.Choice(["angstrom", "voxel"]), + default="angstrom", + help="Unit for --min-skeleton-length: 'angstrom' for Å, 'voxel' for voxels.", +) @add_workers_option @optgroup.group("\nOutput Options", help="Options related to output segmentations.") @add_output_option("segmentation", default_tool="filter-components") @@ -66,6 +79,8 @@ def filter_components( max_size, size_unit, keep_largest, + min_skeleton_length, + length_unit, workers, output_uri, debug, @@ -82,10 +97,17 @@ def filter_components( top of any --min-size/--max-size limits. Run `copick process seg-stats` first to inspect component sizes and choose sensible thresholds. + A multilabel segmentation is filtered label by label and keeps its labels. An instance + segmentation (`?instance=true`) is filtered instance by instance: whole instances are + dropped, and the IDs of the kept ones do not change, so picks and filaments with the same + IDs still match. `--min-skeleton-length` also drops components (or instances) whose + skeleton is shorter than the given length. + URI Format: \b Segmentations: name:user_id/session_id@voxel_spacing + Instance segmentations: append ?instance=true Examples: @@ -107,6 +129,11 @@ def filter_components( # Keep only the single largest connected component copick process filter-components -i "membrane:user1/auto-001@10.0" -o "membrane_main" --keep-largest 1 + \b + # Drop microtubule instances shorter than 100 nm, keeping the IDs of the others + copick process filter-components -i "microtubule:trace/1@10.0?instance=true" -o "microtubule:trace/1-long" \\ + --min-skeleton-length 1000 + See Also: \b @@ -134,6 +161,13 @@ def filter_components( if max_size is not None: max_size = max_size * voxel_volume + if length_unit == "voxel" and min_skeleton_length is not None: + input_params = parse_copick_uri(input_uri, "segmentation") + vs_raw = input_params.get("voxel_spacing") + if vs_raw is None or vs_raw == "*": + raise click.BadParameter("--length-unit voxel requires voxel spacing in the input URI (e.g., @10.0)") + min_skeleton_length = min_skeleton_length * float(vs_raw) + # Create config from URIs with smart defaults try: task_config = create_simple_config( @@ -159,6 +193,8 @@ def filter_components( if keep_largest < 1: raise click.BadParameter("--keep-largest must be >= 1") logger.info(f"Keeping only the {keep_largest} largest component(s)") + if min_skeleton_length is not None: + logger.info(f"Minimum skeleton length: {min_skeleton_length} Å") # Parallel discovery and processing results = filter_components_lazy_batch( @@ -170,6 +206,7 @@ def filter_components( min_size=min_size, max_size=max_size, keep_largest=keep_largest, + min_skeleton_length=min_skeleton_length, ) successful = sum(1 for result in results.values() if result and result.get("processed", 0) > 0) diff --git a/src/copick_utils/cli/seg_stats.py b/src/copick_utils/cli/seg_stats.py index 7bc0860..addbb31 100644 --- a/src/copick_utils/cli/seg_stats.py +++ b/src/copick_utils/cli/seg_stats.py @@ -32,6 +32,12 @@ help="Include background (label=0) connected component analysis in CSV output. " "Background components are never included in plots.", ) +@optgroup.option( + "--skeleton/--no-skeleton", + default=False, + help="Also skeletonize each component (or instance) and report its skeleton length, label radius " + "(median distance from the skeleton to the background) and branch, junction and endpoint counts, in angstroms.", +) @add_workers_option @optgroup.group("\nOutput Options", help="Options related to the output.") @optgroup.option( @@ -55,6 +61,7 @@ def seg_stats( input_uri, connectivity, include_background, + skeleton, workers, output_format, output_path, @@ -68,11 +75,18 @@ def seg_stats( Output can be a CSV file or a histogram plot, and you can choose the voxel connectivity (face, face-edge, or all) used to define components. + An instance segmentation (`?instance=true`) gets one row per instance ID, and a panoptic + segmentation (`?panoptic=true`) one row per (label, instance) segment, plus the connected + components of regions without instances. Without a type flag, the URI matches binary and + multilabel segmentations only. With `--skeleton`, each row also reports the component's + skeleton length and label radius, for example to choose length thresholds for filaments. + URI Format: \b Segmentations: name:user_id/session_id@voxel_spacing Voxel spacing is optional — omit to match all voxel spacings. + Instance or panoptic segmentations: append ?instance=true or ?panoptic=true Examples: @@ -96,6 +110,10 @@ def seg_stats( # Analyze specific runs and export as PNG copick process seg-stats -i "membrane:user1/auto-001@10.0" -f plot -op ./stats.png -r run1 -r run2 + \b + # Per-instance volumes and skeleton lengths of a microtubule instance segmentation + copick process seg-stats -i "microtubule:trace/1@10.0?instance=true" --skeleton -f csv -op ./mt.csv + See Also: \b @@ -121,6 +139,7 @@ def seg_stats( include_background=include_background, run_names=run_names_list, workers=workers, + skeleton=skeleton, ) # Summarize diff --git a/src/copick_utils/cli/separate_components.py b/src/copick_utils/cli/separate_components.py index 3c9da7a..965faf9 100644 --- a/src/copick_utils/cli/separate_components.py +++ b/src/copick_utils/cli/separate_components.py @@ -7,6 +7,7 @@ from copick_utils.cli.util import add_input_option, add_output_option, add_workers_option from copick_utils.util.config_models import create_simple_config +from copick_utils.util.segmentations import type_from_uri @click.command( @@ -62,10 +63,17 @@ def separate_components( For multilabel segmentations the analysis is performed on each label separately. Output segmentations use the `{instance_id}` placeholder for auto-numbering (e.g. `inst-0`, `inst-1`). + Alternatively, write all components into one segmentation: an output URI with + `?instance=true` writes an instance segmentation of a binary input (each component is an + instance, numbered 1, 2, ... by size, largest first), and `?panoptic=true` writes a + panoptic segmentation that keeps every label of a multilabel input and numbers the + components of each label. These outputs need no `{instance_id}` placeholder. + URI Format: \b Segmentations: name:user_id/session_id@voxel_spacing + Instance or panoptic output: append ?instance=true or ?panoptic=true Examples: @@ -83,6 +91,16 @@ def separate_components( copick process separate-components -i "membrane:user1/manual-001@10.0" \\ -o "membrane:components/comp-{instance_id}@10.0" + \b + # All components as one instance segmentation (IDs 1..K by size) + copick process separate-components -i "microtubule:easymode/job006@10.0" --binary \\ + -o "microtubule:components/job006@10.0?instance=true" + + \b + # Components of every label of a multilabel segmentation as one panoptic segmentation + copick process separate-components -i "labels:user1/auto@10.0?multilabel=true" \\ + -o "cell:components/auto@10.0?panoptic=true" + See Also: \b @@ -97,14 +115,16 @@ def separate_components( root = copick.from_file(config) run_names_list = list(run_names) if run_names else None - # Create config from URIs with smart defaults (individual_outputs for {instance_id}) + # One segmentation per component (with {instance_id}), unless the output URI names an instance or panoptic + # segmentation, which holds all components in one store. try: + output_type = type_from_uri(output_uri) task_config = create_simple_config( input_uri=input_uri, input_type="segmentation", output_uri=output_uri, output_type="segmentation", - individual_outputs=True, + individual_outputs=output_type not in ("instance", "panoptic"), command_name="components", ) except ValueError as e: diff --git a/src/copick_utils/cli/split_labels.py b/src/copick_utils/cli/split_labels.py index f770c31..3b8c84d 100644 --- a/src/copick_utils/cli/split_labels.py +++ b/src/copick_utils/cli/split_labels.py @@ -60,10 +60,15 @@ def split( values do not match the config object labels; only the listed values are then split. The input URI must name an exact segmentation (no wildcards) and include a voxel spacing. + A panoptic segmentation (input URI ending in `?panoptic=true`) is split per object: the + voxels with an instance ID become an instance segmentation of that object (keeping the + IDs), and its voxels without one become a binary segmentation. + URI Format: \b Segmentations: name:user_id/session_id@voxel_spacing + Panoptic input: append ?panoptic=true Label-to-Object Mapping: @@ -92,6 +97,10 @@ def split( # Split only specific label values with an explicit name:value map copick process split -i "predictions:model/run-001@10.0" --labels "sample:1,vacuum:2" + \b + # Split a panoptic segmentation into instance and binary segmentations per object + copick process split -i "cell:combine/0@10.0?panoptic=true" + See Also: \b @@ -164,6 +173,7 @@ def split( run_names=run_names_list, workers=workers, labels=labels_map, + panoptic=bool(input_params.get("panoptic")), ) # Aggregate results diff --git a/src/copick_utils/cli/util.py b/src/copick_utils/cli/util.py index aa0b262..d611d1e 100644 --- a/src/copick_utils/cli/util.py +++ b/src/copick_utils/cli/util.py @@ -1,6 +1,6 @@ """CLI utilities for copick-utils commands.""" -from typing import Callable +from typing import Callable, Optional import click from click_option_group import optgroup @@ -340,12 +340,13 @@ def add_input_option(object_type: str, func: click.Command = None, required: boo Add --input/-i option for URI-based input selection. Supports copick URI format with pattern matching: - - Picks/Meshes: object_name:user_id/session_id - - Segmentations: name:user_id/session_id@voxel_spacing?multilabel=true + - Picks/Meshes/Filaments: object_name:user_id/session_id + - Segmentations: name:user_id/session_id@voxel_spacing, with ?multilabel=true, ?instance=true or + ?panoptic=true to select a type (without a flag, binary and multilabel segmentations match) - Tomograms: tomo_type@voxel_spacing Args: - object_type (str): Type of object ('picks', 'mesh', 'segmentation', 'tomogram'). + object_type (str): Type of object ('picks', 'mesh', 'segmentation', 'filaments', 'tomogram'). func (click.Command, optional): The Click command to which the option will be added. required (bool): Whether the option is required. Default is True. Pass False for commands that keep deprecated fallback flags (e.g. --tomo-alg/--voxel-size). @@ -361,12 +362,15 @@ def add_input_option_decorator(_func: click.Command) -> click.Command: "picks": "object_name:user_id/session_id", "mesh": "object_name:user_id/session_id", "segmentation": "name:user_id/session_id@voxel_spacing", + "filaments": "object_name:user_id/session_id", "tomogram": "tomo_type@voxel_spacing", } help_text = ( f"Input {object_type} URI (format: {format_examples.get(object_type, 'URI')}). Supports glob patterns." ) + if object_type == "segmentation": + help_text += " Append ?instance=true or ?panoptic=true to read those segmentation types." opt = optgroup.option( "--input", @@ -389,10 +393,17 @@ def add_output_option( func: click.Command = None, default_tool: str = None, required: bool = True, + flag: str = "--output", + short_flag: Optional[str] = "-o", + param_name: str = "output_uri", + description: Optional[str] = None, ) -> Callable: """ Add --output/-o option for URI-based output specification with smart defaults. + A command with a second output (for example a tracer that writes filaments and an instance segmentation) adds it + with its own ``flag``, ``short_flag`` and ``param_name``; ``required=False`` makes it optional. + Supports copick URI format with smart defaults and pattern matching: - Full format: object_name:user_id/session_id or name:user_id/session_id@voxel_spacing - Smart defaults: Omit components to inherit from input (e.g., just "membrane") @@ -405,9 +416,14 @@ def add_output_option( - voxel_spacing omitted → inherits from input (segmentation only) Args: - object_type (str): Type of object ('picks', 'mesh', 'segmentation'). + object_type (str): Type of object ('picks', 'mesh', 'segmentation', 'filaments', 'tomogram'). func (click.Command, optional): The Click command to which the option will be added. default_tool (str, optional): Default user_id if not specified in URI (deprecated, auto-detected). + required (bool): Whether the option is required. + flag (str): Long option name. Default "--output". + short_flag (str, optional): Short option name, or None. Default "-o". + param_name (str): Name of the command parameter. Default "output_uri". + description (str, optional): Help text to use instead of the generated one. Returns: Callable: The Click command with the output option added. @@ -420,6 +436,7 @@ def add_output_option_decorator(_func: click.Command) -> click.Command: "picks": '"ribosome", "ribosome/my-session", or "/my-session"', "mesh": '"membrane", "membrane/my-session", or "/my-session"', "segmentation": '"membrane", "membrane/my-session", or "/my-session"', + "filaments": '"microtubule", "microtubule/my-session", or "/my-session"', "tomogram": '"wbp@20.0" or "denoised@10.0"', } @@ -432,11 +449,14 @@ def add_output_option_decorator(_func: click.Command) -> click.Command: f"Supports smart defaults (e.g., {shorthand_examples.get(object_type, 'shorthand')}). " f"Full format: object_name:user_id/session_id{voxel_suffix}." ) + if object_type == "segmentation": + help_text += " Append ?instance=true or ?panoptic=true to write those segmentation types." + if description is not None: + help_text = description + names = [flag] + ([short_flag] if short_flag else []) + [param_name] opt = optgroup.option( - "--output", - "-o", - "output_uri", + *names, type=CopickURI(object_type, "output"), required=required, help=help_text, @@ -606,7 +626,9 @@ def add_reference_seg_option_decorator(_func: click.Command) -> click.Command: "ref_seg_uri", type=CopickURI("segmentation", "reference"), required=required, - help="Reference segmentation URI (format: name:user_id/session_id@voxel_spacing). Supports glob patterns.", + help="Reference segmentation URI (format: name:user_id/session_id@voxel_spacing). Supports glob patterns. " + "Append ?instance=true or ?panoptic=true to use an instance or panoptic segmentation (any non-zero " + "voxel of an instance segmentation, or of a panoptic one's label channel, counts).", ) return opt(_func) diff --git a/src/copick_utils/converters/lazy_converter.py b/src/copick_utils/converters/lazy_converter.py index 64439f8..116e7c0 100644 --- a/src/copick_utils/converters/lazy_converter.py +++ b/src/copick_utils/converters/lazy_converter.py @@ -6,6 +6,7 @@ from copick.util.log import get_logger from copick_utils.util.config_models import ReferenceConfig, SelectorConfig, TaskConfig +from copick_utils.util.segmentations import query_kwargs, selection_filters if TYPE_CHECKING: from copick.models import CopickRoot, CopickRun @@ -92,6 +93,16 @@ def create_reference_config( ) +def _accepts_parameter(func: Callable, name: str) -> bool: + """Whether ``func`` names the parameter ``name`` explicitly.""" + import inspect + + try: + return name in inspect.signature(func).parameters + except (TypeError, ValueError): + return False + + def _is_regex_pattern(pattern: str) -> bool: """Check if string is a regex pattern.""" regex_chars = r"[.*+?^${}()|[\]\\" @@ -127,7 +138,7 @@ def discover_tasks_for_run(run: "CopickRun", selector_config: SelectorConfig) -> # Build filter dict based on input type filters = {"pattern_type": pattern_type} - if selector_config.input_type == "picks" or selector_config.input_type == "mesh": + if selector_config.input_type in ("picks", "mesh", "filaments"): filters["object_name"] = selector_config.input_object_name filters["user_id"] = selector_config.input_user_id filters["session_id"] = selector_config.input_session_id @@ -136,6 +147,14 @@ def discover_tasks_for_run(run: "CopickRun", selector_config: SelectorConfig) -> filters["user_id"] = selector_config.input_user_id filters["session_id"] = selector_config.input_session_id filters["voxel_spacing"] = selector_config.voxel_spacing + # Without an instance or panoptic flag on the input URI, only binary and multilabel segmentations match. + filters.update( + selection_filters( + selector_config.input_multilabel, + selector_config.input_instance, + selector_config.input_panoptic, + ), + ) # Find matching input objects using copick's official resolution matching_inputs = get_copick_objects_by_type( @@ -149,12 +168,8 @@ def discover_tasks_for_run(run: "CopickRun", selector_config: SelectorConfig) -> return [] # Generate type-specific input parameter name - if selector_config.input_type == "mesh": - input_param_name = "mesh" - elif selector_config.input_type == "segmentation": - input_param_name = "segmentation" - elif selector_config.input_type == "picks": - input_param_name = "picks" + if selector_config.input_type in ("mesh", "segmentation", "picks", "filaments"): + input_param_name = selector_config.input_type else: input_param_name = "input_object" # fallback @@ -176,6 +191,8 @@ def discover_tasks_for_run(run: "CopickRun", selector_config: SelectorConfig) -> "output_type": selector_config.output_type, "segmentation_name": selector_config.segmentation_name, "voxel_spacing": selector_config.voxel_spacing, + # The segmentation type the output URI names (None: the converter's own choice, as before) + "output_segmentation_type": selector_config.output_segmentation_type, } # Add session ID template for individual outputs @@ -238,6 +255,7 @@ def add_references_to_tasks( user_id=reference_config.user_id, session_id=reference_config.session_id, voxel_size=reference_config.voxel_spacing, + **query_kwargs(reference_config.multilabel, reference_config.instance, reference_config.panoptic), ) ref_key = "reference_segmentation" alt_key = "reference_mesh" @@ -436,12 +454,14 @@ def lazy_conversion_worker( # Single input pattern that expands to N-way union discovered_tasks = discover_tasks_for_run(run, config.selector) - if len(discovered_tasks) < 2: + # A panoptic output can be built from one input plus instance segmentations given separately + minimum = 1 if config.selector.output_segmentation_type == "panoptic" else 2 + if len(discovered_tasks) < minimum: # Not enough matches for union operation return { "processed": 0, "errors": [ - f"Pattern matched {len(discovered_tasks)} segmentation(s) in {run.name}, but union requires at least 2", + f"Pattern matched {len(discovered_tasks)} segmentation(s) in {run.name}, but union requires at least {minimum}", ], } @@ -469,6 +489,7 @@ def lazy_conversion_worker( "session_id": first_task["session_id"], "voxel_spacing": first_task.get("voxel_spacing"), "is_multilabel": False, + "output_segmentation_type": first_task.get("output_segmentation_type"), }, ] @@ -499,6 +520,11 @@ def lazy_conversion_worker( ): task_params.pop(key, None) + # The output URI's segmentation type reaches only converters that ask for it by name, so converters + # that forward **kwargs to helpers are unaffected. + if not _accepts_parameter(converter_func, "output_segmentation_type"): + task_params.pop("output_segmentation_type", None) + task_params.update(converter_kwargs) result = converter_func(**task_params) @@ -519,6 +545,7 @@ def lazy_conversion_worker( or task.get("segmentation") or task.get("mesh") or task.get("picks") + or task.get("filaments") or task.get("segmentation1") or task.get("mesh1") or task.get("picks1") diff --git a/src/copick_utils/io/readers.py b/src/copick_utils/io/readers.py index 27df162..5d62012 100644 --- a/src/copick_utils/io/readers.py +++ b/src/copick_utils/io/readers.py @@ -65,7 +65,17 @@ def tomogram(run, voxel_size: float = 10, algorithm: str = "wbp", raise_error: b return None -def segmentation(run, voxel_spacing: float, name: str, user_id=None, session_id=None, raise_error=False, verbose=True): +def segmentation( + run, + voxel_spacing: float, + name: str, + user_id=None, + session_id=None, + raise_error=False, + verbose=True, + is_instance: bool = False, + is_panoptic: bool = False, +): """Read a segmentation from a copick run as a NumPy array. Resolves the segmentation matching `name` (optionally filtered by `user_id` @@ -82,10 +92,14 @@ def segmentation(run, voxel_spacing: float, name: str, user_id=None, session_id= returning `None`. verbose: Print a diagnostic message listing what is available when no segmentation matches. + is_instance: Read an instance segmentation (voxel value = instance ID). By + default only binary and multilabel segmentations match, even when an + instance or panoptic segmentation shares the name, user and session. + is_panoptic: Read a panoptic segmentation (shape (2, Z, Y, X)). Returns: - The segmentation as a NumPy array of shape (Z, Y, X), or `None` if none - matches and `raise_error` is False. + The segmentation as a NumPy array of shape (Z, Y, X) ((2, Z, Y, X) for a + panoptic one), or `None` if none matches and `raise_error` is False. """ # Construct the Target URI @@ -98,7 +112,15 @@ def segmentation(run, voxel_spacing: float, name: str, user_id=None, session_id= # Try to resolve the segmentation using the Copick URI try: + if is_instance or is_panoptic: + uri += "?instance=true" if is_instance else "?panoptic=true" segs = resolve_copick_objects(uri, run.root, "segmentation", run_name=run.name) + segs = [ + s + for s in segs + if bool(getattr(s, "is_instance", False)) == is_instance + and bool(getattr(s, "is_panoptic", False)) == is_panoptic + ] return segs[0].numpy() except Exception as err: # Force the voxel spacing to be a float diff --git a/src/copick_utils/io/writers.py b/src/copick_utils/io/writers.py index 28bd80d..7138fa2 100644 --- a/src/copick_utils/io/writers.py +++ b/src/copick_utils/io/writers.py @@ -36,6 +36,7 @@ def segmentation( session_id="0", voxel_size=10, multilabel=True, + instance=False, ): """Write a segmentation into a copick run as OME-Zarr. @@ -51,10 +52,20 @@ def segmentation( session_id: Session id for the segmentation. Defaults to `0`. voxel_size: Voxel spacing in angstroms. multilabel: Whether the segmentation holds multiple labels. + instance: Write an instance segmentation (voxel value = instance ID; `name` + is the object). Its values are stored without loss (uint16 or wider) + instead of as `uint8`. """ - # Retrieve or create a segmentation - segmentations = run.get_segmentations(name=name, user_id=user_id, session_id=session_id) + # Retrieve or create a segmentation of the requested type + multilabel = multilabel and not instance + segmentations = run.get_segmentations( + name=name, + user_id=user_id, + session_id=session_id, + is_instance=instance, + is_panoptic=False, + ) # If no segmentation exists or no segmentation at the given voxel size, create a new one if len(segmentations) == 0 or any(seg.voxel_size != voxel_size for seg in segmentations): @@ -64,10 +75,11 @@ def segmentation( session_id=session_id, is_multilabel=multilabel, user_id=user_id, + is_instance=instance, ) else: # Overwrite the current segmentation at the specified voxel size if it exists seg = next(seg for seg in segmentations if seg.voxel_size == voxel_size) # Write the segmentation data - seg.from_numpy(seg_vol, dtype=np.uint8) + seg.from_numpy(seg_vol, dtype=None if instance else np.uint8) diff --git a/src/copick_utils/process/combine_labels.py b/src/copick_utils/process/combine_labels.py index bde770d..d671c1c 100644 --- a/src/copick_utils/process/combine_labels.py +++ b/src/copick_utils/process/combine_labels.py @@ -13,12 +13,81 @@ logger = get_logger(__name__) +def _combine_panoptic( + segmentations: List["CopickSegmentation"], + run: "CopickRun", + object_name: str, + session_id: str, + user_id: str, + instances_uri: Optional[str] = None, +) -> Optional[Tuple["CopickSegmentation", Dict[str, int]]]: + """Combine binary or multilabel segmentations (regions without instances) and per-object instance + segmentations into one panoptic segmentation. Overlaps go to the lowest label; an instance ID survives only + where its object's label won.""" + from copick_utils.util.segmentations import new_segmentation, resolve_segmentations, segmentation_type + + root = run.root + layers = [] # (label value, label mask, instance IDs or None) + for seg in segmentations: + volume = seg.numpy() + if volume is None: + logger.warning(f"Could not load segmentation '{seg.name}', skipping") + continue + if segmentation_type(seg) == "multilabel": + for value in np.unique(volume[volume > 0]): + layers.append((int(value), volume == value, None)) + else: + obj = root.get_object(seg.name) + if obj is None: + logger.warning(f"No pickable object found for '{seg.name}', skipping") + continue + layers.append((int(obj.label), volume > 0, None)) + voxel_size = segmentations[0].voxel_size if segmentations else None + if instances_uri: + uri = instances_uri if "instance=true" in instances_uri else instances_uri + "?instance=true" + for seg in resolve_segmentations(uri, root, run_name=run.name): + voxel_size = voxel_size or seg.voxel_size + obj = root.get_object(seg.name) + if obj is None: + logger.warning(f"No pickable object found for instance segmentation '{seg.name}', skipping") + continue + ids = seg.numpy() + layers.append((int(obj.label), ids > 0, ids)) + if not layers: + logger.error("No segmentations to combine") + return None + + shape = layers[0][1].shape + label_channel = np.zeros(shape, dtype=np.uint32) + instance_channel = np.zeros(shape, dtype=np.uint32) + overlap_count = np.zeros(shape, dtype=np.uint8) + # Paint the highest label first, so the lowest label wins where inputs overlap + for label_value, mask, ids in sorted(layers, key=lambda layer: layer[0], reverse=True): + overlap_count += mask.astype(np.uint8) + label_channel[mask] = label_value + instance_channel[mask] = ids[mask] if ids is not None else 0 + overlapping_voxels = int(np.sum(overlap_count > 1)) + if overlapping_voxels: + logger.warning(f"Detected {overlapping_voxels} overlapping voxels across inputs. Resolved by lowest label.") + + output_seg = new_segmentation(run, voxel_size, object_name, session_id, user_id, "panoptic") + output_seg.from_numpy(np.stack([label_channel, instance_channel])) + stats = { + "labels_combined": len({layer[0] for layer in layers}), + "instances_combined": int(len(np.unique(instance_channel[instance_channel > 0]))), + "overlapping_voxels": overlapping_voxels, + } + return output_seg, stats + + def combine_labels( segmentations: List["CopickSegmentation"], run: "CopickRun", object_name: str, session_id: str, user_id: str, + output_segmentation_type: Optional[str] = None, + instances_uri: Optional[str] = None, **kwargs, ) -> Optional[Tuple["CopickSegmentation", Dict[str, int]]]: """ @@ -28,18 +97,28 @@ def combine_labels( its integer label value. Overlapping regions are resolved by lowest label priority (lowest label value wins). + When the output URI names a panoptic segmentation (``?panoptic=true``), the inputs become its label channel and + the instance segmentations selected by ``instances_uri`` add their objects' labels and instance IDs. + Args: segmentations: List of input CopickSegmentation objects (binary/single-label). run: CopickRun object. object_name: Name for the output multilabel segmentation. session_id: Session ID for the output segmentation. user_id: User ID for the output segmentation. + output_segmentation_type: The segmentation type the output URI names, if any. + instances_uri: Instance segmentations to add to a panoptic output (URI, patterns allowed). **kwargs: Additional keyword arguments from lazy converter. Returns: Tuple of (CopickSegmentation, stats dict) or None if operation failed. """ try: + if output_segmentation_type == "panoptic": + return _combine_panoptic(segmentations, run, object_name, session_id, user_id, instances_uri) + if output_segmentation_type == "instance": + logger.error("combine writes multilabel or panoptic segmentations, not instance ones") + return None if not segmentations: logger.error("No segmentations provided") return None diff --git a/src/copick_utils/process/connected_components.py b/src/copick_utils/process/connected_components.py index ef67260..1b9018f 100644 --- a/src/copick_utils/process/connected_components.py +++ b/src/copick_utils/process/connected_components.py @@ -12,6 +12,8 @@ if TYPE_CHECKING: from copick.models import CopickRun, CopickSegmentation +logger = get_logger(__name__) + def separate_connected_components_3d( volume: np.ndarray, @@ -62,30 +64,23 @@ def separate_connected_components_3d( # Label connected components labeled_volume, num_components = ndimage.label(binary_volume, structure=structure) - - print(f"Found {num_components} connected components") + logger.debug(f"Found {num_components} connected components") + + # Filter by size if specified: one lookup table instead of one pass over the volume per component + if min_size is not None and min_size > 0 and num_components > 0: + sizes = np.bincount(labeled_volume.ravel(), minlength=num_components + 1) + keep = sizes * voxel_spacing**3 >= min_size + keep[0] = False + relabel = np.zeros(num_components + 1, dtype=labeled_volume.dtype) + relabel[keep] = np.arange(1, int(keep.sum()) + 1, dtype=labeled_volume.dtype) + labeled_volume = relabel[labeled_volume] + num_components = int(keep.sum()) + logger.debug(f"After filtering by size (min={min_size} ų): {num_components} components") # Get component properties component_info = {} props = measure.regionprops(labeled_volume) - print(f"Found {len(props)} connected components") - - # Calculate voxel volume in cubic angstroms - voxel_volume = voxel_spacing**3 - - # Filter by size if specified - if min_size is not None and min_size > 0: - for prop in props: - component_volume = prop.area * voxel_volume - if component_volume < min_size: - labeled_volume[labeled_volume == prop.label] = 0 - - # Relabel after filtering - labeled_volume, num_components = ndimage.label(labeled_volume > 0, structure=structure) - props = measure.regionprops(labeled_volume) - print(f"After filtering by size (min={min_size} ų): {num_components} components") - # Store component information for _i, prop in enumerate(props, 1): component_info[prop.label] = { @@ -131,6 +126,35 @@ def print_component_stats(component_info: Dict[int, Dict[str, Any]]) -> None: print(f"{label:<8} {info['volume']:<10} {centroid_str:<25} {info['extent']:<10.3f}") +def components_as_instances( + volume: np.ndarray, + voxel_spacing: float, + connectivity: Union[int, str] = "all", + min_size: Optional[float] = None, +) -> Tuple[np.ndarray, int]: + """ + Label the connected components of a binary volume as instances, numbered 1..K by size (largest first). + + Args: + volume: 3D binary volume (any non-zero voxel is foreground). + voxel_spacing: Voxel spacing in angstroms. + connectivity: "face", "face-edge" or "all" (or 6, 18, 26). + min_size: Minimum component volume in cubic angstroms (ų) to keep (None = keep all). + + Returns: + Tuple of (instance volume, K). The volume holds 0 for background and 1..K for the components. + """ + labeled, n, _ = separate_connected_components_3d(volume, voxel_spacing, connectivity, min_size) + if n == 0: + return np.zeros(volume.shape, dtype=np.uint16), 0 + sizes = np.bincount(labeled.ravel(), minlength=n + 1)[1:] + order = np.argsort(-sizes, kind="stable") # largest first; ties keep their scan order + dtype = np.uint16 if n < 2**16 else np.uint32 + rank = np.zeros(n + 1, dtype=dtype) + rank[order + 1] = np.arange(1, n + 1, dtype=dtype) + return rank[labeled], n + + def separate_segmentation_components( segmentation: "CopickSegmentation", connectivity: Union[int, str] = "all", @@ -160,7 +184,7 @@ def separate_segmentation_components( # Handle deprecated session_id_prefix parameter if session_id_prefix is not None: session_id_template = f"{session_id_prefix}{{instance_id}}" - # Get the segmentation volume + volume = segmentation.numpy() if volume is None: raise ValueError("Could not load segmentation data") @@ -169,90 +193,117 @@ def separate_segmentation_components( voxel_size = segmentation.voxel_size name = segmentation.name - output_segmentations = [] - component_count = 0 - if multilabel: - # Process each label separately unique_labels = np.unique(volume) - unique_labels = unique_labels[unique_labels > 0] # skip background - - print(f"Processing multilabel segmentation with {len(unique_labels)} labels") - - for label_value in unique_labels: - print(f"Processing label {label_value}") - - # Extract binary volume for this label - binary_vol = volume == label_value - - # Separate connected components - labeled_vol, n_components, component_info = separate_connected_components_3d( - binary_vol, - voxel_spacing=voxel_size, - connectivity=connectivity, - min_size=min_size, - ) - - # Extract individual components - individual_components = extract_individual_components(labeled_vol) - - # Create segmentations for each component - for component_vol in individual_components: - session_id = session_id_template.replace("{instance_id}", str(component_count)) - - # Create new segmentation - output_seg = run.new_segmentation( - voxel_size=voxel_size, - name=name, - session_id=session_id, - is_multilabel=False, - user_id=output_user_id, - exist_ok=True, - ) - - # Store the component volume - output_seg.from_numpy(component_vol) - output_segmentations.append(output_seg) - component_count += 1 - + masks = [volume == label_value for label_value in unique_labels[unique_labels > 0]] + logger.info(f"Processing multilabel segmentation with {len(masks)} labels") else: - # Process as binary segmentation - print("Processing binary segmentation") + masks = [volume > 0] - # Separate connected components - labeled_vol, n_components, component_info = separate_connected_components_3d( - volume, + output_segmentations = [] + component_count = 0 + for mask in masks: + labeled_vol, n_components, _ = separate_connected_components_3d( + mask, voxel_spacing=voxel_size, connectivity=connectivity, min_size=min_size, ) - - # Extract individual components - individual_components = extract_individual_components(labeled_vol) - - # Create segmentations for each component - for component_vol in individual_components: - session_id = session_id_template.replace("{instance_id}", str(component_count)) - - # Create new segmentation + # One component volume at a time, so K components never occupy K full volumes at once + for component in range(1, n_components + 1): output_seg = run.new_segmentation( voxel_size=voxel_size, name=name, - session_id=session_id, + session_id=session_id_template.replace("{instance_id}", str(component_count)), is_multilabel=False, user_id=output_user_id, exist_ok=True, ) - - # Store the component volume - output_seg.from_numpy(component_vol) + output_seg.from_numpy((labeled_vol == component).astype(np.uint8)) output_segmentations.append(output_seg) component_count += 1 - print(f"Created {len(output_segmentations)} component segmentations") + logger.info(f"Created {len(output_segmentations)} component segmentations") return output_segmentations +def separate_components_as_instances( + segmentation: "CopickSegmentation", + run: "CopickRun", + object_name: str, + session_id: str, + user_id: str, + output_type: str, + connectivity: Union[int, str] = "all", + min_size: Optional[float] = None, +) -> Tuple["CopickSegmentation", Dict[str, int]]: + """ + Write the connected components of a segmentation as one instance or panoptic segmentation. + + An instance segmentation holds one object, so its input is a binary segmentation (or a multilabel one with a + single label); each component becomes an instance, numbered 1..K by size. A panoptic output keeps every label of + a multilabel input (or the object's label for a binary input) and numbers the components of each label 1..K by + size. + + Args: + segmentation: Input binary or multilabel segmentation. + run: CopickRun object. + object_name: Output name (the object, for an instance segmentation). + session_id: Output session ID. + user_id: Output user ID. + output_type: "instance" or "panoptic". + connectivity: Connectivity for connected components. + min_size: Minimum component volume in cubic angstroms (ų) to keep. + + Returns: + Tuple of (output segmentation, stats dict). + """ + from copick_utils.util.segmentations import new_segmentation, segmentation_type + + volume = segmentation.numpy() + if volume is None: + raise ValueError("Could not load segmentation data") + input_type = segmentation_type(segmentation) + if input_type not in ("binary", "multilabel"): + raise ValueError(f"Components are separated from binary or multilabel segmentations, not {input_type} ones.") + voxel_size = segmentation.voxel_size + labels = np.unique(volume) + labels = labels[labels > 0] + + if output_type == "instance": + if len(labels) > 1: + raise ValueError( + f"An instance segmentation holds one object, but {segmentation.name} has {len(labels)} labels; " + "write a panoptic segmentation (?panoptic=true) instead.", + ) + instances, n = components_as_instances(volume, voxel_size, connectivity, min_size) + output = new_segmentation(run, voxel_size, object_name, session_id, user_id, "instance") + output.from_numpy(instances) + return output, {"components_created": n} + + if output_type != "panoptic": + raise ValueError(f"Unknown output type {output_type!r}.") + if input_type == "binary": + obj = run.root.get_object(segmentation.name) + if obj is None or obj.label is None: + raise ValueError(f"{segmentation.name} is not a pickable object with a label.") + label_of = {int(v): int(obj.label) for v in labels} + else: + label_of = {int(v): int(v) for v in labels} + label_channel = np.zeros(volume.shape, dtype=np.uint32) + instance_channel = np.zeros(volume.shape, dtype=np.uint32) + total = 0 + for value in labels: + instances, n = components_as_instances(volume == value, voxel_size, connectivity, min_size) + kept = instances > 0 + label_channel[kept] = label_of[int(value)] + instance_channel[kept] = instances[kept] + total += n + output = new_segmentation(run, voxel_size, object_name, session_id, user_id, "panoptic") + output.from_numpy(np.stack([label_channel, instance_channel])) + return output, {"components_created": total} + + def separate_components_converter( segmentation: "CopickSegmentation", run: "CopickRun", @@ -262,28 +313,43 @@ def separate_components_converter( connectivity: Union[int, str] = "all", min_size: Optional[float] = None, multilabel: bool = True, + output_segmentation_type: Optional[str] = None, **kwargs, -) -> Optional[Tuple[None, Dict[str, int]]]: +) -> Optional[Tuple[Optional["CopickSegmentation"], Dict[str, int]]]: """ Lazy converter wrapper for separate_segmentation_components. - The session_id from the output URI is used as a template with {instance_id}. + By default each component becomes its own segmentation, and the session_id from the output URI is used as a + template with {instance_id}. When the output URI names an instance (``?instance=true``) or panoptic + (``?panoptic=true``) segmentation, the components are written as the instances of that one segmentation instead. Args: segmentation: Input CopickSegmentation object. run: CopickRun object. - object_name: Unused (components keep the input segmentation name). - session_id: Session ID template (should contain {instance_id}). + object_name: Output name (unused for per-component outputs, which keep the input segmentation name). + session_id: Session ID (a template with {instance_id} for per-component outputs). user_id: User ID for output segmentations. connectivity: Connectivity for connected components. min_size: Minimum component volume in cubic angstroms (ų) to keep. multilabel: Whether to treat input as multilabel segmentation. + output_segmentation_type: The segmentation type the output URI names, if any. **kwargs: Additional keyword arguments from lazy converter. Returns: - Tuple of (None, stats dict) or None if operation failed. + Tuple of (output segmentation or None, stats dict) or None if operation failed. """ try: + if output_segmentation_type in ("instance", "panoptic"): + return separate_components_as_instances( + segmentation, + run, + object_name, + session_id, + user_id, + output_segmentation_type, + connectivity=connectivity, + min_size=min_size, + ) output_segmentations = separate_segmentation_components( segmentation=segmentation, connectivity=connectivity, @@ -292,15 +358,10 @@ def separate_components_converter( output_user_id=user_id, multilabel=multilabel, ) - - stats = { - "components_created": len(output_segmentations), - } - - return None, stats + return None, {"components_created": len(output_segmentations)} except Exception as e: - get_logger(__name__).error(f"Error separating components in {run.name}: {e}") + logger.error(f"Error separating components in {run.name}: {e}") return None diff --git a/src/copick_utils/process/filter_components.py b/src/copick_utils/process/filter_components.py index 7be0e21..381a2cf 100644 --- a/src/copick_utils/process/filter_components.py +++ b/src/copick_utils/process/filter_components.py @@ -14,6 +14,53 @@ logger = get_logger(__name__) +def _keep_mask( + counts: np.ndarray, + voxel_volume: float, + min_size: Optional[float] = None, + max_size: Optional[float] = None, + keep_largest: Optional[int] = None, +) -> np.ndarray: + """Which ids (indices of ``counts``) pass the size filters. Index 0 (background) never passes; ids with no + voxels never pass.""" + volumes = counts * voxel_volume + keep = counts > 0 + keep[0] = False + if min_size is not None: + keep &= volumes >= min_size + if max_size is not None: + keep &= volumes <= max_size + if keep_largest is not None: + present = np.flatnonzero(counts[1:] > 0) + 1 + order = present[np.argsort(counts[present], kind="stable")] + largest = np.zeros(len(counts), dtype=bool) + largest[order[-keep_largest:] if keep_largest > 0 else []] = True + keep &= largest + return keep + + +def _apply_skeleton_length( + labeled: np.ndarray, + keep: np.ndarray, + voxel_spacing: float, + min_skeleton_length: Optional[float], +) -> np.ndarray: + """Also drop the kept ids whose skeleton is shorter than ``min_skeleton_length`` (angstroms).""" + if not min_skeleton_length: + return keep + from scipy.ndimage import find_objects + + from copick_utils.process.skeleton_graph import summarize + + keep = keep.copy() + for value, box in enumerate(find_objects(labeled), start=1): + if box is None or value >= len(keep) or not keep[value]: + continue + if summarize(labeled[box] == value).length * voxel_spacing < min_skeleton_length: + keep[value] = False + return keep + + def _filter_components_by_size( seg: np.ndarray, voxel_spacing: float, @@ -21,6 +68,7 @@ def _filter_components_by_size( min_size: Optional[float] = None, max_size: Optional[float] = None, keep_largest: Optional[int] = None, + min_skeleton_length: Optional[float] = None, ) -> Tuple[np.ndarray, int, int, list]: """ Filter connected components in a segmentation by size. @@ -37,6 +85,8 @@ def _filter_components_by_size( keep_largest: If set, keep only the N largest components by voxel count (None = no limit). Applied in addition to (intersected with) any min/max size filter, so a component is kept only if it both passes the size range and ranks within the N largest. + min_skeleton_length: Minimum skeleton length in angstroms (None = no minimum); drops specks and short + blobs of filamentous structures. Returns: Tuple of (seg_filtered, num_kept, num_removed, component_info) @@ -51,54 +101,44 @@ def _filter_components_by_size( "all": 3, } connectivity_value = connectivity_map.get(connectivity, 3) - struct = generate_binary_structure(seg.ndim, connectivity_value) + labeled_seg, num_components = label(seg, structure=struct) voxel_volume = voxel_spacing**3 - # Use bincount for O(n) counting of all components in one pass - counts = np.bincount(labeled_seg.ravel()) - - # When keeping only the N largest, determine which component ids qualify by rank. - # np.argsort is ascending; take the last N (largest), excluding background id 0. - largest_ids = None - if keep_largest is not None and num_components > 0: - component_ids = np.arange(1, num_components + 1) - order = component_ids[np.argsort(counts[1 : num_components + 1], kind="stable")] - largest_ids = {int(cid) for cid in order[-keep_largest:]} - - seg_filtered = np.zeros_like(seg, dtype=bool) - component_info = [] - num_kept = 0 - num_removed = 0 - - for component_id in range(1, num_components + 1): - component_voxels = int(counts[component_id]) - component_volume = component_voxels * voxel_volume - - passes_filter = True - if min_size is not None and component_volume < min_size: - passes_filter = False - if max_size is not None and component_volume > max_size: - passes_filter = False - if largest_ids is not None and component_id not in largest_ids: - passes_filter = False - - info = { + # One bincount and one lookup table: no pass over the volume per component + counts = np.bincount(labeled_seg.ravel(), minlength=num_components + 1) + keep = _keep_mask(counts, voxel_volume, min_size, max_size, keep_largest) + keep = _apply_skeleton_length(labeled_seg, keep, voxel_spacing, min_skeleton_length) + + component_info = [ + { "component_id": component_id, - "voxels": component_voxels, - "volume": component_volume, - "kept": passes_filter, + "voxels": int(counts[component_id]), + "volume": int(counts[component_id]) * voxel_volume, + "kept": bool(keep[component_id]), } - component_info.append(info) + for component_id in range(1, num_components + 1) + ] + num_kept = int(keep[1:].sum()) + return keep[labeled_seg].astype(np.uint8), num_kept, num_components - num_kept, component_info - if passes_filter: - seg_filtered = np.logical_or(seg_filtered, labeled_seg == component_id) - num_kept += 1 - else: - num_removed += 1 - return seg_filtered.astype(np.uint8), num_kept, num_removed, component_info +def _filter_instances( + instances: np.ndarray, + voxel_spacing: float, + min_size: Optional[float] = None, + max_size: Optional[float] = None, + keep_largest: Optional[int] = None, + min_skeleton_length: Optional[float] = None, +) -> Tuple[np.ndarray, int, int]: + """Drop whole instances of an instance segmentation by volume (and skeleton length); kept IDs are unchanged.""" + counts = np.bincount(instances.ravel()) if instances.size else np.zeros(1, dtype=int) + keep = _keep_mask(counts, voxel_spacing**3, min_size, max_size, keep_largest) + keep = _apply_skeleton_length(instances, keep, voxel_spacing, min_skeleton_length) + present = int((counts[1:] > 0).sum()) + num_kept = int(keep.sum()) + return np.where(keep[instances], instances, 0).astype(instances.dtype), num_kept, present - num_kept def filter_segmentation_components( @@ -111,11 +151,16 @@ def filter_segmentation_components( min_size: Optional[float] = None, max_size: Optional[float] = None, keep_largest: Optional[int] = None, + min_skeleton_length: Optional[float] = None, **kwargs, ) -> Optional[Tuple["CopickSegmentation", Dict[str, int]]]: """ Filter connected components in a segmentation by size. + A binary segmentation is filtered by its connected components, a multilabel one by the components of each label + (labels are kept), and an instance segmentation by whole instances, whose IDs are kept so that picks and + filaments with the same IDs still match. + Args: segmentation: Input CopickSegmentation object. run: CopickRun object. @@ -126,14 +171,16 @@ def filter_segmentation_components( min_size: Minimum component volume in cubic angstroms (ų) to keep. max_size: Maximum component volume in cubic angstroms (ų) to keep. keep_largest: If set, keep only the N largest components by voxel count. + min_skeleton_length: Minimum skeleton length in angstroms to keep. **kwargs: Additional keyword arguments from lazy converter. Returns: Tuple of (CopickSegmentation object, stats dict) or None if operation failed. """ + from copick_utils.util.segmentations import new_segmentation, segmentation_type + try: seg_array = segmentation.numpy() - if seg_array is None: logger.error("Could not load segmentation data") return None @@ -144,31 +191,45 @@ def filter_segmentation_components( # Use actual voxel_size from the segmentation (authoritative source) actual_voxel_spacing = segmentation.voxel_size + seg_type = segmentation_type(segmentation) + filters = { + "min_size": min_size, + "max_size": max_size, + "keep_largest": keep_largest, + "min_skeleton_length": min_skeleton_length, + } - bool_seg = seg_array.astype(bool) - - result_array, num_kept, num_removed, component_info = _filter_components_by_size( - bool_seg, - voxel_spacing=actual_voxel_spacing, - connectivity=connectivity, - min_size=min_size, - max_size=max_size, - keep_largest=keep_largest, - ) - - output_seg = run.new_segmentation( - name=object_name, - user_id=user_id, - session_id=session_id, - is_multilabel=segmentation.is_multilabel, - voxel_size=actual_voxel_spacing, - exist_ok=True, - ) - + if seg_type == "panoptic": + logger.error("filter-components does not filter panoptic segmentations; split them first") + return None + if seg_type == "instance": + result_array, num_kept, num_removed = _filter_instances(seg_array, actual_voxel_spacing, **filters) + elif seg_type == "multilabel": + result_array = np.zeros_like(seg_array) + num_kept = num_removed = 0 + for label_value in np.unique(seg_array[seg_array > 0]): + kept, k, r, _ = _filter_components_by_size( + seg_array == label_value, + voxel_spacing=actual_voxel_spacing, + connectivity=connectivity, + **filters, + ) + result_array[kept > 0] = label_value + num_kept += k + num_removed += r + else: + result_array, num_kept, num_removed, _ = _filter_components_by_size( + seg_array.astype(bool), + voxel_spacing=actual_voxel_spacing, + connectivity=connectivity, + **filters, + ) + + output_seg = new_segmentation(run, actual_voxel_spacing, object_name, session_id, user_id, seg_type) output_seg.from_numpy(result_array) stats = { - "voxels_kept": int(np.sum(result_array)), + "voxels_kept": int(np.count_nonzero(result_array)), "components_kept": num_kept, "components_removed": num_removed, "components_total": num_kept + num_removed, diff --git a/src/copick_utils/process/seg_stats.py b/src/copick_utils/process/seg_stats.py index 462ceba..b37dfc6 100644 --- a/src/copick_utils/process/seg_stats.py +++ b/src/copick_utils/process/seg_stats.py @@ -14,11 +14,60 @@ logger = get_logger(__name__) +def _component_rows( + labeled: np.ndarray, + voxel_spacing: float, + skeleton: bool = False, + ids: Optional[np.ndarray] = None, +) -> List[Dict[str, Any]]: + """ + One row per positive value of a labeled array (a component, or an instance). + + Args: + labeled: Integer array; each positive value is one component or instance. + voxel_spacing: Voxel spacing in angstroms. + skeleton: Also measure each one's skeleton (length, label radius, branches, junctions, endpoints). + ids: The values to report (default: every positive value present). + + Returns: + List of dicts with component_id, volume_voxels, volume_angstroms3 (and skeleton columns). + """ + from scipy.ndimage import find_objects + + counts = np.bincount(labeled.ravel()) if labeled.size else np.zeros(1, dtype=int) + if ids is None: + ids = np.flatnonzero(counts[1:]) + 1 + slices = find_objects(labeled) if skeleton else None + voxel_volume = voxel_spacing**3 + rows = [] + for value in ids: + value = int(value) + voxels = int(counts[value]) if value < len(counts) else 0 + row = {"component_id": value, "volume_voxels": voxels, "volume_angstroms3": voxels * voxel_volume} + if skeleton: + from copick_utils.process.skeleton_graph import summarize + + box = slices[value - 1] if value - 1 < len(slices) else None + summary = summarize(labeled[box] == value) if box is not None else summarize(np.zeros((1, 1, 1), bool)) + row.update( + { + "skeleton_length_angstroms": summary.length * voxel_spacing, + "label_radius_angstroms": None if summary.radius is None else summary.radius * voxel_spacing, + "n_branches": summary.n_branches, + "n_junctions": summary.n_junctions, + "n_endpoints": summary.n_endpoints, + }, + ) + rows.append(row) + return rows + + def _analyze_components_single( seg: np.ndarray, voxel_spacing: float, connectivity: str = "all", include_background: bool = True, + skeleton: bool = False, ) -> List[Dict[str, Any]]: """ Analyze connected components per label in a segmentation. @@ -29,6 +78,7 @@ def _analyze_components_single( connectivity: Connectivity for connected components. "face" = 6-connected, "face-edge" = 18-connected, "all" = 26-connected. include_background: If True, also analyze connected components of the background (label 0). + skeleton: If True, also measure each foreground component's skeleton. Returns: List of dicts, one per component: @@ -41,87 +91,145 @@ def _analyze_components_single( } connectivity_value = connectivity_map.get(connectivity, 3) struct = generate_binary_structure(seg.ndim, connectivity_value) - voxel_volume = voxel_spacing**3 all_unique_labels = np.unique(seg) foreground_labels = all_unique_labels[all_unique_labels != 0] components = [] - for label_value in foreground_labels: - binary_mask = seg == label_value - labeled_array, num_components = label(binary_mask, structure=struct) - - # Use bincount for O(n) counting of all components in one pass - counts = np.bincount(labeled_array.ravel()) - # counts[0] is background within this label mask, skip it - for component_id in range(1, num_components + 1): - component_voxels = int(counts[component_id]) - components.append( - { - "label": int(label_value), - "component_id": component_id, - "volume_voxels": component_voxels, - "volume_angstroms3": component_voxels * voxel_volume, - }, - ) + labeled_array, _ = label(seg == label_value, structure=struct) + for row in _component_rows(labeled_array, voxel_spacing, skeleton=skeleton): + components.append({"label": int(label_value), **row}) # Analyze background connected components (label 0) if include_background and 0 in all_unique_labels: - background_mask = seg == 0 - labeled_bg, num_bg_components = label(background_mask, structure=struct) - - if num_bg_components > 0: - bg_counts = np.bincount(labeled_bg.ravel()) - for component_id in range(1, num_bg_components + 1): - component_voxels = int(bg_counts[component_id]) - components.append( - { - "label": 0, - "component_id": component_id, - "volume_voxels": component_voxels, - "volume_angstroms3": component_voxels * voxel_volume, - }, - ) + labeled_bg, _ = label(seg == 0, structure=struct) + for row in _component_rows(labeled_bg, voxel_spacing): + components.append({"label": 0, **row}) return components +def _analyze_instances( + instances: np.ndarray, + voxel_spacing: float, + label_value: int, + skeleton: bool = False, +) -> List[Dict[str, Any]]: + """One row per instance ID of an instance segmentation (an instance may span several pieces).""" + return [ + {"label": int(label_value), "instance_id": row["component_id"], **row} + for row in _component_rows(instances, voxel_spacing, skeleton=skeleton) + ] + + +def _analyze_panoptic( + panoptic: np.ndarray, + voxel_spacing: float, + connectivity: str = "all", + include_background: bool = True, + skeleton: bool = False, +) -> List[Dict[str, Any]]: + """Rows for a panoptic segmentation: one per (label, instance) segment, and the connected components of each + label's region without instances (instance 0) and of the background.""" + labels, instances = panoptic[0], panoptic[1] + rows = [] + for label_value in np.unique(labels[labels > 0]): + in_label = labels == label_value + ids = np.unique(instances[in_label]) + ids = ids[ids > 0] + if len(ids): + per_label = np.where(in_label, instances, 0) + rows.extend(_analyze_instances(per_label, voxel_spacing, int(label_value), skeleton=skeleton)) + stuff = in_label & (instances == 0) + if stuff.any(): + rows.extend( + _analyze_components_single( + stuff.astype(np.uint8) * int(label_value), + voxel_spacing, + connectivity=connectivity, + include_background=False, + skeleton=skeleton, + ), + ) + if include_background: + rows.extend( + r + for r in _analyze_components_single( + (labels > 0).astype(np.uint8), + voxel_spacing, + connectivity=connectivity, + include_background=True, + ) + if r["label"] == 0 + ) + return rows + + def analyze_segmentation_components( segmentation: "CopickSegmentation", voxel_spacing: float, connectivity: str = "all", include_background: bool = True, + skeleton: bool = False, ) -> Optional[List[Dict[str, Any]]]: """ Analyze connected components in a CopickSegmentation. + Binary and multilabel segmentations give one row per connected component of each label. An instance + segmentation gives one row per instance ID, and a panoptic segmentation one row per (label, instance) segment + plus the connected components of regions without instances. + Args: segmentation: Input CopickSegmentation object. voxel_spacing: Voxel spacing in angstroms. connectivity: Connectivity for connected components. include_background: If True, also analyze background (label 0) components. + skeleton: If True, also measure each component's or instance's skeleton. Returns: List of component dicts, or None if loading failed. """ + from copick_utils.util.segmentations import segmentation_type + try: seg_array = segmentation.numpy() - if seg_array is None: logger.error("Could not load segmentation data") return None - if seg_array.size == 0: logger.error("Empty segmentation data") return None - return _analyze_components_single( - seg_array, - voxel_spacing=voxel_spacing, - connectivity=connectivity, - include_background=include_background, - ) + seg_type = segmentation_type(segmentation) + if seg_type == "instance": + obj = segmentation.run.root.get_object(segmentation.name) + label_value = obj.label if obj is not None else 0 + rows = _analyze_instances(seg_array, voxel_spacing, label_value, skeleton=skeleton) + if include_background: + rows.extend( + r + for r in _analyze_components_single( + (seg_array > 0).astype(np.uint8), + voxel_spacing, + connectivity=connectivity, + include_background=True, + ) + if r["label"] == 0 + ) + elif seg_type == "panoptic": + rows = _analyze_panoptic(seg_array, voxel_spacing, connectivity, include_background, skeleton) + else: + rows = _analyze_components_single( + seg_array, + voxel_spacing=voxel_spacing, + connectivity=connectivity, + include_background=include_background, + skeleton=skeleton, + ) + for row in rows: + row["segmentation_type"] = seg_type + return rows except Exception as e: logger.error(f"Error analyzing segmentation components: {e}") @@ -133,15 +241,17 @@ def _seg_stats_worker( input_uri: str, connectivity: str, include_background: bool = True, + skeleton: bool = False, ) -> Dict[str, Any]: """Worker function for batch segmentation stats. - Uses resolve_copick_objects for proper URI resolution with pattern support. + Uses resolve_copick_objects for proper URI resolution with pattern support. A URI without a type flag selects + binary and multilabel segmentations only; ``?instance=true`` or ``?panoptic=true`` selects those types. """ - from copick.util.uri import resolve_copick_objects + from copick_utils.util.segmentations import resolve_segmentations try: - segmentations = resolve_copick_objects(input_uri, run.root, "segmentation", run_name=run.name) + segmentations = resolve_segmentations(input_uri, run.root, run_name=run.name) if not segmentations: return {"processed": 0, "components": [], "errors": [f"No segmentation found for {run.name}"]} @@ -154,6 +264,7 @@ def _seg_stats_worker( voxel_spacing=segmentation.voxel_size, connectivity=connectivity, include_background=include_background, + skeleton=skeleton, ) if components: @@ -179,6 +290,7 @@ def seg_stats_batch( include_background: bool = True, run_names: Optional[List[str]] = None, workers: int = 8, + skeleton: bool = False, ) -> Dict[str, Any]: """ Batch analyze connected component sizes across multiple runs. @@ -190,6 +302,7 @@ def seg_stats_batch( include_background: If True, also analyze background (label 0) components. run_names: List of run names to process. If None, processes all runs. workers: Number of worker processes. + skeleton: If True, also measure each component's or instance's skeleton. Returns: Dictionary with processing results per run. @@ -207,6 +320,7 @@ def seg_stats_batch( input_uri=input_uri, connectivity=connectivity, include_background=include_background, + skeleton=skeleton, ) return results @@ -233,6 +347,17 @@ def export_stats_csv(results: Dict[str, Any], output_path: str) -> None: output.parent.mkdir(parents=True, exist_ok=True) fieldnames = ["run", "label", "component_id", "volume_voxels", "volume_angstroms3", "voxel_spacing"] + for optional in ( + "segmentation_type", + "instance_id", + "skeleton_length_angstroms", + "label_radius_angstroms", + "n_branches", + "n_junctions", + "n_endpoints", + ): + if any(optional in c for c in all_components): + fieldnames.append(optional) with open(output, "w", newline="") as f: writer = csv.DictWriter(f, fieldnames=fieldnames) diff --git a/src/copick_utils/process/skeleton_graph.py b/src/copick_utils/process/skeleton_graph.py new file mode 100644 index 0000000..dc24999 --- /dev/null +++ b/src/copick_utils/process/skeleton_graph.py @@ -0,0 +1,469 @@ +"""Skeleton graphs: a 3D skeleton's voxels as a graph, split into branches between endpoints and junctions. + +Shared by ``skeletonize`` (spur pruning), ``seg-stats --skeleton`` (skeleton length) and the filament tracer +(``seg2fil``, ``fit-spline``). Coordinates are voxel indices in (z, y, x) order and lengths are in voxels, unless a +function says otherwise. Everything here is numpy and scipy; nothing reads or writes copick objects. + +The graph joins 26-connected skeleton voxels. A thin 26-connected skeleton can still contain staircase corners, where +three voxels are pairwise adjacent; the diagonal edge of such a triangle duplicates the two shorter steps around it and +would make every voxel of a staircase look like a junction. ``skeleton_graph`` therefore drops an edge when a common +neighbour reaches both of its ends by strictly shorter steps. + +A topology (``decompose``) is a list of branches, each an ordered voxel path between two nodes. A node is a junction +(a cluster of adjacent voxels with three or more neighbours), an endpoint (one neighbour, or an isolated voxel), or the +start of a closed loop. +""" + +from dataclasses import dataclass, field +from typing import Dict, List, Optional, Tuple + +import numpy as np +from scipy import sparse +from scipy.sparse.csgraph import connected_components + +#: Half of the 26 neighbour offsets (lexicographically positive), so each edge is generated once. +_OFFSETS = np.array( + [(dz, dy, dx) for dz in (-1, 0, 1) for dy in (-1, 0, 1) for dx in (-1, 0, 1) if (dz, dy, dx) > (0, 0, 0)], +) + + +@dataclass +class SkeletonGraph: + """The voxels of a skeleton and the edges between 26-neighbours. + + Attributes: + coords: (N, 3) voxel indices (z, y, x), in C order. + adjacency: (N, N) symmetric sparse matrix of edge lengths (1, sqrt(2) or sqrt(3)). + """ + + coords: np.ndarray + adjacency: sparse.csr_matrix + + @property + def degree(self) -> np.ndarray: + """The number of neighbours of each voxel.""" + return np.diff(self.adjacency.indptr) + + def neighbours(self, i: int) -> Tuple[np.ndarray, np.ndarray]: + """The neighbours of voxel ``i`` and the lengths of the edges to them.""" + a, b = self.adjacency.indptr[i], self.adjacency.indptr[i + 1] + return self.adjacency.indices[a:b], self.adjacency.data[a:b] + + def subgraph(self, keep: np.ndarray) -> "SkeletonGraph": + """The graph of the kept voxels (a boolean mask over the voxels).""" + keep = np.asarray(keep, dtype=bool) + return SkeletonGraph(self.coords[keep], self.adjacency[keep][:, keep].tocsr()) + + +def skeleton_graph(skeleton: np.ndarray) -> SkeletonGraph: + """The graph of a 3D skeleton's voxels, without the redundant diagonals of staircase corners. + + Args: + skeleton: 3D boolean (or non-zero) array. + + Returns: + The skeleton graph. + """ + coords = np.argwhere(skeleton) + n = len(coords) + if n == 0: + return SkeletonGraph(coords.reshape(0, 3), sparse.csr_matrix((0, 0))) + shape = np.asarray(skeleton.shape) + linear = np.ravel_multi_index(coords.T, skeleton.shape) # ascending: argwhere returns C order + rows, cols, lengths = [], [], [] + for offset in _OFFSETS: + moved = coords + offset + inside = np.all((moved >= 0) & (moved < shape), axis=1) + source = np.flatnonzero(inside) + target_linear = np.ravel_multi_index(moved[inside].T, skeleton.shape) + position = np.minimum(np.searchsorted(linear, target_linear), n - 1) + hit = linear[position] == target_linear + rows.append(source[hit]) + cols.append(position[hit]) + lengths.append(np.full(int(hit.sum()), float(np.linalg.norm(offset)))) + rows, cols, lengths = np.concatenate(rows), np.concatenate(cols), np.concatenate(lengths) + keep = _non_redundant(rows, cols, lengths, n) + rows, cols, lengths = rows[keep], cols[keep], lengths[keep] + adjacency = sparse.coo_matrix( + (np.r_[lengths, lengths], (np.r_[rows, cols], np.r_[cols, rows])), + shape=(n, n), + ).tocsr() + return SkeletonGraph(coords, adjacency) + + +def _non_redundant(rows: np.ndarray, cols: np.ndarray, lengths: np.ndarray, n: int) -> np.ndarray: + """Which edges to keep: an edge is redundant when a common neighbour reaches both its ends by shorter steps.""" + keep = np.ones(len(rows), dtype=bool) + long_edges = np.flatnonzero(lengths > 1.0) + if len(long_edges) == 0: + return keep + full = sparse.coo_matrix( + (np.r_[lengths, lengths], (np.r_[rows, cols], np.r_[cols, rows])), + shape=(n, n), + ).tocsr() + for e in long_edges: + i, j, length = rows[e], cols[e], lengths[e] + ni, li = full.indices[full.indptr[i] : full.indptr[i + 1]], full.data[full.indptr[i] : full.indptr[i + 1]] + nj, lj = full.indices[full.indptr[j] : full.indptr[j + 1]], full.data[full.indptr[j] : full.indptr[j + 1]] + shorter_i = set(ni[li < length].tolist()) + if any(k in shorter_i for k in nj[lj < length].tolist()): + keep[e] = False + return keep + + +@dataclass +class Branch: + """An ordered voxel path between two topology nodes. + + Attributes: + path: Voxel indices into the graph, in order, both end voxels included. + start: Topology node at ``path[0]``. + end: Topology node at ``path[-1]``. + length: Length along the path, in voxels. + """ + + path: np.ndarray + start: int + end: int + length: float + + @property + def closed(self) -> bool: + """Whether the branch starts and ends at the same node.""" + return self.start == self.end + + +@dataclass +class Topology: + """A skeleton as branches between nodes. + + Attributes: + graph: The skeleton graph. + branches: The branches. + node_voxels: Voxel indices of each node (several for a junction cluster). + node_kind: ``"junction"``, ``"endpoint"`` or ``"loop"`` per node. + """ + + graph: SkeletonGraph + branches: List[Branch] = field(default_factory=list) + node_voxels: Dict[int, np.ndarray] = field(default_factory=dict) + node_kind: Dict[int, str] = field(default_factory=dict) + + def incident(self, node: int) -> List[Tuple[int, bool]]: + """The branch ends at a node, as (branch index, True if the branch starts there). A closed branch at the + node appears twice.""" + ends = [] + for b, branch in enumerate(self.branches): + if branch.start == node: + ends.append((b, True)) + if branch.end == node: + ends.append((b, False)) + return ends + + def junctions(self) -> List[int]: + """The junction nodes.""" + return [node for node, kind in self.node_kind.items() if kind == "junction"] + + @property + def length(self) -> float: + """The total length of the branches, in voxels.""" + return float(sum(b.length for b in self.branches)) + + +def decompose(graph: SkeletonGraph, min_loop_length: float = 3.0) -> Topology: + """Split a skeleton graph into branches between endpoints and junction clusters. + + Args: + graph: The skeleton graph. + min_loop_length: A branch that leaves a junction and returns to it is kept only when it is at least this long + (voxels); shorter ones are corners of a junction cluster, not loops. + + Returns: + The topology. + """ + n = len(graph.coords) + topology = Topology(graph) + if n == 0: + return topology + degree = graph.degree + junction_voxels = np.flatnonzero(degree >= 3) + cluster_of = np.full(n, -1) + if len(junction_voxels): + sub = graph.adjacency[junction_voxels][:, junction_voxels] + _, labels = connected_components(sub, directed=False) + cluster_of[junction_voxels] = labels + + node_of = np.full(n, -1) + next_node = 0 + for cluster in range(int(cluster_of.max()) + 1 if len(junction_voxels) else 0): + members = junction_voxels[cluster_of[junction_voxels] == cluster] + node_of[members] = next_node + topology.node_voxels[next_node] = members + topology.node_kind[next_node] = "junction" + next_node += 1 + for i in np.flatnonzero(degree <= 1): + node_of[i] = next_node + topology.node_voxels[next_node] = np.array([i]) + topology.node_kind[next_node] = "endpoint" + next_node += 1 + + visited = np.zeros(n, dtype=bool) + used_edges = set() + for i in np.flatnonzero(node_of >= 0): + visited[i] = True + if degree[i] == 0: + topology.branches.append(Branch(np.array([i]), int(node_of[i]), int(node_of[i]), 0.0)) + continue + neighbours, lengths = graph.neighbours(i) + for j, length in zip(neighbours, lengths): + if node_of[j] >= 0 and node_of[j] == node_of[i]: + continue # an edge inside a junction cluster + if (min(i, j), max(i, j)) in used_edges: + continue + path, total = _walk(graph, i, int(j), float(length), node_of) + used_edges.add((min(path[0], path[1]), max(path[0], path[1]))) + used_edges.add((min(path[-2], path[-1]), max(path[-2], path[-1]))) + visited[path] = True + start, end = int(node_of[path[0]]), int(node_of[path[-1]]) + if start == end and topology.node_kind[start] == "junction" and total < min_loop_length: + continue + topology.branches.append(Branch(np.asarray(path), start, end, total)) + + # Closed loops with no junction or endpoint: every voxel has two neighbours. + for i in np.flatnonzero(~visited): + if visited[i]: + continue + node_of[i] = next_node + topology.node_voxels[next_node] = np.array([i]) + topology.node_kind[next_node] = "loop" + neighbours, lengths = graph.neighbours(i) + path, total = _walk(graph, int(i), int(neighbours[0]), float(lengths[0]), node_of) + visited[path] = True + topology.branches.append(Branch(np.asarray(path), next_node, next_node, total)) + next_node += 1 + return topology + + +def _walk(graph: SkeletonGraph, start: int, step: int, length: float, node_of: np.ndarray) -> Tuple[List[int], float]: + """Follow two-neighbour voxels from ``start`` through ``step`` until a node voxel (or back to ``start``).""" + path = [start, step] + previous, current = start, step + while node_of[current] < 0: + neighbours, lengths = graph.neighbours(current) + others = [(k, w) for k, w in zip(neighbours, lengths) if k != previous] + if not others: + break + k, w = others[0] + previous, current = current, int(k) + path.append(current) + length += float(w) + if current == start: + break + return path, length + + +def splice(topology: Topology) -> Topology: + """Join the two branches at every junction that has exactly two branch ends (a kink, not a junction), and turn + a junction with one branch end into an endpoint.""" + changed = True + while changed: + changed = False + for node in list(topology.junctions()): + ends = topology.incident(node) + if len(ends) == 1: + topology.node_kind[node] = "endpoint" + changed = True + elif len(ends) == 2 and ends[0][0] != ends[1][0]: + (a, a_starts), (b, b_starts) = ends + first, second = topology.branches[a], topology.branches[b] + path_a = first.path[::-1] if a_starts else first.path # ends at the node + path_b = second.path if b_starts else second.path[::-1] # starts at the node + start = first.end if a_starts else first.start + end = second.end if b_starts else second.start + gap = np.linalg.norm(topology.graph.coords[path_a[-1]] - topology.graph.coords[path_b[0]]) + path = np.concatenate([path_a, path_b[1:] if path_a[-1] == path_b[0] else path_b]) + joined = Branch(path, start, end, first.length + second.length + float(gap)) + topology.branches = [br for k, br in enumerate(topology.branches) if k not in (a, b)] + [joined] + del topology.node_kind[node] + del topology.node_voxels[node] + changed = True + break + elif len(ends) == 0: + del topology.node_kind[node] + del topology.node_voxels[node] + changed = True + break + return topology + + +def prune_spurs(topology: Topology, min_length: float) -> Topology: + """Remove terminal branches (one end an endpoint, the other a junction) shorter than ``min_length`` voxels. + + Repeats until no spur is left, splicing the junctions that pruning reduces to two branch ends. A junction whose + branches are all short spurs keeps its longest one. Branches between two endpoints are never removed: they are a + whole piece, not a spur. + + Args: + topology: The topology (modified in place and returned). + min_length: Spurs shorter than this (voxels) are removed. + + Returns: + The pruned topology. + """ + while True: + spurs = [] + for b, branch in enumerate(topology.branches): + kinds = (topology.node_kind.get(branch.start), topology.node_kind.get(branch.end)) + if sorted(kinds) == ["endpoint", "junction"] and branch.length < min_length: + spurs.append(b) + if not spurs: + return topology + spur_set = set(spurs) + for node in topology.junctions(): + ends = topology.incident(node) + if ends and all(b in spur_set for b, _ in ends): + spur_set.discard(max((b for b, _ in ends), key=lambda b: topology.branches[b].length)) + if not spur_set: + return topology + for b in spur_set: + branch = topology.branches[b] + for node in (branch.start, branch.end): + if topology.node_kind.get(node) == "endpoint": + topology.node_kind.pop(node) + topology.node_voxels.pop(node) + topology.branches = [br for k, br in enumerate(topology.branches) if k not in spur_set] + splice(topology) + + +def merge_junctions(topology: Topology, max_length: float) -> Topology: + """Contract branches shorter than ``max_length`` voxels that join two different junctions, so two nearby + junctions (two filaments crossing, seen as two forks joined by a short bridge) become one. + + Args: + topology: The topology (modified in place and returned). + max_length: Bridges up to this length (voxels) are contracted. + + Returns: + The topology. + """ + parent = {node: node for node in topology.node_kind} + + def find(node): + while parent[node] != node: + parent[node] = parent[parent[node]] + node = parent[node] + return node + + bridges = [] + for b, branch in enumerate(topology.branches): + if ( + branch.start != branch.end + and topology.node_kind.get(branch.start) == "junction" + and topology.node_kind.get(branch.end) == "junction" + and branch.length <= max_length + ): + bridges.append(b) + parent[find(branch.start)] = find(branch.end) + if not bridges: + return topology + bridge_set = set(bridges) + merged_voxels: Dict[int, List[np.ndarray]] = {} + for node, voxels in topology.node_voxels.items(): + merged_voxels.setdefault(find(node), []).append(voxels) + for b in bridges: + merged_voxels[find(topology.branches[b].start)].append(topology.branches[b].path) + kept = [] + for b, branch in enumerate(topology.branches): + if b in bridge_set: + continue + branch.start, branch.end = find(branch.start), find(branch.end) + kept.append(branch) + topology.branches = kept + topology.node_kind = {find(node): kind for node, kind in topology.node_kind.items()} + topology.node_voxels = {node: np.unique(np.concatenate(parts)) for node, parts in merged_voxels.items()} + return topology + + +def node_centre(topology: Topology, node: int) -> np.ndarray: + """The mean voxel position (z, y, x) of a node.""" + return topology.graph.coords[topology.node_voxels[node]].mean(axis=0) + + +def pruned_skeleton(skeleton: np.ndarray, min_length: float) -> np.ndarray: + """A skeleton with its spurs shorter than ``min_length`` voxels removed (see ``prune_spurs``). + + Args: + skeleton: 3D boolean (or non-zero) array. + min_length: Spur length threshold, in voxels. + + Returns: + Boolean array of the same shape. + """ + graph = skeleton_graph(skeleton) + topology = prune_spurs(decompose(graph), min_length) + keep = np.zeros(len(graph.coords), dtype=bool) + for branch in topology.branches: + keep[branch.path] = True + for node, voxels in topology.node_voxels.items(): + if topology.node_kind[node] == "junction": + keep[voxels] = True + out = np.zeros(skeleton.shape, dtype=bool) + out[tuple(graph.coords[keep].T)] = True + return out + + +def skeleton_length(skeleton: np.ndarray, voxel_size: float = 1.0) -> float: + """The total length of a skeleton's branches, in the unit of ``voxel_size`` (voxels by default).""" + return decompose(skeleton_graph(skeleton)).length * float(voxel_size) + + +@dataclass +class SkeletonSummary: + """Skeleton measurements of one component or instance. + + Attributes: + length: Total branch length (voxels). + n_branches: Number of branches. + n_junctions: Number of junctions. + n_endpoints: Number of endpoints. + radius: Median distance from the skeleton to the background (voxels), or None without a skeleton. + """ + + length: float + n_branches: int + n_junctions: int + n_endpoints: int + radius: Optional[float] + + +def summarize(mask: np.ndarray, prune_length: float = 0.0) -> SkeletonSummary: + """Skeletonize a binary mask (Lee) and measure its skeleton. + + Args: + mask: 3D boolean array of one component or instance. + prune_length: Remove spurs shorter than this (voxels) before measuring; 0 keeps every branch. + + Returns: + The measurements. + """ + from scipy import ndimage + from skimage.morphology import skeletonize + + padded = np.pad(np.asarray(mask, dtype=bool), 1) + skeleton = skeletonize(padded) + if not skeleton.any(): + return SkeletonSummary(0.0, 0, 0, 0, None) + topology = decompose(skeleton_graph(skeleton)) + if prune_length > 0: + topology = prune_spurs(topology, prune_length) + edt = ndimage.distance_transform_edt(padded) + on_skeleton = topology.graph.coords + radius = float(np.median(edt[tuple(on_skeleton.T)])) if len(on_skeleton) else None + kinds = list(topology.node_kind.values()) + return SkeletonSummary( + length=topology.length, + n_branches=len(topology.branches), + n_junctions=kinds.count("junction"), + n_endpoints=kinds.count("endpoint"), + radius=radius, + ) diff --git a/src/copick_utils/process/split_labels.py b/src/copick_utils/process/split_labels.py index d94174b..1378da8 100644 --- a/src/copick_utils/process/split_labels.py +++ b/src/copick_utils/process/split_labels.py @@ -117,6 +117,68 @@ def split_multilabel_segmentation( return output_segmentations +def split_panoptic_segmentation( + segmentation: "CopickSegmentation", + run: "CopickRun", + output_user_id: str = "split", + labels: Optional[Dict[str, int]] = None, +) -> List["CopickSegmentation"]: + """ + Split a panoptic segmentation into per-object segmentations. + + For each label, the voxels with an instance ID become an instance segmentation named after the object (the IDs + are kept), and the voxels without one (regions such as membranes) become a binary segmentation of that object. + + Args: + segmentation: Input panoptic segmentation. + run: CopickRun object containing the segmentation. + output_user_id: User ID for output segmentations (default: "split"). + labels: Optional explicit ``{name: label_value}`` mapping, as for ``split_multilabel_segmentation``. + + Returns: + List of created CopickSegmentation objects. + """ + from copick_utils.util.segmentations import new_segmentation + + volume = segmentation.numpy() + if volume is None or volume.size == 0: + raise ValueError("Could not load segmentation data") + label_channel, instance_channel = volume[0], volume[1] + root = run.root + if labels: + label_items = [(int(value), name) for name, value in labels.items()] + else: + label_items = [] + for value in np.unique(label_channel[label_channel > 0]): + obj = next((o for o in root.config.pickable_objects if o.label == int(value)), None) + label_items.append((int(value), obj.name if obj is not None else str(int(value)))) + + outputs = [] + for label_value, object_name in label_items: + in_label = label_channel == label_value + things = in_label & (instance_channel > 0) + stuff = in_label & (instance_channel == 0) + for seg_type, mask, data in ( + ("instance", things, np.where(things, instance_channel, 0)), + ("binary", stuff, stuff.astype(np.uint8)), + ): + if not mask.any(): + continue + output = new_segmentation( + run, + segmentation.voxel_size, + object_name, + segmentation.session_id, + output_user_id, + seg_type, + ) + output.from_numpy(data) + outputs.append(output) + if outputs: + logger.info(f"Run '{run.name}': split panoptic segmentation into {len(outputs)} segmentations") + return outputs + + def _split_labels_worker( run: "CopickRun", segmentation_name: str, @@ -125,9 +187,10 @@ def _split_labels_worker( voxel_spacing: float, output_user_id: str, labels: Optional[Dict[str, int]] = None, + panoptic: bool = False, ) -> Dict[str, Any]: """ - Worker function for batch splitting of multilabel segmentations. + Worker function for batch splitting of multilabel (or panoptic) segmentations. Args: run: CopickRun to process @@ -142,6 +205,24 @@ def _split_labels_worker( """ try: # Get the input segmentation + if panoptic: + segmentations = run.get_segmentations( + name=segmentation_name, + user_id=segmentation_user_id, + session_id=segmentation_session_id, + voxel_size=voxel_spacing, + is_panoptic=True, + ) + if not segmentations: + return {"processed": 0, "errors": [f"No panoptic segmentation found for run {run.name}"]} + outputs = split_panoptic_segmentation(segmentations[0], run, output_user_id, labels) + return { + "processed": 1, + "errors": [], + "labels_split": len(outputs), + "object_names": [seg.name for seg in outputs], + } + segmentations = run.get_segmentations( name=segmentation_name, user_id=segmentation_user_id, @@ -195,6 +276,7 @@ def split_labels_batch( run_names: Optional[List[str]] = None, workers: int = 8, labels: Optional[Dict[str, int]] = None, + panoptic: bool = False, ) -> Dict[str, Any]: """ Batch split multilabel segmentations across multiple runs. @@ -210,6 +292,7 @@ def split_labels_batch( workers: Number of worker processes (default: 8) labels: Optional explicit ``{name: label_value}`` mapping for naming outputs (see ``split_multilabel_segmentation``). None = resolve names from config. + panoptic: Split a panoptic segmentation (see ``split_panoptic_segmentation``). Returns: Dictionary with processing results and statistics per run @@ -230,6 +313,7 @@ def split_labels_batch( voxel_spacing=voxel_spacing, output_user_id=output_user_id, labels=labels, + panoptic=panoptic, ) return results diff --git a/src/copick_utils/util/config_models.py b/src/copick_utils/util/config_models.py index 5629f80..07876b1 100644 --- a/src/copick_utils/util/config_models.py +++ b/src/copick_utils/util/config_models.py @@ -4,12 +4,15 @@ from pydantic import BaseModel, Field, field_validator +#: The object types a selector reads or writes. +ObjectType = Literal["picks", "mesh", "segmentation", "filaments"] + class SelectorConfig(BaseModel): """Pydantic model for selector configuration with validation.""" - input_type: Literal["picks", "mesh", "segmentation"] - output_type: Literal["picks", "mesh", "segmentation"] + input_type: ObjectType + output_type: ObjectType input_object_name: str input_user_id: str input_session_id: str @@ -19,6 +22,23 @@ class SelectorConfig(BaseModel): individual_outputs: bool = False segmentation_name: Optional[str] = None voxel_spacing: Optional[float] = None + # Segmentation type flags from the URIs' query parameters (?multilabel=true, ?instance=true, ?panoptic=true). + # An input without an instance or panoptic flag selects binary and multilabel segmentations only. + input_multilabel: Optional[bool] = None + input_instance: Optional[bool] = None + input_panoptic: Optional[bool] = None + output_multilabel: Optional[bool] = None + output_instance: Optional[bool] = None + output_panoptic: Optional[bool] = None + + @property + def output_segmentation_type(self) -> Optional[str]: + """The segmentation type the output URI names (``multilabel``, ``instance``, ``panoptic``), or None.""" + from copick_utils.util.segmentations import type_from_flags + + if self.output_type != "segmentation": + return None + return type_from_flags(self.output_multilabel, self.output_instance, self.output_panoptic) @field_validator("segmentation_name") @classmethod @@ -86,9 +106,9 @@ def validate_output_session_id(cls, v, info): def from_uris( cls, input_uri: str, - input_type: Literal["picks", "mesh", "segmentation"], + input_type: ObjectType, output_uri: str, - output_type: Literal["picks", "mesh", "segmentation"], + output_type: ObjectType, individual_outputs: bool = False, command_name: Optional[str] = None, ) -> "SelectorConfig": @@ -155,9 +175,19 @@ def from_uris( voxel_spacing = None if voxel_spacing == "*" else float(voxel_spacing) config_dict["voxel_spacing"] = voxel_spacing + if input_type == "segmentation": + config_dict.update(_type_flags(input_params, "input")) + if output_type == "segmentation": + config_dict.update(_type_flags(output_params, "output")) + return cls(**config_dict) +def _type_flags(params: Dict[str, Any], prefix: str) -> Dict[str, Optional[bool]]: + """A parsed segmentation URI's type flags as ``_multilabel/instance/panoptic`` fields.""" + return {f"{prefix}_{flag}": params.get(flag) for flag in ("multilabel", "instance", "panoptic")} + + class ReferenceConfig(BaseModel): """Pydantic model for reference discovery configuration.""" @@ -167,6 +197,10 @@ class ReferenceConfig(BaseModel): session_id: Optional[str] = None voxel_spacing: Optional[float] = None tomo_type: Optional[str] = None + # Segmentation references: the URI's type flags (without an instance or panoptic flag, binary or multilabel). + multilabel: Optional[bool] = None + instance: Optional[bool] = None + panoptic: Optional[bool] = None additional_params: Dict[str, Any] = Field(default_factory=dict) @field_validator("voxel_spacing") @@ -255,6 +289,9 @@ def from_uri( user_id=user_id, session_id=session_id, voxel_spacing=voxel_spacing, + multilabel=params.get("multilabel"), + instance=params.get("instance"), + panoptic=params.get("panoptic"), additional_params=additional_params or {}, ) @@ -315,9 +352,9 @@ def validate_reference(cls, v, info): # URI-based convenience functions (simplified interface) def create_simple_config( input_uri: str, - input_type: Literal["picks", "mesh", "segmentation"], + input_type: ObjectType, output_uri: str, - output_type: Literal["picks", "mesh", "segmentation"], + output_type: ObjectType, individual_outputs: bool = False, command_name: Optional[str] = None, ) -> TaskConfig: @@ -357,9 +394,9 @@ def create_simple_config( def create_single_selector_config( input_uri: str, - input_type: Literal["picks", "mesh", "segmentation"], + input_type: ObjectType, output_uri: str, - output_type: Literal["picks", "mesh", "segmentation"], + output_type: ObjectType, command_name: Optional[str] = None, operation: str = "union", ) -> TaskConfig: @@ -416,9 +453,9 @@ def create_single_selector_config( def create_dual_selector_config( input1_uri: str, input2_uri: str, - input_type: Literal["picks", "mesh", "segmentation"], + input_type: ObjectType, output_uri: str, - output_type: Literal["picks", "mesh", "segmentation"], + output_type: ObjectType, pairing_method: str = "index_order", individual_outputs: bool = False, command_name: Optional[str] = None, @@ -490,6 +527,7 @@ def create_dual_selector_config( selector2_dict["segmentation_name"] = seg_name selector2_dict["voxel_spacing"] = voxel_spacing + selector2_dict.update(_type_flags(input2_params, "input")) selector2_config = SelectorConfig(**selector2_dict) @@ -502,9 +540,9 @@ def create_dual_selector_config( def create_multi_selector_config( input_uris: List[str], - input_type: Literal["picks", "mesh", "segmentation"], + input_type: ObjectType, output_uri: str, - output_type: Literal["picks", "mesh", "segmentation"], + output_type: ObjectType, pairing_method: str = "n_way", individual_outputs: bool = False, command_name: Optional[str] = None, @@ -582,6 +620,7 @@ def create_multi_selector_config( selector_dict["segmentation_name"] = object_name selector_dict["voxel_spacing"] = voxel_spacing + selector_dict.update(_type_flags(params, "input")) selectors.append(SelectorConfig(**selector_dict)) @@ -594,9 +633,9 @@ def create_multi_selector_config( def create_reference_config( input_uri: str, - input_type: Literal["picks", "mesh", "segmentation"], + input_type: ObjectType, output_uri: str, - output_type: Literal["picks", "mesh", "segmentation"], + output_type: ObjectType, reference_uri: str, reference_type: Literal["mesh", "segmentation", "tomogram"], additional_params: Optional[Dict[str, Any]] = None, diff --git a/src/copick_utils/util/segmentations.py b/src/copick_utils/util/segmentations.py new file mode 100644 index 0000000..c56745b --- /dev/null +++ b/src/copick_utils/util/segmentations.py @@ -0,0 +1,148 @@ +"""Segmentation types: how a URI's type flags select and create binary, multilabel, instance and panoptic stores. + +copick has four segmentation types. A URI names one with a query flag (``?multilabel=true``, ``?instance=true``, +``?panoptic=true``). copick-utils treats a URI without a flag as binary or multilabel, so a tool that does not ask for +instance or panoptic segmentations never reads one by accident (an instance store can share its name, user and session +with a binary one, and a panoptic store is 4D). +""" + +from typing import TYPE_CHECKING, Any, Dict, Literal, Optional + +if TYPE_CHECKING: + from copick.models import CopickRun, CopickSegmentation + +SegmentationType = Literal["binary", "multilabel", "instance", "panoptic"] + + +def selection_filters( + multilabel: Optional[bool] = None, + instance: Optional[bool] = None, + panoptic: Optional[bool] = None, +) -> Dict[str, Optional[bool]]: + """The type filters for copick's segmentation queries (``get_copick_objects_by_type``), from a URI's flags. + + Args: + multilabel: The URI's multilabel flag, or None. + instance: The URI's instance flag, or None. + panoptic: The URI's panoptic flag, or None. + + Returns: + ``{"multilabel": ..., "instance": ..., "panoptic": ...}``. Without an instance or panoptic flag, both are False. + """ + return {"multilabel": multilabel, "instance": bool(instance), "panoptic": bool(panoptic)} + + +def query_kwargs( + multilabel: Optional[bool] = None, + instance: Optional[bool] = None, + panoptic: Optional[bool] = None, +) -> Dict[str, Optional[bool]]: + """The type keywords for ``CopickRun.get_segmentations``, from a URI's flags (see ``selection_filters``).""" + filters = selection_filters(multilabel, instance, panoptic) + return { + "is_multilabel": filters["multilabel"], + "is_instance": filters["instance"], + "is_panoptic": filters["panoptic"], + } + + +def type_from_flags( + multilabel: Optional[bool] = None, + instance: Optional[bool] = None, + panoptic: Optional[bool] = None, +) -> Optional[SegmentationType]: + """The segmentation type a URI's flags name, or None when it names none. + + Raises: + ValueError: If more than one flag is set. + """ + named = [t for t, flag in (("multilabel", multilabel), ("instance", instance), ("panoptic", panoptic)) if flag] + if len(named) > 1: + raise ValueError(f"A segmentation has one type; the URI names {', '.join(named)}.") + return named[0] if named else None + + +def type_from_uri(uri: str) -> Optional[SegmentationType]: + """The segmentation type a segmentation URI names with its query flags, or None.""" + from copick.util.uri import parse_copick_uri + + params = parse_copick_uri(uri, "segmentation") + return type_from_flags(params.get("multilabel"), params.get("instance"), params.get("panoptic")) + + +def resolve_segmentations(uri: str, root, run_name: Optional[str] = None) -> list: + """The segmentations a URI selects, of the type its flags name (binary or multilabel when it names none). + + Args: + uri: Segmentation URI (patterns allowed). + root: Copick root. + run_name: Restrict to one run. + + Returns: + List of segmentations. + """ + from copick.util.uri import parse_copick_uri, resolve_copick_objects + + params = parse_copick_uri(uri, "segmentation") + wanted = selection_filters(params.get("multilabel"), params.get("instance"), params.get("panoptic")) + found = resolve_copick_objects(uri, root, "segmentation", run_name=run_name) + return [ + s + for s in found + if (wanted["multilabel"] is None or bool(s.is_multilabel) == wanted["multilabel"]) + and bool(getattr(s, "is_instance", False)) == wanted["instance"] + and bool(getattr(s, "is_panoptic", False)) == wanted["panoptic"] + ] + + +def segmentation_type(segmentation: "CopickSegmentation") -> SegmentationType: + """The type of an existing segmentation.""" + if getattr(segmentation, "is_panoptic", False): + return "panoptic" + if getattr(segmentation, "is_instance", False): + return "instance" + return "multilabel" if segmentation.is_multilabel else "binary" + + +def new_segmentation( + run: "CopickRun", + voxel_size: float, + name: str, + session_id: str, + user_id: str, + seg_type: SegmentationType, + exist_ok: bool = True, + **kwargs: Any, +) -> "CopickSegmentation": + """Create (or reuse) a segmentation of the given type. + + Args: + run: Run to write into. + voxel_size: Voxel size in Angstrom. + name: Object name (binary, instance) or descriptive name (multilabel, panoptic). + session_id: Session ID. + user_id: User ID. + seg_type: ``"binary"``, ``"multilabel"``, ``"instance"`` or ``"panoptic"``. + exist_ok: Reuse an existing segmentation with the same key. + **kwargs: Passed to ``CopickRun.new_segmentation``. + + Returns: + The segmentation. + """ + if seg_type not in ("binary", "multilabel", "instance", "panoptic"): + raise ValueError(f"Unknown segmentation type {seg_type!r}.") + type_kwargs = {} + if seg_type == "instance": + type_kwargs["is_instance"] = True + elif seg_type == "panoptic": + type_kwargs["is_panoptic"] = True + return run.new_segmentation( + voxel_size=voxel_size, + name=name, + session_id=session_id, + is_multilabel=seg_type == "multilabel", + user_id=user_id, + exist_ok=exist_ok, + **type_kwargs, + **kwargs, + ) diff --git a/tests/conftest.py b/tests/conftest.py index 8faadb1..807e7c0 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -41,3 +41,9 @@ def copick_root(tmp_path): def run(copick_root): """The fixture project's single run.""" return copick_root.get_run("run1") + + +@pytest.fixture +def config_path(copick_root, tmp_path): + """The fixture project's config file (for CLI tests).""" + return str(tmp_path / "copick_config.json") diff --git a/tests/test_instance_tools.py b/tests/test_instance_tools.py new file mode 100644 index 0000000..619052c --- /dev/null +++ b/tests/test_instance_tools.py @@ -0,0 +1,227 @@ +"""Instance and panoptic segmentations through the component tools, seg-stats, combine and split.""" + +import csv + +import numpy as np +from click.testing import CliRunner +from copick_utils.cli.combine_labels import combine +from copick_utils.cli.filter_components import filter_components +from copick_utils.cli.seg_stats import seg_stats +from copick_utils.cli.separate_components import separate_components +from copick_utils.cli.split_labels import split +from copick_utils.converters.lazy_converter import discover_tasks_for_run +from copick_utils.io import readers +from copick_utils.util.config_models import SelectorConfig + +VS = 10.0 +SHAPE = (20, 40, 60) + + +def _write(run, name, session, array, user="test", **kwargs): + run.new_voxel_spacing(VS, exist_ok=True) + seg = run.new_segmentation(VS, name, session, user_id=user, exist_ok=True, **kwargs) + seg.from_numpy(array) + return seg + + +def _get(run, name, user, session, **kwargs): + run.refresh() + found = run.get_segmentations(name=name, user_id=user, session_id=session, voxel_size=VS, **kwargs) + assert len(found) == 1, found + return found[0] + + +def _invoke(command, args): + result = CliRunner().invoke(command, args, catch_exceptions=False) + assert result.exit_code == 0, result.output + return result + + +def _three_blobs(): + volume = np.zeros(SHAPE, dtype=np.uint8) + volume[2:6, 2:6, 2:6] = 1 # 64 voxels + volume[10:16, 10:16, 10:16] = 1 # 216 voxels + volume[2:4, 30:32, 50:52] = 1 # 8 voxels + return volume + + +def test_separate_components_writes_one_instance_segmentation(run, config_path): + _write(run, "microtubule", "s1", _three_blobs()) + _invoke( + separate_components, + [ + "-c", + config_path, + "-i", + "microtubule:test/s1@10.0", + "--binary", + "-o", + "microtubule:components/s1@10.0?instance=true", + ], + ) + instances = _get(run, "microtubule", "components", "s1", is_instance=True).numpy() + assert sorted(np.unique(instances).tolist()) == [0, 1, 2, 3] + assert (instances == 1).sum() == 216 and (instances == 2).sum() == 64 and (instances == 3).sum() == 8 + + +def test_filter_instances_keeps_ids_and_skeleton_length(run, config_path): + instances = np.zeros(SHAPE, dtype=np.uint16) + instances[2:6, 2:6, 2:6] = 4 # compact blob, 64 voxels + instances[10:13, 10:13, 5:55] = 7 # rod, 450 voxels, skeleton ~45 voxels + instances[2:4, 30:32, 50:52] = 9 # speck, 8 voxels + _write(run, "microtubule", "s2", instances, is_instance=True) + + _invoke( + filter_components, + [ + "-c", + config_path, + "-i", + "microtubule:test/s2@10.0?instance=true", + "-o", + "microtubule:filtered/s2@10.0", + "--min-size", + "20", + "--size-unit", + "voxel", + ], + ) + kept = _get(run, "microtubule", "filtered", "s2", is_instance=True).numpy() + assert sorted(np.unique(kept).tolist()) == [0, 4, 7] + + _invoke( + filter_components, + [ + "-c", + config_path, + "-i", + "microtubule:test/s2@10.0?instance=true", + "-o", + "microtubule:long/s2@10.0", + "--min-skeleton-length", + "200", + ], + ) + long_only = _get(run, "microtubule", "long", "s2", is_instance=True).numpy() + assert sorted(np.unique(long_only).tolist()) == [0, 7] + + +def test_filter_multilabel_keeps_labels(run, config_path): + labels = np.zeros(SHAPE, dtype=np.uint8) + labels[2:6, 2:6, 2:6] = 1 + labels[10:16, 10:16, 10:16] = 3 + labels[2:4, 30:32, 50:52] = 3 # small component of label 3 + _write(run, "labels", "s3", labels, is_multilabel=True) + _invoke( + filter_components, + [ + "-c", + config_path, + "-i", + "labels:test/s3@10.0?multilabel=true", + "-o", + "labels:filtered/s3@10.0", + "--min-size", + "20", + "--size-unit", + "voxel", + ], + ) + out = _get(run, "labels", "filtered", "s3", is_multilabel=True) + assert out.is_multilabel + values = out.numpy() + assert sorted(np.unique(values).tolist()) == [0, 1, 3] + assert values[3, 31, 51] == 0 + + +def test_seg_stats_instance_rows(run, config_path, tmp_path): + instances = np.zeros(SHAPE, dtype=np.uint16) + instances[10:13, 10:13, 5:55] = 7 + instances[2:6, 2:6, 2:6] = 4 + _write(run, "microtubule", "s4", instances, is_instance=True) + out = tmp_path / "stats.csv" + _invoke( + seg_stats, + [ + "-c", + config_path, + "-i", + "microtubule:test/s4@10.0?instance=true", + "--skeleton", + "--no-include-background", + "-f", + "csv", + "-op", + str(out), + ], + ) + rows = {int(r["instance_id"]): r for r in csv.DictReader(out.open())} + assert set(rows) == {4, 7} + assert rows[7]["segmentation_type"] == "instance" and int(rows[7]["label"]) == 1 + assert int(rows[7]["volume_voxels"]) == 450 + assert 400 < float(rows[7]["skeleton_length_angstroms"]) <= 500 + + +def test_untyped_selection_ignores_a_same_key_instance_store(run): + binary = np.zeros(SHAPE, dtype=np.uint8) + binary[1:3, 1:3, 1:3] = 1 + instances = np.zeros(SHAPE, dtype=np.uint16) + instances[5:8, 5:8, 5:8] = 3 + _write(run, "microtubule", "same", binary) + _write(run, "microtubule", "same", instances, is_instance=True) + + selector = SelectorConfig.from_uris( + "microtubule:test/same@10.0", + "segmentation", + "microtubule:out/same@10.0", + "segmentation", + ) + tasks = discover_tasks_for_run(run, selector) + assert len(tasks) == 1 and not tasks[0]["segmentation"].is_instance + + # An output without a voxel spacing inherits the input's spacing and type + selector = SelectorConfig.from_uris( + "microtubule:test/same@10.0?instance=true", + "segmentation", + "microtubule:out/same", + "segmentation", + ) + tasks = discover_tasks_for_run(run, selector) + assert len(tasks) == 1 and tasks[0]["segmentation"].is_instance + assert tasks[0]["output_segmentation_type"] == "instance" + + np.testing.assert_array_equal(readers.segmentation(run, VS, "microtubule", "test", "same"), binary) + np.testing.assert_array_equal( + readers.segmentation(run, VS, "microtubule", "test", "same", is_instance=True), + instances, + ) + + +def test_panoptic_combine_and_split_round_trip(run, config_path): + membrane = np.zeros(SHAPE, dtype=np.uint8) + membrane[0:2] = 1 + instances = np.zeros(SHAPE, dtype=np.uint16) + instances[5:8, 5:8, 5:50] = 2 + instances[12:15, 20:23, 5:50] = 5 + _write(run, "membrane", "p", membrane) + _write(run, "microtubule", "p", instances, is_instance=True) + + _invoke( + combine, + [ + "-c", + config_path, + "-i", + "membrane:test/p@10.0", + "--instances", + "microtubule:test/p@10.0", + "-o", + "cell:combine/p@10.0?panoptic=true", + ], + ) + panoptic = _get(run, "cell", "combine", "p", is_panoptic=True) + assert sorted(panoptic.segments()) == [("membrane", 0), ("microtubule", 2), ("microtubule", 5)] + + _invoke(split, ["-c", config_path, "-i", "cell:combine/p@10.0?panoptic=true", "--output-user-id", "back"]) + np.testing.assert_array_equal(_get(run, "microtubule", "back", "p", is_instance=True).numpy(), instances) + np.testing.assert_array_equal(_get(run, "membrane", "back", "p", is_instance=False).numpy(), membrane) diff --git a/tests/test_skeleton_graph.py b/tests/test_skeleton_graph.py new file mode 100644 index 0000000..e9d3d84 --- /dev/null +++ b/tests/test_skeleton_graph.py @@ -0,0 +1,90 @@ +"""The skeleton graph splits skeletons into branches between endpoints and junctions.""" + +import numpy as np +from copick_utils.process.skeleton_graph import ( + decompose, + merge_junctions, + prune_spurs, + pruned_skeleton, + skeleton_graph, + splice, + summarize, +) + +SHAPE = (20, 40, 40) + + +def _kinds(topology): + return sorted(topology.node_kind.values()) + + +def test_line_and_staircase_are_single_branches(): + line = np.zeros(SHAPE, bool) + line[10, 20, 2:30] = True + topology = decompose(skeleton_graph(line)) + assert len(topology.branches) == 1 and _kinds(topology) == ["endpoint", "endpoint"] + assert topology.length == 27.0 + + # A staircase has pairwise-adjacent corner voxels; their diagonals must not read as junctions. + stairs = np.zeros(SHAPE, bool) + y = x = 5 + for k in range(20): + stairs[10, y, x] = True + if k % 2: + y += 1 + else: + x += 1 + topology = decompose(skeleton_graph(stairs)) + assert len(topology.branches) == 1 and _kinds(topology) == ["endpoint", "endpoint"] + + +def test_junctions_loops_and_isolated_voxels(): + fork = np.zeros(SHAPE, bool) + fork[10, 20, 2:20] = True + for k in range(12): + fork[10, 20 + k, 20 + k] = True + fork[10, 20 - k, 20 + k] = True + topology = splice(decompose(skeleton_graph(fork))) + assert len(topology.branches) == 3 and _kinds(topology).count("junction") == 1 + + ring = np.zeros(SHAPE, bool) + theta = np.linspace(0, 2 * np.pi, 400) + ring[10, np.round(20 + 10 * np.sin(theta)).astype(int), np.round(20 + 10 * np.cos(theta)).astype(int)] = True + from skimage.morphology import skeletonize + + topology = splice(decompose(skeleton_graph(skeletonize(ring)))) + assert len(topology.branches) == 1 and topology.branches[0].closed + + dot = np.zeros(SHAPE, bool) + dot[1, 1, 1] = True + topology = decompose(skeleton_graph(dot)) + assert len(topology.branches) == 1 and _kinds(topology) == ["endpoint"] + + +def test_spurs_are_pruned_and_close_junctions_merge(): + spur = np.zeros(SHAPE, bool) + spur[10, 20, 2:30] = True + spur[10, 21:24, 15] = True + topology = prune_spurs(decompose(skeleton_graph(spur)), min_length=5) + assert len(topology.branches) == 1 and topology.length == 27.0 + assert pruned_skeleton(spur, 5).sum() == 28 + + # Two forks joined by a short bridge, as two crossing filaments often skeletonize + cross = np.zeros(SHAPE, bool) + cross[10, 20, 15:25] = True + for k in range(1, 10): + cross[10, 20 - k, 15 - k] = cross[10, 20 + k, 15 - k] = True + cross[10, 20 - k, 24 + k] = cross[10, 20 + k, 24 + k] = True + topology = splice(decompose(skeleton_graph(cross))) + assert _kinds(topology).count("junction") == 2 + topology = merge_junctions(topology, max_length=12) + assert [len(topology.incident(node)) for node in topology.junctions()] == [4] + + +def test_summary_of_a_rod(): + rod = np.zeros((11, 11, 60), bool) + rod[3:8, 3:8, 10:50] = True + summary = summarize(rod) + assert summary.n_branches == 1 and summary.n_junctions == 0 + assert 30 <= summary.length <= 40 + assert summary.radius == 3.0