Enh: add IncrementalWatershedCut for interactive seeded watershed

[lapertor]

Closes #298

Motivation

labelisation_seeded_watershed recomputes the full watershed from scratch on every seed change, which is expensive in interactive segmentation scenarios. This PR introduces IncrementalWatershedCut, a stateful class that maintains a cached Binary Partition Tree (BPT) and updates only the affected regions when seeds are added or removed.

The implementation follows the algorithm described in:

Q. Lebon, J. Lefevre, J. Cousty, B. Perret.
Interactive Segmentation With Incremental Watershed Cuts, CIARP 2023.
https://hal.science/hal-04069187v1

Algorithm

The key data structures are:

• A canonical BPT built once at construction time via bpt_canonical
• A visitCount array on BPT nodes: tracking how many seed paths pass through each internal node
• A is_cut boolean array on MST edges: an MST edge is a watershed edge when its corresponding BPT node has visitCount >= 2
• A BFS on the MST forest (skipping cut edges) to assign labels to connected components

When a seed is added, the algorithm walks up the BPT from the seed leaf, incrementing visitCount at each node; when it reaches 2, the edge is marked as a cut. When a seed is removed, the walk decrements visitCount; when it drops to 1, the edge is unmarked. A subsequent BFS re-labels the graph. The BFS correctly handles both splits (add seed) and merges (remove seed), which is why a union-find structure is not used here.

Changes

• include/higra/algo/watershed.hpp: new incremental_watershed_cut class
• higra/algo/py_watershed.cpp: pybind11 binding for IncrementalWatershedCut
• higra/algo/watershed.py: Python wrapper class IncrementalWatershedCut
• test/cpp/algo/test_watershed.cpp: 5 new C++ tests
• test/python/test_algo/test_watershed.py: 10 new Python tests in class TestIncrementalWatershed

Usage

import higra as hg
import numpy as np

graph = hg.get_4_adjacency_graph((512, 512))
edge_weights = np.random.rand(graph.num_edges())

iws = hg.IncrementalWatershedCut(graph, edge_weights)  # build BPT once

# Add seeds interactively
iws.add_seeds(np.array([0, 100]), np.array([1, 2]))
labels = iws.get_labeling()  # fast update

# Modify seeds without full recomputation
iws.add_seeds(np.array([2000]), np.array([1]))
labels = iws.get_labeling()

iws.remove_seeds(np.array([100]))
labels = iws.get_labeling()

Performance

Benchmarks on a 500x500 image (20 interactions, 5 seeds each):

Method	Mean per interaction	Total (20 interactions)
labelisation_seeded_watershed	~78 ms	~1575 ms
IncrementalWatershedCut (update only)	~5 ms	~234 ms
Speedup	~15x	~6.7x (init amortized)

On 1000x1000: ~11x speedup per interaction. The speedup is stable across
interactions regardless of the number of accumulated seeds.

Consistency

test_consistency_with_seeded_watershed verifies that IncrementalWatershedCut produces the same labeling as labelisation_seeded_watershed given identical seeds, ensuring backward compatibility.

Benchmark script

"""
Benchmark: IncrementalWatershedCut vs labelisation_seeded_watershed
in an interactive segmentation scenario.

Simulates an interactive session where seeds are added incrementally.
At each interaction, the incremental version only updates its internal state,
while the baseline recomputes the full watershed from scratch.
"""

import higra as hg
import numpy as np
import time


def benchmark_interactive(shape=(500, 500), n_interactions=20, seeds_per_interaction=5,
                          n_runs=3):
    """
    Simulate an interactive session:
    - Start with 2 initial seeds (foreground + background)
    - At each interaction, add a batch of new seeds
    - Compare the time per interaction for both methods
    """
    h, w = shape
    n_vertices = h * w
    graph = hg.get_4_adjacency_graph(shape)

    rng = np.random.RandomState(42)
    edge_weights = rng.rand(graph.num_edges()).astype(np.float64)

    # Pre-generate all seed positions and labels
    all_seed_vertices = rng.choice(n_vertices, size=2 + n_interactions * seeds_per_interaction,
                                   replace=False)
    all_seed_labels = np.ones(len(all_seed_vertices), dtype=np.int64)
    # First seed is foreground (1), second is background (2)
    all_seed_labels[0] = 1
    all_seed_labels[1] = 2
    # Remaining seeds alternate between foreground and background
    for i in range(2, len(all_seed_labels)):
        all_seed_labels[i] = 1 + (i % 2)

    print(f"Shape: {shape}")
    print(f"Num vertices: {n_vertices}")
    print(f"Num edges: {graph.num_edges()}")
    print(f"Interactions: {n_interactions}")
    print(f"Seeds per interaction: {seeds_per_interaction}")
    print(f"Runs per measurement: {n_runs}")
    print()

    # ---- Baseline: labelisation_seeded_watershed from scratch each time ----
    times_baseline = []
    for interaction in range(n_interactions):
        n_seeds = 2 + (interaction + 1) * seeds_per_interaction
        seed_v = all_seed_vertices[:n_seeds]
        seed_l = all_seed_labels[:n_seeds]

        # Build the full seeds array
        vertex_seeds = np.zeros(n_vertices, dtype=np.int64)
        vertex_seeds[seed_v] = seed_l

        elapsed = 0.0
        for _ in range(n_runs):
            t0 = time.perf_counter()
            labels_baseline = hg.labelisation_seeded_watershed(graph, edge_weights, vertex_seeds)
            t1 = time.perf_counter()
            elapsed += (t1 - t0)
        times_baseline.append(elapsed / n_runs)

    # ---- Incremental: IncrementalWatershedCut ----
    times_incremental_total = []  # init + add_seeds + get_labeling
    times_incremental_update = []  # add_seeds + get_labeling only

    t0 = time.perf_counter()
    iws = hg.IncrementalWatershedCut(graph, edge_weights)
    t_init = time.perf_counter() - t0

    # Add initial 2 seeds
    iws.add_seeds(all_seed_vertices[:2], all_seed_labels[:2])
    _ = iws.get_labeling()

    for interaction in range(n_interactions):
        start = 2 + interaction * seeds_per_interaction
        end = start + seeds_per_interaction
        new_v = all_seed_vertices[start:end]
        new_l = all_seed_labels[start:end]

        elapsed = 0.0
        for run in range(n_runs):
            if run == 0:
                # First run: actually modify state
                t0 = time.perf_counter()
                iws.add_seeds(new_v, new_l)
                labels_incr = iws.get_labeling()
                t1 = time.perf_counter()
                elapsed += (t1 - t0)
            else:
                # Subsequent runs: remove then re-add to measure consistently
                iws.remove_seeds(new_v)
                t0 = time.perf_counter()
                iws.add_seeds(new_v, new_l)
                labels_incr = iws.get_labeling()
                t1 = time.perf_counter()
                elapsed += (t1 - t0)

        avg_time = elapsed / n_runs
        times_incremental_update.append(avg_time)
        times_incremental_total.append(avg_time)

    # ---- Report ----
    print("=" * 70)
    print(f"{'Interaction':>12} | {'Baseline (ms)':>14} | {'Incremental (ms)':>17} | {'Speedup':>8}")
    print("-" * 70)
    for i in range(n_interactions):
        t_base = times_baseline[i] * 1000
        t_incr = times_incremental_update[i] * 1000
        speedup = t_base / t_incr if t_incr > 0 else float('inf')
        print(f"{i+1:>12} | {t_base:>14.2f} | {t_incr:>17.2f} | {speedup:>7.1f}x")

    print("-" * 70)
    mean_base = np.mean(times_baseline) * 1000
    mean_incr = np.mean(times_incremental_update) * 1000
    print(f"{'Mean':>12} | {mean_base:>14.2f} | {mean_incr:>17.2f} | {mean_base/mean_incr:>7.1f}x")
    print()
    print(f"Incremental init time: {t_init*1000:.2f} ms (one-time cost)")
    print(f"Baseline total time:   {sum(times_baseline)*1000:.2f} ms")
    print(f"Incremental total time (init + all updates): "
          f"{t_init*1000 + sum(times_incremental_total)*1000:.2f} ms")


def benchmark_scaling():
    """
    Show how the incremental method scales with image size.
    """
    print("\n" + "=" * 70)
    print("SCALING BENCHMARK: per-interaction time vs image size")
    print("=" * 70)

    sizes = [(100, 100), (500, 500), (1000, 1000), (2000, 2000)]
    n_runs = 3

    for shape in sizes:
        h, w = shape
        n_vertices = h * w
        graph = hg.get_4_adjacency_graph(shape)

        rng = np.random.RandomState(42)
        edge_weights = rng.rand(graph.num_edges()).astype(np.float64)

        seed_vertices = rng.choice(n_vertices, size=12, replace=False)
        seed_labels = np.array([1, 2] * 6, dtype=np.int64)

        # Baseline
        vertex_seeds = np.zeros(n_vertices, dtype=np.int64)
        vertex_seeds[seed_vertices] = seed_labels
        elapsed_base = 0.0
        for _ in range(n_runs):
            t0 = time.perf_counter()
            hg.labelisation_seeded_watershed(graph, edge_weights, vertex_seeds)
            elapsed_base += time.perf_counter() - t0
        t_base = (elapsed_base / n_runs) * 1000

        # Incremental: init + add all seeds + get_labeling
        t0 = time.perf_counter()
        iws = hg.IncrementalWatershedCut(graph, edge_weights)
        t_init = (time.perf_counter() - t0) * 1000

        # Add 10 seeds, then measure adding the last 2
        iws.add_seeds(seed_vertices[:10], seed_labels[:10])
        _ = iws.get_labeling()

        elapsed_incr = 0.0
        for run in range(n_runs):
            if run == 0:
                t0 = time.perf_counter()
                iws.add_seeds(seed_vertices[10:], seed_labels[10:])
                _ = iws.get_labeling()
                elapsed_incr += time.perf_counter() - t0
            else:
                iws.remove_seeds(seed_vertices[10:])
                t0 = time.perf_counter()
                iws.add_seeds(seed_vertices[10:], seed_labels[10:])
                _ = iws.get_labeling()
                elapsed_incr += time.perf_counter() - t0
        t_incr = (elapsed_incr / n_runs) * 1000

        speedup = t_base / t_incr if t_incr > 0 else float('inf')
        print(f"{shape[0]:>5}x{shape[1]:<5} | "
              f"Init: {t_init:>8.2f} ms | "
              f"Baseline: {t_base:>8.2f} ms | "
              f"Incremental update: {t_incr:>8.2f} ms | "
              f"Speedup: {speedup:.1f}x")


if __name__ == "__main__":
    benchmark_interactive(shape=(500, 500), n_interactions=20, seeds_per_interaction=5)
    benchmark_interactive(shape=(1000, 1000), n_interactions=20, seeds_per_interaction=5)
    benchmark_scaling()

Benchmark output (Windows, Intel Core i7-14700KF, Python 3.13.12)

Shape: (500, 500)
Num vertices: 250000
Num edges: 499000
Interactions: 20
Seeds per interaction: 5
Runs per measurement: 3

======================================================================
 Interaction |  Baseline (ms) |  Incremental (ms) |  Speedup
----------------------------------------------------------------------
           1 |          83.57 |              5.53 |    15.1x
           2 |          80.64 |              5.20 |    15.5x
           3 |          77.09 |              5.16 |    14.9x
           4 |          86.14 |              5.24 |    16.4x
           5 |          77.70 |              5.14 |    15.1x
           6 |          78.92 |              4.99 |    15.8x
           7 |          77.85 |              5.11 |    15.2x
           8 |          79.57 |              4.87 |    16.3x
           9 |          77.41 |              4.79 |    16.2x
          10 |          78.74 |              5.36 |    14.7x
          11 |          78.33 |              5.07 |    15.4x
          12 |          77.94 |              4.87 |    16.0x
          13 |          80.04 |              4.80 |    16.7x
          14 |          77.10 |              5.04 |    15.3x
          15 |          78.47 |              4.79 |    16.4x
          16 |          76.55 |              4.87 |    15.7x
          17 |          77.44 |              5.20 |    14.9x
          18 |          76.18 |              5.25 |    14.5x
          19 |          77.91 |              4.97 |    15.7x
          20 |          78.03 |              5.05 |    15.4x
----------------------------------------------------------------------
        Mean |          78.78 |              5.07 |    15.6x

Incremental init time: 132.65 ms (one-time cost)
Baseline total time:   1575.63 ms
Incremental total time (init + all updates): 233.96 ms
Shape: (1000, 1000)
Num vertices: 1000000
Num edges: 1998000
Interactions: 20
Seeds per interaction: 5
Runs per measurement: 3

======================================================================
 Interaction |  Baseline (ms) |  Incremental (ms) |  Speedup
----------------------------------------------------------------------
           1 |         387.37 |             35.89 |    10.8x
           2 |         399.85 |             36.06 |    11.1x
           3 |         398.65 |             35.02 |    11.4x
           4 |         396.37 |             36.10 |    11.0x
           5 |         389.45 |             37.91 |    10.3x
           6 |         388.10 |             36.56 |    10.6x
           7 |         383.60 |             35.09 |    10.9x
           8 |         385.06 |             35.98 |    10.7x
           9 |         385.41 |             35.53 |    10.8x
          10 |         409.90 |             37.44 |    10.9x
          11 |         388.49 |             34.13 |    11.4x
          12 |         394.39 |             35.33 |    11.2x
          13 |         403.05 |             34.53 |    11.7x
          14 |         389.26 |             34.39 |    11.3x
          15 |         390.24 |             34.67 |    11.3x
          16 |         390.97 |             33.66 |    11.6x
          17 |         390.84 |             36.26 |    10.8x
          18 |         387.86 |             34.14 |    11.4x
          19 |         386.30 |             35.78 |    10.8x
          20 |         390.45 |             35.16 |    11.1x
----------------------------------------------------------------------
        Mean |         391.78 |             35.48 |    11.0x

Incremental init time: 734.49 ms (one-time cost)
Baseline total time:   7835.60 ms
Incremental total time (init + all updates): 1444.11 ms

======================================================================
SCALING BENCHMARK: per-interaction time vs image size
======================================================================
  100x100   | Init:     3.43 ms | Baseline:     2.51 ms | Incremental update:     0.18 ms | Speedup: 14.3x
  500x500   | Init:   141.11 ms | Baseline:    77.33 ms | Incremental update:     5.49 ms | Speedup: 14.1x
 1000x1000  | Init:   725.87 ms | Baseline:   399.25 ms | Incremental update:    35.18 ms | Speedup: 11.3x
 2000x2000  | Init:  3898.48 ms | Baseline:  2001.76 ms | Incremental update:   179.41 ms | Speedup: 11.2x

[PerretB]

Same here

[PerretB]

This is sub-optimal. When a seed is removed we also know the edge of the mst that is reactivated which can be used to determine the label of the "other side" (except if the other side is also removed in the same batch, in which case relabelling will start somewhere else).

[lapertor]

Pushed a final revision covering all review points, plus a pre-validation fix discovered during testing.

1. add_seeds: two-pass batch handling

Split into Pass 1 (walk up BPT, increment visit_count, mark is_cut on the 0→1→2 transitions) and Pass 2 (one BFS per seed to relabel its component). Pass 2 runs only after the entire batch has updated the cut state, so each seed's BFS operates on the final cut state of the full batch.

2. remove_seeds: two-pass, leverages BPT info

Rewritten to match Algorithm 2 of Lebon et al.:

• Pass 1: for each removed seed, walk up the BPT decrementing visit_count and break on the first 2→1 transition. Record the corresponding MST edge as a de-cut. If the walk-up reaches the root without a 2→1 transition, the seed is "lone" (its component had no other seed). m_is_cut is not modified yet.
• Pass 2a: process each de-cut in batch insertion order. For the MST edge being reactivated, find the surviving label on each side via a bounded BFS (component_seed_label, which respects the current cut state). The side with no surviving seed is repainted with the label from the other side; if both sides are dead, the component goes to background. The BFS runs before flipping m_is_cut, so each side scan stays within its component.
• Pass 2b: lone seeds get their component repainted to background.

3. Batch validation fix

add_seeds and remove_seeds previously mutated m_seed_labels / m_is_cut before validating the entire batch. If an element at position i>0 was invalid (duplicate vertex, out-of-range, seed not found), the exception was raised after elements 0..i-1 had already mutated the object, leaving it in an inconsistent state.
Added Pass 0 pre-validation before Pass 1 in both functions: the entire batch is validated (range, intra-batch duplicates via std::unordered_set<index_t>, seed existence) before any mutation occurs. All checks use hg_assert (no-op in Release, consistent with existing code). No changes to Pass 1/2/2a/2b loop bodies.

4. Documentation

• C++ class: Doxygen Example: block.
• Sphinx: IncrementalWatershedCut added to doc/source/python/watershed.rst (autosummary + autoclass with :members:).
• Added :Complexity: to the Python docstring and C++ Doxygen block. Covers BPT construction, add_seeds, remove_seeds, and get_labeling, based on the current BFS-based implementation.
• Documented the no-relabeling restriction in add_seeds: re-adding a seed on an existing vertex raises an assertion; the documented pattern is remove_seeds(v); add_seeds(v, new_label).

5. Tests

New regression tests added to both Python and C++ suites:

• batch_remove_equals_sequential: batch remove of two seeds should match removing them one at a time.
• batch_remove_both_sides_of_edge: removing the two seeds that sit on each side of the only cut should yield an all-background labeling.
• interactive_churn: alternating add/remove with label changes on a 100×100 grid, comparing the incremental labeling against labelisation_seeded_watershed after every step.
• decuts_ancestor_descendant: 1×8 path with balanced BPT; batch remove_seeds([0, 3]) produces two de-cuts where one BPT node is a descendant of the other.
• add_seed_already_exists_raises / add_seed_label_zero_raises: assert on duplicate seed and zero label.
• add_seeds_batch_duplicate_no_partial_mutation / add_seeds_cross_batch_duplicate_no_partial_mutation: verify invalid add batches do not partially mutate state.
• remove_seeds_batch_duplicate_no_partial_mutation / remove_seeds_cross_batch_duplicate_no_partial_mutation: same for remove batches.
  Suite is now at 26/26 Python tests, 15 C++ test cases (83 assertions).

6. Misc

• component_seed_label uses a generation counter on m_visited so its BFS allocates nothing per call.
• hg_assert added to the unreachable branch of Pass 2a (both sides of a de-cut having different live labels is impossible given the 2→1 invariant).
• CRLF that had slipped into test_watershed.py is now fixed; the file is LF in all commits, and git diff master..HEAD --check is clean.
• All commits carry Signed-off-by.

Let me know if there is something to improve, or if I did something the wrong way, or anything like this!