diff --git a/lib/cluster/tests/conftest.py b/lib/cluster/tests/conftest.py
new file mode 100644
index 00000000000..0626a9d0f0f
--- /dev/null
+++ b/lib/cluster/tests/conftest.py
@@ -0,0 +1,32 @@
+import pytest
+
+import grass.lib.gis as libgis
+
+
+@pytest.fixture
+def libgis_without_session():
+    """Let the clustering functions run without a GRASS session
+
+    Most of the I_cluster_*() functions call G_debug(), which reads the
+    DEBUG variable out of the session's GISRC file. With no session there
+    is no file to read, so the lookup ends in G_fatal_error() and exits
+    the process, taking the pytest run down with it. Clustering itself
+    never touches a project or mapset, so instead of setting up a session
+    for its own sake, GISRC is switched to memory mode: libgis then keeps
+    those variables in memory and never looks for the file. Tools that
+    legitimately run outside a mapset do the same, see
+    general/g.message/main.c and general/g.dirseps/main.c.
+
+    Giving these tests a real session is not an equivalent fix. Reading a
+    GISRC file makes libgis cache that session's project for the rest of
+    the process, so later tests sharing the process would go looking for
+    their own maps in a project that no longer exists.
+
+    The mode is process-global in the C library, so it is restored here.
+    """
+    previous_mode = libgis.G_get_gisrc_mode()
+    libgis.G_set_gisrc_mode(libgis.G_GISRC_MODE_MEMORY)
+    try:
+        yield
+    finally:
+        libgis.G_set_gisrc_mode(previous_mode)
diff --git a/lib/cluster/tests/lib_cluster_algorithm_ctypes_test.py b/lib/cluster/tests/lib_cluster_algorithm_ctypes_test.py
new file mode 100644
index 00000000000..99e890347bc
--- /dev/null
+++ b/lib/cluster/tests/lib_cluster_algorithm_ctypes_test.py
@@ -0,0 +1,261 @@
+"""Tests for the clustering algorithm step ctypes bindings
+
+Most of these functions call G_debug() internally, which needs an
+answer for the session's DEBUG variable even though none of them read or
+write a project or mapset; see conftest.py's libgis_without_session
+fixture for how that is arranged.
+
+The sample data is deliberately two-band rather than one-band: every one
+of these functions loops over bands, and with a single band a loop that
+only ever handled band 0 would still produce correct results.
+"""
+
+from ctypes import byref, c_double, c_int
+
+import pytest
+
+import grass.lib.cluster as libcluster
+
+# Five points whose second band is ten times the first, so that values
+# for the two bands stay easy to check by hand but remain distinguishable
+# from each other.
+SAMPLE_POINTS = [(1.0, 10.0), (2.0, 20.0), (3.0, 30.0), (4.0, 40.0), (5.0, 50.0)]
+
+
+def cluster_with_points(points, nclasses, nbands=2):
+    """A Cluster struct loaded with points and allocated for nclasses
+
+    Mirrors the first few steps of I_cluster_exec(): collect the points,
+    then allocate the per-class arrays that the rest of the algorithm's
+    steps read and write.
+    """
+    cluster = libcluster.Cluster()
+    libcluster.I_cluster_begin(byref(cluster), nbands)
+    for values in points:
+        assert (
+            libcluster.I_cluster_point(byref(cluster), (c_double * nbands)(*values))
+            == 0
+        )
+    cluster.nclasses = nclasses
+    assert libcluster.I_cluster_exec_allocate(byref(cluster)) == 1
+    return cluster
+
+
+def classify(cluster):
+    """Run I_cluster_means() then I_cluster_assign(), as I_cluster_exec()
+    does for its first pass, and return the resulting per-point classes"""
+    libcluster.I_cluster_means(byref(cluster))
+    interrupted = c_int(0)
+    libcluster.I_cluster_assign(byref(cluster), byref(interrupted))
+    class_field = getattr(cluster, "class")
+    return [class_field[p] for p in range(cluster.npoints)]
+
+
+def per_class(array, band, nclasses):
+    """Read one band's row out of a band-by-class array"""
+    return [array[band][c] for c in range(nclasses)]
+
+
+def test_means_spreads_initial_class_means_by_standard_deviation(
+    libgis_without_session,
+) -> None:
+    """I_cluster_means() has no classified points to average yet, so it
+    seeds each class mean at an even spread around the overall mean of
+    each band, scaled by that band's standard deviation"""
+    cluster = cluster_with_points(SAMPLE_POINTS, 3)
+
+    libcluster.I_cluster_means(byref(cluster))
+
+    assert per_class(cluster.mean, 0, 3) == pytest.approx(
+        [1.4188611699158102, 3.0, 4.58113883008419]
+    )
+    assert per_class(cluster.mean, 1, 3) == pytest.approx(
+        [14.188611699158104, 30.0, 45.811388300841898]
+    )
+
+
+def test_assign_puts_each_point_in_its_nearest_class(libgis_without_session) -> None:
+    cluster = cluster_with_points(SAMPLE_POINTS, 3)
+
+    classes = classify(cluster)
+
+    assert classes == [0, 0, 1, 2, 2]
+    assert [cluster.count[c] for c in range(3)] == [2, 1, 2]
+    assert per_class(cluster.sum, 0, 3) == [3.0, 3.0, 9.0]
+    assert per_class(cluster.sum, 1, 3) == [30.0, 30.0, 90.0]
+
+
+def test_reassign_leaves_already_nearest_points_alone(libgis_without_session) -> None:
+    """I_cluster_reassign() is the iterated form of I_cluster_assign():
+    it compares each point against the class sums rather than the seeded
+    means, and reports how many points changed class. Points that are
+    already in their nearest class stay put and are not counted.
+    """
+    cluster = cluster_with_points(SAMPLE_POINTS, 3)
+    classify(cluster)
+    interrupted = c_int(0)
+
+    changes = libcluster.I_cluster_reassign(byref(cluster), byref(interrupted))
+
+    assert changes == 0
+    class_field = getattr(cluster, "class")
+    assert [class_field[p] for p in range(cluster.npoints)] == [0, 0, 1, 2, 2]
+    assert [cluster.count[c] for c in range(3)] == [2, 1, 2]
+
+
+def test_reassign_moves_a_point_to_a_nearer_class(libgis_without_session) -> None:
+    """A point is moved when another class centroid is closer across all
+    bands together, and the counts and sums of both classes are updated
+
+    The single point (1, 9) sits nearer to class 0 on band 0 alone
+    (distance 1 against 4), but nearer to class 1 once band 1 is included
+    (5 against 82), so this also pins down that the comparison really
+    spans every band rather than just the first.
+    """
+    cluster = cluster_with_points([(1.0, 9.0)], 2)
+    # Hand-set the two class centroids: class 0 at (0, 0) and class 1 at
+    # (3, 10), each standing in for 2 points, since I_cluster_reassign()
+    # reads the class sums and counts rather than the means.
+    cluster.count[0] = 2
+    cluster.sum[0][0] = 0.0
+    cluster.sum[1][0] = 0.0
+    cluster.count[1] = 2
+    cluster.sum[0][1] = 6.0
+    cluster.sum[1][1] = 20.0
+    class_field = getattr(cluster, "class")
+    class_field[0] = 0
+    interrupted = c_int(0)
+
+    changes = libcluster.I_cluster_reassign(byref(cluster), byref(interrupted))
+
+    assert changes == 1
+    assert class_field[0] == 1
+    assert [cluster.count[c] for c in range(2)] == [1, 3]
+    assert per_class(cluster.sum, 0, 2) == [-1.0, 7.0]
+    assert per_class(cluster.sum, 1, 2) == [-9.0, 29.0]
+
+
+def test_sum2_accumulates_the_sum_of_squares_per_class(libgis_without_session) -> None:
+    """Unlike GASTATS's variance-style statistics, I_cluster_sum2() stores
+    the plain sum of squares per class; I_cluster_separation() derives
+    variance from it itself"""
+    cluster = cluster_with_points(SAMPLE_POINTS, 3)
+    classify(cluster)
+
+    libcluster.I_cluster_sum2(byref(cluster))
+
+    assert per_class(cluster.sum2, 0, 3) == [5.0, 9.0, 41.0]
+    assert per_class(cluster.sum2, 1, 3) == [500.0, 900.0, 4100.0]
+
+
+def test_separation_returns_the_far_sentinel_for_a_tiny_class(
+    libgis_without_session,
+) -> None:
+    """A class with fewer than 2 points has no variance to compare, so
+    I_cluster_separation() reports it as maximally separated (-1.0)
+    rather than dividing by zero"""
+    cluster = cluster_with_points(SAMPLE_POINTS, 3)
+    classify(cluster)
+    libcluster.I_cluster_sum2(byref(cluster))
+
+    # Class 1 holds only the single point (3, 30), see the assign test.
+    assert libcluster.I_cluster_separation(byref(cluster), 0, 1) == -1.0
+
+
+def test_separation_measures_how_distinguishable_two_classes_are(
+    libgis_without_session,
+) -> None:
+    cluster = cluster_with_points(SAMPLE_POINTS, 3)
+    classify(cluster)
+    libcluster.I_cluster_sum2(byref(cluster))
+
+    # Classes 0 and 2 hold two points each, and the separation sums the
+    # contribution of both bands.
+    separation = libcluster.I_cluster_separation(byref(cluster), 0, 2)
+
+    assert separation == pytest.approx(1.2247448713915889)
+
+
+def test_distinct_finds_the_closest_pair_below_the_separation_threshold(
+    libgis_without_session,
+) -> None:
+    cluster = cluster_with_points(SAMPLE_POINTS, 3)
+    classify(cluster)
+    libcluster.I_cluster_sum2(byref(cluster))
+
+    # A generous threshold finds classes 0 and 2 as the closest pair with
+    # at least 2 points each (class 1 is skipped: it only has one point).
+    assert libcluster.I_cluster_distinct(byref(cluster), 100.0) == 0
+    assert (cluster.merge1, cluster.merge2) == (0, 2)
+
+    # A threshold tighter than any actual separation finds nothing to merge.
+    assert libcluster.I_cluster_distinct(byref(cluster), 0.0001) == 1
+
+
+def test_nclasses_counts_classes_at_or_above_a_minimum_size(
+    libgis_without_session,
+) -> None:
+    cluster = cluster_with_points(SAMPLE_POINTS, 3)
+    classify(cluster)
+
+    # Classes 0 and 2 have 2 points each; class 1 has only 1.
+    assert libcluster.I_cluster_nclasses(byref(cluster), 2) == 2
+
+
+def test_merge_combines_two_classes(libgis_without_session) -> None:
+    cluster = cluster_with_points(SAMPLE_POINTS, 3)
+    classify(cluster)
+    libcluster.I_cluster_sum2(byref(cluster))
+    libcluster.I_cluster_distinct(byref(cluster), 100.0)  # sets merge1/merge2 to 0, 2
+
+    libcluster.I_cluster_merge(byref(cluster))
+
+    class_field = getattr(cluster, "class")
+    # Every point that was in class 2 now reports class 0; class 2's count
+    # and per-band sums move into class 0 and are zeroed at their old slot.
+    assert [class_field[p] for p in range(cluster.npoints)] == [0, 0, 1, 0, 0]
+    assert [cluster.count[c] for c in range(3)] == [4, 1, 0]
+    assert per_class(cluster.sum, 0, 3) == [12.0, 3.0, 0.0]
+    assert per_class(cluster.sum, 1, 3) == [120.0, 30.0, 0.0]
+
+
+def test_reclass_compacts_out_a_class_below_the_minimum_size(
+    libgis_without_session,
+) -> None:
+    """I_cluster_reclass() removes classes below minsize and shifts the
+    remaining ones down to fill the resulting gaps, remapping every
+    point's class along the way"""
+    cluster = cluster_with_points([(0.0, 0.0)] * 3, 3)
+    # Hand-set an empty class 1 between two populated ones, rather than
+    # reaching this state through classification, to isolate reclass()
+    # from the rest of the algorithm.
+    for class_index, (count, band0, band1) in enumerate(
+        [(5, 50.0, 500.0), (0, 0.0, 0.0), (3, 30.0, 300.0)]
+    ):
+        cluster.count[class_index] = count
+        cluster.sum[0][class_index] = band0
+        cluster.sum[1][class_index] = band1
+    class_field = getattr(cluster, "class")
+    class_field[0] = 0
+    class_field[1] = 0
+    class_field[2] = 2
+
+    ret = libcluster.I_cluster_reclass(byref(cluster), 1)
+
+    assert ret == 0
+    assert cluster.nclasses == 2
+    assert [class_field[p] for p in range(3)] == [0, 0, 1]
+    assert [cluster.count[c] for c in range(2)] == [5, 3]
+    assert per_class(cluster.sum, 0, 2) == [50.0, 30.0]
+    assert per_class(cluster.sum, 1, 2) == [500.0, 300.0]
+
+
+def test_reclass_returns_1_when_every_class_already_meets_the_minimum(
+    libgis_without_session,
+) -> None:
+    cluster = cluster_with_points([(0.0, 0.0)], 2)
+    cluster.count[0] = 5
+    cluster.count[1] = 3
+
+    assert libcluster.I_cluster_reclass(byref(cluster), 1) == 1
+    assert cluster.nclasses == 2
diff --git a/lib/cluster/tests/lib_cluster_exec_ctypes_test.py b/lib/cluster/tests/lib_cluster_exec_ctypes_test.py
new file mode 100644
index 00000000000..9de83ae3f4b
--- /dev/null
+++ b/lib/cluster/tests/lib_cluster_exec_ctypes_test.py
@@ -0,0 +1,369 @@
+"""Tests for the full clustering workflow ctypes binding
+
+I_cluster_exec() drives the whole algorithm end to end: initial means,
+assignment, iterative re-assignment, merging classes that are not
+distinct enough, and finally building the output signatures. It calls
+G_debug() internally, so it needs the same libgis_without_session
+fixture as lib_cluster_algorithm_ctypes_test.py; see conftest.py for why.
+
+Before any of that, the function quietly replaces out-of-range arguments
+with defaults of its own, and most of the tests below pin those down.
+None of the corrected values can be read back off the struct, so each is
+checked by the effect it has on the result, alongside a contrasting
+in-range value that produces a different one.
+
+The sample data is two-band for the reason given in
+lib_cluster_algorithm_ctypes_test.py: a single band cannot tell a
+correct per-band loop apart from one that only ever handles band 0.
+"""
+
+from contextlib import contextmanager
+from ctypes import CFUNCTYPE, byref, c_double, c_int
+
+import pytest
+
+import grass.lib.cluster as libcluster
+import grass.lib.gis as libgis
+
+# Two groups far enough apart that the clustering result is unambiguous,
+# with the second band on a different scale from the first.
+GROUP_A = [(0.0, 0.0), (1.0, 1.0), (2.0, 2.0), (-1.0, -1.0), (0.5, 0.5)]
+GROUP_B = [
+    (100.0, 200.0),
+    (101.0, 201.0),
+    (99.0, 199.0),
+    (100.5, 200.5),
+    (99.5, 199.5),
+]
+
+# One tight group plus a single distant point. The outlier ends up alone
+# in a class of its own, which is what the final minimum-size pass then
+# has to decide about.
+GROUP_WITH_OUTLIER = [
+    (0.0, 0.0),
+    (1.0, 1.0),
+    (2.0, 2.0),
+    (3.0, 3.0),
+    (4.0, 4.0),
+    (100.0, 100.0),
+]
+
+# A continuum rather than distinct groups. Where the class boundaries
+# fall is arbitrary enough that the first re-assignment pass still moves
+# a point, leaving the run short of a 98% stable result on that pass.
+EVENLY_SPACED = [(float(i), float(i) * 2.0) for i in range(12)]
+
+# Eight of the nine points of a 3x3 grid. Classes carved out of it
+# overlap heavily, so the separation between them comes to about 0.35:
+# below the 0.5 default, but above a smaller threshold passed explicitly.
+UNIT_GRID = [
+    (0.0, 0.0),
+    (1.0, 0.0),
+    (2.0, 0.0),
+    (0.0, 1.0),
+    (1.0, 1.0),
+    (2.0, 1.0),
+    (0.0, 2.0),
+    (1.0, 2.0),
+]
+
+# I_cluster_exec() reports progress through this callback, which is also
+# the only way to reach into the middle of a run. Phase 1 is the initial
+# means, phase 2 the initial assignment, phase 3 each re-assignment pass,
+# and phase 4 a merge of two classes that were not distinct enough.
+CHECKPOINT = libcluster.I_cluster_exec.argtypes[6]
+
+# The default handler for G_warning() looks up GISBASE to decide where to
+# log the message, which is a fatal error when it is not set, and on some
+# platforms (observed: glibc, not musl) that kills the process silently
+# rather than just losing the message. Routing messages to a handler of
+# our own avoids that lookup entirely, following the same
+# G_set_error_routine() pattern used in gui/wxpython/vdigit/wxdisplay.py
+# and gui/wxpython/nviz/wxnviz.py, and in lib/arraystats' tests.
+ERROR_ROUTINE = CFUNCTYPE(libgis.UNCHECKED(c_int), libgis.String, c_int)
+
+
+@contextmanager
+def swallow_grass_messages():
+    @ERROR_ROUTINE
+    def handler(message, is_fatal):
+        return 1  # non-zero: caller (this test) has handled it
+
+    libgis.G_set_error_routine(handler)
+    try:
+        yield
+    finally:
+        libgis.G_unset_error_routine()
+
+
+def cluster_with_points(points, nbands=2):
+    cluster = libcluster.Cluster()
+    libcluster.I_cluster_begin(byref(cluster), nbands)
+    for values in points:
+        assert (
+            libcluster.I_cluster_point(byref(cluster), (c_double * nbands)(*values))
+            == 0
+        )
+    return cluster
+
+
+def exec_with(
+    points,
+    maxclass=2,
+    iterations=20,
+    convergence=98.0,
+    separation=0.5,
+    min_class_size=2,
+    interrupt_at_phase=None,
+):
+    """Cluster points with I_cluster_exec(), recording the phases it reports
+
+    Every parameter defaults to an in-range value, so a test only has to
+    name the one it is about. Passing interrupt_at_phase raises the
+    interrupt flag from inside the checkpoint callback, which is the only
+    point at which a caller could realistically raise it mid-run.
+
+    Returns the return value, the Cluster struct and the phases seen.
+    """
+    cluster = cluster_with_points(points)
+    phases = []
+    interrupted = c_int(0)
+
+    @CHECKPOINT
+    def checkpoint(_cluster, phase):
+        phases.append(phase)
+        if phase == interrupt_at_phase:
+            interrupted.value = 1
+        return 0
+
+    ret = libcluster.I_cluster_exec(
+        byref(cluster),
+        maxclass,
+        iterations,
+        convergence,
+        separation,
+        min_class_size,
+        checkpoint,
+        byref(interrupted),
+    )
+    return ret, cluster, phases
+
+
+def classes_of(cluster):
+    """Each point's class, with -1 for points whose class was eliminated"""
+    class_field = getattr(cluster, "class")
+    return [class_field[p] for p in range(cluster.npoints)]
+
+
+def counts_of(cluster):
+    return [cluster.count[c] for c in range(cluster.nclasses)]
+
+
+def band_mean(points, band):
+    return sum(point[band] for point in points) / len(points)
+
+
+def test_exec_separates_two_well_separated_groups(libgis_without_session) -> None:
+    """The checkpoint callback lets a caller (e.g. i.cluster's progress
+    reporting) observe each phase as it runs; it is exercised here as
+    well as the final result"""
+    ret, cluster, phases = exec_with(GROUP_A + GROUP_B)
+
+    assert ret == 0
+    assert cluster.nclasses == 2
+    classes = classes_of(cluster)
+    # The two groups are far enough apart that every point in GROUP_A must
+    # land in one class and every point in GROUP_B in the other, whichever
+    # of the two class numbers the algorithm happens to assign to each.
+    assert len(set(classes[: len(GROUP_A)])) == 1
+    assert len(set(classes[len(GROUP_A) :])) == 1
+    assert classes[0] != classes[-1]
+    # No phase 4, since the two groups are already distinct on the first
+    # pass and so nothing is ever merged.
+    assert phases == [1, 2, 3]
+
+    # The signature of each class is the per-band mean of its own group.
+    assert cluster.S.nsigs == 2
+    signatures = sorted(
+        (cluster.S.sig[c].mean[0], cluster.S.sig[c].mean[1]) for c in range(2)
+    )
+    assert signatures == pytest.approx(
+        [
+            (band_mean(GROUP_A, 0), band_mean(GROUP_A, 1)),
+            (band_mean(GROUP_B, 0), band_mean(GROUP_B, 1)),
+        ]
+    )
+    for c in range(2):
+        assert cluster.S.sig[c].npoints == len(GROUP_A)
+        assert cluster.S.sig[c].status == 1
+
+
+def test_exec_rejects_too_few_points(libgis_without_session) -> None:
+    """I_cluster_exec() checks the point count before doing anything else,
+    so a null checkpoint callback is safe here: it is never reached
+
+    This is the one case in this file that reaches I_cluster_exec()'s
+    G_warning() call, so it needs swallow_grass_messages(); see that
+    function's comment for why.
+    """
+    cluster = cluster_with_points([(1.0, 1.0)])
+    interrupted = c_int(0)
+    no_checkpoint = CHECKPOINT()  # a null function pointer
+
+    with swallow_grass_messages():
+        ret = libcluster.I_cluster_exec(
+            byref(cluster), 2, 20, 98.0, 0.5, 2, no_checkpoint, byref(interrupted)
+        )
+
+    assert ret == 1
+
+
+def test_exec_turns_a_negative_maxclass_into_a_single_class(
+    libgis_without_session,
+) -> None:
+    """A negative class count becomes 1 rather than being rejected, so
+    the two groups come back merged into one class"""
+    ret, cluster, _ = exec_with(GROUP_A + GROUP_B, maxclass=-1)
+
+    assert ret == 0
+    assert cluster.nclasses == 1
+    assert classes_of(cluster) == [0] * len(GROUP_A + GROUP_B)
+    assert counts_of(cluster) == [len(GROUP_A + GROUP_B)]
+    assert cluster.S.nsigs == 1
+
+
+def test_exec_raises_a_non_positive_min_class_size_to_seventeen(
+    libgis_without_session,
+) -> None:
+    """A minimum class size of 0 is replaced by 17, not treated as "no
+    minimum"
+
+    Both groups here hold five points, so the substituted minimum wipes
+    out every class and leaves no signatures at all, while a minimum of 5
+    keeps both. That difference is the only way to see which value was
+    actually used, since the corrected one is never stored.
+    """
+    ret, cluster, _ = exec_with(GROUP_A + GROUP_B, min_class_size=0)
+
+    assert ret == 0
+    assert cluster.nclasses == 0
+    assert classes_of(cluster) == [-1] * len(GROUP_A + GROUP_B)
+    assert cluster.S.nsigs == 0
+
+    _, kept, _ = exec_with(GROUP_A + GROUP_B, min_class_size=5)
+    assert kept.nclasses == 2
+    assert counts_of(kept) == [5, 5]
+
+
+def test_exec_raises_a_min_class_size_of_one_to_two(libgis_without_session) -> None:
+    """A minimum class size of 1 is raised to 2, so a class holding a
+    single point is still discarded
+
+    Taken literally, a minimum of 1 would keep the outlier's one-point
+    class; it is dropped instead, which is what shows the value was
+    raised.
+    """
+    ret, cluster, _ = exec_with(GROUP_WITH_OUTLIER, min_class_size=1)
+
+    assert ret == 0
+    assert cluster.nclasses == 1
+    assert counts_of(cluster) == [5]
+    assert classes_of(cluster) == [0, 0, 0, 0, 0, -1]
+
+
+def test_exec_replaces_non_positive_iterations_with_twenty(
+    libgis_without_session,
+) -> None:
+    """Asking for no iterations runs twenty of them rather than none
+
+    The convergence threshold is set beyond 100% so that it can never be
+    met, which leaves the iteration limit as the only thing that can end
+    the loop and makes the substituted value countable.
+    """
+    ret, cluster, phases = exec_with(GROUP_A + GROUP_B, iterations=0, convergence=200.0)
+
+    assert ret == 0
+    assert phases.count(3) == 20
+    assert cluster.iteration == 20
+
+    _, capped, capped_phases = exec_with(
+        GROUP_A + GROUP_B, iterations=3, convergence=200.0
+    )
+    assert capped_phases.count(3) == 3
+    assert capped.iteration == 3
+
+
+def test_exec_replaces_a_non_positive_convergence_with_ninety_eight(
+    libgis_without_session,
+) -> None:
+    """A convergence of 0 becomes 98%, so the run keeps iterating rather
+    than stopping as soon as any point at all is stable
+
+    The first re-assignment pass over this data leaves it short of 98%
+    stable, so the substituted threshold forces a second pass, while a
+    threshold the first pass already clears stops after one.
+    """
+    ret, cluster, phases = exec_with(EVENLY_SPACED, maxclass=4, convergence=0.0)
+
+    assert ret == 0
+    assert phases.count(3) == 2
+    assert cluster.iteration == 2
+
+    _, explicit, explicit_phases = exec_with(
+        EVENLY_SPACED, maxclass=4, convergence=98.0
+    )
+    assert explicit_phases.count(3) == 2
+    assert explicit.iteration == 2
+
+    _, lenient, lenient_phases = exec_with(EVENLY_SPACED, maxclass=4, convergence=50.0)
+    assert lenient_phases.count(3) == 1
+    assert lenient.iteration == 1
+
+
+def test_exec_replaces_a_negative_separation_with_half(
+    libgis_without_session,
+) -> None:
+    """A negative separation becomes 0.5, which is enough to merge the
+    two classes this grid is split into
+
+    The merge is also the only case in this file that reaches phase 4.
+    The two classes are about 0.35 apart, so the substituted 0.5 merges
+    them exactly as an explicit 0.5 does, while a smaller threshold
+    passed on purpose leaves them alone.
+    """
+    ret, cluster, phases = exec_with(UNIT_GRID, separation=-1.0)
+
+    assert ret == 0
+    assert phases == [1, 2, 3, 4, 3]
+    assert cluster.nclasses == 1
+
+    _, explicit, explicit_phases = exec_with(UNIT_GRID, separation=0.5)
+    assert explicit_phases == [1, 2, 3, 4, 3]
+    assert explicit.nclasses == 1
+
+    _, distinct, distinct_phases = exec_with(UNIT_GRID, separation=0.1)
+    assert distinct_phases == [1, 2, 3]
+    assert distinct.nclasses == 2
+
+
+def test_exec_reports_an_interrupt_during_the_initial_assignment(
+    libgis_without_session,
+) -> None:
+    """Raising the flag from the phase 1 checkpoint stops the run inside
+    I_cluster_assign(), before any point has been assigned"""
+    ret, _, phases = exec_with(GROUP_A + GROUP_B, interrupt_at_phase=1)
+
+    assert ret == -2
+    assert phases == [1]
+
+
+def test_exec_reports_an_interrupt_before_the_reassignment_loop(
+    libgis_without_session,
+) -> None:
+    """Raising the flag from the phase 2 checkpoint stops the run at the
+    top of the re-assignment loop, which is a separate check from the one
+    the initial assignment goes through"""
+    ret, _, phases = exec_with(GROUP_A + GROUP_B, interrupt_at_phase=2)
+
+    assert ret == -2
+    assert phases == [1, 2]
diff --git a/lib/cluster/tests/lib_cluster_points_ctypes_test.py b/lib/cluster/tests/lib_cluster_points_ctypes_test.py
new file mode 100644
index 00000000000..2b5e0dea22a
--- /dev/null
+++ b/lib/cluster/tests/lib_cluster_points_ctypes_test.py
@@ -0,0 +1,110 @@
+"""Tests for the point-collection ctypes bindings
+
+I_cluster_begin(), I_cluster_point() and the point-set variants only
+allocate and fill the Cluster struct's own arrays; none of them call
+G_debug(), so unlike most of the clustering algorithm itself they need
+nothing from libgis' session environment and no fixture (see
+conftest.py and lib_cluster_algorithm_ctypes_test.py).
+"""
+
+from ctypes import byref, c_double
+
+import grass.lib.cluster as libcluster
+import grass.lib.raster as libraster
+
+
+def new_cluster(nbands):
+    """A freshly initialized Cluster struct for the given band count"""
+    cluster = libcluster.Cluster()
+    assert libcluster.I_cluster_begin(byref(cluster), nbands) == 0
+    return cluster
+
+
+def point(*values):
+    """Build a ctypes double array from per-band values"""
+    return (c_double * len(values))(*values)
+
+
+def null_point(nbands):
+    """A point whose bands are all set to the DCELL NULL value"""
+    values = point(*([0.0] * nbands))
+    libraster.Rast_set_d_null_value(values, nbands)
+    return values
+
+
+def test_begin_rejects_a_non_positive_band_count() -> None:
+    cluster = libcluster.Cluster()
+    assert libcluster.I_cluster_begin(byref(cluster), 0) == 1
+
+
+def test_point_accumulates_values_and_running_sums() -> None:
+    cluster = new_cluster(2)
+
+    assert libcluster.I_cluster_point(byref(cluster), point(1.0, 2.0)) == 0
+    assert libcluster.I_cluster_point(byref(cluster), point(3.0, 4.0)) == 0
+
+    assert cluster.npoints == 2
+    assert [cluster.points[0][p] for p in range(2)] == [1.0, 3.0]
+    assert [cluster.points[1][p] for p in range(2)] == [2.0, 4.0]
+    # band_sum/band_sum2 are running totals kept for I_cluster_means(), not
+    # recomputed from the points array, so they are checked here too.
+    assert [cluster.band_sum[b] for b in range(2)] == [4.0, 6.0]
+    assert [cluster.band_sum2[b] for b in range(2)] == [10.0, 20.0]
+
+
+def test_point_rejects_a_point_with_a_null_band() -> None:
+    """A point with any band set to NULL is rejected outright, not just
+    skipped in that one band, since a partial pixel cannot be clustered"""
+    cluster = new_cluster(2)
+
+    assert libcluster.I_cluster_point(byref(cluster), null_point(2)) == 1
+    assert cluster.npoints == 0
+
+
+def test_point_set_batch_workflow_drops_all_zero_points() -> None:
+    """I_cluster_begin_point_set()/I_cluster_point_part()/
+    I_cluster_end_point_set() is the batch equivalent of repeated
+    I_cluster_point() calls, used when bands are filled one at a time
+    rather than one point at a time. A point whose every band is exactly
+    zero (as opposed to NULL) is silently dropped, matching the "point
+    contains no data" convention used elsewhere in this library.
+    """
+    cluster = new_cluster(2)
+    values = [(1.0, 10.0), (0.0, 0.0), (2.0, 20.0)]
+
+    assert libcluster.I_cluster_begin_point_set(byref(cluster), len(values)) == 0
+    for n, bands in enumerate(values):
+        for band, value in enumerate(bands):
+            assert (
+                libcluster.I_cluster_point_part(
+                    byref(cluster), c_double(value), band, n
+                )
+                == 0
+            )
+    kept = libcluster.I_cluster_end_point_set(byref(cluster), len(values))
+
+    assert kept == 2
+    assert cluster.npoints == 2
+    assert [cluster.points[0][p] for p in range(2)] == [1.0, 2.0]
+    assert [cluster.points[1][p] for p in range(2)] == [10.0, 20.0]
+
+
+def test_point_part_rejects_a_null_value() -> None:
+    cluster = new_cluster(1)
+    libcluster.I_cluster_begin_point_set(byref(cluster), 1)
+
+    null_value = c_double()
+    libraster.Rast_set_d_null_value(byref(null_value), 1)
+    assert libcluster.I_cluster_point_part(byref(cluster), null_value, 0, 0) == 1
+
+
+def test_clear_resets_the_structure() -> None:
+    cluster = new_cluster(2)
+    libcluster.I_cluster_point(byref(cluster), point(1.0, 2.0))
+
+    assert libcluster.I_cluster_clear(byref(cluster)) == 0
+
+    assert cluster.nbands == 0
+    assert not cluster.points
+    assert not cluster.band_sum
+    assert not cluster.band_sum2