diff --git a/README.md b/README.md
index 4816d19..e3b68fc 100644
--- a/README.md
+++ b/README.md
@@ -17,12 +17,57 @@ It provides a **clean Python API** backed by **custom CUDA kernels**, enabling h
 
 ---
 
-## 📦 Installation
+## 📦 Local build and development install
 
-Clone and build locally:
+The instructions in this section are for **building TorchMorph from this repository** and using it locally in a development environment. TorchMorph builds a custom CUDA extension during install, which means your PyTorch runtime and CUDA compiler (`nvcc`) must be compatible.
 
+### Recommended workflow (Stable & Decoupled)
+
+1. Create and activate a fresh conda environment:
+   ```bash
+   conda create -n torchmorph python=3.11 -y
+   conda activate torchmorph
+   ```
+
+2. Install the stable PyTorch release (e.g., CUDA 12.4):
+   ```bash
+   pip install torch --index-url https://download.pytorch.org/whl/cu124
+   ```
+
+3. Ensure you have a compatible `nvcc` compiler.
+   ```bash
+   conda install -c nvidia cuda-nvcc=12.4 -y
+   ```
+
+4. Install TorchMorph's dependencies and build the extension:
+   ```bash
+   pip install -r "requirements-dev.txt"
+   pip install --no-build-isolation -e .
+   ```
+
+### Verify the environment before building
+
+Check your PyTorch and CUDA compiler versions:
 ```bash
-git clone https://github.com/torchmorph/torchmorph.git
-cd torchmorph
-pip install -e .
+python -c "import torch; print(f'PyTorch: {torch.__version__} | PyTorch CUDA: {torch.version.cuda}')"
+nvcc --version
+```
+As long as your `nvcc` version matchs your PyTorch CUDA version (e.g., `nvcc 12.4` with PyTorch `cu124`), the extension will compile successfully.
+
+### Minimal validation
+
+After installation succeeds, verify import and a simple CUDA kernel call:
+
+```bash
+import torch
+import torchmorph as tm
+
+print(torch.__version__, torch.version.cuda, torch.cuda.is_available())
+
+if torch.cuda.is_available():
+    x = torch.tensor([[1.0, 2.0], [3.0, 4.0]], device="cuda")
+    y = tm.add(x, 1.5)
+    print(y)
+else:
+    print("CUDA not available; install verified for import only.")
 ```
diff --git a/benchmark/binary.py b/benchmark/binary.py
new file mode 100644
index 0000000..319df3c
--- /dev/null
+++ b/benchmark/binary.py
@@ -0,0 +1,317 @@
+import scipy.ndimage as ndi
+import torch
+import torch.utils.benchmark as benchmark
+from prettytable import PrettyTable
+
+import torchmorph as tm
+
+image_size = [64, 128, 256, 1024]
+batch_size = [1, 2, 4, 16]
+MIN_RUN_TIME = 1.0
+
+
+def bench_binary_erosion():
+    print("\n============================================")
+    print(" Benchmark: binary erosion ")
+    print("============================================")
+
+    for batch in batch_size:
+        table = PrettyTable()
+        table.field_names = [
+            "Size",
+            "SciPy(ms)",
+            "Torch 1x(ms)",
+            "Torch batch(ms)",
+            "Speedup 1x",
+            "Speedup batch",
+        ]
+        for column in table.field_names:
+            table.align[column] = "r"
+
+        for size in image_size:
+            x = (torch.randn(batch, 1, size, size, device="cuda") > 0).to(torch.float32)
+            # scipy data
+            x_np_list = [x[i, 0].detach().cpu().numpy() for i in range(batch)]
+            # torch 1x data
+            x_one = [x[i : i + 1] for i in range(batch)]
+
+            t_scipy = benchmark.Timer(
+                stmt="[scipy_erosion(data) for data in x_np_list]",
+                globals={
+                    "scipy_erosion": ndi.binary_erosion,
+                    "x_np_list": x_np_list,
+                },
+            ).blocked_autorange(min_run_time=MIN_RUN_TIME)
+
+            # warmup
+            for data in x_one:
+                tm.binary_erosion(data)
+            torch.cuda.synchronize()
+
+            t_torch_1x = benchmark.Timer(
+                stmt="[tm_erosion(data) for data in x_one]",
+                globals={
+                    "tm_erosion": tm.binary_erosion,
+                    "x_one": x_one,
+                },
+            ).blocked_autorange(min_run_time=MIN_RUN_TIME)
+
+            t_torch_batch = benchmark.Timer(
+                stmt="tm_erosion(x)",
+                globals={
+                    "tm_erosion": tm.binary_erosion,
+                    "x": x,
+                },
+            ).blocked_autorange(min_run_time=MIN_RUN_TIME)
+
+            scipy_ms = t_scipy.median * 1e3 / batch
+            torch_1x_ms = t_torch_1x.median * 1e3 / batch
+            torch_batch_ms = t_torch_batch.median * 1e3 / batch
+
+            speedup1x = scipy_ms / torch_1x_ms
+            speedupbatch = scipy_ms / torch_batch_ms
+
+            table.add_row(
+                [
+                    size,
+                    f"{scipy_ms:.3f}",
+                    f"{torch_1x_ms:.3f}",
+                    f"{torch_batch_ms:.3f}",
+                    f"{speedup1x:.1f}x",
+                    f"{speedupbatch:.1f}x",
+                ]
+            )
+        print(f"\n=========== Batch size : {batch} ===========")
+        print(table)
+
+
+def bench_binary_dilation():
+    print("\n============================================")
+    print(" Benchmark: binary dilation ")
+    print("============================================")
+
+    for batch in batch_size:
+        table = PrettyTable()
+        table.field_names = [
+            "Size",
+            "SciPy(ms)",
+            "Torch 1x(ms)",
+            "Torch batch(ms)",
+            "Speedup 1x",
+            "Speedup batch",
+        ]
+        for column in table.field_names:
+            table.align[column] = "r"
+
+        for size in image_size:
+            x = (torch.randn(batch, 1, size, size, device="cuda") > 0).to(torch.float32)
+            # scipy data
+            x_np_list = [x[i, 0].detach().cpu().numpy() for i in range(batch)]
+            # torch 1x data
+            x_one = [x[i : i + 1] for i in range(batch)]
+
+            t_scipy = benchmark.Timer(
+                stmt="[scipy_dilation(data) for data in x_np_list]",
+                globals={
+                    "scipy_dilation": ndi.binary_dilation,
+                    "x_np_list": x_np_list,
+                },
+            ).blocked_autorange(min_run_time=MIN_RUN_TIME)
+
+            # warmup
+            for data in x_one:
+                tm.binary_dilation(data)
+            torch.cuda.synchronize()
+
+            t_torch_1x = benchmark.Timer(
+                stmt="[tm_dilation(data) for data in x_one]",
+                globals={
+                    "tm_dilation": tm.binary_dilation,
+                    "x_one": x_one,
+                },
+            ).blocked_autorange(min_run_time=MIN_RUN_TIME)
+
+            t_torch_batch = benchmark.Timer(
+                stmt="tm_dilation(x)",
+                globals={
+                    "tm_dilation": tm.binary_dilation,
+                    "x": x,
+                },
+            ).blocked_autorange(min_run_time=MIN_RUN_TIME)
+
+            scipy_ms = t_scipy.median * 1e3 / batch
+            torch_1x_ms = t_torch_1x.median * 1e3 / batch
+            torch_batch_ms = t_torch_batch.median * 1e3 / batch
+
+            speedup1x = scipy_ms / torch_1x_ms
+            speedupbatch = scipy_ms / torch_batch_ms
+
+            table.add_row(
+                [
+                    size,
+                    f"{scipy_ms:.3f}",
+                    f"{torch_1x_ms:.3f}",
+                    f"{torch_batch_ms:.3f}",
+                    f"{speedup1x:.1f}x",
+                    f"{speedupbatch:.1f}x",
+                ]
+            )
+        print(f"\n=========== Batch size : {batch} ===========")
+        print(table)
+
+
+def bench_binary_opening():
+    print("\n============================================")
+    print(" Benchmark: binary opening ")
+    print("============================================")
+
+    for batch in batch_size:
+        table = PrettyTable()
+        table.field_names = [
+            "Size",
+            "SciPy(ms)",
+            "Torch 1x(ms)",
+            "Torch batch(ms)",
+            "Speedup 1x",
+            "Speedup batch",
+        ]
+        for column in table.field_names:
+            table.align[column] = "r"
+
+        for size in image_size:
+            x = (torch.randn(batch, 1, size, size, device="cuda") > 0).to(torch.float32)
+            # scipy data
+            x_np_list = [x[i, 0].detach().cpu().numpy() for i in range(batch)]
+            # torch 1x data
+            x_one = [x[i : i + 1] for i in range(batch)]
+
+            t_scipy = benchmark.Timer(
+                stmt="[scipy_opening(data) for data in x_np_list]",
+                globals={
+                    "scipy_opening": ndi.binary_opening,
+                    "x_np_list": x_np_list,
+                },
+            ).blocked_autorange(min_run_time=MIN_RUN_TIME)
+
+            # warmup
+            for data in x_one:
+                tm.binary_opening(data)
+
+            torch.cuda.synchronize()
+            t_torch_1x = benchmark.Timer(
+                stmt="[tm_opening(data) for data in x_one]",
+                globals={
+                    "tm_opening": tm.binary_opening,
+                    "x_one": x_one,
+                },
+            ).blocked_autorange(min_run_time=MIN_RUN_TIME)
+
+            t_torch_batch = benchmark.Timer(
+                stmt="tm_opening(x)",
+                globals={
+                    "tm_opening": tm.binary_opening,
+                    "x": x,
+                },
+            ).blocked_autorange(min_run_time=MIN_RUN_TIME)
+
+            scipy_ms = t_scipy.median * 1e3 / batch
+            torch_1x_ms = t_torch_1x.median * 1e3 / batch
+            torch_batch_ms = t_torch_batch.median * 1e3 / batch
+
+            speedup1x = scipy_ms / torch_1x_ms
+            speedupbatch = scipy_ms / torch_batch_ms
+
+            table.add_row(
+                [
+                    size,
+                    f"{scipy_ms:.3f}",
+                    f"{torch_1x_ms:.3f}",
+                    f"{torch_batch_ms:.3f}",
+                    f"{speedup1x:.1f}x",
+                    f"{speedupbatch:.1f}x",
+                ]
+            )
+        print(f"\n=========== Batch size : {batch} ===========")
+        print(table)
+
+
+def bench_binary_closing():
+    print("\n============================================")
+    print(" Benchmark: binary closing ")
+    print("============================================")
+
+    for batch in batch_size:
+        table = PrettyTable()
+        table.field_names = [
+            "Size",
+            "SciPy(ms)",
+            "Torch 1x(ms)",
+            "Torch batch(ms)",
+            "Speedup 1x",
+            "Speedup batch",
+        ]
+        for column in table.field_names:
+            table.align[column] = "r"
+
+        for size in image_size:
+            x = (torch.randn(batch, 1, size, size, device="cuda") > 0).to(torch.float32)
+            # scipy data
+            x_np_list = [x[i, 0].detach().cpu().numpy() for i in range(batch)]
+            # torch 1x data
+            x_one = [x[i : i + 1] for i in range(batch)]
+
+            t_scipy = benchmark.Timer(
+                stmt="[scipy_closing(data) for data in x_np_list]",
+                globals={
+                    "scipy_closing": ndi.binary_closing,
+                    "x_np_list": x_np_list,
+                },
+            ).blocked_autorange(min_run_time=MIN_RUN_TIME)
+
+            # warmup
+            for data in x_one:
+                tm.binary_closing(data)
+
+            torch.cuda.synchronize()
+            t_torch_1x = benchmark.Timer(
+                stmt="[tm_closing(data) for data in x_one]",
+                globals={
+                    "tm_closing": tm.binary_closing,
+                    "x_one": x_one,
+                },
+            ).blocked_autorange(min_run_time=MIN_RUN_TIME)
+
+            t_torch_batch = benchmark.Timer(
+                stmt="tm_closing(x)",
+                globals={
+                    "tm_closing": tm.binary_closing,
+                    "x": x,
+                },
+            ).blocked_autorange(min_run_time=MIN_RUN_TIME)
+
+            scipy_ms = t_scipy.median * 1e3 / batch
+            torch_1x_ms = t_torch_1x.median * 1e3 / batch
+            torch_batch_ms = t_torch_batch.median * 1e3 / batch
+
+            speedup1x = scipy_ms / torch_1x_ms
+            speedupbatch = scipy_ms / torch_batch_ms
+
+            table.add_row(
+                [
+                    size,
+                    f"{scipy_ms:.3f}",
+                    f"{torch_1x_ms:.3f}",
+                    f"{torch_batch_ms:.3f}",
+                    f"{speedup1x:.1f}x",
+                    f"{speedupbatch:.1f}x",
+                ]
+            )
+        print(f"\n=========== Batch size : {batch} ===========")
+        print(table)
+
+
+if __name__ == "__main__":
+    bench_binary_erosion()
+    bench_binary_dilation()
+    bench_binary_opening()
+    bench_binary_closing()
diff --git a/benchmark/dilation_erosion.py b/benchmark/dilation_erosion.py
deleted file mode 100644
index afe7e8e..0000000
--- a/benchmark/dilation_erosion.py
+++ /dev/null
@@ -1,91 +0,0 @@
-import scipy.ndimage as ndi
-import torch
-import torch.utils.benchmark as benchmark
-from prettytable import PrettyTable
-
-import torchmorph as tm
-
-sizes = [64, 128, 256, 512]
-batches = [1, 4, 8, 16]
-dtype = torch.float32
-device = "cuda"
-MIN_RUN = 1.0  # seconds per measurement
-
-torch.set_num_threads(torch.get_num_threads())
-
-
-def bench_single_op(op_name):
-    """
-    op_name: "dilation" or "erosion"
-    """
-
-    scipy_op = ndi.binary_dilation if op_name == "dilation" else ndi.binary_erosion
-    torch_op = tm.binary_dilation if op_name == "dilation" else tm.binary_erosion
-
-    print("\n==============================")
-    print(f"   Benchmark: Binary {op_name}")
-    print("==============================")
-
-    for B in batches:
-        table = PrettyTable()
-        table.field_names = [
-            "Size",
-            "SciPy (ms/img)",
-            "Torch 1× (ms/img)",
-            "Torch batch (ms/img)",
-            "Speedup 1×",
-            "Speedup batch",
-        ]
-        for c in table.field_names:
-            table.align[c] = "r"
-
-        for s in sizes:
-            # Generate binary input
-            x = (torch.randn(B, 1, s, s, device=device) > 0).to(dtype)
-            x_np_list = [x[i, 0].detach().cpu().numpy() for i in range(B)]
-            x_imgs = [x[i : i + 1] for i in range(B)]  # (1, 1, H, W)
-            # SciPy (CPU, one-by-one)
-            stmt_scipy = "out = [scipy_op(arr) for arr in x_np_list]"
-            t_scipy = benchmark.Timer(
-                stmt=stmt_scipy,
-                globals={"x_np_list": x_np_list, "scipy_op": scipy_op},
-            ).blocked_autorange(min_run_time=MIN_RUN)
-            scipy_ms = (t_scipy.median * 1e3) / B
-
-            # Torch CUDA (one-by-one)
-            stmt_torch1 = "out = [torch_op(img) for img in x_imgs]"
-            t_torch1 = benchmark.Timer(
-                stmt=stmt_torch1,
-                globals={"x_imgs": x_imgs, "torch_op": torch_op},
-            ).blocked_autorange(min_run_time=MIN_RUN)
-            torch1_ms = (t_torch1.median * 1e3) / B
-
-            # Torch CUDA (batched)
-            t_batch = benchmark.Timer(
-                stmt="torch_op(x)",
-                globals={"x": x, "torch_op": torch_op},
-            ).blocked_autorange(min_run_time=MIN_RUN)
-            torchB_ms = (t_batch.median * 1e3) / B
-
-            # Speedups
-            speed1 = scipy_ms / torch1_ms
-            speedB = scipy_ms / torchB_ms
-
-            table.add_row(
-                [
-                    s,
-                    f"{scipy_ms:.3f}",
-                    f"{torch1_ms:.3f}",
-                    f"{torchB_ms:.3f}",
-                    f"{speed1:.1f}×",
-                    f"{speedB:.1f}×",
-                ]
-            )
-
-        print(f"\n=== Batch Size: {B} ===")
-        print(table)
-
-
-print("Loaded from:", tm.__file__)
-bench_single_op("dilation")
-bench_single_op("erosion")
diff --git a/benchmark/generate.py b/benchmark/generate.py
new file mode 100644
index 0000000..aebddd8
--- /dev/null
+++ b/benchmark/generate.py
@@ -0,0 +1,59 @@
+import scipy.ndimage as ndi
+import torch.utils.benchmark as benchmark
+from prettytable import PrettyTable
+
+import torchmorph as tm
+
+DATA_CASES = [(1, 1), (2, 1), (4, 1), (8, 1), (10, 1), (10, 4), (10, 10), (12, 1)]
+MIN_RUN = 1.0
+
+
+def bench_generate_binary_structure():
+    print("\n========================================")
+    print("\n Benchmark: generate binary structure ")
+    print("\n========================================")
+
+    table = PrettyTable()
+    table.field_names = ["Rank", "Connectivity", "SciPy(ms)", "Torch(ms)", "Speedup"]
+    for column in table.field_names:
+        table.align[column] = "r"
+
+    for rank, connectivity in DATA_CASES:
+        t_scipy = benchmark.Timer(
+            stmt="generate_binary_structure(rank, connectivity)",
+            globals={
+                "generate_binary_structure": ndi.generate_binary_structure,
+                "rank": rank,
+                "connectivity": connectivity,
+            },
+        ).blocked_autorange(min_run_time=MIN_RUN)
+
+        t_torch = benchmark.Timer(
+            stmt="generate_binary_structure(rank, connectivity)",
+            globals={
+                "generate_binary_structure": tm.generate_binary_structure,
+                "rank": rank,
+                "connectivity": connectivity,
+            },
+        ).blocked_autorange(min_run_time=MIN_RUN)
+
+        scipy_ms = t_scipy.median * 1e3
+        torch_ms = t_torch.median * 1e3
+        speedup = scipy_ms / torch_ms
+
+        table.add_row(
+            [
+                rank,
+                connectivity,
+                f"{scipy_ms:.3f}",
+                f"{torch_ms:.3f}",
+                f"{speedup:.1f}x",
+            ]
+        )
+
+    print("Load from:", tm.__file__)
+    print(table)
+
+
+if __name__ == "__main__":
+    bench_generate_binary_structure()
diff --git a/test/test_add.py b/test/test_add.py
index 5f647ff..76de045 100644
--- a/test/test_add.py
+++ b/test/test_add.py
@@ -18,4 +18,4 @@ def test_add():
     torch.testing.assert_close(y, expected)
     assert y.device.type == "cuda"
     assert y.shape == x.shape
-    print("tm.bar test passed ✅")
+    print("tm.add test passed ✅")
diff --git a/test/test_binary.py b/test/test_binary.py
new file mode 100644
index 0000000..c961b45
--- /dev/null
+++ b/test/test_binary.py
@@ -0,0 +1,326 @@
+import numpy as np  # noqa: F401
+import pytest
+import torch
+from scipy.ndimage import (
+    binary_closing,
+    binary_dilation,
+    binary_erosion,
+    binary_opening,
+    generate_binary_structure,
+)
+
+import torchmorph as tm
+
+pytestmark = pytest.mark.skipif(
+    not torch.cuda.is_available(), reason="CUDA is required for torchmorph tests"
+)
+
+
+def apply_scipy_to_batch(np_input, scipy_func, **kwargs):
+    """
+    Apply a SciPy operator to each (B, C) sample independently.
+    The input is expected to have shape (B, C, Spatial...). The batch and
+    channel dimensions are flattened, `scipy_func` is applied to each spatial
+    sample separately, and the results are reshaped back to the original
+    layout.
+    Args:
+        np_input: Input array in (B, C, Spatial...) format.
+        scipy_func: SciPy function applied to each spatial sample.
+        **kwargs: Additional keyword arguments forwarded to `scipy_func`.
+    Returns:
+        np.ndarray: Output array with the same shape as `np_input`.
+    """
+    input_shape = np_input.shape
+    batch_shape = input_shape[:2]
+    spatial_shape = input_shape[2:]
+    batch_size = int(np.prod(batch_shape))
+    flatten_input = np_input.reshape(batch_size, *spatial_shape)
+
+    mask = kwargs.pop("mask", None)
+    output = kwargs.pop("output", None)
+    if mask is not None and mask.shape != spatial_shape:
+        raise ValueError(f"mask shape must be {spatial_shape}, got {mask.shape}")
+    if output is not None and output.shape != spatial_shape:
+        raise ValueError(f"output shape must be {spatial_shape}, got {output.shape}")
+
+    results = []
+    for sample in flatten_input:
+        sample_kwargs = dict(kwargs)
+        if mask is not None:
+            sample_kwargs["mask"] = mask
+        if output is not None:
+            sample_output = np.empty_like(output)
+            sample_kwargs["output"] = sample_output
+        result = scipy_func(sample, **sample_kwargs)
+        results.append(np.asarray(result).copy())
+    return np.stack(results, axis=0).reshape(*batch_shape, *spatial_shape)
+
+
+def batch_scipy(
+    np_input: np.ndarray,
+    scipy_func: callable,
+    structure: np.ndarray | None = None,
+    iterations: int = 1,
+    mask: np.ndarray | None = None,
+    output: np.ndarray | None = None,
+    border_value: bool = False,
+    origin: int | tuple[int, ...] = 0,
+) -> np.ndarray:
+    return apply_scipy_to_batch(
+        np_input,
+        scipy_func,
+        structure=structure,
+        iterations=iterations,
+        mask=mask,
+        output=output,
+        border_value=border_value,
+        origin=origin,
+    )
+
+
+# case
+case_2d = np.array(
+    [
+        [
+            [
+                [0, 1, 0],
+                [1, 1, 1],
+                [0, 1, 0],
+            ]
+        ]
+    ],
+    dtype=np.float32,
+)
+case_3d = np.zeros((2, 1, 5, 5, 5), dtype=bool)
+case_3d[0, 0, 2:3, 2:4, 1:2] = True
+case_4d = np.zeros((2, 1, 4, 4, 4, 4), dtype=bool)
+case_4d[0, 0, 1:3, 1:2, 2:3] = True
+
+# structure
+structure_2d = np.array([[0, 1, 0], [1, 1, 1], [0, 1, 0]])
+
+structure_3d_1 = generate_binary_structure(rank=3, connectivity=2)
+structure_3d_2 = generate_binary_structure(rank=3, connectivity=3)
+
+structure_4d = generate_binary_structure(rank=4, connectivity=4)
+
+# mask
+mask_2d_np = np.array([[1, 1, 0], [1, 1, 1], [1, 1, 0]], dtype=bool)
+mask_2d_tm = torch.tensor([[[[1, 1, 0], [1, 1, 1], [1, 1, 0]]]], dtype=bool)
+
+mask_3d_np = np.zeros((5, 5, 5), dtype=bool)
+mask_3d_np[0:2, 0:1, 0:3] = True
+mask_3d_tm = torch.zeros((2, 1, 5, 5, 5), dtype=bool)
+mask_3d_tm[:, :, 0:2, 0:1, 0:3] = True
+
+mask_4d_np = np.zeros((4, 4, 4, 4), dtype=bool)
+mask_4d_np[1:2, 0:1, 1:3] = True
+mask_4d_tm = torch.zeros((2, 1, 4, 4, 4, 4), dtype=bool)
+mask_4d_tm[:, :, 1:2, 0:1, 1:3] = True
+
+# output
+output_2d_np = np.empty([3, 3])
+output_2d_tm = torch.empty([1, 1, 3, 3])
+
+output_3d_np = np.empty([5, 5, 5])
+output_3d_tm = torch.empty([2, 1, 5, 5, 5])
+
+output_4d_np = np.empty([4, 4, 4, 4])
+output_4d_tm = torch.empty([2, 1, 4, 4, 4, 4])
+
+
+@pytest.mark.parametrize(
+    ("np_input, scipy_func, structure, iterations, origin, border_value"),
+    [
+        pytest.param(case_2d, binary_erosion, None, 1, 0, False, id="er_2D_basic"),
+        pytest.param(case_2d, binary_erosion, structure_2d, 1, 0, False, id="er_2D_structure"),
+        pytest.param(case_2d, binary_erosion, None, 2, 0, False, id="er_2D_2iterations"),
+        pytest.param(case_2d, binary_erosion, None, -1, 0, False, id="er_2D_-1iterations"),
+        pytest.param(case_2d, binary_erosion, None, 1, 1, False, id="er_2D_1origin_False"),
+        pytest.param(case_2d, binary_erosion, None, 1, 1, True, id="er_2D_1origin_True"),
+        pytest.param(case_3d, binary_erosion, None, 1, 0, False, id="er_3D_basic"),
+        pytest.param(case_3d, binary_erosion, structure_3d_1, 1, 0, False, id="er_3D_2Dstructure1"),
+        pytest.param(case_3d, binary_erosion, structure_3d_2, 1, 0, False, id="er_3D_3Dstructure2"),
+        pytest.param(case_3d, binary_erosion, None, 2, 0, False, id="er_3D_2iterations"),
+        pytest.param(case_3d, binary_erosion, None, -1, 0, False, id="er_3D_-1iterations"),
+        pytest.param(case_3d, binary_erosion, None, 1, 1, False, id="er_3D_1origin_False"),
+        pytest.param(case_3d, binary_erosion, None, 1, 1, True, id="er_3D_1origin_True"),
+        pytest.param(case_4d, binary_erosion, None, 1, 0, False, id="er_4D_basic"),
+        pytest.param(case_4d, binary_erosion, structure_4d, 1, 0, False, id="er_4D_4Dstructure"),
+        pytest.param(case_4d, binary_erosion, None, 2, 0, False, id="er_4D_2iterations"),
+        pytest.param(case_4d, binary_erosion, None, -1, 0, False, id="er_4D_-1iterations"),
+        pytest.param(case_4d, binary_erosion, None, 1, 1, False, id="er_4D_1origin_False"),
+        pytest.param(case_4d, binary_erosion, None, 1, 1, True, id="er_4D_1origin_True"),
+        pytest.param(case_2d, binary_dilation, None, 1, 0, False, id="di_2D_basic"),
+        pytest.param(case_2d, binary_dilation, structure_2d, 1, 0, False, id="di_2D_structure"),
+        pytest.param(case_2d, binary_dilation, None, 2, 0, False, id="di_2D_2iterations"),
+        pytest.param(case_2d, binary_dilation, None, -1, 0, False, id="di_2D_-1iterations"),
+        pytest.param(case_2d, binary_dilation, None, 1, 1, False, id="di_2D_1origin_False"),
+        pytest.param(case_2d, binary_dilation, None, 1, 1, True, id="di_2D_1origin_True"),
+        pytest.param(case_3d, binary_dilation, None, 1, 0, False, id="di_3D_basic"),
+        pytest.param(
+            case_3d, binary_dilation, structure_3d_1, 1, 0, False, id="di_3D_2Dstructure1"
+        ),
+        pytest.param(
+            case_3d, binary_dilation, structure_3d_2, 1, 0, False, id="di_3D_3Dstructure2"
+        ),
+        pytest.param(case_3d, binary_dilation, None, 2, 0, False, id="di_3D_2iterations"),
+        pytest.param(case_3d, binary_dilation, None, -1, 0, False, id="di_3D_-1iterations"),
+        pytest.param(case_3d, binary_dilation, None, 1, 1, False, id="di_3D_1origin_False"),
+        pytest.param(case_3d, binary_dilation, None, 1, 1, True, id="di_3D_1origin_True"),
+        pytest.param(case_4d, binary_dilation, None, 1, 0, False, id="di_4D_basic"),
+        pytest.param(case_4d, binary_dilation, structure_4d, 1, 0, False, id="di_4D_4Dstructure"),
+        pytest.param(case_4d, binary_dilation, None, 2, 0, False, id="di_4D_2iterations"),
+        pytest.param(case_4d, binary_dilation, None, -1, 0, False, id="di_4D_-1iterations"),
+        pytest.param(case_4d, binary_dilation, None, 1, 1, False, id="di_4D_1origin_False"),
+        pytest.param(case_4d, binary_dilation, None, 1, 1, True, id="di_4D_1origin_True"),
+        pytest.param(case_2d, binary_opening, None, 1, 0, False, id="op_2D_basic"),
+        pytest.param(case_2d, binary_opening, structure_2d, 1, 0, False, id="op_2D_structure"),
+        pytest.param(case_2d, binary_opening, None, 2, 0, False, id="op_2D_2iterations"),
+        pytest.param(case_2d, binary_opening, None, -1, 0, False, id="op_2D_-1iterations"),
+        pytest.param(case_2d, binary_opening, None, 1, 1, False, id="op_2D_1origin_False"),
+        pytest.param(case_2d, binary_opening, None, 1, 1, True, id="op_2D_1origin_True"),
+        pytest.param(case_3d, binary_opening, None, 1, 0, False, id="op_3D_basic"),
+        pytest.param(case_3d, binary_opening, structure_3d_1, 1, 0, False, id="op_3D_2Dstructure1"),
+        pytest.param(case_3d, binary_opening, structure_3d_2, 1, 0, False, id="op_3D_3Dstructure2"),
+        pytest.param(case_3d, binary_opening, None, 2, 0, False, id="op_3D_2iterations"),
+        pytest.param(case_3d, binary_opening, None, -1, 0, False, id="op_3D_-1iterations"),
+        pytest.param(case_3d, binary_opening, None, 1, 1, False, id="op_3D_1origin_False"),
+        pytest.param(case_3d, binary_opening, None, 1, 1, True, id="op_3D_1origin_True"),
+        pytest.param(case_4d, binary_opening, None, 1, 0, False, id="op_4D_basic"),
+        pytest.param(case_4d, binary_opening, structure_4d, 1, 0, False, id="op_4D_4Dstructure"),
+        pytest.param(case_4d, binary_opening, None, 2, 0, False, id="op_4D_2iterations"),
+        pytest.param(case_4d, binary_opening, None, -1, 0, False, id="op_4D_-1iterations"),
+        pytest.param(case_4d, binary_opening, None, 1, 1, False, id="op_4D_1origin_False"),
+        pytest.param(case_4d, binary_opening, None, 1, 1, True, id="op_4D_1origin_True"),
+        pytest.param(case_2d, binary_closing, None, 1, 0, False, id="cl_2D_basic"),
+        pytest.param(case_2d, binary_closing, structure_2d, 1, 0, False, id="cl_2D_structure"),
+        pytest.param(case_2d, binary_closing, None, 2, 0, False, id="cl_2D_2iterations"),
+        pytest.param(case_2d, binary_closing, None, -1, 0, False, id="cl_2D_-1iterations"),
+        pytest.param(case_2d, binary_closing, None, 1, 1, False, id="cl_2D_1origin_False"),
+        pytest.param(case_2d, binary_closing, None, 1, 1, True, id="cl_2D_1origin_True"),
+        pytest.param(case_3d, binary_closing, None, 1, 0, False, id="cl_3D_basic"),
+        pytest.param(case_3d, binary_closing, structure_3d_1, 1, 0, False, id="cl_3D_2Dstructure1"),
+        pytest.param(case_3d, binary_closing, structure_3d_2, 1, 0, False, id="cl_3D_3Dstructure2"),
+        pytest.param(case_3d, binary_closing, None, 2, 0, False, id="cl_3D_2iterations"),
+        pytest.param(case_3d, binary_closing, None, -1, 0, False, id="cl_3D_-1iterations"),
+        pytest.param(case_3d, binary_closing, None, 1, 1, False, id="cl_3D_1origin_False"),
+        pytest.param(case_3d, binary_closing, None, 1, 1, True, id="cl_3D_1origin_True"),
+        pytest.param(case_4d, binary_closing, None, 1, 0, False, id="cl_4D_basic"),
+        pytest.param(case_4d, binary_closing, structure_4d, 1, 0, False, id="cl_4D_4Dstructure"),
+        pytest.param(case_4d, binary_closing, None, 2, 0, False, id="cl_4D_2iterations"),
+        pytest.param(case_4d, binary_closing, None, -1, 0, False, id="cl_4D_-1iterations"),
+        pytest.param(case_4d, binary_closing, None, 1, 1, False, id="cl_4D_1origin_False"),
+        pytest.param(case_4d, binary_closing, None, 1, 1, True, id="cl_4D_1origin_True"),
+    ],
+)
+def test_binary_basic(
+    np_input,
+    scipy_func,
+    structure,
+    iterations,
+    origin,
+    border_value,
+):
+    x = torch.as_tensor(np_input, dtype=torch.float32)
+    if structure is not None:
+        structure_cuda = torch.as_tensor(structure, dtype=torch.float32)
+    else:
+        structure_cuda = None
+    tm_func = getattr(tm, scipy_func.__name__)
+    actual = tm_func(
+        x,
+        structure=structure_cuda,
+        iterations=iterations,
+        origin=origin,
+        border_value=border_value,
+    )
+    expected_np = batch_scipy(
+        np_input,
+        scipy_func,
+        structure=structure,
+        iterations=iterations,
+        origin=origin,
+        border_value=border_value,
+    )
+    expected = torch.as_tensor(expected_np)
+    torch.testing.assert_close(actual, expected)
+
+
+@pytest.mark.parametrize(
+    ("input_np, scipy_func, mask_np, mask_tm"),
+    [
+        pytest.param(case_2d, binary_erosion, mask_2d_np, mask_2d_tm, id="er_2d_mask"),
+        pytest.param(case_2d, binary_dilation, mask_2d_np, mask_2d_tm, id="di_2d_mask"),
+        pytest.param(case_2d, binary_opening, mask_2d_np, mask_2d_tm, id="op_2d_mask"),
+        pytest.param(case_2d, binary_closing, mask_2d_np, mask_2d_tm, id="cl_2d_mask"),
+        pytest.param(case_3d, binary_erosion, mask_3d_np, mask_3d_tm, id="er_3d_mask"),
+        pytest.param(case_3d, binary_dilation, mask_3d_np, mask_3d_tm, id="di_3d_mask"),
+        pytest.param(case_3d, binary_opening, mask_3d_np, mask_3d_tm, id="op_3d_mask"),
+        pytest.param(case_3d, binary_closing, mask_3d_np, mask_3d_tm, id="cl_3d_mask"),
+        pytest.param(case_4d, binary_erosion, mask_4d_np, mask_4d_tm, id="er_4d_mask"),
+        pytest.param(case_4d, binary_dilation, mask_4d_np, mask_4d_tm, id="di_4d_mask"),
+        pytest.param(case_4d, binary_opening, mask_4d_np, mask_4d_tm, id="op_4d_mask"),
+        pytest.param(case_4d, binary_closing, mask_4d_np, mask_4d_tm, id="cl_4d_mask"),
+    ],
+)
+def test_binary_mask(
+    input_np,
+    scipy_func,
+    mask_np,
+    mask_tm,
+):
+    x = torch.as_tensor(input_np, dtype=torch.float32)
+    tm_func = getattr(tm, scipy_func.__name__)
+    actual = tm_func(
+        x,
+        mask=mask_tm,
+    )
+
+    expected_np = batch_scipy(
+        input_np,
+        scipy_func,
+        mask=mask_np,
+    )
+    expected = torch.as_tensor(expected_np)
+    torch.testing.assert_close(actual, expected)
+
+
+@pytest.mark.parametrize(
+    "input_np,scipy_func, output_np, output_tm",
+    [
+        pytest.param(case_2d, binary_erosion, output_2d_np, output_2d_tm, id='er_output_2d'),
+        pytest.param(case_2d, binary_dilation, output_2d_np, output_2d_tm, id='di_output_2d'),
+        pytest.param(case_2d, binary_opening, output_2d_np, output_2d_tm, id='op_output_2d'),
+        pytest.param(case_2d, binary_closing, output_2d_np, output_2d_tm, id='cl_output_2d'),
+        pytest.param(case_3d, binary_erosion, output_3d_np, output_3d_tm, id='er_output_3d'),
+        pytest.param(case_3d, binary_dilation, output_3d_np, output_3d_tm, id='di_output_3d'),
+        pytest.param(case_3d, binary_opening, output_3d_np, output_3d_tm, id='op_output_3d'),
+        pytest.param(case_3d, binary_closing, output_3d_np, output_3d_tm, id='cl_output_3d'),
+        pytest.param(case_4d, binary_erosion, output_4d_np, output_4d_tm, id='er_output_4d'),
+        pytest.param(case_4d, binary_dilation, output_4d_np, output_4d_tm, id='di_output_4d'),
+        pytest.param(case_4d, binary_opening, output_4d_np, output_4d_tm, id='op_output_4d'),
+        pytest.param(case_4d, binary_closing, output_4d_np, output_4d_tm, id='cl_output_4d'),
+    ],
+)
+def test_binary_output(
+    input_np,
+    scipy_func,
+    output_np,
+    output_tm,
+):
+    x = torch.as_tensor(input_np, dtype=torch.float32)
+    tm_func = getattr(tm, scipy_func.__name__)
+
+    actual = tm_func(
+        x,
+        output=output_tm,
+    )
+
+    expected_np = batch_scipy(
+        input_np,
+        scipy_func,
+        output=output_np,
+    )
+    expected = torch.as_tensor(expected_np, dtype=torch.float32)
+    torch.testing.assert_close(actual, expected)
diff --git a/test/test_structure.py b/test/test_structure.py
new file mode 100644
index 0000000..55726ce
--- /dev/null
+++ b/test/test_structure.py
@@ -0,0 +1,28 @@
+import numpy as np  # noqa: F401
+import pytest
+import torch
+from scipy.ndimage import generate_binary_structure as scipy_generate_binary_structure
+
+import torchmorph as tm  # noqa: F401
+
+
+@pytest.mark.parametrize(
+    "rank, connectivity",
+    [
+        pytest.param(1, 1, id="1D"),
+        pytest.param(2, 1, id="2D1C"),
+        pytest.param(2, 2, id="2D2C"),
+        pytest.param(3, 1, id="3D1C"),
+        pytest.param(3, 3, id="3D3C"),
+        pytest.param(4, 4, id="4D4C"),
+        pytest.param(7, 2, id="7D2C"),
+    ],
+)
+def test_generate_binary_structure(rank, connectivity):
+    expected = scipy_generate_binary_structure(rank, connectivity)
+    actual = tm.generate_binary_structure(rank, connectivity)
+
+    np.testing.assert_array_equal(actual, expected)
+
+    expected_tensor = torch.as_tensor(expected)
+    torch.testing.assert_close(actual, expected_tensor)
diff --git a/torchmorph/__init__.py b/torchmorph/__init__.py
index 70b9520..be99a6e 100644
--- a/torchmorph/__init__.py
+++ b/torchmorph/__init__.py
@@ -1,16 +1,25 @@
 from .add import add
-from .dilation_erosion import binary_dilation, binary_erosion
 from .distance_transform import (
     brute_force_distance_transform,
     chamfer_distance_transform,
     euclidean_distance_transform,
 )
+from .morphology import (
+    binary_closing,
+    binary_dilation,
+    binary_erosion,
+    binary_opening,
+    generate_binary_structure,
+)
 
 __all__ = [
     "add",
     "euclidean_distance_transform",
     "chamfer_distance_transform",
     "brute_force_distance_transform",
-    "binary_dilation",
+    "generate_binary_structure",
     "binary_erosion",
+    "binary_dilation",
+    "binary_opening",
+    "binary_closing",
 ]
diff --git a/torchmorph/dilation_erosion.py b/torchmorph/dilation_erosion.py
deleted file mode 100644
index 4e35d1b..0000000
--- a/torchmorph/dilation_erosion.py
+++ /dev/null
@@ -1,321 +0,0 @@
-from typing import Optional, Sequence, Tuple, Union
-
-import torch
-import torch.nn.functional as F
-
-
-def _to_bool_tensor(x: torch.Tensor) -> torch.Tensor:
-    """
-    Convert an input value into a boolean PyTorch tensor.
-
-    This helper function ensures that the input is represented as a
-    `torch.bool` tensor, which is the internal format required by
-    binary morphological operations (e.g., dilation/erosion).
-
-    Behavior:
-    - If `x` is not already a tensor, it is converted using `torch.tensor(x)`.
-    - Non-zero values become `True`; zero values become `False`.
-
-    Args:
-        x (torch.Tensor or array-like):
-            Input data. May be a Python list, scalar, NumPy array, or torch.Tensor.
-
-    Returns:
-        torch.Tensor (dtype=torch.bool):
-            Boolean tensor where each element is `True` if corresponding input value
-            is non-zero, otherwise `False`.
-
-    Examples:
-        >>> _to_bool_tensor([0, 1, 2])
-        tensor([False, True, True])
-
-        >>> _to_bool_tensor(torch.tensor([3.0, 0.0]))
-        tensor([True, False])
-    """
-    # If x is not a tensor yet (e.g., list, numpy array, int, float), convert to tensor.
-    if not torch.is_tensor(x):
-        x = torch.tensor(x)
-
-    # Convert input tensor into boolean by checking non-zero status.
-    # Non-zero -> True, zero -> False.
-    return x != 0
-
-
-def _normalize_structure(structure: Optional[torch.Tensor], ndim: int) -> torch.Tensor:
-    """
-    Normalize a structuring element into a boolean tensor with the correct
-    number of spatial dimensions.
-
-    This utility function standardizes user-provided structuring elements
-    for binary morphological operations (e.g., dilation and erosion).
-
-    Behavior:
-    1. If `structure` is None, a default full-connectivity structuring
-       element of shape (3, 3, ..., 3) with `ndim` dimensions is created.
-       This matches the default behavior of scipy.ndimage morphology.
-    2. If `structure` is provided, it is converted into a boolean tensor,
-       where non-zero values are treated as True.
-    3. The dimensionality of the structuring element is strictly checked
-       to ensure it matches the spatial dimensionality of the input.
-       A mismatch indicates an invalid morphological definition and
-       raises a ValueError.
-
-    Args:
-        structure (Optional[torch.Tensor]):
-            Structuring element defining the neighborhood for morphology.
-            If None, a full (3,) * ndim boolean structure is used.
-        ndim (int):
-            Number of spatial dimensions of the input (e.g., 2 for H×W,
-            3 for D×H×W). Batch and channel dimensions are excluded.
-
-    Returns:
-        torch.Tensor (dtype=torch.bool):
-            An `ndim`-dimensional boolean tensor representing the normalized
-            structuring element.
-
-    Raises:
-        ValueError:
-            If the provided structuring element does not have exactly
-            `ndim` dimensions.
-
-    Notes:
-        - This function does not enforce any particular kernel size other
-          than dimensionality; arbitrary shapes are allowed.
-        - Channel and batch dimensions are intentionally not supported
-          for structuring elements, as morphology is defined purely in
-          spatial dimensions.
-
-    Examples:
-        >>> _normalize_structure(None, ndim=2)
-        tensor([[True, True, True],
-                [True, True, True],
-                [True, True, True]])
-
-        >>> _normalize_structure([[0, 1, 0],
-        ...                       [1, 1, 1],
-        ...                       [0, 1, 0]], ndim=2)
-        tensor([[False, True, False],
-                [ True, True, True],
-                [False, True, False]])
-    """
-    # Case 1: No structuring element provided by the user.
-    # Use a default full-connectivity neighborhood of size 3 in each
-    # spatial dimension (e.g., 3×3 for 2D, 3×3×3 for 3D).
-    if structure is None:
-        shape = (3,) * ndim
-        return torch.ones(shape, dtype=torch.bool)
-
-    # Case 2: A structuring element is provided.
-    # Convert it to a boolean tensor so that non-zero values indicate
-    # active neighbors and zero values are ignored.
-    st = _to_bool_tensor(structure)
-
-    # Validate dimensionality: the structuring element must have the same
-    # number of dimensions as the spatial dimensions of the input tensor.
-    if st.ndim != ndim:
-        raise ValueError(f"structure must be {ndim}-D (got {st.ndim}-D)")
-
-    # Return the normalized boolean structuring element.
-    return st
-
-
-def _origin_to_tuple(
-    origin: Union[int, Sequence[int], Tuple[int, ...]], ndim: int
-) -> Tuple[int, ...]:
-    """
-    Normalize the `origin` argument into an ndim-length tuple.
-
-    The origin defines the anchor point of the structuring element,
-    consistent with SciPy's definition.
-
-    Args:
-        origin (int or sequence of int):
-            If an int is given, it is broadcast to all spatial dimensions.
-            If a sequence is given, its length must match `ndim`.
-        ndim (int):
-            Number of spatial dimensions.
-
-    Returns:
-        Tuple[int, ...]:
-            Origin offset per spatial dimension.
-    """
-    # If a scalar is given, replicate it across all dimensions.
-    if isinstance(origin, int):
-        return tuple([origin] * ndim)
-
-    # Otherwise, ensure it is a tuple with correct dimensionality.
-    origin = tuple(origin)
-    if len(origin) != ndim:
-        raise ValueError("origin must match spatial ndim")
-
-    return origin
-
-
-def _pad_for_kernel(
-    kernel_shape: Sequence[int], origin: Sequence[int]
-) -> Tuple[Tuple[int, int], ...]:
-    """
-    Compute per-dimension padding sizes required to keep output shape
-    identical to input shape after convolution.
-
-    This takes into account the kernel size and the origin offset.
-
-    Returns:
-        Tuple of (pad_before, pad_after) for each spatial dimension.
-    """
-    pads = []
-    for k, o in zip(kernel_shape, origin):
-        # Default symmetric padding would be k//2,
-        # but origin shifts the effective center.
-        pad_before = k // 2 - o
-        pad_after = k - 1 - pad_before
-
-        # Padding must be non-negative.
-        pad_before = max(pad_before, 0)
-        pad_after = max(pad_after, 0)
-
-        pads.append((pad_before, pad_after))
-    return tuple(pads)
-
-
-def _make_padding_tuple_for_Fpad(pads: Tuple[Tuple[int, int], ...]) -> Tuple[int, ...]:
-    """
-    Convert per-dimension padding into the flattened format required
-    by torch.nn.functional.pad.
-
-    PyTorch expects padding in reverse order:
-        (pad_last_dim_left, pad_last_dim_right, ..., pad_first_dim_left, pad_first_dim_right)
-    """
-    flat = []
-    for pb, pa in reversed(pads):
-        flat.append(pb)
-        flat.append(pa)
-    return tuple(flat)
-
-
-def _conv_nd(x: torch.Tensor, kernel: torch.Tensor, ndim: int) -> torch.Tensor:
-    """
-    Dispatch N-dimensional convolution based on spatial dimensionality.
-
-    Args:
-        x (torch.Tensor):
-            Input tensor of shape (B*C, 1, *spatial_dims)
-        kernel (torch.Tensor):
-            Structuring element kernel.
-        ndim (int):
-            Number of spatial dimensions (1, 2, or 3).
-
-    Returns:
-        torch.Tensor:
-            Convolution result.
-    """
-    # Convert kernel into convolution weight:
-    # shape -> (out_channels=1, in_channels=1, *kernel_shape)
-    weight = kernel.to(dtype=x.dtype, device=x.device).unsqueeze(0).unsqueeze(0)
-
-    if ndim == 1:
-        return F.conv1d(x, weight)
-    elif ndim == 2:
-        return F.conv2d(x, weight)
-    elif ndim == 3:
-        return F.conv3d(x, weight)
-    else:
-        raise NotImplementedError("Only supports 1D/2D/3D")
-
-
-def _morph_op(
-    input_tensor: torch.Tensor,
-    structure: Optional[torch.Tensor],
-    iterations: int,
-    origin: Union[int, Sequence[int]],
-    border_value: int,
-    mode: str,
-) -> torch.Tensor:
-    """
-    Core implementation of binary dilation and erosion using convolution.
-
-    This function supports batch and channel dimensions by flattening
-    (B, C) into a single dimension and applying morphology independently
-    per channel.
-
-    Args:
-        input_tensor (torch.Tensor):
-            Input binary tensor.
-        structure (Optional[torch.Tensor]):
-            Structuring element.
-        iterations (int):
-            Number of times to apply the operation.
-        origin:
-            Origin offset of the structuring element.
-        border_value (int):
-            Value used for padding outside image boundaries.
-        mode (str):
-            Either 'dilation' or 'erosion'.
-
-    Returns:
-        torch.Tensor (dtype=torch.bool):
-            Output binary tensor.
-    """
-    if mode not in ('dilation', 'erosion'):
-        raise ValueError("mode must be 'dilation' or 'erosion'")
-
-    x = input_tensor
-    if not torch.is_tensor(x):
-        x = torch.tensor(x)
-
-    # Convert input to boolean (binary morphology).
-    x_bool = x != 0
-    # Supported input shapes:
-    # (H,W), (C,H,W), (B,C,H,W), (B,C,D,H,W)
-    full_ndim = x_bool.ndim
-
-    if full_ndim < 2:
-        raise NotImplementedError("Need at least 2D (H,W)")
-    if full_ndim > 5:
-        raise NotImplementedError("Only supports up to 5D (B,C,D,H,W)")
-    spatial_ndim = full_ndim - 2  # remove (B,C)
-    if not (1 <= spatial_ndim <= 3):
-        raise NotImplementedError("Supports 1D/2D/3D spatial dims")
-
-    B, C = x_bool.shape[0], x_bool.shape[1]
-    spatial_shape = x_bool.shape[2:]
-    st = _normalize_structure(structure, spatial_ndim)
-    origin_t = _origin_to_tuple(origin, spatial_ndim)
-
-    k_sum = st.sum().item()
-    kernel = st.to(torch.float32)
-
-    # Apply origin shift by rolling kernel.
-    for axis, o in enumerate(origin_t):
-        if o != 0:
-            kernel = torch.roll(kernel, shifts=-o, dims=axis)
-    pads = _pad_for_kernel(kernel.shape, origin_t)
-    pad_tuple = _make_padding_tuple_for_Fpad(pads)
-
-    cur = x_bool.to(torch.float32)
-
-    # Flatten (B,C) -> (B*C,1)
-    cur = cur.view(B * C, 1, *spatial_shape)
-    for _ in range(max(1, iterations)):
-        x_pad = F.pad(cur, pad_tuple, value=float(border_value))
-        conv_res = _conv_nd(x_pad, kernel, spatial_ndim)
-
-        if mode == 'dilation':
-            # Any overlap -> True
-            cur = (conv_res > 0).to(torch.float32)
-        else:
-            # Full overlap -> True
-            if k_sum == 0:
-                cur = torch.ones_like(cur)
-            else:
-                cur = (conv_res >= (k_sum - 1e-6)).to(torch.float32)
-    out = cur.view(B, C, *spatial_shape)
-    return out.to(torch.bool)
-
-
-def binary_dilation(input_tensor, structure=None, iterations=1, origin=0, border_value=0):
-    return _morph_op(input_tensor, structure, iterations, origin, border_value, mode="dilation")
-
-
-def binary_erosion(input_tensor, structure=None, iterations=1, origin=0, border_value=0):
-    return _morph_op(input_tensor, structure, iterations, origin, border_value, mode="erosion")
diff --git a/torchmorph/morphology/__init__.py b/torchmorph/morphology/__init__.py
new file mode 100644
index 0000000..0efa28d
--- /dev/null
+++ b/torchmorph/morphology/__init__.py
@@ -0,0 +1,10 @@
+from .binary import binary_closing, binary_dilation, binary_erosion, binary_opening
+from .structure import generate_binary_structure
+
+__all__ = [
+    "generate_binary_structure",
+    "binary_erosion",
+    "binary_dilation",
+    "binary_opening",
+    "binary_closing",
+]
diff --git a/torchmorph/morphology/_convnd.py b/torchmorph/morphology/_convnd.py
new file mode 100644
index 0000000..b1ee854
--- /dev/null
+++ b/torchmorph/morphology/_convnd.py
@@ -0,0 +1,60 @@
+import torch
+from torch import Tensor
+
+_NATIVE_CONVS = {
+    1: torch.nn.functional.conv1d,
+    2: torch.nn.functional.conv2d,
+    3: torch.nn.functional.conv3d,
+}
+
+
+def conv_nd(x: Tensor, weight: Tensor) -> Tensor:
+    """Apply correlation on tensors shaped as ``(N, C, Spatial...)``.
+
+    This is a minimal private helper for ``torchmorph.morphology.binary``.
+    Inputs are expected to already be padded, use channel-first layout, and
+    convolve over the trailing spatial dimensions with stride=1, dilation=1,
+    and groups=1.
+    """
+
+    num_spatial = x.ndim - 2
+    if num_spatial < 1:
+        raise ValueError(f"expected at least 1 spatial dim, got input ndim={x.ndim}")
+    if weight.ndim != x.ndim:
+        raise ValueError(
+            f"expected weight ndim {x.ndim} to match input ndim {x.ndim}, got {weight.ndim}"
+        )
+
+    return _conv_core(x, weight)
+
+
+def _conv_core(x: Tensor, weight: Tensor) -> Tensor:
+    num_spatial = x.ndim - 2
+    if num_spatial in _NATIVE_CONVS:
+        return _NATIVE_CONVS[num_spatial](x, weight, bias=None, stride=1, padding=0, dilation=1)
+    return _conv_recursive(x, weight)
+
+
+def _conv_recursive(x: Tensor, weight: Tensor) -> Tensor:
+    kernel_size = weight.shape[2]
+    input_size = x.shape[2]
+    output_size = input_size - kernel_size + 1
+    batch_size = x.shape[0]
+
+    accumulated = None
+    for kernel_index in range(kernel_size):
+        x_slice = x[:, :, kernel_index : kernel_index + output_size]
+        flattened = x_slice.moveaxis(2, 1).reshape(
+            batch_size * output_size,
+            x.shape[1],
+            *x_slice.shape[3:],
+        )
+
+        partial = _conv_core(flattened, weight[:, :, kernel_index])
+        partial = partial.reshape(batch_size, output_size, partial.shape[1], *partial.shape[2:])
+        partial = partial.moveaxis(2, 1)
+        accumulated = partial if accumulated is None else accumulated + partial
+
+    if accumulated is None:
+        raise RuntimeError("recursive convolution produced no output")
+    return accumulated
diff --git a/torchmorph/morphology/binary.py b/torchmorph/morphology/binary.py
new file mode 100644
index 0000000..49ec8b9
--- /dev/null
+++ b/torchmorph/morphology/binary.py
@@ -0,0 +1,222 @@
+import torch
+import torch.nn.functional as F
+from torch import Tensor
+
+from ._convnd import conv_nd
+from .structure import generate_binary_structure
+
+
+def _prepare_origin(origin: int | tuple[int, ...], ndim=int) -> tuple[int, ...]:
+    """change the origin into tuple"""
+    if isinstance(origin, int):
+        return (origin,) * ndim
+    origin = tuple(origin)
+
+    if (len(origin)) != ndim:
+        raise ValueError(f"origin dimension is not {ndim}, got {len(origin)}")
+
+    return origin
+
+
+def _extend_pad(kernel_shape: torch.Size, origin: tuple[int, ...]) -> list[int]:
+    """extend the padlist for kernel"""
+    pad = []
+    for dim in range(len(kernel_shape) - 1, -1, -1):
+        center = kernel_shape[dim] // 2
+        pad_before = center + origin[dim]
+        pad_after = kernel_shape[dim] - 1 - pad_before
+        pad.extend([pad_before, pad_after])
+    return pad
+
+
+def _flip_structure(structure: Tensor) -> Tensor:
+    dim = tuple(range(structure.ndim))
+    return torch.flip(structure, dim)
+
+
+def _binary_morphology(
+    input: Tensor,
+    structure: Tensor | None,
+    iterations: int,
+    mask: Tensor | None,
+    output: Tensor | None,
+    border_value: bool,
+    origin: int | tuple[int, ...],
+    *,
+    mode: str,
+) -> Tensor:
+    iterations_flag = iterations < 1
+
+    spatial_ndim = input.ndim - 2
+
+    if structure is None:
+        structure = generate_binary_structure(spatial_ndim, 1)
+
+    batch, channels = input.shape[:2]
+    spatial_shape = input.shape[2:]
+
+    x = (input != 0).to(dtype=torch.float32).reshape(batch * channels, 1, *spatial_shape)
+    structure = (structure != 0).to(device=input.device, dtype=torch.float32)
+    if mode == "dilation":
+        structure = _flip_structure(structure)
+    kernel = structure.unsqueeze(0).unsqueeze(0)
+    kernel_sum = kernel.sum()
+
+    origin = _prepare_origin(origin, spatial_ndim)
+    if mode == "dilation":
+        origin = tuple(-value for value in origin)
+    pad = _extend_pad(structure.shape, origin)
+    pad_value = float(bool(border_value))
+
+    if mask is not None:
+        mask_flat = mask.to(dtype=torch.bool).reshape(batch * channels, 1, *spatial_shape)
+        input_flat = (
+            (input != 0).to(dtype=torch.float32).reshape(batch * channels, 1, *spatial_shape)
+        )
+    else:
+        mask_flat = None
+        input_flat = None
+
+    if iterations_flag:
+        old = None
+        while True:
+            x_padded = F.pad(x, pad, value=pad_value)
+            conv = conv_nd(x_padded, kernel)
+            if mode == "erosion":
+                x = (conv == kernel_sum).to(dtype=torch.float32)
+            else:
+                x = (conv > 0).to(dtype=torch.float32)
+            if mask_flat is not None:
+                x = torch.where(mask_flat, x, input_flat)
+
+            if old is not None and torch.equal(x, old):
+                break
+
+            old = x.clone()
+    else:
+        for _ in range(iterations):
+            x_padded = F.pad(x, pad, value=pad_value)
+            conv = conv_nd(x_padded, kernel)
+            if mode == "erosion":
+                x = (conv == kernel_sum).to(dtype=torch.float32)
+            else:
+                x = (conv > 0).to(dtype=torch.float32)
+
+            if mask_flat is not None:
+                x = torch.where(mask_flat, x, input_flat)
+
+    result = x.reshape(batch, channels, *spatial_shape).to(dtype=torch.bool)
+    if output is not None:
+        output.copy_(result)
+        return output
+    return result
+
+
+def binary_erosion(
+    input: Tensor,
+    structure: Tensor | None = None,
+    iterations: int = 1,
+    mask: Tensor | None = None,
+    output: Tensor | None = None,
+    border_value: bool = False,
+    origin: int | tuple[int, ...] = 0,
+) -> Tensor:
+    """
+    N-dimensional binary erosion for `(B, C, Spatial...)` tensors.
+
+    For a single image or volume, add batch and channel dimensions first.
+    """
+    return _binary_morphology(
+        input=input,
+        structure=structure,
+        iterations=iterations,
+        mask=mask,
+        output=output,
+        border_value=border_value,
+        origin=origin,
+        mode="erosion",
+    )
+
+
+def binary_dilation(
+    input: Tensor,
+    structure: Tensor | None = None,
+    iterations: int = 1,
+    mask: Tensor | None = None,
+    output: Tensor | None = None,
+    border_value: bool = False,
+    origin: int | tuple[int, ...] = 0,
+) -> Tensor:
+    """binary dilation for `(B, C, Spatial...)` tensors."""
+    return _binary_morphology(
+        input,
+        structure,
+        iterations,
+        mask,
+        output,
+        border_value,
+        origin,
+        mode="dilation",
+    )
+
+
+def binary_opening(
+    input: Tensor,
+    structure: Tensor | None = None,
+    iterations: int = 1,
+    mask: Tensor | None = None,
+    output: Tensor | None = None,
+    border_value: bool = False,
+    origin: int | tuple[int, ...] = 0,
+) -> Tensor:
+    """binary opening for '(B, C, ...)' Tensors ."""
+    x = binary_erosion(
+        input,
+        structure,
+        iterations,
+        mask,
+        output,
+        border_value,
+        origin,
+    )
+    x = binary_dilation(
+        x,
+        structure,
+        iterations,
+        mask,
+        output,
+        border_value,
+        origin,
+    )
+    return x
+
+
+def binary_closing(
+    input: Tensor,
+    structure: Tensor | None = None,
+    iterations: int = 1,
+    mask: Tensor | None = None,
+    output: Tensor | None = None,
+    border_value: bool = False,
+    origin: int | tuple[int, ...] = 0,
+) -> Tensor:
+    """binary closing for (B, C, ...) Tensors."""
+    x = binary_dilation(
+        input,
+        structure,
+        iterations,
+        mask,
+        output,
+        border_value,
+        origin,
+    )
+    x = binary_erosion(
+        x,
+        structure,
+        iterations,
+        mask,
+        output,
+        border_value,
+        origin,
+    )
+    return x
diff --git a/torchmorph/morphology/structure.py b/torchmorph/morphology/structure.py
new file mode 100644
index 0000000..497c42a
--- /dev/null
+++ b/torchmorph/morphology/structure.py
@@ -0,0 +1,16 @@
+import torch
+from torch import Tensor
+
+
+def generate_binary_structure(rank: int, connectivity: int) -> Tensor:
+    """N-D generate binary structure"""
+    if connectivity < 1 or connectivity > rank:
+        raise ValueError(f"connectivity must be in [1, rank], got {connectivity}")
+
+    if rank < 1:
+        raise ValueError(f"rank must be >= 1, got {rank}")
+
+    axes = [torch.tensor([-1, 0, 1]) for _ in range(rank)]
+    grids = torch.meshgrid(*axes, indexing="ij")
+    offsets = torch.stack(grids, dim=0)
+    return (offsets != 0).sum(dim=0) <= connectivity