diff --git a/GWYDDION_FEATURE_GAP.md b/GWYDDION_FEATURE_GAP.md
index 965d3db..e708040 100644
--- a/GWYDDION_FEATURE_GAP.md
+++ b/GWYDDION_FEATURE_GAP.md
@@ -29,7 +29,7 @@ Reference for future implementation. Grouped by value to typical SPM workflows.
 |---|---------|---------------|-------------|
 | 15 | Correlation / Pattern Matching | crosscor.c, maskcor.c | Find repeated features or align images via cross-correlation. |
 | 16 | Slope Distribution | slope_dist.c | Angular histogram of surface slopes. Characterizes surface texture directionality. |
-| 17 | Grain Filtering | grain_filter.c | Remove particles by size, height, or border contact. Refine grain masks post-detection. |
+| 17 | Grain Filtering | grain_filter.c | Remove grains by size, height, or border contact. Refine grain masks post-detection. |
 | 18 | Field Arithmetic | arithmetic.c | Add/subtract/multiply/divide two DATA_FIELDs. Useful for difference maps, normalization. |
 | 19 | Spot Removal | spotremove.c | Interpolate over selected point defects (dust, spikes). |
 | 20 | Tip Modeling / Deconvolution | tip_blind.c, tip_model.c | Estimate tip shape from image, deconvolve to recover true surface. |
@@ -94,8 +94,8 @@ For reference, these Gwyddion equivalents are already covered:
 | Threshold Mask | mask | threshold.c, otsu_threshold.c |
 | Mask Morphology | mask | mask_morph.c (erode, dilate, open, close) |
 | Mask Invert | mask | — |
-| Mask Combine | mask | — (boolean AND, OR, XOR, subtract) |
+| Mask Operations | mask | — (boolean logic on two masks: AND, OR, XOR, NAND, NOR, XNOR, implication, etc.) |
 | Grain Distance Transform | mask | mask_edt.c |
-| Watershed Segmentation | particles | grain_wshed.c |
-| Particle Analysis | particles | grain_stat.c |
+| Watershed Segmentation | grains | grain_wshed.c |
+| Grain Analysis | grains | grain_stat.c |
 | Preview / 3D View / Print Table | display | Presentation, 3D view |
diff --git a/backend/data_types.py b/backend/data_types.py
index 95ae9c4..e199387 100644
--- a/backend/data_types.py
+++ b/backend/data_types.py
@@ -32,7 +32,7 @@ PREVIEW_MARKUP_REFERENCE_DIM = 512
 
 
 class RecordTable(list):
-    """Tabular rows with a shared schema, e.g. particle statistics."""
+    """Tabular rows with a shared schema, e.g. grain statistics."""
 
 
 class MeasureTable(list):
diff --git a/backend/execution_context.py b/backend/execution_context.py
index ca19cfb..5a5a287 100644
--- a/backend/execution_context.py
+++ b/backend/execution_context.py
@@ -2,6 +2,7 @@ from __future__ import annotations
 
 from contextlib import contextmanager
 from contextvars import ContextVar
+import inspect
 from typing import Any, Callable
 
 Callback = Callable[[str, Any], None]
@@ -12,6 +13,15 @@ _callbacks_var: ContextVar[dict[str, Callback | None]] = ContextVar(
 )
 _node_id_var: ContextVar[str | None] = ContextVar("argonode_execution_node_id", default=None)
 
+_LEGACY_CALLBACK_ATTRS = {
+    "preview": "_broadcast_fn",
+    "table": "_broadcast_table_fn",
+    "mesh": "_broadcast_mesh_fn",
+    "overlay": "_broadcast_overlay_fn",
+    "value": "_broadcast_value_fn",
+    "warning": "_broadcast_warning_fn",
+}
+
 
 @contextmanager
 def execution_callbacks(
@@ -50,12 +60,42 @@ def current_node_id() -> str | None:
     return _node_id_var.get()
 
 
+def _legacy_emit(kind: str, payload: Any) -> bool:
+    callback_attr = _LEGACY_CALLBACK_ATTRS.get(kind)
+    if not callback_attr:
+        return False
+
+    frame = inspect.currentframe()
+    try:
+        frame = frame.f_back
+        while frame is not None:
+            for owner_name in ("self", "cls"):
+                owner = frame.f_locals.get(owner_name)
+                if owner is None:
+                    continue
+
+                candidate = owner if isinstance(owner, type) else owner.__class__
+                callback = getattr(candidate, callback_attr, None)
+                node_id = getattr(candidate, "_current_node_id", None)
+                if callback is not None and node_id:
+                    callback(str(node_id), payload)
+                    return True
+            frame = frame.f_back
+    finally:
+        del frame
+
+    return False
+
+
 def _emit(kind: str, payload: Any) -> None:
     callbacks = _callbacks_var.get()
     callback = callbacks.get(kind)
     node_id = current_node_id()
     if callback is not None and node_id:
         callback(node_id, payload)
+        return
+
+    _legacy_emit(kind, payload)
 
 
 def emit_preview(payload: Any) -> None:
diff --git a/backend/node_menu.py b/backend/node_menu.py
index e718e73..d12c127 100644
--- a/backend/node_menu.py
+++ b/backend/node_menu.py
@@ -25,11 +25,11 @@ MENU_LAYOUT: dict[str, list[str]] = {
     ],
     "Output": [
         "PreviewImage",
+        "ValueDisplay",
+        "View3D",
         "Save",
         "SaveImage",
-        "View3D",
         "PrintTable",
-        "ValueDisplay",
     ],
     "Overlay": [
         "Markup",
@@ -37,10 +37,10 @@ MENU_LAYOUT: dict[str, list[str]] = {
         "AngleMeasure",
     ],
     "Modify": [
-        "ColormapAdjust",
         "CropResizeField",
         "RotateField",
         "FlipField",
+        "ColormapAdjust",
     ],
     "Filter": [
         "GaussianFilter",
@@ -51,24 +51,24 @@ MENU_LAYOUT: dict[str, list[str]] = {
         "ScarRemoval",
     ],
     "Frequency": [
-        "FFT2D",
         "PSDF",
+        "FFT2D",
         "InverseFFT2D",
     ],
     "Flatten": [
-        "PlaneLevelField",
-        "FacetLevelField",
-        "PolyLevelField",
         "FixZero",
         "LineCorrection",
+        "PlaneLevelField",
+        "PolyLevelField",
+        "FacetLevelField",
     ],
     "Measure": [
         "CrossSection",
-        "Curvature",
         "Histogram",
+        "Cursors",
+        "Curvature",
         "FractalDimension",
         "ACF",
-        "Cursors",
         "Statistics",
         "Stats",
     ],
@@ -77,12 +77,12 @@ MENU_LAYOUT: dict[str, list[str]] = {
         "ThresholdMask",
         "MaskMorphology",
         "MaskInvert",
-        "MaskCombine",
-        "GrainDistanceTransform",
+        "MaskOperations",
     ],
-    "Particles": [
+    "Grains": [
+        "GrainAnalysis",
+        "GrainDistanceTransform",
         "WatershedSegmentation",
-        "ParticleAnalysis",
     ],
 }
 
diff --git a/backend/nodes/__init__.py b/backend/nodes/__init__.py
index 1e1d1a7..2016cb5 100644
--- a/backend/nodes/__init__.py
+++ b/backend/nodes/__init__.py
@@ -32,7 +32,7 @@ from backend.nodes import (
     threshold_mask,
     mask_morphology,
     mask_invert,
-    mask_combine,
+    mask_operations,
     grain_distance_transform,
     # Correction
     scar_removal,
@@ -59,9 +59,5 @@ from backend.nodes import (
     cross_section,
     stats,
     watershed_segmentation,
+    grain_analysis,
 )
-
-try:
-    from backend.nodes import particle_analysis
-except ImportError:
-    pass
diff --git a/backend/nodes/particle_analysis.py b/backend/nodes/grain_analysis.py
similarity index 76%
rename from backend/nodes/particle_analysis.py
rename to backend/nodes/grain_analysis.py
index 75f26aa..f2d2b7c 100644
--- a/backend/nodes/particle_analysis.py
+++ b/backend/nodes/grain_analysis.py
@@ -5,8 +5,8 @@ from backend.data_types import DataField, RecordTable
 from backend.nodes.helpers import _square_unit
 
 
-@register_node(display_name="Particle Analysis")
-class ParticleAnalysis:
+@register_node(display_name="Grain Analysis")
+class GrainAnalysis:
     @classmethod
     def INPUT_TYPES(cls):
         return {
@@ -18,44 +18,44 @@ class ParticleAnalysis:
         }
 
     RETURN_TYPES = ("RECORD_TABLE",)
-    RETURN_NAMES = ("particle_stats",)
+    RETURN_NAMES = ("grain_stats",)
     FUNCTION = "process"
 
     DESCRIPTION = (
-        "Label connected particle regions in a binary mask and compute per-particle "
+        "Label connected grain regions in a binary mask and compute per-grain "
         "statistics: area, equivalent diameter, mean/max height, bounding box. "
-        "Equivalent to gwy_data_field_particles_get_values."
+        "Equivalent to Gwyddion's grain statistics tools."
     )
 
     def process(self, field: DataField, mask: np.ndarray, min_size: int) -> tuple:
         from scipy.ndimage import label
 
         binary = (mask > 127).astype(np.int32)
-        labeled, n_particles = label(binary)
+        labeled, n_grains = label(binary)
 
         pixel_area = field.dx * field.dy
         xy_unit = str(field.si_unit_xy or "").strip()
         z_unit = str(field.si_unit_z or "").strip()
 
         rows = RecordTable()
-        for pid in range(1, n_particles + 1):
-            particle_pixels = labeled == pid
-            area_px = int(particle_pixels.sum())
+        for gid in range(1, n_grains + 1):
+            grain_pixels = labeled == gid
+            area_px = int(grain_pixels.sum())
             if area_px < min_size:
                 continue
 
             area_m2 = area_px * pixel_area
             equiv_diam = float(2.0 * np.sqrt(area_m2 / np.pi))
 
-            heights = field.data[particle_pixels]
+            heights = field.data[grain_pixels]
             mean_h = float(heights.mean())
             max_h = float(heights.max())
 
-            ys, xs = np.where(particle_pixels)
+            ys, xs = np.where(grain_pixels)
             bbox = f"({int(xs.min())},{int(ys.min())})-({int(xs.max())},{int(ys.max())})"
 
             rows.append({
-                "particle_id":  pid,
+                "grain_id":     gid,
                 "area_px":      area_px,
                 "area_px_unit": _square_unit("px"),
                 "area_m2":      area_m2,
diff --git a/backend/nodes/mask_combine.py b/backend/nodes/mask_combine.py
deleted file mode 100644
index ab25489..0000000
--- a/backend/nodes/mask_combine.py
+++ /dev/null
@@ -1,61 +0,0 @@
-from __future__ import annotations
-import numpy as np
-from backend.node_registry import register_node
-from backend.execution_context import emit_preview
-from backend.data_types import DataField, encode_preview
-from backend.nodes.helpers import _mask_overlay
-
-
-@register_node(display_name="Mask Combine")
-class MaskCombine:
-    _CUSTOM_PREVIEW = True
-
-    @classmethod
-    def INPUT_TYPES(cls):
-        return {
-            "required": {
-                "mask_a": ("IMAGE",),
-                "mask_b": ("IMAGE",),
-                "operation": (["and", "or", "xor", "subtract"],),
-            },
-            "optional": {
-                "field": ("DATA_FIELD",),
-            }
-        }
-
-    RETURN_TYPES = ("IMAGE",)
-    RETURN_NAMES = ("mask",)
-    FUNCTION = "process"
-
-    DESCRIPTION = (
-        "Combine two binary masks with a boolean operation. "
-        "AND keeps overlap, OR merges, XOR keeps non-overlapping regions, "
-        "subtract removes mask_b from mask_a."
-    )
-
-    _broadcast_fn = None
-    _current_node_id: str = ""
-
-    def process(self, mask_a: np.ndarray, mask_b: np.ndarray, operation: str,
-                field: DataField | None = None) -> tuple:
-        a = mask_a > 127
-        b = mask_b > 127
-
-        if operation == "and":
-            result = a & b
-        elif operation == "or":
-            result = a | b
-        elif operation == "xor":
-            result = a ^ b
-        elif operation == "subtract":
-            result = a & ~b
-        else:
-            raise ValueError(f"Unknown mask operation: {operation}")
-
-        out = result.astype(np.uint8) * 255
-
-        if field is not None:
-            overlay = _mask_overlay(field, out)
-            emit_preview(encode_preview(overlay))
-
-        return (out,)
diff --git a/backend/nodes/mask_operations.py b/backend/nodes/mask_operations.py
new file mode 100644
index 0000000..1e7ed01
--- /dev/null
+++ b/backend/nodes/mask_operations.py
@@ -0,0 +1,65 @@
+from __future__ import annotations
+import numpy as np
+from backend.node_registry import register_node
+
+
+_MASK_BOOLEAN_OPERATIONS = {
+    "and": lambda a, b: a & b,
+    "or": lambda a, b: a | b,
+    "xor": lambda a, b: a ^ b,
+    "xnor": lambda a, b: ~(a ^ b),
+    "nand": lambda a, b: ~(a & b),
+    "nor": lambda a, b: ~(a | b),
+    "a_minus_b": lambda a, b: a & ~b,
+    "b_minus_a": lambda a, b: b & ~a,
+    "a": lambda a, b: a,
+    "b": lambda a, b: b,
+    "not_a": lambda a, b: ~a,
+    "not_b": lambda a, b: ~b,
+    "a_implies_b": lambda a, b: ~a | b,
+    "b_implies_a": lambda a, b: ~b | a,
+    "false": lambda a, b: np.zeros_like(a, dtype=bool),
+    "true": lambda a, b: np.ones_like(a, dtype=bool),
+}
+
+
+@register_node(display_name="Mask Operations")
+class MaskOperations:
+    _CUSTOM_PREVIEW = True
+
+    @classmethod
+    def INPUT_TYPES(cls):
+        return {
+            "required": {
+                "mask_a": ("IMAGE",),
+                "mask_b": ("IMAGE",),
+                "operation": (list(_MASK_BOOLEAN_OPERATIONS.keys()),),
+            },
+        }
+
+    RETURN_TYPES = ("IMAGE",)
+    RETURN_NAMES = ("mask",)
+    FUNCTION = "process"
+
+    DESCRIPTION = (
+        "Apply boolean logic to two binary masks. Includes AND, OR, XOR, NAND, NOR, "
+        "XNOR, directional subtraction, implication, pass-through, and constant true/false outputs."
+    )
+
+    def process(
+        self,
+        mask_a: np.ndarray,
+        mask_b: np.ndarray,
+        operation: str,
+    ) -> tuple:
+        a = mask_a > 127
+        b = mask_b > 127
+
+        op = _MASK_BOOLEAN_OPERATIONS.get(operation)
+        if op is None:
+            raise ValueError(f"Unknown mask operation: {operation}")
+        result = op(a, b)
+
+        out = result.astype(np.uint8) * 255
+
+        return (out,)
diff --git a/backend/nodes/plane_level_field.py b/backend/nodes/plane_level_field.py
index 6bf8f10..bb8c15f 100644
--- a/backend/nodes/plane_level_field.py
+++ b/backend/nodes/plane_level_field.py
@@ -4,6 +4,49 @@ from backend.node_registry import register_node
 from backend.data_types import DataField
 
 
+def _normalize_mask(mask: np.ndarray | None, shape: tuple[int, int]) -> np.ndarray | None:
+    if mask is None:
+        return None
+
+    mask_array = np.asarray(mask)
+    if mask_array.shape[:2] != shape:
+        raise ValueError(f"Mask shape {mask_array.shape} does not match field shape {shape}.")
+    return mask_array > 127
+
+
+def _fit_plane(
+    data: np.ndarray,
+    mask: np.ndarray | None,
+    masking: str,
+) -> tuple[float, float, float, np.ndarray, np.ndarray]:
+    yres, xres = data.shape
+    x = np.linspace(0.0, 1.0, xres)
+    y = np.linspace(0.0, 1.0, yres)
+    xx, yy = np.meshgrid(x, y)
+
+    if mask is None or masking == "ignore":
+        valid = np.ones(data.shape, dtype=bool)
+    elif masking == "include":
+        valid = mask
+    elif masking == "exclude":
+        valid = ~mask
+    else:
+        raise ValueError(f"Unknown masking mode: {masking}")
+
+    if np.count_nonzero(valid) < 3:
+        raise ValueError("Plane Level requires at least three usable pixels for fitting.")
+
+    A = np.column_stack([
+        np.ones(int(np.count_nonzero(valid)), dtype=np.float64),
+        xx[valid].ravel(),
+        yy[valid].ravel(),
+    ])
+    z = data[valid].ravel()
+    coeffs, _, _, _ = np.linalg.lstsq(A, z, rcond=None)
+    pa, pbx, pby = coeffs
+    return float(pa), float(pbx), float(pby), xx, yy
+
+
 @register_node(display_name="Plane Level")
 class PlaneLevelField:
     @classmethod
@@ -11,7 +54,11 @@ class PlaneLevelField:
         return {
             "required": {
                 "field": ("DATA_FIELD",),
-            }
+                "masking": (["ignore", "include", "exclude"], {"default": "ignore"}),
+            },
+            "optional": {
+                "mask": ("IMAGE",),
+            },
         }
 
     RETURN_TYPES = ("DATA_FIELD",)
@@ -19,27 +66,19 @@ class PlaneLevelField:
     FUNCTION = "process"
 
     DESCRIPTION = (
-        "Fit and subtract a least-squares plane from the data. "
-        "Equivalent to gwy_data_field_fit_plane + gwy_data_field_plane_level."
+        "Fit and subtract a least-squares plane from the data. Supports include/exclude mask fitting "
+        "for flattening around features, similar to masked plane fitting workflows in Gwyddion."
     )
 
-    def process(self, field: DataField) -> tuple:
+    def process(
+        self,
+        field: DataField,
+        masking: str = "ignore",
+        mask: np.ndarray | None = None,
+    ) -> tuple:
         data = field.data.copy()
-        yres, xres = data.shape
-
-        x = np.linspace(0.0, 1.0, xres)
-        y = np.linspace(0.0, 1.0, yres)
-        xx, yy = np.meshgrid(x, y)
-
-        A = np.column_stack([
-            np.ones(xres * yres),
-            xx.ravel(),
-            yy.ravel(),
-        ])
-        z = data.ravel()
-
-        coeffs, _, _, _ = np.linalg.lstsq(A, z, rcond=None)
-        pa, pbx, pby = coeffs
+        mask_array = _normalize_mask(mask, data.shape)
+        pa, pbx, pby, xx, yy = _fit_plane(data, mask_array, masking)
 
         plane = (pa + pbx * xx + pby * yy)
         return (field.replace(data=data - plane),)
diff --git a/frontend/src/constants.js b/frontend/src/constants.js
index baef17a..719c7f7 100644
--- a/frontend/src/constants.js
+++ b/frontend/src/constants.js
@@ -38,7 +38,7 @@ export const CAT_COLORS = {
   modify:    '#0f766e',
   level:     '#1b5e20',
   analysis:  '#4a148c',
-  particles: '#bf360c',
+  Grains:    '#bf360c',
   display:   '#212121',
 };
 
diff --git a/frontend/tests/valueFormatting.test.mjs b/frontend/tests/valueFormatting.test.mjs
index a970e84..d2f3847 100644
--- a/frontend/tests/valueFormatting.test.mjs
+++ b/frontend/tests/valueFormatting.test.mjs
@@ -11,7 +11,7 @@ import {
 test('getTableColumns hides companion record-table unit columns', () => {
   const columns = getTableColumns([
     {
-      particle_id: 1,
+      grain_id: 1,
       area_m2: 2.5e-12,
       area_m2_unit: 'm^2',
       mean_height: 1.5e-9,
@@ -20,7 +20,7 @@ test('getTableColumns hides companion record-table unit columns', () => {
     },
   ]);
 
-  assert.deepEqual(columns, ['particle_id', 'area_m2', 'mean_height', 'bbox']);
+  assert.deepEqual(columns, ['grain_id', 'area_m2', 'mean_height', 'bbox']);
 });
 
 test('formatTableRowCell appends companion units inline for record-table values', () => {
diff --git a/pytest.ini b/pytest.ini
index e64a3bf..7f49248 100644
--- a/pytest.ini
+++ b/pytest.ini
@@ -1,2 +1,3 @@
 [pytest]
-norecursedirs = .git .venv .pytest_cache frontend/node_modules frontend/dist pytest-cache-files-*
+testpaths = tests
+norecursedirs = .git .venv .pytest_cache frontend/node_modules frontend/dist pytest-cache-files-* desktop-dist desktop-build
diff --git a/tests/test_grains.py b/tests/test_grains.py
index ade075d..8aca46f 100644
--- a/tests/test_grains.py
+++ b/tests/test_grains.py
@@ -1,12 +1,12 @@
 """
-Thorough tests for the particles/particle analysis pipeline:
-    ThresholdMask → GrainAnalysis
+Thorough tests for the grain-analysis pipeline:
+    ThresholdMask -> GrainAnalysis
 
 Covers synthetic geometry (known answers), the demo nanoparticles image,
 edge cases, and physical-unit correctness.
 
 Run from project root:
-    .venv/bin/python -m tests.test_particles
+    .venv/bin/python -m tests.test_grains
 """
 
 import sys
@@ -100,7 +100,7 @@ def test_threshold_full_mask():
 def test_single_circle_area():
     """A single filled circle — verify pixel count and physical area."""
     print("=== Test: Single circle area ===")
-    from backend.nodes.particle_analysis import GrainAnalysis
+    from backend.nodes.grain_analysis import GrainAnalysis
     node = GrainAnalysis()
 
     N = 200
@@ -118,35 +118,35 @@ def test_single_circle_area():
     field = make_field(data, xreal=XREAL, yreal=XREAL)
     table, = node.process(field, mask=mask, min_size=1)
 
-    assert len(table) == 1, f"Expected 1 particles, got {len(table)}"
-    particles = table[0]
+    assert len(table) == 1, f"Expected 1 grain, got {len(table)}"
+    grain = table[0]
 
     # Pixel area of a discrete circle: should be close to π r²
     expected_px = np.pi * r ** 2
-    assert abs(particles["area_px"] - expected_px) / expected_px < 0.02, \
-        f"area_px={particles['area_px']}, expected≈{expected_px:.0f}"
+    assert abs(grain["area_px"] - expected_px) / expected_px < 0.02, \
+        f"area_px={grain['area_px']}, expected≈{expected_px:.0f}"
 
     # Physical area
     pixel_area = (XREAL / N) ** 2
-    expected_m2 = particles["area_px"] * pixel_area
-    assert abs(particles["area_m2"] - expected_m2) < 1e-20, \
-        f"area_m2 mismatch: {particles['area_m2']} vs {expected_m2}"
+    expected_m2 = grain["area_px"] * pixel_area
+    assert abs(grain["area_m2"] - expected_m2) < 1e-20, \
+        f"area_m2 mismatch: {grain['area_m2']} vs {expected_m2}"
 
     # Equivalent diameter should be close to 2r in physical units
     expected_diam = 2 * r * (XREAL / N)
-    assert abs(particles["equiv_diam_m"] - expected_diam) / expected_diam < 0.02, \
-        f"equiv_diam={particles['equiv_diam_m']:.3e}, expected≈{expected_diam:.3e}"
+    assert abs(grain["equiv_diam_m"] - expected_diam) / expected_diam < 0.02, \
+        f"equiv_diam={grain['equiv_diam_m']:.3e}, expected≈{expected_diam:.3e}"
 
     # Heights
-    assert abs(particles["mean_height"] - 5.0) < 1e-10
-    assert abs(particles["max_height"] - 5.0) < 1e-10
+    assert abs(grain["mean_height"] - 5.0) < 1e-10
+    assert abs(grain["max_height"] - 5.0) < 1e-10
     print("  PASS\n")
 
 
-def test_multiple_particles_separation():
-    """Three well-separated particles of different sizes — check each is reported."""
-    print("=== Test: Multiple particles separation ===")
-    from backend.nodes.particle_analysis import GrainAnalysis
+def test_multiple_grains_separation():
+    """Three well-separated grains of different sizes — check each is reported."""
+    print("=== Test: Multiple grains separation ===")
+    from backend.nodes.grain_analysis import GrainAnalysis
     node = GrainAnalysis()
 
     N = 128
@@ -168,7 +168,7 @@ def test_multiple_particles_separation():
     field = make_field(data)
     table, = node.process(field, mask=mask, min_size=1)
 
-    assert len(table) == 3, f"Expected 3 particles, got {len(table)}"
+    assert len(table) == 3, f"Expected 3 grains, got {len(table)}"
 
     table.sort(key=lambda r: r["area_px"], reverse=True)
     assert table[0]["area_px"] == 400   # 20×20
@@ -182,24 +182,24 @@ def test_multiple_particles_separation():
 
 
 def test_min_size_filtering():
-    """min_size should exclude particles smaller than the threshold."""
+    """min_size should exclude grains smaller than the threshold."""
     print("=== Test: min_size filtering ===")
-    from backend.nodes.particle_analysis import GrainAnalysis
+    from backend.nodes.grain_analysis import GrainAnalysis
     node = GrainAnalysis()
 
     N = 64
     data = np.zeros((N, N))
     mask = np.zeros((N, N), dtype=np.uint8)
 
-    # Large particles: 15×15 = 225 px
+    # Large grains: 15x15 = 225 px
     data[5:20, 5:20] = 1.0
     mask[5:20, 5:20] = 255
 
-    # Medium particles: 8×8 = 64 px
+    # Medium grains: 8x8 = 64 px
     data[30:38, 30:38] = 1.0
     mask[30:38, 30:38] = 255
 
-    # Tiny particles: 3×3 = 9 px
+    # Tiny grains: 3x3 = 9 px
     data[50:53, 50:53] = 1.0
     mask[50:53, 50:53] = 255
 
@@ -224,16 +224,16 @@ def test_min_size_filtering():
     print("  PASS\n")
 
 
-def test_particles_bounding_box():
-    """Bounding box should match the particles extents."""
+def test_grains_bounding_box():
+    """Bounding box should match the grain extents."""
     print("=== Test: Grain bounding box ===")
-    from backend.nodes.particle_analysis import GrainAnalysis
+    from backend.nodes.grain_analysis import GrainAnalysis
     node = GrainAnalysis()
 
     N = 64
     data = np.zeros((N, N))
     mask = np.zeros((N, N), dtype=np.uint8)
-    # Place a particles at rows 20:35, cols 10:45
+    # Place a grain at rows 20:35, cols 10:45
     data[20:35, 10:45] = 2.0
     mask[20:35, 10:45] = 255
 
@@ -247,10 +247,10 @@ def test_particles_bounding_box():
     print("  PASS\n")
 
 
-def test_empty_mask_produces_no_particles():
-    """An all-zero mask should yield zero particles."""
-    print("=== Test: Empty mask → no particles ===")
-    from backend.nodes.particle_analysis import GrainAnalysis
+def test_empty_mask_produces_no_grains():
+    """An all-zero mask should yield zero grains."""
+    print("=== Test: Empty mask -> no grains ===")
+    from backend.nodes.grain_analysis import GrainAnalysis
     node = GrainAnalysis()
 
     field = make_field(np.ones((64, 64)))
@@ -261,10 +261,10 @@ def test_empty_mask_produces_no_particles():
     print("  PASS\n")
 
 
-def test_particles_at_image_edge():
-    """A particles touching the image border should still be detected."""
+def test_grains_at_image_edge():
+    """A grain touching the image border should still be detected."""
     print("=== Test: Grain at image edge ===")
-    from backend.nodes.particle_analysis import GrainAnalysis
+    from backend.nodes.grain_analysis import GrainAnalysis
     node = GrainAnalysis()
 
     N = 64
@@ -282,11 +282,11 @@ def test_particles_at_image_edge():
     print("  PASS\n")
 
 
-def test_adjacent_particles_connectivity():
-    """Two diagonally-touching blocks should be separate particles
+def test_adjacent_grains_connectivity():
+    """Two diagonally-touching blocks should be separate grains
     (scipy.ndimage.label uses 4-connectivity by default)."""
-    print("=== Test: Diagonal adjacency → separate particles ===")
-    from backend.nodes.particle_analysis import GrainAnalysis
+    print("=== Test: Diagonal adjacency -> separate grains ===")
+    from backend.nodes.grain_analysis import GrainAnalysis
     node = GrainAnalysis()
 
     N = 32
@@ -305,7 +305,7 @@ def test_adjacent_particles_connectivity():
     table, = node.process(field, mask=mask, min_size=1)
     # Default label() uses structure that connects diagonals? Let's verify.
     # scipy.ndimage.label default is cross-shaped (no diagonals) for 2D
-    assert len(table) == 2, f"Expected 2 separate particles, got {len(table)}"
+    assert len(table) == 2, f"Expected 2 separate grains, got {len(table)}"
     print("  PASS\n")
 
 
@@ -315,17 +315,17 @@ def test_adjacent_particles_connectivity():
 
 def test_pipeline_synthetic():
     """Full pipeline on a synthetic image with known geometry."""
-    print("=== Test: Full pipeline on synthetic particles ===")
+    print("=== Test: Full pipeline on synthetic grains ===")
     from backend.nodes.threshold_mask import ThresholdMask
-    from backend.nodes.particle_analysis import GrainAnalysis
+    from backend.nodes.grain_analysis import GrainAnalysis
 
     N = 200
     XREAL = 10e-6  # 10 µm
     rng = np.random.default_rng(99)
 
-    # Background at 0 with small noise, particles as raised bumps
+    # Background at 0 with small noise, grains as raised bumps
     bg = rng.normal(0, 0.1, (N, N))
-    particles = np.zeros((N, N))
+    grains = np.zeros((N, N))
 
     yy, xx = np.mgrid[0:N, 0:N]
 
@@ -338,32 +338,32 @@ def test_pipeline_synthetic():
     ]
     for cx, cy, r, h in specs:
         inside = ((xx - cx) ** 2 + (yy - cy) ** 2) <= r ** 2
-        particles[inside] = h
+        grains[inside] = h
 
-    data = bg + particles
+    data = bg + grains
     field = make_field(data, xreal=XREAL, yreal=XREAL)
 
     # Step 1: threshold
     thresh = ThresholdMask()
     mask, = thresh.process(field, method="absolute", threshold=1.0, direction="above")
 
-    # Particles are well above noise, so mask should capture all 5
-    assert mask.max() == 255, "No particles detected"
+    # Grains are well above noise, so mask should capture all 5
+    assert mask.max() == 255, "No grains detected"
 
-    # Step 2: particles analysis
+    # Step 2: grain analysis
     ga = GrainAnalysis()
     table, = ga.process(field, mask=mask, min_size=5)
 
-    assert len(table) == 5, f"Expected 5 particles, got {len(table)}"
+    assert len(table) == 5, f"Expected 5 grains, got {len(table)}"
 
     # Verify that detected areas are in the right ballpark
     table.sort(key=lambda r: r["area_px"], reverse=True)
     expected_areas = sorted([np.pi * r ** 2 for _, _, r, _ in specs], reverse=True)
 
-    for particles, expected_px in zip(table, expected_areas):
-        ratio = particles["area_px"] / expected_px
+    for grain, expected_px in zip(table, expected_areas):
+        ratio = grain["area_px"] / expected_px
         assert 0.85 < ratio < 1.15, \
-            f"particles area_px={particles['area_px']}, expected≈{expected_px:.0f}, ratio={ratio:.2f}"
+            f"grain area_px={grain['area_px']}, expected≈{expected_px:.0f}, ratio={ratio:.2f}"
 
     print("  PASS\n")
 
@@ -373,7 +373,7 @@ def test_pipeline_demo_image():
     print("=== Test: Full pipeline on demo nanoparticles.npy ===")
     from pathlib import Path
     from backend.nodes.threshold_mask import ThresholdMask
-    from backend.nodes.particle_analysis import GrainAnalysis
+    from backend.nodes.grain_analysis import GrainAnalysis
     from backend.runtime_paths import demo_dir
 
     npy_path = demo_dir() / "nanoparticles.npy"
@@ -385,31 +385,31 @@ def test_pipeline_demo_image():
     # The demo image is a 5 µm × 5 µm scan
     field = make_field(data, xreal=5e-6, yreal=5e-6)
 
-    # Threshold to find particles (they are raised above background)
+    # Threshold to find grains (they are raised above background)
     thresh = ThresholdMask()
     mask, = thresh.process(field, method="otsu", threshold=0.0, direction="above")
 
-    # Should detect particles
-    assert mask.max() == 255, "No particles found in demo image"
+    # Should detect grains
+    assert mask.max() == 255, "No grains found in demo image"
     particle_fraction = (mask == 255).sum() / mask.size
     assert 0.01 < particle_fraction < 0.5, \
         f"Suspicious particle fraction: {particle_fraction:.3f}"
-    print(f"  Mask: {particle_fraction*100:.1f}% of pixels are particles")
+    print(f"  Mask: {particle_fraction*100:.1f}% of pixels are grains")
 
     # Grain analysis
     ga = GrainAnalysis()
     table, = ga.process(field, mask=mask, min_size=20)
 
-    assert len(table) > 0, "No particles detected"
-    print(f"  Found {len(table)} particles (min_size=20)")
+    assert len(table) > 0, "No grains detected"
+    print(f"  Found {len(table)} grains (min_size=20)")
 
-    # Sanity checks on particles properties
-    for particles in table:
-        assert particles["area_px"] >= 20
-        assert particles["area_m2"] > 0
-        assert particles["equiv_diam_m"] > 0
-        assert particles["max_height"] >= particles["mean_height"]
-        assert particles["mean_height"] > 0
+    # Sanity checks on grain properties
+    for grain in table:
+        assert grain["area_px"] >= 20
+        assert grain["area_m2"] > 0
+        assert grain["equiv_diam_m"] > 0
+        assert grain["max_height"] >= grain["mean_height"]
+        assert grain["mean_height"] > 0
 
     # Physical size sanity: equivalent diameters should be in the nm–µm range
     diams_nm = [g["equiv_diam_m"] * 1e9 for g in table]
@@ -433,15 +433,15 @@ if __name__ == "__main__":
 
     # GrainAnalysis
     test_single_circle_area()
-    test_multiple_particles_separation()
+    test_multiple_grains_separation()
     test_min_size_filtering()
-    test_particles_bounding_box()
-    test_empty_mask_produces_no_particles()
-    test_particles_at_image_edge()
-    test_adjacent_particles_connectivity()
+    test_grains_bounding_box()
+    test_empty_mask_produces_no_grains()
+    test_grains_at_image_edge()
+    test_adjacent_grains_connectivity()
 
     # End-to-end pipeline
     test_pipeline_synthetic()
     test_pipeline_demo_image()
 
-    print("All particles tests passed!")
+    print("All grain tests passed!")
diff --git a/tests/test_nodes.py b/tests/test_nodes.py
index 6980828..46f7e27 100644
--- a/tests/test_nodes.py
+++ b/tests/test_nodes.py
@@ -523,6 +523,31 @@ def test_plane_level():
     # The signal should remain (correlation with original sine)
     corr = np.corrcoef(result.data.ravel(), signal.ravel())[0, 1]
     assert corr > 0.98, f"Signal correlation after leveling: {corr}"
+
+    yy_px, xx_px = np.mgrid[0:N, 0:N]
+
+    def fit_pixel_plane(data_in: np.ndarray, region: np.ndarray) -> tuple[float, float, float]:
+        A = np.column_stack([
+            np.ones(int(np.count_nonzero(region)), dtype=np.float64),
+            xx_px[region].astype(np.float64),
+            yy_px[region].astype(np.float64),
+        ])
+        coeffs, _, _, _ = np.linalg.lstsq(A, data_in[region].ravel().astype(np.float64), rcond=None)
+        return float(coeffs[0]), float(coeffs[1]), float(coeffs[2])
+
+    mask = np.zeros((N, N), dtype=np.uint8)
+    mask[20:44, 22:46] = 255
+    feature = np.zeros((N, N), dtype=np.float64)
+    feature[mask > 0] = 35.0
+    masked_field = make_field(data=100 * x + 50 * y + feature)
+
+    unmasked, = node.process(masked_field)
+    masked, = node.process(masked_field, masking="exclude", mask=mask)
+
+    outside = mask == 0
+    _, unmasked_bx, unmasked_by = fit_pixel_plane(unmasked.data, outside)
+    _, masked_bx, masked_by = fit_pixel_plane(masked.data, outside)
+    assert np.hypot(masked_bx, masked_by) < np.hypot(unmasked_bx, unmasked_by) * 1e-3
     print("  PASS\n")
 
 
@@ -1261,10 +1286,10 @@ def test_mask_invert():
     print("  PASS\n")
 
 
-def test_mask_combine():
-    print("=== Test: MaskCombine ===")
-    from backend.nodes.mask_combine import MaskCombine
-    node = MaskCombine()
+def test_mask_operations():
+    print("=== Test: MaskOperations ===")
+    from backend.nodes.mask_operations import MaskOperations
+    node = MaskOperations()
 
     # Two overlapping squares
     a = np.zeros((64, 64), dtype=np.uint8)
@@ -1291,12 +1316,26 @@ def test_mask_combine():
     assert result_xor[35, 35] == 255  # b-only
     assert result_xor[25, 25] == 0    # overlap excluded
 
-    # Subtract — a minus b
-    result_sub, = node.process(a, b, operation="subtract")
+    # A minus B
+    result_sub, = node.process(a, b, operation="a_minus_b")
     assert result_sub[15, 15] == 255  # a-only kept
     assert result_sub[25, 25] == 0    # overlap removed
     assert result_sub[35, 35] == 0    # b-only not included
 
+    # NAND — everything except overlap
+    result_nand, = node.process(a, b, operation="nand")
+    assert result_nand[15, 15] == 255
+    assert result_nand[35, 35] == 255
+    assert result_nand[25, 25] == 0
+    assert result_nand[5, 5] == 255
+
+    # XNOR — overlap plus shared background
+    result_xnor, = node.process(a, b, operation="xnor")
+    assert result_xnor[25, 25] == 255
+    assert result_xnor[5, 5] == 255
+    assert result_xnor[15, 15] == 0
+    assert result_xnor[35, 35] == 0
+
     print("  PASS\n")
 
 
@@ -1347,17 +1386,17 @@ def test_draw_mask():
     print("  PASS\n")
 
 
-def test_particle_analysis():
-    print("=== Test: ParticleAnalysis ===")
-    from backend.nodes.particle_analysis import ParticleAnalysis
-    node = ParticleAnalysis()
+def test_grain_analysis():
+    print("=== Test: GrainAnalysis ===")
+    from backend.nodes.grain_analysis import GrainAnalysis
+    node = GrainAnalysis()
 
-    # Create a field with two distinct particles
+    # Create a field with two distinct grains
     N = 64
     data = np.zeros((N, N))
-    # Particle 1: 10x10 block at top-left with height 5
+    # Grain 1: 10x10 block at top-left with height 5
     data[5:15, 5:15] = 5.0
-    # Particle 2: 8x8 block at bottom-right with height 3
+    # Grain 2: 8x8 block at bottom-right with height 3
     data[45:53, 45:53] = 3.0
     field = make_field(data=data, xreal=1e-6, yreal=1e-6)
 
@@ -1367,7 +1406,7 @@ def test_particle_analysis():
     mask[45:53, 45:53] = 255
 
     table, = node.process(field, mask=mask, min_size=10)
-    assert len(table) == 2, f"Expected 2 particles, got {len(table)}"
+    assert len(table) == 2, f"Expected 2 grains, got {len(table)}"
 
     # Sort by area descending
     table.sort(key=lambda r: r["area_px"], reverse=True)
@@ -1381,7 +1420,7 @@ def test_particle_analysis():
     assert table[0]["mean_height_unit"] == "m"
     assert table[0]["max_height_unit"] == "m"
 
-    # min_size filtering: only keep particles >= 80 px
+    # min_size filtering: only keep grains >= 80 px
     table_filtered, = node.process(field, mask=mask, min_size=80)
     assert len(table_filtered) == 1
     assert table_filtered[0]["area_px"] == 100
@@ -3140,11 +3179,11 @@ if __name__ == "__main__":
     test_threshold_mask()
     test_mask_morphology()
     test_mask_invert()
-    test_mask_combine()
+    test_mask_operations()
     test_draw_mask()
 
     # Grains
-    test_particle_analysis()
+    test_grain_analysis()
 
     # I/O
     test_load_file()