Fix one-pixel shift when smooth_edge_size is even

cv2.morphologyEx(MORPH_OPEN) applies a single anchor to both the erosion
and the dilation it performs. That is only correct when the anchor lands
on the structuring element's centre of symmetry: true for odd kernel
sizes, false for even ones. For even smooth_edge_size the result was the
correct opening translated one pixel down and to the right, and was not
anti-extensive -- smoothing could add pixels to a class rather than only
removing them. It was most visible at smooth_edge_size=2, where the
opening is otherwise a no-op on a solid blob so the shift was the only
effect.

Run the erosion and dilation separately, anchoring the dilation at the
reflection of the erosion's anchor. For odd sizes the two anchors
coincide and output is bit-identical to before, including border
handling; the Landsat regression fixture (smooth_edge_size=3) is
unchanged.

The dilation writes back into the erosion's buffer, as morphologyEx does
internally. Allocating a second full-size buffer instead costs more than
the morphology itself at small kernels -- 2.6x on a 4000x4000 mask at
smooth_edge_size=2, where allocation dominates the trivial compute. As
written, runtime and peak memory match the previous implementation.

This changes output for even smooth_edge_size, including the default
of 2.

Tests: pin the opening against an independent implementation of the
definition, exact at borders as well as in the interior, and cover the
even sizes the parameter sweeps previously skipped (SWEEP_SMOOTH was
[0, 1, 3]).

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

Author: wrignj08

CHANGELOG.md

@@ -4,6 +4,21 @@ All notable changes to MultiClean are documented here.
 
 ## [Unreleased]
 
+### Fixed
+- Edge smoothing no longer translates the output by one pixel down and to the
+  right when `smooth_edge_size` is even. `cv2.morphologyEx(MORPH_OPEN)` applies
+  a single anchor to both the erosion and the dilation, which is only correct
+  when that anchor coincides with the structuring element's centre of symmetry
+  — true for odd kernel sizes, false for even ones. The result was a shifted
+  opening that could also add pixels to a class rather than only removing them.
+  The erosion and dilation are now run separately with the dilation anchored at
+  the reflection of the erosion's anchor, giving a true opening at every size.
+
+  **This changes output for even `smooth_edge_size`, including the default of
+  `2`.** Output for odd values is bit-identical to previous releases. If you
+  have cached results produced with an even `smooth_edge_size`, regenerate
+  them or expect a one-pixel offset against new output.
+
 ## [0.4.0] - 2026-07-28
 
 ### Changed

multiclean/utils.py

@@ -93,12 +93,39 @@ def smooth_edges_to_codes(
     kernel = create_circle_kernel(smooth_edge_size)
     codes = np.zeros(array.shape, dtype=code_dtype)
 
+    # A true opening dilates with the *reflected* structuring element. cv2's
+    # ``morphologyEx(MORPH_OPEN)`` reuses one anchor for both passes, which is
+    # self-correcting only when the anchor sits on the element's centre of
+    # symmetry -- true for odd ``smooth_edge_size``, false for even, where it
+    # translated the result down and right by one pixel (and so was not even
+    # anti-extensive). Anchoring the erosion at cv2's default and the dilation
+    # at its reflection removes the shift; for odd sizes the two anchors
+    # coincide and output is bit-identical to the previous implementation.
+    erode_anchor = smooth_edge_size // 2
+    dilate_anchor = smooth_edge_size - 1 - erode_anchor
+    kernel_reflected = np.ascontiguousarray(kernel[::-1, ::-1])
+
     def _opened_for_class(cv_) -> Tuple[object, np.ndarray]:
         # bool storage is 1 byte/element so ``.view(np.uint8)`` is a zero-
         # copy reinterpretation -- avoids the bool→uint8 astype copy.
         class_mask_u8 = (array == cv_).view(np.uint8)
-        opened_u8 = cv2.morphologyEx(
-            class_mask_u8, cv2.MORPH_OPEN, kernel, iterations=1
+        eroded_u8 = cv2.erode(
+            class_mask_u8,
+            kernel,
+            anchor=(erode_anchor, erode_anchor),
+            iterations=1,
+        )
+        # Dilate back into the erosion's buffer, which is what
+        # ``morphologyEx(MORPH_OPEN)`` does internally. Allocating a second
+        # full-size buffer instead costs more than the morphology itself at
+        # small kernel sizes. ``eroded_u8`` is local to this call, so this
+        # stays safe under the thread pool below.
+        opened_u8 = cv2.dilate(
+            eroded_u8,
+            kernel_reflected,
+            dst=eroded_u8,
+            anchor=(dilate_anchor, dilate_anchor),
+            iterations=1,
         )
         return cv_, opened_u8.view(bool)
 

tests/test_multiclean.py

@@ -4,6 +4,7 @@
 import pytest
 
 from multiclean import clean_array
+from multiclean.utils import create_circle_kernel, smooth_edges_to_codes
 
 TEST_DATA_DIR = Path(__file__).resolve().parent / "data"
 LANDSAT_INPUT = TEST_DATA_DIR / "Landsat cloud and cloud shadow.tif"
@@ -100,6 +101,239 @@ def test_smoothing_removes_single_pixel_when_enabled():
     assert out[2, 2] == 0
 
 
+def _centred_blob(size: int = 24, half_width: int = 6) -> np.ndarray:
+    """Solid square of class 1 centred in a background-0 array.
+
+    The array is symmetric under a 180-degree rotation, which the smoothing
+    kernels are too, so any correct opening must preserve that symmetry.
+    """
+    arr = np.zeros((size, size), dtype=np.uint8)
+    centre = (size - 1) / 2
+    lo, hi = int(centre - half_width + 1), int(centre + half_width + 1)
+    arr[lo:hi, lo:hi] = 1
+    return arr
+
+
+@pytest.mark.parametrize("smooth_edge_size", [1, 2, 3, 4, 5, 6, 7, 8])
+def test_smoothing_does_not_translate_blobs(smooth_edge_size):
+    # Regression: cv2's MORPH_OPEN applies one anchor to both the erosion and
+    # the dilation, which is only correct when the structuring element's centre
+    # of symmetry lands on that anchor. For even kernel sizes it does not, and
+    # the whole blob came back shifted one pixel down and right -- most obvious
+    # at smooth_edge_size=2, where the opening should otherwise be a no-op on a
+    # solid blob.
+    arr = _centred_blob()
+    out = clean_array(
+        arr,
+        class_values=[0, 1],
+        smooth_edge_size=smooth_edge_size,
+        min_island_size=0,
+        connectivity=4,
+        max_workers=1,
+    )
+
+    # The blob is far larger than any kernel here, so it must survive with its
+    # bounding box unmoved (a circular opening rounds the corners off a square
+    # but leaves the edge midpoints, so the extent is unchanged).
+    ys, xs = np.where(out == 1)
+    in_ys, in_xs = np.where(arr == 1)
+    assert (ys.min(), ys.max()) == (in_ys.min(), in_ys.max())
+    assert (xs.min(), xs.max()) == (in_xs.min(), in_xs.max())
+
+
+@pytest.mark.parametrize("smooth_edge_size", [1, 2, 3, 4, 5, 6, 7, 8])
+def test_smoothing_preserves_symmetry(smooth_edge_size):
+    # Sharper form of the same regression: the kernels are all symmetric under
+    # a 180-degree rotation, so a symmetric input must smooth to a symmetric
+    # result. A one-pixel translation in either axis breaks this even when the
+    # bounding box survives.
+    #
+    # Asserted against the smoothing stage rather than clean_array because the
+    # nearest-neighbour fill that follows it breaks ties between equidistant
+    # source pixels arbitrarily, which is asymmetric by design.
+    arr = _centred_blob()
+    codes, _ = smooth_edges_to_codes(
+        arr,
+        smooth_edge_size=smooth_edge_size,
+        target_class_values=[1],
+        background_class_values=[0],
+        all_class_values=[0, 1],
+        max_workers=1,
+    )
+    assert np.array_equal(codes, codes[::-1, ::-1])
+
+
+@pytest.mark.parametrize("smooth_edge_size", [1, 2, 3, 4, 5, 6, 7, 8])
+def test_smoothing_never_grows_a_class(smooth_edge_size):
+    # An opening is anti-extensive: it can only remove pixels from a class,
+    # never add them. The even-size anchor bug broke this -- the shifted blob
+    # covered pixels that were background in the input. Checked against a
+    # ragged multiclass array so it is not just the solid-blob case.
+    arr = np.zeros((40, 40), dtype=np.uint8)
+    arr[5:20, 5:20] = 1
+    arr[22:36, 22:36] = 1
+    arr[10:14, 24:30] = 2
+    arr[30, 3] = 1  # thin spur that smoothing is expected to erase
+
+    out = clean_array(
+        arr,
+        class_values=[1, 2],
+        smooth_edge_size=smooth_edge_size,
+        min_island_size=0,
+        connectivity=4,
+        max_workers=1,
+    )
+
+    # Every pixel that came out as a smoothed class must have held that class
+    # on input; fill may reassign a pixel to another class, but smoothing must
+    # not extend one beyond its original footprint.
+    for cv in (1, 2):
+        grew = (out == cv) & (arr != cv)
+        # Fill can only draw from surviving neighbours, so any growth here is
+        # the smoothing step inventing coverage.
+        assert not grew.any(), f"class {cv} grew by {int(grew.sum())} pixels"
+
+
+def _reference_opening(mask: np.ndarray, kernel: np.ndarray, anchor: int) -> np.ndarray:
+    """Morphological opening straight from the definition, for small arrays.
+
+    ``A opened by B`` is the union of every translate of ``B`` that fits
+    entirely inside ``A``. Written out as an explicit slide so it shares no
+    machinery with the cv2 erode/dilate pair under test.
+
+    Two details model cv2's finite-image behaviour, so this is exact at the
+    borders and not just in the interior:
+
+    * Pixels outside the image read as foreground, matching the border value
+      cv2 uses for erosion -- content is not eaten away merely because the
+      image ends.
+    * Translate positions are restricted to the image domain, because that is
+      the domain cv2's intermediate erosion is defined on.
+
+    The second point is why ``anchor`` has to be named: it fixes where a
+    translate sits relative to the position that must stay in-domain. Away
+    from the border the choice cannot matter (a translated structuring element
+    is still the same set of pixels), and the interior test below asserts
+    exactly that. Within one kernel width of the border it does matter, so the
+    border test passes cv2's own anchor.
+    """
+    ks = kernel.shape[0]
+    height, width = mask.shape
+    offsets = [(i - anchor, j - anchor) for i, j in np.argwhere(kernel > 0)]
+
+    # Foreground-padded view, wide enough that no translate can run off it.
+    extended = np.ones((height + 2 * ks, width + 2 * ks), dtype=np.uint8)
+    extended[ks : ks + height, ks : ks + width] = mask
+
+    out = np.zeros_like(mask)
+    for row in range(height):
+        for col in range(width):
+            if all(extended[row + ks + dr, col + ks + dc] for dr, dc in offsets):
+                for dr, dc in offsets:
+                    r, c = row + dr, col + dc
+                    if 0 <= r < height and 0 <= c < width:
+                        out[r, c] = 1
+    return out
+
+
+def _smoothed_mask(arr: np.ndarray, smooth_edge_size: int) -> np.ndarray:
+    """Run the smoothing stage alone and return class 1's mask."""
+    codes, code_to_value = smooth_edges_to_codes(
+        arr,
+        smooth_edge_size=smooth_edge_size,
+        target_class_values=[1],
+        background_class_values=[],
+        all_class_values=[0, 1],
+        max_workers=1,
+    )
+    return (code_to_value[codes] == 1).astype(np.uint8)
+
+
+@pytest.mark.parametrize("smooth_edge_size", [1, 2, 3, 4, 5, 6, 7, 8])
+def test_smoothing_matches_reference_opening(smooth_edge_size):
+    # The other smoothing tests assert properties (no shift, symmetric,
+    # anti-extensive). Properties alone cannot distinguish a correct opening
+    # from a differently-wrong one -- silently rounding even kernel sizes up to
+    # odd, for instance, satisfies every one of them while changing how much
+    # smoothing the caller actually asked for. This pins the exact result
+    # against an independent implementation of the definition instead.
+    rng = np.random.default_rng(7)
+    kernel = create_circle_kernel(smooth_edge_size)
+
+    for _ in range(10):
+        arr = np.zeros((34, 34), dtype=np.uint8)
+        arr[9:25, 9:25] = rng.random((16, 16)) > 0.35  # wide zero margin
+
+        expected = _reference_opening(arr, kernel, anchor=smooth_edge_size // 2)
+        assert np.array_equal(_smoothed_mask(arr, smooth_edge_size), expected)
+
+        # Away from the border the reference must not depend on how the
+        # structuring element is anchored. This is the no-shift property
+        # restated at the definition level, and it keeps the assertion above
+        # from silently inheriting the implementation's anchor convention.
+        for alt_anchor in (0, smooth_edge_size - 1):
+            assert np.array_equal(
+                _reference_opening(arr, kernel, anchor=alt_anchor), expected
+            )
+
+
+@pytest.mark.parametrize("smooth_edge_size", [1, 2, 3, 4, 5, 6, 7, 8])
+def test_smoothing_matches_reference_opening_at_borders(smooth_edge_size):
+    # Same equivalence, but with content running flush to all four edges, where
+    # cv2's border convention decides the answer. Untested until now: the
+    # interior test deliberately keeps a margin so border handling cannot
+    # affect it, which left the edges of every real raster unpinned.
+    rng = np.random.default_rng(11)
+    kernel = create_circle_kernel(smooth_edge_size)
+
+    for _ in range(6):
+        arr = (rng.random((20, 20)) > 0.35).astype(np.uint8)
+        expected = _reference_opening(arr, kernel, anchor=smooth_edge_size // 2)
+        assert np.array_equal(_smoothed_mask(arr, smooth_edge_size), expected)
+
+
+@pytest.mark.parametrize("smooth_edge_size", [1, 2, 3, 4, 5, 6, 7, 8])
+def test_smoothing_does_not_erode_content_at_the_image_edge(smooth_edge_size):
+    # The human-readable half of the border contract: a band running the full
+    # width of the image, flush against the top, left and right edges, must
+    # come through an opening completely intact. If the erosion treated
+    # out-of-image pixels as background it would chew a kernel-wide bite out of
+    # all three edges, which on tiled processing would show up as seams along
+    # every tile boundary.
+    #
+    # A band rather than a square block: its only boundary is the straight
+    # edge along the bottom, which an opening preserves exactly. A block would
+    # additionally have an interior corner, and a circular kernel rounds those
+    # off by design -- correct behaviour that has nothing to do with borders.
+    arr = np.zeros((24, 24), dtype=np.uint8)
+    arr[:12, :] = 1
+
+    assert np.array_equal(_smoothed_mask(arr, smooth_edge_size), arr)
+
+
+@pytest.mark.parametrize("smooth_edge_size", [0, 2, 3])
+def test_output_is_independent_of_max_workers(smooth_edge_size):
+    # Smoothing runs one class per thread and the opening now dilates back into
+    # the erosion's own buffer. That is safe because the buffer is created per
+    # call, but hoisting it out to avoid a per-class allocation would be an easy
+    # and plausible "optimisation" -- and would corrupt results only under
+    # concurrency, which every other test pins to a single worker count.
+    rng = np.random.default_rng(3)
+    arr = rng.integers(0, 6, size=(64, 64), dtype=np.uint8)
+    arr[10:30, 10:30] = 2  # solid regions so smoothing has real work to do
+    arr[35:60, 35:60] = 4
+
+    kwargs = dict(
+        class_values=[1, 2, 3, 4, 5],
+        smooth_edge_size=smooth_edge_size,
+        min_island_size=6,
+        connectivity=4,
+    )
+    baseline = clean_array(arr, max_workers=1, **kwargs)
+    for workers in (2, 4, 8):
+        assert np.array_equal(clean_array(arr, max_workers=workers, **kwargs), baseline)
+
+
 def test_island_threshold_strictness_preserves_area_equal_to_threshold():
     # 2-pixel island (area = 2) should be preserved when min_island_size = 2
     arr = np.zeros((5, 5), dtype=np.int32)
@@ -183,7 +417,7 @@ def test_class_values_absent_from_array_are_ignored():
     arr = np.full((64, 64), 254, dtype=np.uint8)
     arr[10:50, 10:50] = 1  # only class 1 (and background 254) present
 
-    for smooth_edge_size in (0, 3):
+    for smooth_edge_size in (0, 2, 3):
         kwargs = dict(
             smooth_edge_size=smooth_edge_size,
             min_island_size=10,
@@ -206,7 +440,10 @@ def test_class_values_absent_from_array_are_ignored():
 # any single branch still fails the suite.
 
 SWEEP_DTYPES = [np.uint8, np.int16, np.int32, np.float32]
-SWEEP_SMOOTH = [0, 1, 3]
+# Even kernel sizes are included deliberately: the anchor bug that translated
+# smoothed blobs by one pixel only ever fired for even ``smooth_edge_size``,
+# and the sweep's previous [0, 1, 3] never touched that half of the space.
+SWEEP_SMOOTH = [0, 1, 2, 3, 4]
 SWEEP_ISLAND = [0, 25]