From 5de93e6c4d0ee183b2b8070c48e074d5ca4d2bab Mon Sep 17 00:00:00 2001 From: matei jordache Date: Sat, 4 Apr 2026 00:25:53 -0700 Subject: [PATCH] low pri features --- GWYDDION_FEATURE_GAP.md | 38 +-- backend/node_menu.py | 15 + backend/nodes/calibration.py | 143 ++++++++ backend/nodes/displacement_field.py | 95 ++++++ backend/nodes/distribution_coercion.py | 81 +++++ backend/nodes/dwt_anisotropy.py | 198 +++++++++++ backend/nodes/grain_visualization.py | 238 +++++++++++++ backend/nodes/logistic_classification.py | 229 +++++++++++++ backend/nodes/mark_disconnected.py | 71 ++++ backend/nodes/mask_noisify.py | 84 +++++ backend/nodes/mask_shift.py | 98 ++++++ backend/nodes/neural_classification.py | 204 ++++++++++++ backend/nodes/pixel_classification.py | 209 ++++++++++++ backend/nodes/presentation_ops.py | 88 +++++ backend/nodes/psf_estimation.py | 177 ++++++++++ backend/nodes/super_resolution.py | 126 +++++++ backend/nodes/tip_shape_estimate.py | 351 ++++++++++++++++++++ docs/nodes/Calibration.md | 37 +++ docs/nodes/DWT Anisotropy.md | 32 ++ docs/nodes/Displacement Field.md | 34 ++ docs/nodes/Distribution Coercion.md | 31 ++ docs/nodes/Grain Visualization.md | 34 ++ docs/nodes/Logistic Classification.md | 37 +++ docs/nodes/Mark Disconnected.md | 30 ++ docs/nodes/Mask Noisify.md | 32 ++ docs/nodes/Mask Shift.md | 31 ++ docs/nodes/Neural Classification.md | 35 ++ docs/nodes/PSF Estimation.md | 35 ++ docs/nodes/Pixel Classification.md | 34 ++ docs/nodes/Presentation Ops.md | 31 ++ docs/nodes/Super Resolution.md | 32 ++ docs/nodes/Tip Shape Estimate.md | 33 ++ tests/node_tests/calibration.py | 77 +++++ tests/node_tests/displacement_field.py | 41 +++ tests/node_tests/distribution_coercion.py | 61 ++++ tests/node_tests/dwt_anisotropy.py | 64 ++++ tests/node_tests/grain_visualization.py | 64 ++++ tests/node_tests/logistic_classification.py | 111 +++++++ tests/node_tests/mark_disconnected.py | 43 +++ tests/node_tests/mask_noisify.py | 50 +++ tests/node_tests/mask_shift.py | 74 +++++ tests/node_tests/neural_classification.py | 72 ++++ tests/node_tests/pixel_classification.py | 49 +++ tests/node_tests/presentation_ops.py | 43 +++ tests/node_tests/psf_estimation.py | 56 ++++ tests/node_tests/super_resolution.py | 51 +++ tests/node_tests/tip_shape_estimate.py | 86 +++++ 47 files changed, 3866 insertions(+), 19 deletions(-) create mode 100644 backend/nodes/calibration.py create mode 100644 backend/nodes/displacement_field.py create mode 100644 backend/nodes/distribution_coercion.py create mode 100644 backend/nodes/dwt_anisotropy.py create mode 100644 backend/nodes/grain_visualization.py create mode 100644 backend/nodes/logistic_classification.py create mode 100644 backend/nodes/mark_disconnected.py create mode 100644 backend/nodes/mask_noisify.py create mode 100644 backend/nodes/mask_shift.py create mode 100644 backend/nodes/neural_classification.py create mode 100644 backend/nodes/pixel_classification.py create mode 100644 backend/nodes/presentation_ops.py create mode 100644 backend/nodes/psf_estimation.py create mode 100644 backend/nodes/super_resolution.py create mode 100644 backend/nodes/tip_shape_estimate.py create mode 100644 docs/nodes/Calibration.md create mode 100644 docs/nodes/DWT Anisotropy.md create mode 100644 docs/nodes/Displacement Field.md create mode 100644 docs/nodes/Distribution Coercion.md create mode 100644 docs/nodes/Grain Visualization.md create mode 100644 docs/nodes/Logistic Classification.md create mode 100644 docs/nodes/Mark Disconnected.md create mode 100644 docs/nodes/Mask Noisify.md create mode 100644 docs/nodes/Mask Shift.md create mode 100644 docs/nodes/Neural Classification.md create mode 100644 docs/nodes/PSF Estimation.md create mode 100644 docs/nodes/Pixel Classification.md create mode 100644 docs/nodes/Presentation Ops.md create mode 100644 docs/nodes/Super Resolution.md create mode 100644 docs/nodes/Tip Shape Estimate.md create mode 100644 tests/node_tests/calibration.py create mode 100644 tests/node_tests/displacement_field.py create mode 100644 tests/node_tests/distribution_coercion.py create mode 100644 tests/node_tests/dwt_anisotropy.py create mode 100644 tests/node_tests/grain_visualization.py create mode 100644 tests/node_tests/logistic_classification.py create mode 100644 tests/node_tests/mark_disconnected.py create mode 100644 tests/node_tests/mask_noisify.py create mode 100644 tests/node_tests/mask_shift.py create mode 100644 tests/node_tests/neural_classification.py create mode 100644 tests/node_tests/pixel_classification.py create mode 100644 tests/node_tests/presentation_ops.py create mode 100644 tests/node_tests/psf_estimation.py create mode 100644 tests/node_tests/super_resolution.py create mode 100644 tests/node_tests/tip_shape_estimate.py diff --git a/GWYDDION_FEATURE_GAP.md b/GWYDDION_FEATURE_GAP.md index 8476d14..d4aad33 100644 --- a/GWYDDION_FEATURE_GAP.md +++ b/GWYDDION_FEATURE_GAP.md @@ -107,23 +107,23 @@ All features from the original gap analysis are implemented: ### Lower Priority — Specialized or niche -| # | Feature | Gwyddion Source | Description | -|---|---------|---------------|-------------| -| 77 | Mark Disconnected Regions | mark_disconn.c | Mask topologically isolated surface regions using threshold and radius criteria. | -| 78 | Mask Shift | mask_shift.c | Translate mask by pixel offset in any direction. | -| 79 | Mask Noisify | mask_noisify.c | Add random perturbation to mask boundaries. Useful for testing mask sensitivity. | -| 80 | DWT Anisotropy | dwtanisotropy.c | Quantify surface anisotropy using discrete wavelet transform decomposition. | -| 81 | Displacement Field | displfield.c | Distort images using displacement fields (Gaussian, tear, image-based). Simulates scanning artifacts. | -| 82 | Pixel Classification | classify.c | Classify pixels into categories using decision trees on height, slope, and curvature criteria. | -| 83 | Neural Network Classification | neural.c | Train and apply neural networks for pixel-level feature classification. | -| 84 | Logistic Classification | logistic.c | Classify features using logistic regression on Gaussian derivative features. | -| 85 | Super-Resolution | superresolution.c | Combine multiple aligned low-resolution scans to produce a higher-resolution image. | -| 86 | PSF Estimation | psf.c, psf-fit.c | Estimate and fit point spread functions from image features for deconvolution. | -| 87 | Tip Shape from Features | tipshape.c | Estimate SPM tip shape from known calibration feature convolutions. | -| 88 | Presentation Ops | presentationops.c | Manage presentation overlays (extract, attach, remove presentation layers). | -| 89 | Calibration Coefficients | calcoefs_*.c, calibrate.c | Load, create, and apply lateral/height calibration corrections. | -| 90 | Distribution Coercion | coerce.c | Transform data distribution to match target (uniform, Gaussian, custom). | -| 91 | Grain Selection Visualization | grain_makesel.c | Visualize grains as discs, circles, or bounding boxes for selection. | +| # | Feature | Gwyddion Source | tono Node | Status | +|---|---------|---------------|-----------|--------| +| 77 | Mark Disconnected Regions | mark_disconn.c | MarkDisconnected | **DONE** | +| 78 | Mask Shift | mask_shift.c | MaskShift | **DONE** | +| 79 | Mask Noisify | mask_noisify.c | MaskNoisify | **DONE** | +| 80 | DWT Anisotropy | dwtanisotropy.c | DWTAnisotropy | **DONE** | +| 81 | Displacement Field | displfield.c | DisplacementField | **DONE** | +| 82 | Pixel Classification | classify.c | PixelClassification | **DONE** | +| 83 | Neural Network Classification | neural.c | NeuralClassification | **DONE** | +| 84 | Logistic Classification | logistic.c | LogisticClassification | **DONE** | +| 85 | Super-Resolution | superresolution.c | SuperResolution | **DONE** | +| 86 | PSF Estimation | psf.c, psf-fit.c | PSFEstimation | **DONE** | +| 87 | Tip Shape from Features | tipshape.c | TipShapeEstimate | **DONE** | +| 88 | Presentation Ops | presentationops.c | PresentationOps | **DONE** | +| 89 | Calibration Coefficients | calcoefs_*.c, calibrate.c | Calibration | **DONE** | +| 90 | Distribution Coercion | coerce.c | DistributionCoercion | **DONE** | +| 91 | Grain Selection Visualization | grain_makesel.c | GrainVisualization | **DONE** | ### Synthesis — Additional surface generation patterns @@ -181,8 +181,8 @@ Gwyddion supports 155+ file format modules. tono currently handles a smaller set | High Value (41–51) | 11 | **All 11 done** | | Medium Value (52–70) | 19 | **All 19 done** | | SPM Mode-Specific (71–76) | 6 | **All 6 done** | -| Lower Priority (77–91) | 15 | Pending | +| Lower Priority (77–91) | 15 | **All 15 done** | | Synthesis Patterns (92–113) | 22 | **All 22 done** | | File Formats | 10+ | Pending | -**97 of 98 tracked features implemented.** 15 remaining gaps identified (lower priority + file formats). +**112 of 113 tracked features implemented.** Only file format support gaps remain. diff --git a/backend/node_menu.py b/backend/node_menu.py index 0884f31..0566bb8 100644 --- a/backend/node_menu.py +++ b/backend/node_menu.py @@ -35,6 +35,7 @@ MENU_LAYOUT: dict[str, list[str]] = { "Save", "SaveImage", "Shade", + "PresentationOps", ], "Overlay": [ "Markup", @@ -56,6 +57,7 @@ MENU_LAYOUT: dict[str, list[str]] = { "PixelBinning", "ExtendPad", "FieldArithmetic", + "DisplacementField", ], "Level & Correct": [ "FixZero", @@ -72,6 +74,8 @@ MENU_LAYOUT: dict[str, list[str]] = { "ScanLineReorder", "Tilt", "WrapValue", + "DistributionCoercion", + "Calibration", ], "Filter": [ "GaussianFilter", @@ -98,6 +102,7 @@ MENU_LAYOUT: dict[str, list[str]] = { "LogPolarPSDF", "FrequencySplit", "CrossCorrelate", + "SuperResolution", ], "Measure": [ "CrossSection", @@ -116,6 +121,7 @@ MENU_LAYOUT: dict[str, list[str]] = { "MultipleProfiles", "StraightenPath", "RelateFields", + "DWTAnisotropy", ], "Detect": [ "FeatureDetection", @@ -130,6 +136,9 @@ MENU_LAYOUT: dict[str, list[str]] = { "LateralForceSim", "SEMSimulation", "SMMAnalysis", + "PixelClassification", + "NeuralClassification", + "LogisticClassification", ], "Mask": [ "DrawMask", @@ -139,6 +148,9 @@ MENU_LAYOUT: dict[str, list[str]] = { "MaskMorphology", "MaskInvert", "MaskOperations", + "MarkDisconnected", + "MaskShift", + "MaskNoisify", ], "Grains": [ "GrainDistanceTransform", @@ -150,11 +162,14 @@ MENU_LAYOUT: dict[str, list[str]] = { "LevelGrains", "GrainEdge", "GrainCross", + "GrainVisualization", ], "Tip": [ "TipModel", "TipDeconvolution", "BlindTipEstimate", + "TipShapeEstimate", + "PSFEstimation", ], } diff --git a/backend/nodes/calibration.py b/backend/nodes/calibration.py new file mode 100644 index 0000000..cb2486a --- /dev/null +++ b/backend/nodes/calibration.py @@ -0,0 +1,143 @@ +"""Calibration — apply lateral and height calibration corrections.""" + +from __future__ import annotations + +import numpy as np + +from backend.node_registry import register_node +from backend.data_types import DataField + + +@register_node(display_name="Calibration") +class Calibration: + @classmethod + def INPUT_TYPES(cls): + return { + "required": { + "field": ("DATA_FIELD",), + "xy_mode": (["keep", "set_size", "scale"], {"default": "keep"}), + "z_mode": (["keep", "set_range", "scale", "offset"], {"default": "keep"}), + "xreal_new": ("FLOAT", { + "default": 1e-6, + "min": 1e-12, + "max": 1.0, + "step": 1e-9, + "show_when_widget_value": {"xy_mode": ["set_size"]}, + }), + "yreal_new": ("FLOAT", { + "default": 1e-6, + "min": 1e-12, + "max": 1.0, + "step": 1e-9, + "show_when_widget_value": {"xy_mode": ["set_size"]}, + }), + "xy_scale": ("FLOAT", { + "default": 1.0, + "min": 0.001, + "max": 1000.0, + "step": 0.001, + "show_when_widget_value": {"xy_mode": ["scale"]}, + }), + "z_min": ("FLOAT", { + "default": 0.0, + "min": -1e-3, + "max": 1e-3, + "step": 1e-12, + "show_when_widget_value": {"z_mode": ["set_range"]}, + }), + "z_max": ("FLOAT", { + "default": 1e-9, + "min": -1e-3, + "max": 1e-3, + "step": 1e-12, + "show_when_widget_value": {"z_mode": ["set_range"]}, + }), + "z_scale": ("FLOAT", { + "default": 1.0, + "min": 0.001, + "max": 1000.0, + "step": 0.001, + "show_when_widget_value": {"z_mode": ["scale"]}, + }), + "z_offset": ("FLOAT", { + "default": 0.0, + "min": -1e-3, + "max": 1e-3, + "step": 1e-12, + "show_when_widget_value": {"z_mode": ["offset"]}, + }), + "xy_unit": ("STRING", {"default": ""}), + "z_unit": ("STRING", {"default": ""}), + } + } + + OUTPUTS = ( + ('DATA_FIELD', 'result'), + ) + FUNCTION = "process" + + DESCRIPTION = ( + "Apply lateral and height calibration corrections to a DATA_FIELD. " + "Lateral mode can set explicit physical size or scale uniformly. " + "Height mode can rescale to a target range, multiply by a factor, " + "or add a constant offset. Optionally override the XY or Z unit strings. " + "Equivalent to Gwyddion's calibrate functionality." + ) + + def process( + self, + field: DataField, + xy_mode: str, + z_mode: str, + xreal_new: float, + yreal_new: float, + xy_scale: float, + z_min: float, + z_max: float, + z_scale: float, + z_offset: float, + xy_unit: str, + z_unit: str, + ) -> tuple: + data = np.asarray(field.data, dtype=np.float64).copy() + xreal = float(field.xreal) + yreal = float(field.yreal) + si_unit_xy = field.si_unit_xy + si_unit_z = field.si_unit_z + + # --- lateral calibration --- + if xy_mode == "set_size": + xreal = float(xreal_new) + yreal = float(yreal_new) + elif xy_mode == "scale": + xreal *= float(xy_scale) + yreal *= float(xy_scale) + # "keep" → no change + + # --- height calibration --- + if z_mode == "set_range": + cur_min = float(data.min()) + cur_max = float(data.max()) + if cur_max > cur_min: + data = float(z_min) + (data - cur_min) * (float(z_max) - float(z_min)) / (cur_max - cur_min) + else: + data[:] = float(z_min) + elif z_mode == "scale": + data *= float(z_scale) + elif z_mode == "offset": + data += float(z_offset) + # "keep" → no change + + # --- unit overrides --- + if xy_unit: + si_unit_xy = xy_unit + if z_unit: + si_unit_z = z_unit + + return (field.replace( + data=data, + xreal=xreal, + yreal=yreal, + si_unit_xy=si_unit_xy, + si_unit_z=si_unit_z, + ),) diff --git a/backend/nodes/displacement_field.py b/backend/nodes/displacement_field.py new file mode 100644 index 0000000..9cbe329 --- /dev/null +++ b/backend/nodes/displacement_field.py @@ -0,0 +1,95 @@ +"""Displacement field — distort images using displacement fields.""" + +from __future__ import annotations + +import numpy as np +from scipy.ndimage import gaussian_filter, gaussian_filter1d, map_coordinates + +from backend.node_registry import register_node +from backend.data_types import DataField + + +@register_node(display_name="Displacement Field") +class DisplacementField: + @classmethod + def INPUT_TYPES(cls): + return { + "required": { + "field": ("DATA_FIELD",), + "method": (["gaussian_1d", "gaussian_2d", "tear"], {"default": "gaussian_1d"}), + "sigma": ("FLOAT", {"default": 5.0, "min": 0.1, "max": 100.0, "step": 0.1}), + "tau": ("FLOAT", {"default": 20.0, "min": 1.0, "max": 500.0, "step": 1.0}), + "density": ("FLOAT", { + "default": 0.02, + "min": 0.001, + "max": 0.25, + "step": 0.001, + "show_when_widget_value": {"method": ["tear"]}, + }), + "seed": ("INT", {"default": 42, "min": 0, "max": 999999}), + } + } + + OUTPUTS = ( + ('DATA_FIELD', 'result'), + ) + FUNCTION = "process" + + DESCRIPTION = ( + "Distort an image using synthetic displacement fields. " + "Supports 1D Gaussian (row-correlated), 2D Gaussian (fully correlated), " + "and tear (random horizontal tear lines) distortion modes. " + "Equivalent to Gwyddion's displfield.c module." + ) + + def process( + self, + field: DataField, + method: str, + sigma: float, + tau: float, + density: float, + seed: int, + ) -> tuple: + data = np.asarray(field.data, dtype=np.float64) + yres, xres = data.shape + rng = np.random.default_rng(seed) + + y_grid, x_grid = np.mgrid[:yres, :xres].astype(np.float64) + + if method == "gaussian_1d": + dx_1d = gaussian_filter1d(rng.standard_normal(xres), tau) * sigma + dx = np.tile(dx_1d, (yres, 1)) + dy = np.zeros_like(dx) + + elif method == "gaussian_2d": + dx = gaussian_filter(rng.standard_normal((yres, xres)), tau) * sigma + dy = gaussian_filter(rng.standard_normal((yres, xres)), tau) * sigma + + elif method == "tear": + dx = np.zeros((yres, xres), dtype=np.float64) + dy = np.zeros((yres, xres), dtype=np.float64) + + # Select tear rows based on density + tear_mask = rng.random(yres) < density + tear_rows = np.nonzero(tear_mask)[0] + + for row in tear_rows: + offset = rng.standard_normal() * sigma + # Apply offset that decays exponentially away from the tear line + for i in range(yres): + dist = abs(i - row) + dx[i] += offset * np.exp(-dist / max(tau, 1.0)) + + # Smooth the displacement to avoid sharp edges + for i in range(yres): + dx[i] = gaussian_filter1d(dx[i], tau / 4.0) + + else: + raise ValueError(f"Unknown method: {method!r}") + + coords_y = y_grid + dy + coords_x = x_grid + dx + result = map_coordinates(data, [coords_y, coords_x], mode='reflect', order=3) + + return (field.replace(data=result),) diff --git a/backend/nodes/distribution_coercion.py b/backend/nodes/distribution_coercion.py new file mode 100644 index 0000000..bcf3250 --- /dev/null +++ b/backend/nodes/distribution_coercion.py @@ -0,0 +1,81 @@ +"""Distribution coercion — transform data to match a target distribution.""" + +from __future__ import annotations + +import numpy as np +from math import ceil +from scipy.stats import norm + +from backend.node_registry import register_node +from backend.data_types import DataField + + +def _coerce_block(data: np.ndarray, distribution: str, n_levels: int) -> np.ndarray: + """Coerce a flat or 2-D block to the target distribution, returning same shape.""" + shape = data.shape + flat = data.ravel().astype(np.float64) + n_pixels = flat.size + if n_pixels == 0: + return data.copy() + + indices = np.argsort(flat, kind="mergesort") + + if distribution == "uniform": + target = np.linspace(float(flat.min()), float(flat.max()), n_pixels) + elif distribution == "gaussian": + eps = 0.5 / n_pixels + quantiles = np.linspace(eps, 1.0 - eps, n_pixels) + target = norm.ppf(quantiles) * float(flat.std()) + float(flat.mean()) + elif distribution == "levels": + n_levels = max(2, int(n_levels)) + level_values = np.linspace(float(flat.min()), float(flat.max()), n_levels) + target = np.repeat(level_values, ceil(n_pixels / n_levels))[:n_pixels] + else: + raise ValueError(f"Unknown distribution: {distribution}") + + result = np.empty_like(flat) + result[indices] = target + return result.reshape(shape) + + +@register_node(display_name="Distribution Coercion") +class DistributionCoercion: + @classmethod + def INPUT_TYPES(cls): + return { + "required": { + "field": ("DATA_FIELD",), + "distribution": (["uniform", "gaussian", "levels"], {"default": "uniform"}), + "n_levels": ("INT", { + "default": 4, + "min": 2, + "max": 1000, + "show_when_widget_value": {"distribution": ["levels"]}, + }), + "processing": (["field", "rows"], {"default": "field"}), + } + } + + OUTPUTS = ( + ('DATA_FIELD', 'result'), + ) + FUNCTION = "process" + + DESCRIPTION = ( + "Transform pixel values so their distribution matches a target shape " + "(uniform, Gaussian, or discrete levels) using rank-based reassignment. " + "Equivalent to Gwyddion's coerce.c module." + ) + + def process(self, field: DataField, distribution: str, n_levels: int, + processing: str) -> tuple: + data = np.asarray(field.data, dtype=np.float64) + + if processing == "rows": + result = np.empty_like(data) + for i in range(data.shape[0]): + result[i] = _coerce_block(data[i], distribution, n_levels) + else: + result = _coerce_block(data, distribution, n_levels) + + return (field.replace(data=result),) diff --git a/backend/nodes/dwt_anisotropy.py b/backend/nodes/dwt_anisotropy.py new file mode 100644 index 0000000..5671dae --- /dev/null +++ b/backend/nodes/dwt_anisotropy.py @@ -0,0 +1,198 @@ +"""DWT anisotropy — quantify surface anisotropy using wavelet decomposition.""" + +from __future__ import annotations + +import numpy as np + +from backend.data_types import DataField, RecordTable +from backend.node_registry import register_node + + +def _next_power_of_2(n: int) -> int: + """Return the smallest power of 2 >= n.""" + p = 1 + while p < n: + p <<= 1 + return p + + +def _pad_to_pow2(data: np.ndarray) -> np.ndarray: + """Pad *data* to the next power-of-2 dimensions using edge values.""" + rows, cols = data.shape + new_rows = _next_power_of_2(rows) + new_cols = _next_power_of_2(cols) + if new_rows == rows and new_cols == cols: + return data.copy() + out = np.zeros((new_rows, new_cols), dtype=np.float64) + out[:rows, :cols] = data + # edge-pad + if new_cols > cols: + out[:rows, cols:] = data[:, -1:] + if new_rows > rows: + out[rows:, :cols] = data[-1:, :] + if new_rows > rows and new_cols > cols: + out[rows:, cols:] = data[-1, -1] + return out + + +def _haar_decompose_2d(data: np.ndarray) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: + """ + One level of 2-D Haar wavelet decomposition. + + Returns (LL, LH, HL, HH) each of shape (rows//2, cols//2). + LL = (a+b+c+d)/2 approximation + LH = (a+b-c-d)/2 horizontal detail (captures vertical features) + HL = (a-b+c-d)/2 vertical detail (captures horizontal features) + HH = (a-b-c+d)/2 diagonal detail + """ + rows, cols = data.shape + a = data[0::2, 0::2] # top-left + b = data[0::2, 1::2] # top-right + c = data[1::2, 0::2] # bottom-left + d = data[1::2, 1::2] # bottom-right + + ll = (a + b + c + d) / 2.0 + lh = (a + b - c - d) / 2.0 + hl = (a - b + c - d) / 2.0 + hh = (a - b - c + d) / 2.0 + return ll, lh, hl, hh + + +def _compute_dwt_anisotropy( + data: np.ndarray, + n_levels: int, +) -> tuple[list[float], list[float], list[float], list[np.ndarray]]: + """ + Multi-level 2-D Haar decomposition with per-level energy ratios. + + Returns + ------- + x_energies : per-level sum(HL**2) + y_energies : per-level sum(LH**2) + ratios : per-level x_energy / y_energy + ratio_maps : per-level ratio arrays (at decomposition resolution) + """ + padded = _pad_to_pow2(np.asarray(data, dtype=np.float64)) + current = padded + + x_energies: list[float] = [] + y_energies: list[float] = [] + ratios: list[float] = [] + ratio_maps: list[np.ndarray] = [] + + for _ in range(n_levels): + if current.shape[0] < 2 or current.shape[1] < 2: + break + ll, lh, hl, hh = _haar_decompose_2d(current) + + x_energy = float(np.sum(hl ** 2)) + y_energy = float(np.sum(lh ** 2)) + ratio = x_energy / (y_energy + 1e-30) + + x_energies.append(x_energy) + y_energies.append(y_energy) + ratios.append(ratio) + + # Build a per-pixel ratio map at this level's resolution + hl_sq = hl ** 2 + lh_sq = lh ** 2 + level_map = hl_sq / (lh_sq + 1e-30) + ratio_maps.append(level_map) + + current = ll + + return x_energies, y_energies, ratios, ratio_maps + + +def _build_anisotropy_map( + ratio_maps: list[np.ndarray], + orig_rows: int, + orig_cols: int, +) -> np.ndarray: + """ + Combine per-level ratio maps into a single anisotropy map at + the original field resolution by upsampling and averaging. + """ + if not ratio_maps: + return np.ones((orig_rows, orig_cols), dtype=np.float64) + + from scipy.ndimage import zoom + + target = np.zeros((orig_rows, orig_cols), dtype=np.float64) + for level_map in ratio_maps: + zy = orig_rows / level_map.shape[0] + zx = orig_cols / level_map.shape[1] + upsampled = zoom(level_map, (zy, zx), order=1) + # zoom may produce shape off by one — trim or pad + upsampled = upsampled[:orig_rows, :orig_cols] + if upsampled.shape[0] < orig_rows or upsampled.shape[1] < orig_cols: + tmp = np.ones((orig_rows, orig_cols), dtype=np.float64) + tmp[: upsampled.shape[0], : upsampled.shape[1]] = upsampled + upsampled = tmp + target += upsampled + + target /= len(ratio_maps) + return target + + +@register_node(display_name="DWT Anisotropy") +class DWTAnisotropy: + @classmethod + def INPUT_TYPES(cls): + return { + "required": { + "field": ("DATA_FIELD",), + "n_levels": ( + "INT", + {"default": 4, "min": 1, "max": 10}, + ), + "ratio_threshold": ( + "FLOAT", + {"default": 0.2, "min": 0.001, "max": 10.0, "step": 0.01}, + ), + } + } + + OUTPUTS = ( + ('DATA_FIELD', 'anisotropy_map'), + ('RECORD_TABLE', 'statistics'), + ) + FUNCTION = "process" + + DESCRIPTION = ( + "Quantify surface anisotropy using a multi-level 2-D Haar wavelet decomposition. " + "At each level, horizontal (HL) and vertical (LH) detail energies are compared to " + "produce an X/Y energy ratio. Ratio > 1 indicates more horizontal features; " + "ratio < 1 indicates more vertical features. Equivalent to Gwyddion's dwtanisotropy.c." + ) + + def process( + self, + field: DataField, + n_levels: int, + ratio_threshold: float, + ) -> tuple: + data = np.asarray(field.data, dtype=np.float64) + orig_rows, orig_cols = data.shape + + x_energies, y_energies, ratios, ratio_maps = _compute_dwt_anisotropy( + data, int(n_levels), + ) + + # Build per-pixel anisotropy map + aniso_map = _build_anisotropy_map(ratio_maps, orig_rows, orig_cols) + aniso_field = field.replace(data=aniso_map) + + # Build statistics table + rows = [] + for i, (xe, ye, r) in enumerate(zip(x_energies, y_energies, ratios)): + rows.append({ + "level": i + 1, + "x_energy": float(xe), + "y_energy": float(ye), + "ratio": float(r), + "anisotropic": abs(r - 1.0) > float(ratio_threshold), + }) + stats = RecordTable(rows) + + return (aniso_field, stats) diff --git a/backend/nodes/grain_visualization.py b/backend/nodes/grain_visualization.py new file mode 100644 index 0000000..f2e7c6d --- /dev/null +++ b/backend/nodes/grain_visualization.py @@ -0,0 +1,238 @@ +"""Grain visualization — visualize grains as geometric shapes.""" + +from __future__ import annotations + +import numpy as np +from scipy.ndimage import label, find_objects, distance_transform_edt + +from backend.node_registry import register_node +from backend.data_types import DataField +from backend.nodes.helpers import mask_to_bool, bool_to_mask + + +def _grain_centroid(grain_mask: np.ndarray, slc: tuple[slice, slice]) -> tuple[float, float]: + """Return (cy, cx) centroid of a grain within its bounding slice.""" + ys, xs = np.where(grain_mask[slc]) + cy = float(ys.mean()) + slc[0].start + cx = float(xs.mean()) + slc[1].start + return cy, cx + + +def _grain_inscribed_radius(grain_mask: np.ndarray, slc: tuple[slice, slice]) -> float: + """Return the inscribed disc radius for a grain region.""" + region = grain_mask[slc] + if not np.any(region): + return 0.0 + dt = distance_transform_edt(region) + return float(dt.max()) + + +def _grain_inertia(grain_mask: np.ndarray, slc: tuple[slice, slice]) -> tuple[float, float, float]: + """Return (semi_major, semi_minor, angle_rad) from the inertia tensor.""" + ys, xs = np.where(grain_mask[slc]) + cy_local = ys.mean() + cx_local = xs.mean() + dy = ys - cy_local + dx = xs - cx_local + n = len(ys) + if n < 2: + return 1.0, 1.0, 0.0 + + Ixx = np.sum(dy * dy) / n + Iyy = np.sum(dx * dx) / n + Ixy = -np.sum(dx * dy) / n + + # Eigenvalues of the 2x2 inertia tensor + mean_I = (Ixx + Iyy) / 2.0 + diff_I = (Ixx - Iyy) / 2.0 + discriminant = max(0.0, diff_I ** 2 + Ixy ** 2) + sqrt_disc = np.sqrt(discriminant) + + lambda1 = mean_I + sqrt_disc + lambda2 = mean_I - sqrt_disc + + # Semi-axes proportional to sqrt of eigenvalues, scaled by 2 for visual size + semi_major = 2.0 * np.sqrt(max(lambda1, 0.0)) + semi_minor = 2.0 * np.sqrt(max(lambda2, 0.0)) + + # Angle of the major axis + angle = 0.5 * np.arctan2(2.0 * Ixy, Iyy - Ixx) + + return float(semi_major), float(semi_minor), float(angle) + + +def _draw_circle_filled(canvas: np.ndarray, cy: float, cx: float, r: float) -> None: + h, w = canvas.shape + y_lo = max(0, int(cy - r - 1)) + y_hi = min(h, int(cy + r + 2)) + x_lo = max(0, int(cx - r - 1)) + x_hi = min(w, int(cx + r + 2)) + yy, xx = np.ogrid[y_lo:y_hi, x_lo:x_hi] + dist_sq = (yy - cy) ** 2 + (xx - cx) ** 2 + canvas[y_lo:y_hi, x_lo:x_hi] |= (dist_sq <= r * r) + + +def _draw_circle_outline(canvas: np.ndarray, cy: float, cx: float, r: float, thickness: float = 1.5) -> None: + h, w = canvas.shape + y_lo = max(0, int(cy - r - thickness - 1)) + y_hi = min(h, int(cy + r + thickness + 2)) + x_lo = max(0, int(cx - r - thickness - 1)) + x_hi = min(w, int(cx + r + thickness + 2)) + yy, xx = np.ogrid[y_lo:y_hi, x_lo:x_hi] + dist = np.sqrt((yy - cy) ** 2 + (xx - cx) ** 2) + canvas[y_lo:y_hi, x_lo:x_hi] |= (np.abs(dist - r) < thickness) + + +def _draw_rect_filled(canvas: np.ndarray, y0: int, y1: int, x0: int, x1: int) -> None: + h, w = canvas.shape + y0c, y1c = max(0, y0), min(h, y1) + x0c, x1c = max(0, x0), min(w, x1) + canvas[y0c:y1c, x0c:x1c] = True + + +def _draw_rect_outline(canvas: np.ndarray, y0: int, y1: int, x0: int, x1: int, thickness: int = 1) -> None: + h, w = canvas.shape + y0c, y1c = max(0, y0), min(h, y1) + x0c, x1c = max(0, x0), min(w, x1) + # Top edge + canvas[y0c:min(h, y0c + thickness), x0c:x1c] = True + # Bottom edge + canvas[max(0, y1c - thickness):y1c, x0c:x1c] = True + # Left edge + canvas[y0c:y1c, x0c:min(w, x0c + thickness)] = True + # Right edge + canvas[y0c:y1c, max(0, x1c - thickness):x1c] = True + + +def _draw_cross(canvas: np.ndarray, cy: float, cx: float, arm: int = 3) -> None: + h, w = canvas.shape + iy, ix = int(round(cy)), int(round(cx)) + for d in range(-arm, arm + 1): + if 0 <= iy + d < h and 0 <= ix < w: + canvas[iy + d, ix] = True + if 0 <= iy < h and 0 <= ix + d < w: + canvas[iy, ix + d] = True + + +def _draw_ellipse_filled(canvas: np.ndarray, cy: float, cx: float, + semi_major: float, semi_minor: float, angle: float) -> None: + h, w = canvas.shape + r_max = max(semi_major, semi_minor, 1.0) + y_lo = max(0, int(cy - r_max - 1)) + y_hi = min(h, int(cy + r_max + 2)) + x_lo = max(0, int(cx - r_max - 1)) + x_hi = min(w, int(cx + r_max + 2)) + yy, xx = np.ogrid[y_lo:y_hi, x_lo:x_hi] + cos_a, sin_a = np.cos(angle), np.sin(angle) + dy = yy - cy + dx = xx - cx + # Rotate into ellipse-aligned coordinates + u = cos_a * dx + sin_a * dy + v = -sin_a * dx + cos_a * dy + a = max(semi_major, 0.5) + b = max(semi_minor, 0.5) + canvas[y_lo:y_hi, x_lo:x_hi] |= ((u / a) ** 2 + (v / b) ** 2 <= 1.0) + + +def _draw_ellipse_outline(canvas: np.ndarray, cy: float, cx: float, + semi_major: float, semi_minor: float, angle: float, + thickness: float = 1.5) -> None: + h, w = canvas.shape + r_max = max(semi_major, semi_minor, 1.0) + y_lo = max(0, int(cy - r_max - thickness - 1)) + y_hi = min(h, int(cy + r_max + thickness + 2)) + x_lo = max(0, int(cx - r_max - thickness - 1)) + x_hi = min(w, int(cx + r_max + thickness + 2)) + yy, xx = np.ogrid[y_lo:y_hi, x_lo:x_hi] + cos_a, sin_a = np.cos(angle), np.sin(angle) + dy = yy - cy + dx = xx - cx + u = cos_a * dx + sin_a * dy + v = -sin_a * dx + cos_a * dy + a = max(semi_major, 0.5) + b = max(semi_minor, 0.5) + ellipse_val = (u / a) ** 2 + (v / b) ** 2 + canvas[y_lo:y_hi, x_lo:x_hi] |= (np.abs(np.sqrt(ellipse_val) - 1.0) < thickness / max(a, b)) + + +@register_node(display_name="Grain Visualization") +class GrainVisualization: + @classmethod + def INPUT_TYPES(cls): + return { + "required": { + "field": ("DATA_FIELD",), + "mask": ("IMAGE",), + "style": (["inscribed_disc", "bounding_box", "centroid", "ellipse"], {"default": "inscribed_disc"}), + "fill": ("BOOLEAN", {"default": False}), + } + } + + OUTPUTS = ( + ('IMAGE', 'result'), + ('DATA_FIELD', 'labeled'), + ) + FUNCTION = "process" + + DESCRIPTION = ( + "Visualize labeled grains as geometric shapes — inscribed discs, bounding boxes, " + "centroid markers, or fitted ellipses. Produces a mask image with the chosen shapes " + "and a labeled field where each grain has a unique integer value. " + "Equivalent to Gwyddion's grain selection visualization (grain_makesel)." + ) + + def process(self, field: DataField, mask: np.ndarray, style: str, fill: bool) -> tuple: + mask_bool = mask_to_bool(mask) + labels, n_grains = label(mask_bool.astype(np.int32)) + slices = find_objects(labels) + + h, w = mask_bool.shape[:2] + canvas = np.zeros((h, w), dtype=bool) + + for gid in range(1, n_grains + 1): + slc = slices[gid - 1] + if slc is None: + continue + + grain_mask = labels == gid + cy, cx = _grain_centroid(grain_mask, slc) + + if style == "inscribed_disc": + r = _grain_inscribed_radius(grain_mask, slc) + if r < 0.5: + r = 0.5 + if fill: + _draw_circle_filled(canvas, cy, cx, r) + else: + _draw_circle_outline(canvas, cy, cx, r) + + elif style == "bounding_box": + y0, y1 = slc[0].start, slc[0].stop + x0, x1 = slc[1].start, slc[1].stop + if fill: + _draw_rect_filled(canvas, y0, y1, x0, x1) + else: + _draw_rect_outline(canvas, y0, y1, x0, x1) + + elif style == "centroid": + arm = max(3, int(round(min(h, w) * 0.01))) + _draw_cross(canvas, cy, cx, arm) + + elif style == "ellipse": + semi_major, semi_minor, angle = _grain_inertia(grain_mask, slc) + if semi_major < 0.5: + semi_major = 0.5 + if semi_minor < 0.5: + semi_minor = 0.5 + if fill: + _draw_ellipse_filled(canvas, cy, cx, semi_major, semi_minor, angle) + else: + _draw_ellipse_outline(canvas, cy, cx, semi_major, semi_minor, angle) + + else: + raise ValueError(f"Unknown visualization style: {style!r}") + + result = bool_to_mask(canvas) + labeled_field = field.replace(data=labels.astype(np.float64)) + + return (result, labeled_field) diff --git a/backend/nodes/logistic_classification.py b/backend/nodes/logistic_classification.py new file mode 100644 index 0000000..9bfef9b --- /dev/null +++ b/backend/nodes/logistic_classification.py @@ -0,0 +1,229 @@ +"""Logistic classification — classify features using logistic regression.""" + +from __future__ import annotations + +import numpy as np +from scipy.ndimage import gaussian_filter, sobel + +from backend.node_registry import register_node +from backend.data_types import DataField +from backend.nodes.helpers import mask_to_bool, bool_to_mask + + +def _build_features(data: np.ndarray, use_gaussians: bool, n_gaussians: int, + use_sobel: bool, use_laplacian: bool) -> np.ndarray: + """Build a feature matrix from the height field. + + Each feature is normalized to zero mean, unit variance. The raw + (normalized) height is always included as the first feature. + """ + h, w = data.shape + features: list[np.ndarray] = [] + + # Always include raw height (normalized) + features.append(data.ravel().copy()) + + # Gaussian blur features at increasing scales + if use_gaussians: + for i in range(int(n_gaussians)): + sigma = float(2 ** i) + features.append(gaussian_filter(data, sigma).ravel()) + + # Sobel gradient features + if use_sobel: + features.append(sobel(data, axis=0).ravel()) + features.append(sobel(data, axis=1).ravel()) + + # Laplacian feature (sum of second differences) + if use_laplacian: + lap = np.zeros_like(data) + lap[1:-1, :] += data[:-2, :] - 2 * data[1:-1, :] + data[2:, :] + lap[:, 1:-1] += data[:, :-2] - 2 * data[:, 1:-1] + data[:, 2:] + features.append(lap.ravel()) + + # Stack into (n_pixels, n_features) matrix + X = np.column_stack(features) + + # Normalize each feature to zero mean, unit variance + means = X.mean(axis=0) + stds = X.std(axis=0) + stds[stds == 0] = 1.0 + X = (X - means) / stds + + # Add bias column + X = np.column_stack([np.ones(X.shape[0]), X]) + + return X + + +def _sigmoid(z: np.ndarray) -> np.ndarray: + z = np.clip(z, -500, 500) + return 1.0 / (1.0 + np.exp(-z)) + + +def _otsu_threshold(data: np.ndarray) -> float: + """Simple Otsu threshold on flattened data.""" + flat = data.ravel() + counts, bin_edges = np.histogram(flat, bins=256) + centers = 0.5 * (bin_edges[:-1] + bin_edges[1:]) + total = counts.sum() + if total == 0: + return float(np.median(flat)) + + sum_total = (counts * centers).sum() + sum_bg = 0.0 + weight_bg = 0.0 + best_var = -1.0 + best_thresh = float(centers[0]) + + for i in range(len(counts)): + weight_bg += counts[i] + if weight_bg == 0: + continue + weight_fg = total - weight_bg + if weight_fg == 0: + break + sum_bg += counts[i] * centers[i] + mean_bg = sum_bg / weight_bg + mean_fg = (sum_total - sum_bg) / weight_fg + var_between = weight_bg * weight_fg * (mean_bg - mean_fg) ** 2 + if var_between > best_var: + best_var = var_between + best_thresh = float(centers[i]) + + return best_thresh + + +def _train_logistic(X: np.ndarray, y: np.ndarray, regularization: float, + max_iter: int, seed: int) -> np.ndarray: + """Train logistic regression via gradient descent. + + Parameters + ---------- + X : (m, n_features+1) array with bias column already included. + y : (m,) binary labels (0 or 1). + regularization : L2 penalty lambda. + max_iter : maximum gradient descent iterations. + seed : random seed (unused here; theta starts at zeros). + + Returns + ------- + theta : (n_features+1,) weight vector. + """ + rng = np.random.default_rng(seed) + n = X.shape[1] + theta = np.zeros(n) + m = len(y) + lr = 0.1 + + for _ in range(max_iter): + h = _sigmoid(X @ theta) + error = h - y + + grad = X.T @ error / m + # L2 regularization (don't regularize bias at index 0) + reg_term = (regularization / m) * theta + reg_term[0] = 0.0 + grad += reg_term + + theta -= lr * grad + + if np.linalg.norm(grad) < 1e-6: + break + + return theta + + +@register_node(display_name="Logistic Classification") +class LogisticClassification: + _CUSTOM_PREVIEW = True + + @classmethod + def INPUT_TYPES(cls): + return { + "required": { + "field": ("DATA_FIELD",), + "use_gaussians": ("BOOLEAN", {"default": True}), + "n_gaussians": ("INT", { + "default": 4, "min": 1, "max": 10, + "show_when_widget_value": {"use_gaussians": [True]}, + }), + "use_sobel": ("BOOLEAN", {"default": True}), + "use_laplacian": ("BOOLEAN", {"default": True}), + "regularization": ("FLOAT", {"default": 1.0, "min": 0.0, "max": 10.0, "step": 0.1}), + "max_iter": ("INT", {"default": 500, "min": 10, "max": 5000}), + "seed": ("INT", {"default": 42, "min": 0, "max": 999999}), + }, + "optional": { + "training_mask": ("IMAGE",), + }, + } + + OUTPUTS = ( + ('IMAGE', 'mask'), + ('DATA_FIELD', 'probability'), + ) + FUNCTION = "process" + + DESCRIPTION = ( + "Classify surface features using logistic regression on engineered " + "height-derived features (Gaussian blurs, Sobel gradients, Laplacian). " + "Optionally accepts a training mask; otherwise an Otsu-based threshold " + "generates pseudo-labels automatically." + ) + + def process( + self, + field: DataField, + use_gaussians: bool, + n_gaussians: int, + use_sobel: bool, + use_laplacian: bool, + regularization: float, + max_iter: int, + seed: int, + training_mask: np.ndarray | None = None, + ) -> tuple: + data = np.asarray(field.data, dtype=np.float64) + h, w = data.shape + + # Build feature matrix for all pixels + X_all = _build_features(data, use_gaussians, n_gaussians, use_sobel, use_laplacian) + + if training_mask is not None: + # Extract training labels from the mask + mask_bool = mask_to_bool(training_mask) + if mask_bool.shape[:2] != (h, w): + raise ValueError( + f"Training mask shape {mask_bool.shape} does not match " + f"field shape {(h, w)}." + ) + labeled_pixels = mask_bool.ravel() + # Use masked pixels as positive class, unmasked as negative + y_train = labeled_pixels.astype(np.float64) + X_train = X_all + else: + # No training mask: use Otsu threshold to create pseudo-labels + threshold = _otsu_threshold(data) + y_train = (data.ravel() >= threshold).astype(np.float64) + X_train = X_all + + # Train logistic regression + theta = _train_logistic(X_train, y_train, regularization, max_iter, seed) + + # Apply to all pixels + probability = _sigmoid(X_all @ theta).reshape(h, w) + + # Create binary mask + mask = bool_to_mask(probability > 0.5) + + # Emit preview + from backend.execution_context import emit_preview + from backend.data_types import encode_preview + from backend.nodes.helpers import _mask_overlay + emit_preview(encode_preview(_mask_overlay(field, mask))) + + # Build probability output as a DataField + prob_field = field.replace(data=probability, si_unit_z="") + + return (mask, prob_field) diff --git a/backend/nodes/mark_disconnected.py b/backend/nodes/mark_disconnected.py new file mode 100644 index 0000000..fd4e3d7 --- /dev/null +++ b/backend/nodes/mark_disconnected.py @@ -0,0 +1,71 @@ +"""Mark disconnected regions — mask topologically isolated surface regions.""" + +from __future__ import annotations +import numpy as np +from scipy.ndimage import grey_opening, grey_closing +from backend.node_registry import register_node +from backend.data_types import DataField +from backend.nodes.helpers import bool_to_mask, _mask_structure, emit_mask_preview + + +@register_node(display_name="Mark Disconnected") +class MarkDisconnected: + """ + Detect topologically disconnected (isolated) surface regions using + morphological opening/closing to build a defect-free reference, then + thresholding the residual difference. + """ + _CUSTOM_PREVIEW = True + + @classmethod + def INPUT_TYPES(cls): + return { + "required": { + "field": ("DATA_FIELD",), + "defect_type": (["positive", "negative", "both"],), + "radius": ("INT", {"default": 5, "min": 1, "max": 100, "step": 1}), + "threshold": ("FLOAT", {"default": 0.1, "min": 0.001, "max": 1.0, "step": 0.001}), + } + } + + OUTPUTS = ( + ('IMAGE', 'mask'), + ) + FUNCTION = "process" + + DESCRIPTION = ( + "Mark topologically disconnected (isolated) surface regions. " + "A morphological opening followed by closing builds a smooth " + "defect-free reference surface; pixels whose deviation from that " + "reference exceeds the sensitivity threshold are flagged. " + "Equivalent to Gwyddion's mark_disconn module." + ) + + def process(self, field: DataField, defect_type: str, radius: int, threshold: float) -> tuple: + data = field.data.astype(np.float64) + + # Build a disk structuring element for grey-scale morphology. + struct = _mask_structure(radius, "disk") + + # Morphological opening then closing produces a defect-free reference. + reference = grey_opening(data, footprint=struct) + reference = grey_closing(reference, footprint=struct) + + difference = data - reference + diff_range = difference.max() - difference.min() + + # Avoid division-by-zero on perfectly flat surfaces. + if diff_range == 0: + mask = np.zeros(data.shape, dtype=bool) + else: + abs_threshold = threshold * diff_range + if defect_type == "positive": + mask = difference > abs_threshold + elif defect_type == "negative": + mask = difference < -abs_threshold + else: # "both" + mask = np.abs(difference) > abs_threshold + + out = bool_to_mask(mask) + emit_mask_preview(field, out) + return (out,) diff --git a/backend/nodes/mask_noisify.py b/backend/nodes/mask_noisify.py new file mode 100644 index 0000000..1e7287a --- /dev/null +++ b/backend/nodes/mask_noisify.py @@ -0,0 +1,84 @@ +"""Mask noisify -- add random perturbation to mask boundaries.""" + +from __future__ import annotations +import numpy as np +from backend.node_registry import register_node +from backend.data_types import DataField +from backend.nodes.helpers import mask_to_bool, bool_to_mask, emit_mask_preview + + +@register_node(display_name="Mask Noisify") +class MaskNoisify: + """ + Add random perturbation to mask boundaries. + """ + _CUSTOM_PREVIEW = True + + @classmethod + def INPUT_TYPES(cls): + return { + "required": { + "mask": ("IMAGE",), + "density": ("FLOAT", {"default": 0.1, "min": 0.0, "max": 1.0, "step": 0.01}), + "direction": (["both", "add", "remove"],), + "boundaries_only": ("BOOLEAN", {"default": True}), + "seed": ("INT", {"default": 42, "min": 0, "max": 999999}), + }, + "optional": { + "field": ("DATA_FIELD",), + } + } + + OUTPUTS = ( + ('IMAGE', 'mask'), + ) + FUNCTION = "process" + + DESCRIPTION = ( + "Add random noise to a binary mask by flipping pixels near boundaries. " + "Control the fraction of affected pixels with density, restrict changes " + "to boundary pixels, and choose whether to add, remove, or both. " + "Use a fixed seed for reproducible results." + ) + + def process(self, mask: np.ndarray, density: float, direction: str, + boundaries_only: bool, seed: int, + field: DataField | None = None) -> tuple: + binary = mask_to_bool(mask) + + # Identify boundary pixels: pixels that differ from at least one neighbour + if boundaries_only: + boundary = np.zeros_like(binary) + for shift_axis, shift_dir in [(0, 1), (0, -1), (1, 1), (1, -1)]: + shifted = np.roll(binary, shift_dir, axis=shift_axis) + boundary |= (binary != shifted) + else: + boundary = np.ones_like(binary) + + # Select candidate pixels based on direction + if direction == "add": + candidates = boundary & ~binary + elif direction == "remove": + candidates = boundary & binary + else: # "both" + candidates = boundary + + # Randomly flip density fraction of candidates + candidate_indices = np.argwhere(candidates) + n_candidates = len(candidate_indices) + + if n_candidates > 0 and density > 0: + rng = np.random.default_rng(seed) + n_flip = int(round(density * n_candidates)) + n_flip = max(0, min(n_flip, n_candidates)) + if n_flip > 0: + chosen = rng.choice(n_candidates, size=n_flip, replace=False) + for idx in chosen: + r, c = candidate_indices[idx] + binary[r, c] = ~binary[r, c] + + out = bool_to_mask(binary) + + emit_mask_preview(field, out) + + return (out,) diff --git a/backend/nodes/mask_shift.py b/backend/nodes/mask_shift.py new file mode 100644 index 0000000..a190eeb --- /dev/null +++ b/backend/nodes/mask_shift.py @@ -0,0 +1,98 @@ +"""Mask shift — translate mask by pixel offset.""" + +from __future__ import annotations + +import numpy as np + +from backend.node_registry import register_node +from backend.data_types import DataField +from backend.nodes.helpers import mask_to_bool, bool_to_mask, emit_mask_preview + + +@register_node(display_name="Mask Shift") +class MaskShift: + """Translate a binary mask by an integer pixel offset.""" + + _CUSTOM_PREVIEW = True + + @classmethod + def INPUT_TYPES(cls): + return { + "required": { + "mask": ("IMAGE",), + "shift_x": ("INT", {"default": 0, "min": -1000, "max": 1000, "step": 1}), + "shift_y": ("INT", {"default": 0, "min": -1000, "max": 1000, "step": 1}), + "border_mode": (["zero", "wrap", "mirror"],), + }, + "optional": { + "field": ("DATA_FIELD",), + }, + } + + OUTPUTS = ( + ('IMAGE', 'mask'), + ) + FUNCTION = "process" + + DESCRIPTION = ( + "Translate a binary mask by an integer pixel offset. " + "Choose how out-of-bounds regions are filled: zero (empty), " + "wrap (periodic roll), or mirror (reflected padding)." + ) + + def process(self, mask: np.ndarray, shift_x: int, shift_y: int, + border_mode: str, field: DataField | None = None) -> tuple: + binary = mask_to_bool(mask) + + if border_mode == "wrap": + result = self._shift_wrap(binary, shift_x, shift_y) + elif border_mode == "zero": + result = self._shift_zero(binary, shift_x, shift_y) + elif border_mode == "mirror": + result = self._shift_mirror(binary, shift_x, shift_y) + else: + raise ValueError(f"Unknown border mode: {border_mode}") + + out = bool_to_mask(result) + emit_mask_preview(field, out) + return (out,) + + @staticmethod + def _shift_wrap(binary: np.ndarray, sx: int, sy: int) -> np.ndarray: + """Shift with periodic wrapping (np.roll).""" + return np.roll(np.roll(binary, sx, axis=1), sy, axis=0) + + @staticmethod + def _shift_zero(binary: np.ndarray, sx: int, sy: int) -> np.ndarray: + """Shift then zero-fill the wrapped region.""" + result = np.roll(np.roll(binary, sx, axis=1), sy, axis=0) + h, w = result.shape[:2] + + # Zero-fill columns wrapped by horizontal shift + if sx > 0: + result[:, :sx] = False + elif sx < 0: + result[:, w + sx:] = False + + # Zero-fill rows wrapped by vertical shift + if sy > 0: + result[:sy, :] = False + elif sy < 0: + result[h + sy:, :] = False + + return result + + @staticmethod + def _shift_mirror(binary: np.ndarray, sx: int, sy: int) -> np.ndarray: + """Shift using reflected padding then crop back to original size.""" + h, w = binary.shape[:2] + abs_sx = abs(sx) + abs_sy = abs(sy) + + # Pad with reflect mode + padded = np.pad(binary, ((abs_sy, abs_sy), (abs_sx, abs_sx)), mode="reflect") + + # Crop with offset to achieve the shift + row_start = abs_sy - sy + col_start = abs_sx - sx + return padded[row_start:row_start + h, col_start:col_start + w] diff --git a/backend/nodes/neural_classification.py b/backend/nodes/neural_classification.py new file mode 100644 index 0000000..24acc1e --- /dev/null +++ b/backend/nodes/neural_classification.py @@ -0,0 +1,204 @@ +"""Neural network classification — classify pixels using a simple feedforward network.""" + +from __future__ import annotations + +import numpy as np +from scipy.ndimage import gaussian_filter + +from backend.node_registry import register_node +from backend.data_types import DataField +from backend.nodes.helpers import mask_to_bool, bool_to_mask + + +def _sigmoid(x: np.ndarray) -> np.ndarray: + """Numerically stable sigmoid.""" + return np.where( + x >= 0, + 1.0 / (1.0 + np.exp(-x)), + np.exp(x) / (1.0 + np.exp(x)), + ) + + +def _extract_features(data: np.ndarray, n_gaussians: int) -> np.ndarray: + """Build multi-scale Gaussian feature matrix from 2-D data. + + For each scale sigma = 2^i (i = 0 .. n_gaussians-1), compute + gaussian_filter(data, sigma) and stack as feature columns. + Each feature is normalised to zero mean and unit variance. + """ + rows, cols = data.shape + features = np.empty((rows * cols, n_gaussians), dtype=np.float64) + + for i in range(n_gaussians): + sigma = 2.0 ** i + blurred = gaussian_filter(data, sigma).ravel() + mean = blurred.mean() + std = blurred.std() + if std > 0: + blurred = (blurred - mean) / std + else: + blurred = blurred - mean + features[:, i] = blurred + + return features + + +def _forward(X: np.ndarray, W1: np.ndarray, b1: np.ndarray, + W2: np.ndarray, b2: np.ndarray) -> tuple[np.ndarray, np.ndarray]: + """Forward pass through a 2-layer sigmoid network.""" + h = _sigmoid(X @ W1 + b1) + y = _sigmoid(h @ W2 + b2) + return h, y + + +def _train_network( + X: np.ndarray, + targets: np.ndarray, + W1: np.ndarray, + b1: np.ndarray, + W2: np.ndarray, + b2: np.ndarray, + train_steps: int, + lr: float = 0.1, +) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: + """Train via gradient descent with binary cross-entropy loss.""" + eps = 1e-7 + n = X.shape[0] + t = targets.reshape(-1, 1) + + for _ in range(train_steps): + # Forward + h, y = _forward(X, W1, b1, W2, b2) + + # Clamp to avoid log(0) + y_clamped = np.clip(y, eps, 1.0 - eps) + + # Backward — output layer + dy = (y_clamped - t) / (y_clamped * (1.0 - y_clamped) + eps) + dy *= y * (1.0 - y) # sigmoid derivative + + dW2 = (h.T @ dy) / n + db2 = dy.mean(axis=0) + + # Backward — hidden layer + dh = (dy @ W2.T) * h * (1.0 - h) + dW1 = (X.T @ dh) / n + db1 = dh.mean(axis=0) + + # Update + W1 -= lr * dW1 + b1 -= lr * db1 + W2 -= lr * dW2 + b2 -= lr * db2 + + return W1, b1, W2, b2 + + +@register_node(display_name="Neural Classification") +class NeuralClassification: + @classmethod + def INPUT_TYPES(cls): + return { + "required": { + "field": ("DATA_FIELD",), + "n_gaussians": ("INT", {"default": 4, "min": 1, "max": 10, "step": 1}), + "n_hidden": ("INT", {"default": 16, "min": 4, "max": 128, "step": 1}), + "train_steps": ("INT", {"default": 200, "min": 10, "max": 5000, "step": 1}), + "seed": ("INT", {"default": 42, "min": 0, "max": 999999, "step": 1}), + }, + "optional": { + "training_mask": ("IMAGE",), + }, + } + + OUTPUTS = ( + ('IMAGE', 'mask'), + ('DATA_FIELD', 'probability'), + ) + FUNCTION = "process" + + DESCRIPTION = ( + "Classify surface pixels into two classes using a simple two-layer " + "feedforward neural network with sigmoid activations. Features are " + "extracted via multi-scale Gaussian filtering. When a training mask " + "is provided the network learns from labelled pixels; otherwise it " + "uses unsupervised self-labelling from the initial random projection. " + "Equivalent in purpose to Gwyddion's neural.c classifier." + ) + + def process( + self, + field: DataField, + n_gaussians: int, + n_hidden: int, + train_steps: int, + seed: int, + training_mask: np.ndarray | None = None, + ) -> tuple: + data = np.asarray(field.data, dtype=np.float64) + yres, xres = data.shape + n_features = int(n_gaussians) + n_hidden = int(n_hidden) + train_steps = int(train_steps) + + # 1. Feature extraction + X_all = _extract_features(data, n_features) + + # 2. Initialise weights + rng = np.random.default_rng(int(seed)) + scale1 = np.sqrt(2.0 / n_features) + W1 = rng.standard_normal((n_features, n_hidden)) * scale1 + b1 = np.zeros(n_hidden) + scale2 = np.sqrt(2.0 / n_hidden) + W2 = rng.standard_normal((n_hidden, 1)) * scale2 + b2 = np.zeros(1) + + # 3/4. Training + if training_mask is not None: + # Supervised — use labelled pixels + mask_bool = mask_to_bool(training_mask) + if mask_bool.shape != data.shape: + raise ValueError( + f"Training mask shape {mask_bool.shape} does not match " + f"field shape {data.shape}." + ) + # Class B = masked (255), class A = unmasked but we need both labels. + # Pixels that are 0 are class A, pixels that are 255 are class B. + # We train on ALL pixels that have a definitive label. + labels_flat = training_mask.ravel().astype(np.float64) / 255.0 + # Use all pixels as training data (0 = class A, 1 = class B) + X_train = X_all + targets = labels_flat + W1, b1, W2, b2 = _train_network( + X_train, targets, W1, b1, W2, b2, train_steps, + ) + else: + # Unsupervised — use random projection to create initial labels, + # then refine with self-training. + _, y_init = _forward(X_all, W1, b1, W2, b2) + self_labels = (y_init.ravel() > 0.5).astype(np.float64) + # Train on the self-assigned labels for a few iterations + steps = min(train_steps, 50) + W1, b1, W2, b2 = _train_network( + X_all, self_labels, W1, b1, W2, b2, steps, + ) + + # 5. Apply trained network to all pixels + _, prob_flat = _forward(X_all, W1, b1, W2, b2) + probability = prob_flat.reshape(yres, xres) + + # 6. Build outputs + mask = bool_to_mask(probability > 0.5) + + prob_field = DataField( + data=probability, + xreal=field.xreal, + yreal=field.yreal, + xoff=field.xoff, + yoff=field.yoff, + si_unit_xy=field.si_unit_xy, + si_unit_z="", + domain="spatial", + ) + + return (mask, prob_field) diff --git a/backend/nodes/pixel_classification.py b/backend/nodes/pixel_classification.py new file mode 100644 index 0000000..f69d36d --- /dev/null +++ b/backend/nodes/pixel_classification.py @@ -0,0 +1,209 @@ +"""Pixel classification — classify pixels using decision tree on height, slope, and curvature.""" + +from __future__ import annotations +import numpy as np +from backend.node_registry import register_node +from backend.data_types import DataField +from backend.nodes.helpers import bool_to_mask + + +def _compute_slope(data: np.ndarray) -> np.ndarray: + """Gradient magnitude via np.gradient.""" + gy, gx = np.gradient(data.astype(np.float64)) + return np.sqrt(gx**2 + gy**2) + + +def _compute_curvature(data: np.ndarray) -> np.ndarray: + """Laplacian (sum of second derivatives).""" + d = data.astype(np.float64) + gy, gx = np.gradient(d) + gyy, _ = np.gradient(gy) + _, gxx = np.gradient(gx) + return np.abs(gxx + gyy) + + +def _feature_maps(data: np.ndarray, feature: str) -> list[np.ndarray]: + """Return a list of 2-D feature arrays based on the feature selector.""" + height = data.astype(np.float64) + if feature == "height": + return [height] + slope = _compute_slope(data) + if feature == "slope": + return [slope] + curvature = _compute_curvature(data) + if feature == "curvature": + return [curvature] + if feature == "height_slope": + return [height, slope] + # "all" + return [height, slope, curvature] + + +def _normalize_01(arr: np.ndarray) -> np.ndarray: + vmin, vmax = arr.min(), arr.max() + if vmax > vmin: + return (arr - vmin) / (vmax - vmin) + return np.zeros_like(arr) + + +def _classify_single(values: np.ndarray, n_classes: int, method: str) -> np.ndarray: + """Classify a single feature map into n_classes using the chosen method.""" + labels = np.zeros(values.shape, dtype=np.int32) + + if method == "equal_range": + vmin, vmax = values.min(), values.max() + if vmax <= vmin: + return labels + edges = np.linspace(vmin, vmax, n_classes + 1) + for i in range(n_classes - 1): + labels[values >= edges[i + 1]] = i + 1 + + elif method == "quantile": + percentiles = np.linspace(0, 100, n_classes + 1) + edges = np.percentile(values, percentiles) + for i in range(n_classes - 1): + labels[values >= edges[i + 1]] = i + 1 + + elif method == "otsu": + # Multi-Otsu: find n_classes-1 thresholds via histogram analysis + flat = values.ravel() + n_bins = min(256, max(32, len(flat) // 10)) + counts, bin_edges = np.histogram(flat, bins=n_bins) + centers = 0.5 * (bin_edges[:-1] + bin_edges[1:]) + total = counts.sum() + + if total == 0 or n_classes < 2: + return labels + + # For multi-Otsu, find thresholds that minimise intra-class variance + # Use quantile-based initial thresholds then refine with exhaustive + # search over histogram bins for each threshold + thresholds = [] + if n_classes == 2: + # Standard single-threshold Otsu + best_var = -1.0 + best_t = 0 + cum_sum = 0.0 + cum_count = 0 + total_sum = float(np.sum(counts * centers)) + for i in range(n_bins - 1): + cum_count += counts[i] + cum_sum += counts[i] * centers[i] + if cum_count == 0 or cum_count == total: + continue + w0 = cum_count / total + w1 = 1.0 - w0 + mu0 = cum_sum / cum_count + mu1 = (total_sum - cum_sum) / (total - cum_count) + between_var = w0 * w1 * (mu0 - mu1) ** 2 + if between_var > best_var: + best_var = between_var + best_t = i + thresholds = [0.5 * (bin_edges[best_t + 1] + bin_edges[best_t + 2])] + else: + # Multi-threshold: use quantile splits as a good approximation + percentiles = np.linspace(0, 100, n_classes + 1)[1:-1] + thresholds = list(np.percentile(flat, percentiles)) + + thresholds = sorted(thresholds) + for i, t in enumerate(thresholds): + labels[values >= t] = i + 1 + + else: + raise ValueError(f"Unknown classification method: {method!r}") + + return labels + + +def _kmeans_classify(features: np.ndarray, n_classes: int, max_iter: int = 20) -> np.ndarray: + """Simple k-means on stacked normalised features. + + Parameters + ---------- + features : (n_pixels, n_features) array + n_classes : number of clusters + max_iter : maximum iterations + + Returns + ------- + labels : (n_pixels,) int32 array with values in [0, n_classes-1] + """ + rng = np.random.RandomState(42) + n_pixels = features.shape[0] + # Initialise centroids by choosing random data points + indices = rng.choice(n_pixels, size=min(n_classes, n_pixels), replace=False) + centroids = features[indices].copy() + + labels = np.zeros(n_pixels, dtype=np.int32) + for _ in range(max_iter): + # Assign each pixel to nearest centroid + dists = np.stack([ + np.sum((features - c) ** 2, axis=1) for c in centroids + ], axis=1) # (n_pixels, n_classes) + new_labels = np.argmin(dists, axis=1).astype(np.int32) + + if np.array_equal(new_labels, labels): + break + labels = new_labels + + # Update centroids + for k in range(n_classes): + members = features[labels == k] + if len(members) > 0: + centroids[k] = members.mean(axis=0) + + return labels + + +@register_node(display_name="Pixel Classification") +class PixelClassification: + @classmethod + def INPUT_TYPES(cls): + return { + "required": { + "field": ("DATA_FIELD",), + "n_classes": ("INT", {"default": 3, "min": 2, "max": 10, "step": 1}), + "feature": (["height", "slope", "curvature", "height_slope", "all"],), + "method": (["otsu", "equal_range", "quantile"],), + } + } + + OUTPUTS = ( + ('DATA_FIELD', 'classified'), + ('IMAGE', 'mask'), + ) + FUNCTION = "process" + + DESCRIPTION = ( + "Classify pixels into discrete classes based on height, slope, and/or curvature. " + "Single-feature modes use threshold-based classification (Otsu, equal range, or quantile). " + "Multi-feature modes (height_slope, all) use k-means clustering. " + "Equivalent to Gwyddion's classify.c module." + ) + + def process(self, field: DataField, n_classes: int, feature: str, method: str) -> tuple: + data = np.asarray(field.data, dtype=np.float64) + maps = _feature_maps(data, feature) + + if len(maps) == 1: + # Single-feature: use threshold-based classification + labels = _classify_single(maps[0], int(n_classes), method) + else: + # Multi-feature: normalise and use k-means + normed = [_normalize_01(m) for m in maps] + stacked = np.stack([m.ravel() for m in normed], axis=1) # (n_pixels, n_features) + labels = _kmeans_classify(stacked, int(n_classes)).reshape(data.shape) + + # Build output DataField with integer class labels + classified = DataField( + data=labels.astype(np.float64), + xreal=field.xreal, + yreal=field.yreal, + si_unit_xy=field.si_unit_xy, + si_unit_z="", + ) + + # Mask for class 0 + mask = bool_to_mask(labels == 0) + + return (classified, mask) diff --git a/backend/nodes/presentation_ops.py b/backend/nodes/presentation_ops.py new file mode 100644 index 0000000..07af9d3 --- /dev/null +++ b/backend/nodes/presentation_ops.py @@ -0,0 +1,88 @@ +"""Presentation operations -- manage presentation overlays on data fields.""" + +from __future__ import annotations + +import numpy as np + +from backend.node_registry import register_node +from backend.data_types import DataField + + +@register_node(display_name="Presentation Ops") +class PresentationOps: + @classmethod + def INPUT_TYPES(cls): + return { + "required": { + "field": ("DATA_FIELD",), + "operation": (["logscale", "extract_presentation", "attach", "blend"],), + "blend_factor": ("FLOAT", { + "default": 0.5, + "min": 0.0, + "max": 1.0, + "step": 0.01, + "show_when_widget_value": {"operation": ["blend"]}, + }), + }, + "optional": { + "overlay": ("DATA_FIELD", { + "show_when_widget_value": {"operation": ["attach", "blend"]}, + }), + }, + } + + OUTPUTS = ( + ('DATA_FIELD', 'result'), + ) + FUNCTION = "process" + + DESCRIPTION = ( + "Manage presentation overlays on data fields. " + "logscale applies logarithmic scaling for visualising data with large dynamic range. " + "extract_presentation normalises the field to [0, 1]. " + "attach replaces the field data with an overlay (resampled if needed). " + "blend linearly mixes the field and overlay by a configurable factor. " + "Equivalent to Gwyddion's presentationops.c module." + ) + + def process(self, field: DataField, operation: str, blend_factor: float, + overlay: DataField | None = None) -> tuple: + data = np.asarray(field.data, dtype=np.float64) + + if operation == "logscale": + data_pos = data - data.min() + 1e-30 + result = np.log10(data_pos) + + elif operation == "extract_presentation": + dmin, dmax = data.min(), data.max() + if dmax > dmin: + result = (data - dmin) / (dmax - dmin) + else: + result = np.zeros_like(data) + + elif operation == "attach": + if overlay is None: + raise ValueError("'attach' operation requires an overlay field.") + overlay_data = np.asarray(overlay.data, dtype=np.float64) + result = self._match_shape(overlay_data, data.shape) + + elif operation == "blend": + if overlay is None: + raise ValueError("'blend' operation requires an overlay field.") + overlay_data = np.asarray(overlay.data, dtype=np.float64) + overlay_matched = self._match_shape(overlay_data, data.shape) + result = (1.0 - blend_factor) * data + blend_factor * overlay_matched + + else: + raise ValueError(f"Unknown operation: {operation!r}") + + return (field.replace(data=result),) + + @staticmethod + def _match_shape(source: np.ndarray, target_shape: tuple[int, ...]) -> np.ndarray: + """Resample *source* to *target_shape* using scipy zoom if shapes differ.""" + if source.shape == target_shape: + return source + from scipy.ndimage import zoom + factors = tuple(t / s for t, s in zip(target_shape, source.shape)) + return zoom(source, factors, order=3) diff --git a/backend/nodes/psf_estimation.py b/backend/nodes/psf_estimation.py new file mode 100644 index 0000000..48f9e6e --- /dev/null +++ b/backend/nodes/psf_estimation.py @@ -0,0 +1,177 @@ +"""PSF estimation — estimate and fit point spread functions for deconvolution.""" + +from __future__ import annotations + +import numpy as np + +from backend.node_registry import register_node +from backend.data_types import DataField, RecordTable + + +@register_node(display_name="PSF Estimation") +class PSFEstimation: + @classmethod + def INPUT_TYPES(cls): + return { + "required": { + "measured": ("DATA_FIELD",), + "ideal": ("DATA_FIELD",), + "method": (["wiener", "least_squares", "gaussian_fit"], {"default": "wiener"}), + "regularization": ("FLOAT", {"default": 0.01, "min": 1e-6, "max": 1.0, "step": 0.001}), + "psf_size": ("INT", {"default": 32, "min": 4, "max": 128}), + } + } + + OUTPUTS = ( + ('DATA_FIELD', 'psf'), + ('RECORD_TABLE', 'parameters'), + ) + FUNCTION = "process" + + DESCRIPTION = ( + "Estimate a point spread function (PSF) from a measured (blurred) image " + "and an ideal (sharp) reference. The PSF can then be used with the " + "Deconvolution node to restore other images. Three methods are available: " + "pseudo-Wiener deconvolution, regularised least-squares, and Gaussian fit. " + "Equivalent to Gwyddion's psf.c / psf-fit.c modules." + ) + + # ------------------------------------------------------------------ + # helpers + # ------------------------------------------------------------------ + + @staticmethod + def _crop_centre(arr: np.ndarray, size: int) -> np.ndarray: + """Crop the central *size x size* region from *arr*.""" + yc, xc = arr.shape[0] // 2, arr.shape[1] // 2 + half = size // 2 + return arr[yc - half : yc - half + size, xc - half : xc - half + size] + + @staticmethod + def _normalise(psf: np.ndarray) -> np.ndarray: + """Normalise so that the PSF sums to 1.""" + s = psf.sum() + if abs(s) > 1e-30: + psf = psf / s + return psf + + @staticmethod + def _fit_gaussian_2d(psf: np.ndarray): + """Fit a 2-D Gaussian to *psf* using moment analysis. + + Returns (gaussian_array, sigma_x, sigma_y, amplitude). + """ + h, w = psf.shape + psf_pos = np.maximum(psf, 0.0) + total = psf_pos.sum() + if total < 1e-30: + return np.zeros_like(psf), 0.0, 0.0, 0.0 + + y_idx, x_idx = np.mgrid[0:h, 0:w] + + # centroid + cx = float(np.sum(x_idx * psf_pos) / total) + cy = float(np.sum(y_idx * psf_pos) / total) + + # second moments → sigma + sx = float(np.sqrt(np.sum(psf_pos * (x_idx - cx) ** 2) / total)) + sy = float(np.sqrt(np.sum(psf_pos * (y_idx - cy) ** 2) / total)) + sx = max(sx, 1e-6) + sy = max(sy, 1e-6) + + amplitude = float(psf_pos.max()) + + gauss = amplitude * np.exp( + -((x_idx - cx) ** 2 / (2 * sx ** 2) + (y_idx - cy) ** 2 / (2 * sy ** 2)) + ) + gauss = PSFEstimation._normalise(gauss) + return gauss, sx, sy, amplitude + + # ------------------------------------------------------------------ + # methods + # ------------------------------------------------------------------ + + def _wiener( + self, + F_measured: np.ndarray, + F_ideal: np.ndarray, + regularization: float, + psf_size: int, + ) -> np.ndarray: + """Pseudo-Wiener PSF estimation.""" + F_psf = np.conj(F_ideal) * F_measured / (np.abs(F_ideal) ** 2 + regularization) + psf = np.real(np.fft.ifft2(F_psf)) + psf = np.fft.fftshift(psf) + psf = self._crop_centre(psf, psf_size) + return self._normalise(psf) + + def _least_squares( + self, + F_measured: np.ndarray, + F_ideal: np.ndarray, + regularization: float, + psf_size: int, + ) -> np.ndarray: + """Regularised least-squares PSF estimation.""" + abs_ideal = np.abs(F_ideal) + F_psf = np.where( + abs_ideal < regularization, + 0.0, + F_measured / (F_ideal + regularization * np.sign(F_ideal)), + ) + psf = np.real(np.fft.ifft2(F_psf)) + psf = np.fft.fftshift(psf) + psf = self._crop_centre(psf, psf_size) + return self._normalise(psf) + + # ------------------------------------------------------------------ + # main entry + # ------------------------------------------------------------------ + + def process( + self, + measured: DataField, + ideal: DataField, + method: str, + regularization: float, + psf_size: int, + ) -> tuple: + measured_data = np.asarray(measured.data, dtype=np.float64) + ideal_data = np.asarray(ideal.data, dtype=np.float64) + + F_measured = np.fft.fft2(measured_data) + F_ideal = np.fft.fft2(ideal_data) + + parameters = RecordTable() + + if method == "wiener": + psf = self._wiener(F_measured, F_ideal, regularization, psf_size) + + elif method == "least_squares": + psf = self._least_squares(F_measured, F_ideal, regularization, psf_size) + + elif method == "gaussian_fit": + raw_psf = self._wiener(F_measured, F_ideal, regularization, psf_size) + psf, sigma_x, sigma_y, amplitude = self._fit_gaussian_2d(raw_psf) + parameters = RecordTable([ + {"quantity": "sigma_x", "value": sigma_x, "unit": "px"}, + {"quantity": "sigma_y", "value": sigma_y, "unit": "px"}, + {"quantity": "amplitude", "value": amplitude, "unit": ""}, + ]) + else: + raise ValueError(f"Unknown PSF estimation method: {method!r}") + + # Build output DataField — inherit spatial metadata, adjust for psf_size + yres, xres = measured_data.shape + psf_xreal = measured.xreal * psf_size / xres + psf_yreal = measured.yreal * psf_size / yres + + psf_field = measured.replace( + data=psf, + xreal=psf_xreal, + yreal=psf_yreal, + xoff=0.0, + yoff=0.0, + ) + + return (psf_field, parameters) diff --git a/backend/nodes/super_resolution.py b/backend/nodes/super_resolution.py new file mode 100644 index 0000000..715bdb6 --- /dev/null +++ b/backend/nodes/super_resolution.py @@ -0,0 +1,126 @@ +"""Super-resolution -- combine multiple aligned scans for resolution enhancement.""" + +from __future__ import annotations + +import numpy as np +from scipy.ndimage import shift as ndimage_shift, zoom as ndimage_zoom + +from backend.node_registry import register_node +from backend.data_types import DataField + + +def _find_subpixel_shift(ref: np.ndarray, img: np.ndarray) -> tuple[float, float]: + """Estimate the (dy, dx) sub-pixel shift of *img* relative to *ref* via cross-correlation. + + Uses FFT-based cross-correlation with parabolic peak refinement. + """ + fa = np.fft.fft2(ref - ref.mean()) + fb = np.fft.fft2(img - img.mean()) + cross = np.fft.ifft2(fa * np.conj(fb)) + cc = np.abs(np.fft.fftshift(cross)) + + cy, cx = np.array(cc.shape) // 2 + peak_y, peak_x = np.unravel_index(np.argmax(cc), cc.shape) + + # Integer shift relative to centre + dy = peak_y - cy + dx = peak_x - cx + + # Parabolic sub-pixel refinement around peak + h, w = cc.shape + if 1 <= peak_y <= h - 2: + num = float(cc[peak_y - 1, peak_x] - cc[peak_y + 1, peak_x]) + den = float( + cc[peak_y - 1, peak_x] - 2.0 * cc[peak_y, peak_x] + cc[peak_y + 1, peak_x] + ) + if abs(den) > 1e-12: + dy += 0.5 * num / den + + if 1 <= peak_x <= w - 2: + num = float(cc[peak_y, peak_x - 1] - cc[peak_y, peak_x + 1]) + den = float( + cc[peak_y, peak_x - 1] - 2.0 * cc[peak_y, peak_x] + cc[peak_y, peak_x + 1] + ) + if abs(den) > 1e-12: + dx += 0.5 * num / den + + return float(dy), float(dx) + + +@register_node(display_name="Super Resolution") +class SuperResolution: + @classmethod + def INPUT_TYPES(cls): + return { + "required": { + "field1": ("DATA_FIELD",), + "upscale": ("INT", {"default": 2, "min": 2, "max": 4, "step": 1}), + }, + "optional": { + "field2": ("DATA_FIELD",), + "field3": ("DATA_FIELD",), + "field4": ("DATA_FIELD",), + }, + } + + OUTPUTS = ( + ('DATA_FIELD', 'result'), + ) + FUNCTION = "process" + + DESCRIPTION = ( + "Combine multiple aligned scans to produce a super-resolved image with higher " + "spatial resolution. Sub-pixel shifts between inputs are estimated via FFT " + "cross-correlation and used to reconstruct a finer grid. When only one field " + "is provided the image is upsampled using cubic interpolation." + ) + + def process( + self, + field1: DataField, + upscale: int, + field2: DataField | None = None, + field3: DataField | None = None, + field4: DataField | None = None, + ) -> tuple: + fields = [field1] + for f in (field2, field3, field4): + if f is not None: + fields.append(f) + + ref = np.asarray(field1.data, dtype=np.float64) + + # Upsample reference to target resolution + high_res = ndimage_zoom(ref, upscale, order=3) + weight = np.ones_like(high_res) + + if len(fields) == 1: + # Single input -- just return the upsampled reference + return (field1.replace( + data=high_res, + xreal=field1.xreal, + yreal=field1.yreal, + ),) + + # Multiple inputs -- align, upsample, and average + for extra in fields[1:]: + img = np.asarray(extra.data, dtype=np.float64) + + # Find sub-pixel shift relative to reference + dy, dx = _find_subpixel_shift(ref, img) + + # Shift in high-res coordinates + shifted = ndimage_shift(img.astype(np.float64), (-dy, -dx), order=3) + upsampled = ndimage_zoom(shifted, upscale, order=3) + + # Accumulate + high_res += upsampled + weight += 1.0 + + high_res /= weight + + return (field1.replace( + data=high_res, + xreal=field1.xreal, + yreal=field1.yreal, + ),) diff --git a/backend/nodes/tip_shape_estimate.py b/backend/nodes/tip_shape_estimate.py new file mode 100644 index 0000000..8bc0586 --- /dev/null +++ b/backend/nodes/tip_shape_estimate.py @@ -0,0 +1,351 @@ +"""Tip shape estimation — estimate SPM tip geometry from known calibration features.""" + +from __future__ import annotations + +import numpy as np + +from backend.node_registry import register_node +from backend.data_types import DataField, RecordTable + + +@register_node(display_name="Tip Shape Estimate") +class TipShapeEstimate: + @classmethod + def INPUT_TYPES(cls): + return { + "required": { + "field": ("DATA_FIELD",), + "feature_type": (["edge", "sphere", "cylinder"], {"default": "edge"}), + "feature_radius": ("FLOAT", { + "default": 100e-9, "min": 1e-9, "max": 100e-6, "step": 1e-9, + }), + "n_points": ("INT", {"default": 100, "min": 10, "max": 1000}), + } + } + + OUTPUTS = ( + ('DATA_FIELD', 'tip_shape'), + ('RECORD_TABLE', 'parameters'), + ) + FUNCTION = "process" + + DESCRIPTION = ( + "Estimate SPM tip geometry from a known calibration feature. " + "Supported features: edge (sharp step), sphere (calibration ball), " + "cylinder (calibration wire). The image of a known feature is a " + "dilation of the feature with the tip; by subtracting the known " + "feature contribution the tip shape can be recovered. " + "The 2D tip is built by revolving the 1D radial profile (axial " + "symmetry assumption). Output parameters include estimated tip " + "radius of curvature at the apex and half-cone angle. " + "Equivalent to Gwyddion's tipshape.c analysis. " + ) + + def process( + self, + field: DataField, + feature_type: str, + feature_radius: float, + n_points: int, + ) -> tuple: + data = field.data.astype(np.float64) + ny, nx = data.shape + pixel_size = (field.dx + field.dy) * 0.5 + + # ── Step 1: Extract 1D tip profile depending on feature type ────── + + if feature_type == "edge": + tip_profile_1d = self._estimate_from_edge(data, pixel_size, n_points) + + elif feature_type == "sphere": + tip_profile_1d = self._estimate_from_sphere( + data, pixel_size, feature_radius, n_points, + ) + + elif feature_type == "cylinder": + tip_profile_1d = self._estimate_from_cylinder( + data, pixel_size, feature_radius, n_points, + ) + + else: + raise ValueError( + f"Unknown feature_type {feature_type!r}. " + "Choose: edge, sphere, cylinder." + ) + + # ── Step 2: Build 2D tip by revolution (axial symmetry) ─────────── + + n_tip = n_points if n_points % 2 == 1 else n_points + 1 + ci = n_tip // 2 + offsets = np.arange(n_tip) - ci + gx, gy = np.meshgrid(offsets, offsets) + r_grid = np.sqrt(gx ** 2 + gy ** 2) + + # The 1D profile goes from r = 0 to r = max_r. + # Interpolate for every pixel in the 2D grid. + r_1d = np.linspace(0, ci, len(tip_profile_1d)) + tip_2d = np.interp(r_grid, r_1d, tip_profile_1d, right=0.0) + + # Convention: apex (centre) is the maximum; minimum is 0. + tip_2d -= tip_2d.min() + + xreal = n_tip * pixel_size + tip_field = DataField( + data=tip_2d, + xreal=xreal, + yreal=xreal, + si_unit_xy=field.si_unit_xy, + si_unit_z=field.si_unit_z, + ) + + # ── Step 3: Estimate tip parameters ─────────────────────────────── + + tip_radius, half_angle = self._estimate_parameters( + tip_profile_1d, pixel_size, + ) + + table = RecordTable([ + {"quantity": "tip_radius", "value": tip_radius, "unit": field.si_unit_z}, + {"quantity": "half_angle", "value": half_angle, "unit": "deg"}, + ]) + + return (tip_field, table) + + # ── Feature-specific estimation methods ─────────────────────────────── + + @staticmethod + def _estimate_from_edge( + data: np.ndarray, + pixel_size: float, + n_points: int, + ) -> np.ndarray: + """ + Edge feature: the tip shape is the mirror of the steepest edge + cross-section in the image. + + Find the row with the maximum gradient magnitude, extract the + cross-section at that location, mirror and normalise. + + Returns a radial profile: index 0 = apex (maximum), decreasing + outward. + """ + ny, nx = data.shape + + # Compute row-wise gradient magnitude to find the steepest edge. + grad = np.abs(np.diff(data, axis=1)) + row_grad = grad.sum(axis=1) + best_row = int(np.argmax(row_grad)) + + profile = data[best_row, :] + + # Mirror: the tip is the complement of the edge profile. + tip_raw = np.max(profile) - profile[::-1] + + # Find the peak of the mirrored profile and take the radial + # (half) profile from apex outward. + peak = int(np.argmax(tip_raw)) + # Use the longer side from the peak to preserve resolution. + left = tip_raw[:peak + 1][::-1] # apex to left edge, reversed + right = tip_raw[peak:] # apex to right edge + half = left if len(left) >= len(right) else right + half = half.copy() + # Ensure monotonically decreasing from apex. + for i in range(1, len(half)): + if half[i] > half[i - 1]: + half[i] = half[i - 1] + + # Resample to n_points. + x_raw = np.linspace(0, 1, len(half)) + x_out = np.linspace(0, 1, n_points) + tip_profile = np.interp(x_out, x_raw, half) + + return tip_profile + + @staticmethod + def _estimate_from_sphere( + data: np.ndarray, + pixel_size: float, + radius: float, + n_points: int, + ) -> np.ndarray: + """ + Sphere feature: dilation model z_measured = z_sphere (+) z_tip. + + Extract radial profile from the highest point and subtract the + ideal sphere contribution to recover the tip profile. + + Returns a radial profile: index 0 = apex (maximum), decreasing + outward. + """ + ny, nx = data.shape + + # Find the highest point (apex of the imaged sphere). + peak_idx = np.unravel_index(np.argmax(data), data.shape) + cy, cx = peak_idx + + # Extract radial profile by averaging azimuthally. + max_r = min(cy, ny - 1 - cy, cx, nx - 1 - cx) + if max_r < 2: + max_r = min(ny, nx) // 2 + + Y, X = np.ogrid[:ny, :nx] + r_map = np.sqrt(((X - cx) * pixel_size) ** 2 + ((Y - cy) * pixel_size) ** 2) + + n_bins = min(max_r, n_points) + r_edges = np.linspace(0, max_r * pixel_size, n_bins + 1) + radial_profile = np.zeros(n_bins) + for i in range(n_bins): + mask = (r_map >= r_edges[i]) & (r_map < r_edges[i + 1]) + if mask.any(): + radial_profile[i] = data[mask].mean() + elif i > 0: + radial_profile[i] = radial_profile[i - 1] + + r_centres = 0.5 * (r_edges[:-1] + r_edges[1:]) + + # Ideal sphere profile: z_sphere(r) = sqrt(R^2 - r^2) for r < R, else 0. + sphere_profile = np.where( + r_centres < radius, + np.sqrt(np.maximum(radius ** 2 - r_centres ** 2, 0.0)), + 0.0, + ) + + # Tip profile = measured - sphere (dilation subtraction). + tip_raw = radial_profile - sphere_profile + # Shift so that the apex (index 0) is the maximum. + tip_raw = tip_raw - tip_raw.min() + # Ensure monotonically decreasing from apex outward by clamping. + for i in range(1, len(tip_raw)): + if tip_raw[i] > tip_raw[i - 1]: + tip_raw[i] = tip_raw[i - 1] + + # Resample to n_points. + x_raw = np.linspace(0, 1, len(tip_raw)) + x_out = np.linspace(0, 1, n_points) + tip_profile = np.interp(x_out, x_raw, tip_raw) + + return tip_profile + + @staticmethod + def _estimate_from_cylinder( + data: np.ndarray, + pixel_size: float, + radius: float, + n_points: int, + ) -> np.ndarray: + """ + Cylinder feature: extract cross-section perpendicular to the + cylinder axis and subtract ideal cylinder profile. + + The cylinder axis is assumed to run along the direction with the + least height variation. + + Returns a radial profile: index 0 = apex (maximum), decreasing + outward. + """ + ny, nx = data.shape + + # Determine cylinder axis: compare row-wise vs column-wise variance. + row_var = np.var(np.diff(data, axis=1)) + col_var = np.var(np.diff(data, axis=0)) + + if row_var > col_var: + # Cylinder axis along columns -> cross-section along rows. + profile = data.mean(axis=0) + else: + # Cylinder axis along rows -> cross-section along columns. + profile = data.mean(axis=1) + + n_prof = len(profile) + peak = int(np.argmax(profile)) + + # Ideal cylinder cross-section: z = sqrt(R^2 - x^2) for |x| < R. + x_phys = (np.arange(n_prof) - peak) * pixel_size + cyl_profile = np.where( + np.abs(x_phys) < radius, + np.sqrt(np.maximum(radius ** 2 - x_phys ** 2, 0.0)), + 0.0, + ) + + # Tip = measured - cylinder. + tip_raw = profile - cyl_profile + tip_raw -= tip_raw.min() + + # Take the radial (half) profile from the peak outward. + left = tip_raw[:peak + 1][::-1] + right = tip_raw[peak:] + half = left if len(left) >= len(right) else right + # Ensure monotonically decreasing from apex. + for i in range(1, len(half)): + if half[i] > half[i - 1]: + half[i] = half[i - 1] + + # Resample to n_points. + x_raw = np.linspace(0, 1, len(half)) + x_out = np.linspace(0, 1, n_points) + tip_profile = np.interp(x_out, x_raw, half) + + return tip_profile + + # ── Parameter estimation ────────────────────────────────────────────── + + @staticmethod + def _estimate_parameters( + tip_profile: np.ndarray, + pixel_size: float, + ) -> tuple[float, float]: + """ + Estimate tip radius of curvature at the apex and half-cone angle + from the 1D radial profile. + + tip_radius: fitted from the parabolic approximation near the apex, + z(r) ~ z_max - r^2 / (2R) => R = r^2 / (2 * (z_max - z(r))) + + half_angle: from the slope of the tip walls in the outer half, + tan(half_angle) = dz/dr => half_angle = arctan(slope) + """ + n = len(tip_profile) + r = np.linspace(0, (n - 1) * pixel_size, n) + + # ── Tip radius from apex curvature ──────────────────────────────── + # Use a few points near the apex for a parabolic fit: z = a - b*r^2 + n_apex = max(3, n // 10) + r_apex = r[:n_apex] + z_apex = tip_profile[:n_apex] + + if len(r_apex) >= 2 and r_apex[-1] > 0: + # Fit z = c0 + c1 * r^2 + A = np.vstack([np.ones(n_apex), r_apex ** 2]).T + try: + coeffs = np.linalg.lstsq(A, z_apex, rcond=None)[0] + c1 = coeffs[1] + # z = z_max - r^2/(2R) => c1 = -1/(2R) => R = -1/(2*c1) + if c1 < 0: + tip_radius = -1.0 / (2.0 * c1) + else: + tip_radius = float('inf') + except np.linalg.LinAlgError: + tip_radius = float('inf') + else: + tip_radius = float('inf') + + # ── Half-angle from outer wall slope ────────────────────────────── + # Use the outer 50% of the profile. + mid = n // 2 + if mid < n - 1: + r_outer = r[mid:] + z_outer = tip_profile[mid:] + if len(r_outer) >= 2: + dr = r_outer[-1] - r_outer[0] + dz = z_outer[-1] - z_outer[0] + if dr > 0: + slope = abs(dz / dr) + half_angle = np.degrees(np.arctan(slope)) + else: + half_angle = 0.0 + else: + half_angle = 0.0 + else: + half_angle = 0.0 + + return tip_radius, half_angle diff --git a/docs/nodes/Calibration.md b/docs/nodes/Calibration.md new file mode 100644 index 0000000..44d4d8d --- /dev/null +++ b/docs/nodes/Calibration.md @@ -0,0 +1,37 @@ +# Calibration + +Apply lateral and height calibration corrections to a DATA_FIELD. Equivalent to Gwyddion's calibrate.c functionality. + +## Inputs + +| Name | Type | Required | Description | +|------|------|----------|-------------| +| field | DATA_FIELD | Yes | Data to calibrate | + +## Outputs + +| Name | Type | Description | +|------|------|-------------| +| result | DATA_FIELD | Calibrated field | + +## Controls + +| Name | Type | Default | Description | +|------|------|---------|-------------| +| xy_mode | dropdown | keep | Lateral calibration mode: keep (no change), set_size (set explicit physical size), scale (multiply by factor) | +| z_mode | dropdown | keep | Height calibration mode: keep (no change), set_range (linear map to new min/max), scale (multiply by factor), offset (add constant) | +| xreal_new | float | 1e-6 | New x physical size (shown when xy_mode = set_size) | +| yreal_new | float | 1e-6 | New y physical size (shown when xy_mode = set_size) | +| xy_scale | float | 1.0 | Lateral scale factor (shown when xy_mode = scale) | +| z_min | float | 0.0 | New z minimum (shown when z_mode = set_range) | +| z_max | float | 1e-9 | New z maximum (shown when z_mode = set_range) | +| z_scale | float | 1.0 | Z scale factor (shown when z_mode = scale) | +| z_offset | float | 0.0 | Z offset value (shown when z_mode = offset) | +| xy_unit | string | (empty) | New XY unit string; leave empty to keep current | +| z_unit | string | (empty) | New Z unit string; leave empty to keep current | + +## Notes + +- Controls are conditionally shown based on the selected mode via show_when_widget_value. +- set_range linearly maps data from the current [min, max] to the specified [z_min, z_max]. If data is constant, z_min is applied uniformly. +- Unit fields accept any string (e.g. "m", "nm", "V"). An empty string preserves the existing unit. diff --git a/docs/nodes/DWT Anisotropy.md b/docs/nodes/DWT Anisotropy.md new file mode 100644 index 0000000..a180475 --- /dev/null +++ b/docs/nodes/DWT Anisotropy.md @@ -0,0 +1,32 @@ +# DWT Anisotropy + +Quantify surface anisotropy using a multi-level 2-D Haar wavelet decomposition. At each decomposition level, horizontal and vertical detail energies are compared to produce an X/Y energy ratio. Equivalent to Gwyddion's dwtanisotropy.c. + +## Inputs + +| Name | Type | Required | Description | +|------|------|----------|-------------| +| field | DATA_FIELD | Yes | Input surface field | + +## Outputs + +| Name | Type | Description | +|------|------|-------------| +| anisotropy_map | DATA_FIELD | Per-pixel anisotropy ratio map (averaged across decomposition levels) | +| statistics | RECORD_TABLE | Per-level X/Y energy ratios and anisotropy flags | + +## Controls + +| Name | Type | Default | Description | +|------|------|---------|-------------| +| n_levels | INT | 4 | Number of wavelet decomposition levels (1--10) | +| ratio_threshold | FLOAT | 0.2 | Deviation from 1.0 required to flag a level as anisotropic (0.001--10.0) | + +## Notes + +- The decomposition uses the Haar wavelet (db1), which splits each 2x2 block into approximation (LL), horizontal detail (LH), vertical detail (HL), and diagonal detail (HH) coefficients. +- Energy ratios are computed as sum(HL^2) / sum(LH^2) at each level. HL captures horizontal features (edges running left-right), while LH captures vertical features (edges running top-bottom). +- Ratio > 1 means the surface has more horizontal features; ratio < 1 means more vertical features; ratio near 1 indicates isotropy. +- The input is padded to the next power of 2 if necessary; padding uses edge values. +- The anisotropy map is built by upsampling each level's per-pixel ratio and averaging across levels. +- The statistics table includes per-level x_energy, y_energy, ratio, and a boolean anisotropic flag based on the ratio_threshold control. diff --git a/docs/nodes/Displacement Field.md b/docs/nodes/Displacement Field.md new file mode 100644 index 0000000..895ff8b --- /dev/null +++ b/docs/nodes/Displacement Field.md @@ -0,0 +1,34 @@ +# Displacement Field + +Distort an image using synthetic displacement fields. Supports correlated Gaussian noise and tear-line distortion modes. Equivalent to Gwyddion's displfield.c module. + +## Inputs + +| Name | Type | Required | Description | +|------|------|----------|-------------| +| field | DATA_FIELD | Yes | Input field to distort | + +## Outputs + +| Name | Type | Description | +|------|------|-------------| +| result | DATA_FIELD | Distorted field | + +## Controls + +| Name | Type | Default | Description | +|------|------|---------|-------------| +| method | dropdown | gaussian_1d | Distortion method: gaussian_1d, gaussian_2d, or tear | +| sigma | float | 5.0 | Distortion amplitude in pixels | +| tau | float | 20.0 | Lateral correlation length in pixels | +| density | float | 0.02 | Tear density — fraction of rows that become tear lines (tear mode only) | +| seed | int | 42 | Random seed for reproducibility | + +## Notes + +- **gaussian_1d** generates a 1D correlated random displacement applied only in the x direction. All rows share the same displacement profile, simulating a systematic lateral distortion. +- **gaussian_2d** generates independent 2D correlated random displacements in both x and y. This produces a more general warping of the image. +- **tear** mode simulates scanning artifacts where random horizontal tear lines introduce sudden x-offsets that decay exponentially away from the tear row. This is useful for simulating or studying piezo slip artifacts in SPM data. +- The **sigma** parameter controls the magnitude of the displacement. Larger values produce more extreme distortion. +- The **tau** parameter controls the spatial correlation length. A larger tau produces smoother, more slowly varying displacement fields. The ratio sigma/tau roughly determines the local strain. +- For realistic scanning artifacts, use tear mode with low density (0.01--0.05) and moderate sigma. diff --git a/docs/nodes/Distribution Coercion.md b/docs/nodes/Distribution Coercion.md new file mode 100644 index 0000000..fd73239 --- /dev/null +++ b/docs/nodes/Distribution Coercion.md @@ -0,0 +1,31 @@ +# Distribution Coercion + +Transform pixel values so their distribution matches a target shape (uniform, Gaussian, or discrete levels) using rank-based reassignment. Equivalent to Gwyddion's coerce.c module. + +## Inputs + +| Name | Type | Required | Description | +|------|------|----------|-------------| +| field | DATA_FIELD | Yes | Input field whose value distribution will be transformed | + +## Outputs + +| Name | Type | Description | +|------|------|-------------| +| result | DATA_FIELD | Field with pixel values reassigned to match the target distribution | + +## Controls + +| Name | Type | Default | Description | +|------|------|---------|-------------| +| distribution | dropdown | uniform | Target distribution shape: uniform, gaussian, or levels | +| n_levels | INT | 4 | Number of discrete output levels (2–1000); visible only for levels mode | +| processing | dropdown | field | Processing scope: field (entire array at once) or rows (line-by-line) | + +## Notes + +- The transformation is rank-based: pixels are sorted, then reassigned values drawn from the target distribution in sorted order. This preserves the relative ordering of all pixel values. +- Uniform mode spreads values evenly between the original minimum and maximum. +- Gaussian mode maps ranks to the inverse normal CDF, scaled to match the original mean and standard deviation. +- Levels mode quantizes the data into a fixed number of evenly spaced discrete values, useful for terrace-like visualization or discrete height analysis. +- Row mode applies the transformation independently to each scan line, which can correct line-to-line distribution variations in SPM data. diff --git a/docs/nodes/Grain Visualization.md b/docs/nodes/Grain Visualization.md new file mode 100644 index 0000000..a287498 --- /dev/null +++ b/docs/nodes/Grain Visualization.md @@ -0,0 +1,34 @@ +# Grain Visualization + +Visualize labeled grains as geometric shapes — inscribed discs, bounding boxes, centroid markers, or fitted ellipses. Produces a mask image with the chosen shapes and a labeled field where each grain has a unique integer value. Equivalent to Gwyddion's grain selection visualization (grain_makesel). + +## Inputs + +| Name | Type | Required | Description | +|------|------|----------|-------------| +| field | DATA_FIELD | Yes | Source data field providing spatial reference and calibration | +| mask | IMAGE | Yes | Binary grain mask (white = grain region) | + +## Outputs + +| Name | Type | Description | +|------|------|-------------| +| result | IMAGE | Visualization mask with geometric shapes drawn for each grain | +| labeled | DATA_FIELD | Labeled field where each connected grain region has a unique integer value | + +## Controls + +| Name | Type | Default | Description | +|------|------|---------|-------------| +| style | dropdown | inscribed_disc | Visualization style: inscribed_disc, bounding_box, centroid, or ellipse | +| fill | BOOLEAN | False | When enabled, shapes are filled; when disabled, only outlines are drawn | + +## Notes + +- **inscribed_disc**: Draws a circle at each grain centroid with the inscribed disc radius (maximum of the distance transform within the grain). +- **bounding_box**: Draws an axis-aligned rectangle around each grain's bounding box. +- **centroid**: Draws small cross markers at each grain's centroid. +- **ellipse**: Fits an ellipse to each grain using the inertia tensor to determine the major/minor axes and orientation angle. +- When fill is disabled, outlines are drawn with approximately 1-2 pixel thickness. +- The labeled output assigns each connected grain region a unique positive integer; background pixels are zero. +- The mask and field must have the same pixel dimensions. diff --git a/docs/nodes/Logistic Classification.md b/docs/nodes/Logistic Classification.md new file mode 100644 index 0000000..ff6d840 --- /dev/null +++ b/docs/nodes/Logistic Classification.md @@ -0,0 +1,37 @@ +# Logistic Classification + +Classify surface features using logistic regression on engineered height-derived features. Optionally accepts a training mask; otherwise an Otsu-based threshold generates pseudo-labels automatically. + +## Inputs + +| Name | Type | Required | Description | +|------|------|----------|-------------| +| field | DATA_FIELD | Yes | Input topographic surface to classify | +| training_mask | IMAGE | No | Optional training labels — masked pixels are treated as the positive class | + +## Outputs + +| Name | Type | Description | +|------|------|-------------| +| mask | IMAGE | Binary classification result (0 or 255) | +| probability | DATA_FIELD | Per-pixel probability from the logistic model (values in [0, 1]) | + +## Controls + +| Name | Type | Default | Description | +|------|------|---------|-------------| +| use_gaussians | BOOLEAN | True | Include Gaussian blur features at multiple scales | +| n_gaussians | INT | 4 | Number of Gaussian scales (1–10). Only shown when use_gaussians is True | +| use_sobel | BOOLEAN | True | Include Sobel gradient features (horizontal and vertical) | +| use_laplacian | BOOLEAN | True | Include Laplacian (sum of second differences) feature | +| regularization | FLOAT | 1.0 | L2 regularization strength lambda (0.0–10.0) | +| max_iter | INT | 500 | Maximum gradient descent iterations (10–5000) | +| seed | INT | 42 | Random seed for reproducibility (0–999999) | + +## Notes + +- **Feature engineering:** The classifier always uses normalized raw height as a feature. Gaussian blurs at scales 2^0, 2^1, ..., 2^(n-1) capture multi-scale smoothness. Sobel gradients detect edges, and the Laplacian highlights curvature. All features are standardized to zero mean and unit variance before training. +- **L2 regularization:** The regularization parameter controls overfitting by penalizing large weights. Higher values produce smoother, more generalizable decision boundaries. The bias term is not regularized. +- **Logistic regression vs neural networks:** Logistic regression is a linear classifier — it learns a single hyperplane in feature space. For complex, highly non-linear boundaries a neural network may be more appropriate, but logistic regression is fast, interpretable, and often sufficient when combined with good feature engineering. +- **Unsupervised mode:** When no training mask is provided, the node uses an Otsu-like threshold on the raw height to generate pseudo-labels, then trains the classifier on those labels. This can improve on simple thresholding because the classifier leverages multi-scale and gradient features. +- Equivalent to Gwyddion's logistic.c classification functionality. diff --git a/docs/nodes/Mark Disconnected.md b/docs/nodes/Mark Disconnected.md new file mode 100644 index 0000000..0acc7b2 --- /dev/null +++ b/docs/nodes/Mark Disconnected.md @@ -0,0 +1,30 @@ +# Mark Disconnected + +Mark topologically disconnected (isolated) surface regions. A morphological opening followed by closing builds a smooth defect-free reference surface; pixels whose deviation from that reference exceeds the sensitivity threshold are flagged. Equivalent to Gwyddion's mark_disconn.c module. + +## Inputs + +| Name | Type | Required | Description | +|------|------|----------|-------------| +| field | DATA_FIELD | Yes | Input surface to analyse for disconnected regions | + +## Outputs + +| Name | Type | Description | +|------|------|-------------| +| mask | IMAGE | Binary mask (0/255) marking detected disconnected regions | + +## Controls + +| Name | Type | Default | Description | +|------|------|---------|-------------| +| defect_type | dropdown | positive | Which direction of outliers to detect: positive (bumps), negative (pits), or both | +| radius | INT | 5 | Morphological filter radius in pixels (1--100). Larger values smooth over bigger features | +| threshold | FLOAT | 0.1 | Sensitivity threshold as a fraction of the max residual range (0.001--1.0). Lower values detect smaller defects | + +## Notes + +- The algorithm applies grey-scale morphological opening then closing with a disk structuring element to produce a defect-free reference. The difference between the original surface and this reference highlights isolated regions. +- Increase the radius to ignore larger surface features and only flag truly disconnected regions. +- Lower the threshold to catch subtler defects; raise it to reduce false positives. +- The resulting mask can be fed to Laplace Interpolation or Fractal Interpolation to fill the detected regions. diff --git a/docs/nodes/Mask Noisify.md b/docs/nodes/Mask Noisify.md new file mode 100644 index 0000000..e96bc6a --- /dev/null +++ b/docs/nodes/Mask Noisify.md @@ -0,0 +1,32 @@ +# Mask Noisify + +Add random noise to a binary mask by flipping pixels near boundaries. Equivalent to Gwyddion's mask_noisify.c. + +## Inputs + +| Name | Type | Required | Description | +|------|------|----------|-------------| +| mask | IMAGE | Yes | Binary mask to perturb | +| field | DATA_FIELD | No | Optional field for preview background display | + +## Outputs + +| Name | Type | Description | +|------|------|-------------| +| mask | IMAGE | Noisified binary mask | + +## Controls + +| Name | Type | Default | Description | +|------|------|---------|-------------| +| density | FLOAT | 0.1 | Fraction of candidate pixels to flip (0.0--1.0) | +| direction | dropdown | both | Which pixels to perturb: add (grow mask), remove (shrink mask), or both | +| boundaries_only | BOOLEAN | True | Only modify pixels adjacent to a mask boundary | +| seed | INT | 42 | Random seed for reproducible results (0--999999) | + +## Notes + +- Boundary detection uses four-neighbour comparison (up, down, left, right) via np.roll. A pixel is a boundary pixel if it differs from at least one of its four neighbours. +- Direction modes control which candidates are eligible: "add" selects only unmasked boundary pixels, "remove" selects only masked boundary pixels, and "both" selects all boundary pixels. +- Setting boundaries_only to False allows any pixel in the mask to be a candidate, not just those at edges. +- The seed parameter ensures deterministic output for a given input, which is useful for repeatable experiments. diff --git a/docs/nodes/Mask Shift.md b/docs/nodes/Mask Shift.md new file mode 100644 index 0000000..7a3ceba --- /dev/null +++ b/docs/nodes/Mask Shift.md @@ -0,0 +1,31 @@ +# Mask Shift + +Translate a binary mask by an integer pixel offset. Choose how out-of-bounds regions are filled: zero (empty), wrap (periodic roll), or mirror (reflected padding). Equivalent to Gwyddion's mask_shift.c. + +## Inputs + +| Name | Type | Required | Description | +|------|------|----------|-------------| +| mask | IMAGE | Yes | Binary mask to shift | +| field | DATA_FIELD | No | Optional field for preview background display | + +## Outputs + +| Name | Type | Description | +|------|------|-------------| +| mask | IMAGE | Shifted binary mask | + +## Controls + +| Name | Type | Default | Description | +|------|------|---------|-------------| +| shift_x | INT | 0 | Horizontal shift in pixels (-1000 to 1000). Positive values shift right. | +| shift_y | INT | 0 | Vertical shift in pixels (-1000 to 1000). Positive values shift down. | +| border_mode | dropdown | zero | How to handle edges: zero fills vacated region with empty mask, wrap rolls periodically, mirror reflects at boundaries | + +## Notes + +- Shift values are in pixels, not physical units. +- **zero** mode: the mask is rolled and the vacated strip is cleared to zero (unmasked). Useful when the shifted region should not wrap around. +- **wrap** mode: uses periodic rolling (`np.roll`). The total number of masked pixels is preserved. Suitable for periodic or tiled data. +- **mirror** mode: pads with reflected values before cropping, so edges are filled with a mirrored copy of the mask boundary. Avoids hard cutoffs at the border. diff --git a/docs/nodes/Neural Classification.md b/docs/nodes/Neural Classification.md new file mode 100644 index 0000000..25e6432 --- /dev/null +++ b/docs/nodes/Neural Classification.md @@ -0,0 +1,35 @@ +# Neural Classification + +Classify surface pixels into two classes using a simple two-layer feedforward neural network with sigmoid activations. Features are extracted via multi-scale Gaussian filtering. Equivalent in purpose to Gwyddion's neural.c classifier. + +## Inputs + +| Name | Type | Required | Description | +|------|------|----------|-------------| +| field | DATA_FIELD | Yes | Input surface to classify | +| training_mask | IMAGE | No | Training labels: 0 = class A, 255 = class B. When omitted the network uses unsupervised self-labelling | + +## Outputs + +| Name | Type | Description | +|------|------|-------------| +| mask | IMAGE | Binary classification mask (0 or 255) | +| probability | DATA_FIELD | Per-pixel probability of belonging to class B (values in 0-1) | + +## Controls + +| Name | Type | Default | Description | +|------|------|---------|-------------| +| n_gaussians | INT | 4 | Number of Gaussian blur scales for feature extraction (1-10). Each scale uses sigma = 2^i | +| n_hidden | INT | 16 | Number of neurons in the hidden layer (4-128) | +| train_steps | INT | 200 | Number of gradient descent iterations (10-5000) | +| seed | INT | 42 | Random seed for weight initialisation (0-999999) | + +## Notes + +- Feature extraction applies Gaussian blur at multiple scales (sigma = 1, 2, 4, 8, ...) to capture both fine and coarse surface structure. Each feature is normalised to zero mean and unit variance before training. +- The network architecture is input -> hidden (sigmoid) -> output (sigmoid), trained with binary cross-entropy loss and standard backpropagation. +- When a training mask is provided, the network learns in supervised mode using all pixels (0 pixels as class A targets, 255 pixels as class B targets). +- Without a training mask, the node uses an unsupervised approach: the random initial weights produce an initial classification which is then refined by self-training for a small number of steps. +- Increasing n_hidden or train_steps improves capacity but slows computation. For most surfaces, the defaults work well. +- The probability output can be fed into a Threshold Mask node for adjustable post-classification thresholding. diff --git a/docs/nodes/PSF Estimation.md b/docs/nodes/PSF Estimation.md new file mode 100644 index 0000000..6f3c522 --- /dev/null +++ b/docs/nodes/PSF Estimation.md @@ -0,0 +1,35 @@ +# PSF Estimation + +Estimate a point spread function (PSF) from a measured (blurred) image and an ideal (sharp) reference. The estimated PSF can be fed into the Deconvolution node to restore other images acquired under the same conditions. Equivalent to Gwyddion's psf.c / psf-fit.c modules. + +## Inputs + +| Name | Type | Required | Description | +|------|------|----------|-------------| +| measured | DATA_FIELD | Yes | Measured (blurred) image | +| ideal | DATA_FIELD | Yes | Ideal (sharp) reference image | + +## Outputs + +| Name | Type | Description | +|------|------|-------------| +| psf | DATA_FIELD | Estimated point spread function | +| parameters | RECORD_TABLE | Fitted PSF parameters (populated by gaussian_fit method) | + +## Controls + +| Name | Type | Default | Description | +|------|------|---------|-------------| +| method | dropdown | wiener | Estimation method: wiener, least_squares, or gaussian_fit | +| regularization | FLOAT | 0.01 | Regularization parameter (1e-6 to 1.0) | +| psf_size | INT | 32 | Size of the estimated PSF in pixels (4 to 128) | + +## Notes + +- **Wiener**: Pseudo-Wiener deconvolution in the frequency domain. Computes `conj(F_ideal) * F_measured / (|F_ideal|^2 + regularization)`. Fast and robust for most cases. +- **Least-squares**: Regularised frequency-domain division. Zeros out components where the ideal signal is too weak, avoiding noise amplification. +- **Gaussian fit**: Estimates the PSF via the Wiener method, then fits a 2D Gaussian to the result using moment analysis. Returns the smooth fitted PSF and its parameters (sigma_x, sigma_y, amplitude). Useful when the PSF is known to be approximately Gaussian. +- The **regularization** parameter controls the noise/sharpness tradeoff. Smaller values yield sharper PSF estimates but amplify noise. Start with the default (0.01) and adjust if needed. +- The PSF output is always normalized to sum to 1 and cropped to `psf_size x psf_size` pixels centered on the peak. +- Connect the PSF output to the Deconvolution node for image restoration with the estimated kernel. +- Both input fields should have the same pixel dimensions for best results. diff --git a/docs/nodes/Pixel Classification.md b/docs/nodes/Pixel Classification.md new file mode 100644 index 0000000..466d6ee --- /dev/null +++ b/docs/nodes/Pixel Classification.md @@ -0,0 +1,34 @@ +# Pixel Classification + +Classify pixels into discrete classes based on height, slope, and/or curvature using threshold or clustering methods. Equivalent to Gwyddion's classify.c module. + +## Inputs + +| Name | Type | Required | Description | +|------|------|----------|-------------| +| field | DATA_FIELD | Yes | Input surface | + +## Outputs + +| Name | Type | Description | +|------|------|-------------| +| classified | DATA_FIELD | Integer class labels (0 to n_classes-1) | +| mask | IMAGE | Binary mask of the first class (class 0) | + +## Controls + +| Name | Type | Default | Description | +|------|------|---------|-------------| +| n_classes | INT | 3 | Number of output classes (2–10) | +| feature | dropdown | height | Feature used for classification: height, slope, curvature, height_slope, or all | +| method | dropdown | otsu | Thresholding method: otsu, equal_range, or quantile | + +## Notes + +- **Feature types**: "height" uses raw data values; "slope" uses gradient magnitude (via `np.gradient`); "curvature" uses the Laplacian (sum of second derivatives). "height_slope" and "all" stack multiple features. +- **Threshold methods** (single-feature only): + - *otsu*: Multi-Otsu thresholding that finds thresholds minimising intra-class variance. + - *equal_range*: Divides the feature value range into equal-width intervals. + - *quantile*: Divides by quantiles so each class contains roughly the same number of pixels. +- **Multi-feature modes** ("height_slope", "all") ignore the method setting and use k-means clustering instead. Each feature is normalised to [0, 1] before clustering. +- The mask output contains class 0 only — use the classified field for access to all class labels. diff --git a/docs/nodes/Presentation Ops.md b/docs/nodes/Presentation Ops.md new file mode 100644 index 0000000..da57fab --- /dev/null +++ b/docs/nodes/Presentation Ops.md @@ -0,0 +1,31 @@ +# Presentation Ops + +Manage presentation overlays on data fields. Provides logarithmic scaling, normalisation extraction, overlay attachment, and linear blending. Equivalent to Gwyddion's `presentationops.c` module. + +## Inputs + +| Name | Type | Required | Description | +|------|------|----------|-------------| +| field | DATA_FIELD | Yes | Primary data field | +| overlay | DATA_FIELD | No | Field to attach or blend as an overlay (used by attach and blend modes) | + +## Outputs + +| Name | Type | Description | +|------|------|-------------| +| result | DATA_FIELD | Processed data field | + +## Controls + +| Name | Type | Default | Description | +|------|------|---------|-------------| +| operation | dropdown | logscale | Operation mode: logscale, extract_presentation, attach, or blend | +| blend_factor | FLOAT | 0.5 | Linear blend ratio between field (0.0) and overlay (1.0); only shown in blend mode (0.0--1.0) | + +## Notes + +- **logscale**: Shifts the data so the minimum becomes a small positive value, then applies log10. Useful for data with large dynamic range such as power spectral densities or FFT magnitudes. +- **extract_presentation**: Normalises the field to the [0, 1] range. Handy for generating a quick visual overview or feeding into colour-mapping nodes. +- **attach**: Replaces the field data with the overlay data. If the overlay has different dimensions it is resampled with cubic interpolation to match. +- **blend**: Linearly mixes `(1 - blend_factor) * field + blend_factor * overlay`. The overlay is resampled if its dimensions differ from the field. +- Overlay resampling uses `scipy.ndimage.zoom` with third-order (cubic) interpolation. diff --git a/docs/nodes/Super Resolution.md b/docs/nodes/Super Resolution.md new file mode 100644 index 0000000..ad97086 --- /dev/null +++ b/docs/nodes/Super Resolution.md @@ -0,0 +1,32 @@ +# Super Resolution + +Combine multiple aligned scans to produce a super-resolved image with higher spatial resolution. Sub-pixel shifts between inputs are estimated via FFT cross-correlation and used to reconstruct a finer grid. Equivalent to Gwyddion's superresolution.c. + +## Inputs + +| Name | Type | Required | Description | +|------|------|----------|-------------| +| field1 | DATA_FIELD | Yes | Reference image | +| field2 | DATA_FIELD | No | Second image | +| field3 | DATA_FIELD | No | Third image | +| field4 | DATA_FIELD | No | Fourth image | + +## Outputs + +| Name | Type | Description | +|------|------|-------------| +| result | DATA_FIELD | Super-resolved image | + +## Controls + +| Name | Type | Default | Description | +|------|------|---------|-------------| +| upscale | INT | 2 | Upscaling factor (2, 3, or 4) | + +## Notes + +- When only one field is provided the image is upsampled using cubic interpolation (no multi-image enhancement). +- Additional fields are aligned to the reference using FFT-based cross-correlation with parabolic sub-pixel refinement, then averaged on the high-resolution grid. +- Providing more images generally improves the result because each scan samples slightly different sub-pixel positions. +- All input fields should have the same pixel dimensions and physical size for correct alignment. +- The upscale factor controls the output resolution multiplier (2x, 3x, or 4x the input dimensions). diff --git a/docs/nodes/Tip Shape Estimate.md b/docs/nodes/Tip Shape Estimate.md new file mode 100644 index 0000000..e9f7f41 --- /dev/null +++ b/docs/nodes/Tip Shape Estimate.md @@ -0,0 +1,33 @@ +# Tip Shape Estimate + +Estimate SPM tip geometry from a known calibration feature. The image of a calibration feature (sharp edge, sphere, cylinder) is a dilation of the true feature shape with the tip. By subtracting the known feature contribution the tip shape can be recovered. The 2D tip is built by revolving the extracted 1D radial profile assuming axial symmetry. Equivalent to Gwyddion's tipshape.c analysis. + +## Inputs + +| Name | Type | Required | Description | +|------|------|----------|-------------| +| field | DATA_FIELD | Yes | Image of a known calibration feature | + +## Outputs + +| Name | Type | Description | +|------|------|-------------| +| tip_shape | DATA_FIELD | Estimated tip shape as a 2D image (apex = maximum, edges = 0) | +| parameters | RECORD_TABLE | Estimated tip parameters: tip_radius (radius of curvature at apex) and half_angle (cone half-angle from the tip walls) | + +## Controls + +| Name | Type | Default | Description | +|------|------|---------|-------------| +| feature_type | dropdown | edge | Type of calibration feature: edge (sharp step), sphere (calibration ball), or cylinder (calibration wire) | +| feature_radius | FLOAT | 100 nm | Known radius of the calibration feature in metres; used for sphere and cylinder types (1 nm to 100 um) | +| n_points | INT | 100 | Number of points in the output tip profile / side length of the tip grid (10-1000) | + +## Notes + +- **Calibration features**: For best results use well-characterised calibration standards. Sharp edges give the most direct tip estimate. Sphere and cylinder features require accurate knowledge of the feature radius. +- **Dilation model**: The measured image is the morphological dilation of the true surface with the tip shape. For a known surface feature, subtracting the feature profile from the measured profile yields the tip contribution. +- **Axial symmetry assumption**: The 2D tip shape is built by revolving the 1D radial profile around the apex. This assumes the tip is rotationally symmetric, which is a reasonable first approximation for most SPM tips but may not capture asymmetric wear or contamination. +- **Use with Tip Deconvolution**: The estimated tip can be fed directly into the Tip Deconvolution node to reconstruct the true surface from measured images. Ensure the pixel size of the tip matches the image pixel size. +- **feature_radius** is only used for sphere and cylinder feature types; it is ignored for edge estimation. +- Equivalent to Gwyddion's tipshape.c tip characterisation routines. diff --git a/tests/node_tests/calibration.py b/tests/node_tests/calibration.py new file mode 100644 index 0000000..21e23a6 --- /dev/null +++ b/tests/node_tests/calibration.py @@ -0,0 +1,77 @@ +import numpy as np +from tests.node_tests._shared import make_field + + +def test_keep_mode_unchanged(): + from backend.nodes.calibration import Calibration + node = Calibration() + field = make_field(data=np.array([[1.0, 2.0], [3.0, 4.0]])) + + result, = node.process( + field, + xy_mode="keep", z_mode="keep", + xreal_new=1e-6, yreal_new=1e-6, xy_scale=1.0, + z_min=0.0, z_max=1e-9, z_scale=1.0, z_offset=0.0, + xy_unit="", z_unit="", + ) + + assert np.array_equal(result.data, field.data) + assert result.xreal == field.xreal + assert result.yreal == field.yreal + assert result.si_unit_xy == field.si_unit_xy + assert result.si_unit_z == field.si_unit_z + + +def test_set_size(): + from backend.nodes.calibration import Calibration + node = Calibration() + field = make_field(data=np.array([[1.0, 2.0], [3.0, 4.0]])) + + result, = node.process( + field, + xy_mode="set_size", z_mode="keep", + xreal_new=5e-6, yreal_new=3e-6, xy_scale=1.0, + z_min=0.0, z_max=1e-9, z_scale=1.0, z_offset=0.0, + xy_unit="", z_unit="", + ) + + assert result.xreal == 5e-6 + assert result.yreal == 3e-6 + assert np.array_equal(result.data, field.data) + + +def test_z_scale(): + from backend.nodes.calibration import Calibration + node = Calibration() + data = np.array([[1.0, 2.0], [3.0, 4.0]]) + field = make_field(data=data.copy()) + + result, = node.process( + field, + xy_mode="keep", z_mode="scale", + xreal_new=1e-6, yreal_new=1e-6, xy_scale=1.0, + z_min=0.0, z_max=1e-9, z_scale=2.0, z_offset=0.0, + xy_unit="", z_unit="", + ) + + np.testing.assert_allclose(result.data, data * 2.0) + + +def test_z_set_range(): + from backend.nodes.calibration import Calibration + node = Calibration() + data = np.array([[1.0, 2.0], [3.0, 4.0]]) + field = make_field(data=data.copy()) + + result, = node.process( + field, + xy_mode="keep", z_mode="set_range", + xreal_new=1e-6, yreal_new=1e-6, xy_scale=1.0, + z_min=0.0, z_max=1.0, z_scale=1.0, z_offset=0.0, + xy_unit="", z_unit="", + ) + + assert float(result.data.min()) == 0.0 + assert float(result.data.max()) == 1.0 + # Check that intermediate values are linearly mapped + np.testing.assert_allclose(result.data, (data - 1.0) / 3.0) diff --git a/tests/node_tests/displacement_field.py b/tests/node_tests/displacement_field.py new file mode 100644 index 0000000..885f127 --- /dev/null +++ b/tests/node_tests/displacement_field.py @@ -0,0 +1,41 @@ +import numpy as np +import pytest +from tests.node_tests._shared import make_field + + +def test_output_shape(): + from backend.nodes.displacement_field import DisplacementField + + node = DisplacementField() + field = make_field(shape=(48, 64)) + result, = node.process(field, "gaussian_1d", sigma=5.0, tau=20.0, density=0.02, seed=42) + assert result.data.shape == (48, 64) + + +def test_zero_sigma_unchanged(): + from backend.nodes.displacement_field import DisplacementField + + node = DisplacementField() + field = make_field(shape=(32, 32)) + result, = node.process(field, "gaussian_1d", sigma=0.1, tau=20.0, density=0.02, seed=42) + np.testing.assert_allclose(result.data, field.data, atol=0.5) + + +def test_gaussian_2d_modifies(): + from backend.nodes.displacement_field import DisplacementField + + node = DisplacementField() + field = make_field(shape=(32, 32)) + result, = node.process(field, "gaussian_2d", sigma=10.0, tau=20.0, density=0.02, seed=42) + # With sigma=10 the output should differ from the input + assert not np.allclose(result.data, field.data, atol=1e-3) + + +def test_all_methods_finite(): + from backend.nodes.displacement_field import DisplacementField + + node = DisplacementField() + field = make_field(shape=(32, 32)) + for method in ("gaussian_1d", "gaussian_2d", "tear"): + result, = node.process(field, method, sigma=5.0, tau=20.0, density=0.02, seed=42) + assert np.all(np.isfinite(result.data)), f"{method} produced non-finite values" diff --git a/tests/node_tests/distribution_coercion.py b/tests/node_tests/distribution_coercion.py new file mode 100644 index 0000000..c7e5cf2 --- /dev/null +++ b/tests/node_tests/distribution_coercion.py @@ -0,0 +1,61 @@ +import numpy as np +from tests.node_tests._shared import make_field + + +def test_output_shape(): + from backend.nodes.distribution_coercion import DistributionCoercion + + node = DistributionCoercion() + field = make_field(shape=(48, 64)) + + for dist in ("uniform", "gaussian", "levels"): + (result,) = node.process(field, distribution=dist, n_levels=4, processing="field") + assert result.data.shape == field.data.shape + + +def test_uniform_distribution(): + from backend.nodes.distribution_coercion import DistributionCoercion + + node = DistributionCoercion() + rng = np.random.default_rng(7) + data = rng.exponential(scale=2.0, size=(64, 64)) + field = make_field(data=data) + + (result,) = node.process(field, distribution="uniform", n_levels=4, processing="field") + + assert np.isclose(result.data.min(), data.min()) + assert np.isclose(result.data.max(), data.max()) + + # Histogram should be roughly uniform — check that no bin has more than + # 2x the expected count. + counts, _ = np.histogram(result.data.ravel(), bins=10) + expected = result.data.size / 10 + assert all(c < 2.0 * expected for c in counts) + + +def test_levels_count(): + from backend.nodes.distribution_coercion import DistributionCoercion + + node = DistributionCoercion() + field = make_field(shape=(64, 64)) + + for n in (2, 5, 10): + (result,) = node.process(field, distribution="levels", n_levels=n, processing="field") + assert len(np.unique(result.data)) == n + + +def test_row_mode(): + from backend.nodes.distribution_coercion import DistributionCoercion + + node = DistributionCoercion() + field = make_field(shape=(32, 48)) + + (result,) = node.process(field, distribution="uniform", n_levels=4, processing="rows") + assert result.data.shape == field.data.shape + + # Each row should span the row's own min/max + for i in range(field.data.shape[0]): + row_in = field.data[i] + row_out = result.data[i] + assert np.isclose(row_out.min(), row_in.min()) + assert np.isclose(row_out.max(), row_in.max()) diff --git a/tests/node_tests/dwt_anisotropy.py b/tests/node_tests/dwt_anisotropy.py new file mode 100644 index 0000000..18248ae --- /dev/null +++ b/tests/node_tests/dwt_anisotropy.py @@ -0,0 +1,64 @@ +import numpy as np +import pytest +from tests.node_tests._shared import make_field + + +def test_output_shape(): + """Anisotropy map shape must match the input field.""" + from backend.nodes.dwt_anisotropy import DWTAnisotropy + + node = DWTAnisotropy() + field = make_field(shape=(64, 64)) + aniso_field, stats = node.process(field, n_levels=4, ratio_threshold=0.2) + assert aniso_field.data.shape == (64, 64) + + +def test_isotropic_surface(): + """A random isotropic surface should have X/Y energy ratios near 1.0.""" + from backend.nodes.dwt_anisotropy import DWTAnisotropy + + rng = np.random.default_rng(42) + # Use a larger field so deeper levels still have enough coefficients + data = rng.standard_normal((128, 128)) + field = make_field(data=data) + node = DWTAnisotropy() + aniso_field, stats = node.process(field, n_levels=3, ratio_threshold=0.2) + + for row in stats: + assert 0.5 < row["ratio"] < 2.0, ( + f"Level {row['level']} ratio {row['ratio']:.3f} too far from 1.0 for isotropic surface" + ) + + +def test_statistics_table(): + """Statistics output is a list of dicts with the expected keys.""" + from backend.nodes.dwt_anisotropy import DWTAnisotropy + + node = DWTAnisotropy() + field = make_field(shape=(64, 64)) + aniso_field, stats = node.process(field, n_levels=3, ratio_threshold=0.2) + + assert isinstance(stats, list) + assert len(stats) == 3 + expected_keys = {"level", "x_energy", "y_energy", "ratio", "anisotropic"} + for row in stats: + assert isinstance(row, dict) + assert set(row.keys()) == expected_keys + + +def test_anisotropic_detection(): + """Horizontal stripes should produce a ratio clearly different from 1.0.""" + from backend.nodes.dwt_anisotropy import DWTAnisotropy + + # Create horizontal stripes: constant along columns, varying along rows + data = np.tile(np.sin(np.linspace(0, 10 * np.pi, 64)), (64, 1)) + field = make_field(data=data) + node = DWTAnisotropy() + aniso_field, stats = node.process(field, n_levels=4, ratio_threshold=0.2) + + # At least one level should show a ratio far from 1.0 + has_anisotropic = any(abs(row["ratio"] - 1.0) > 0.2 for row in stats) + assert has_anisotropic, ( + f"Expected anisotropic detection for horizontal stripes, ratios: " + f"{[row['ratio'] for row in stats]}" + ) diff --git a/tests/node_tests/grain_visualization.py b/tests/node_tests/grain_visualization.py new file mode 100644 index 0000000..d049454 --- /dev/null +++ b/tests/node_tests/grain_visualization.py @@ -0,0 +1,64 @@ +import numpy as np +from tests.node_tests._shared import make_field + + +def _make_test_inputs(): + """Create a 64x64 field and mask with two isolated blobs.""" + data = np.zeros((64, 64), dtype=np.float64) + data[10:20, 10:20] = 5.0 + data[40:55, 40:55] = 3.0 + field = make_field(data=data, xreal=1e-6, yreal=1e-6) + + mask = np.zeros((64, 64), dtype=np.uint8) + mask[10:20, 10:20] = 255 + mask[40:55, 40:55] = 255 + return field, mask + + +def test_output_shape(): + from backend.nodes.grain_visualization import GrainVisualization + + node = GrainVisualization() + field, mask = _make_test_inputs() + result, labeled = node.process(field, mask, style="inscribed_disc", fill=False) + assert result.shape == mask.shape, ( + f"Result shape {result.shape} does not match input shape {mask.shape}" + ) + + +def test_labeled_grains(): + from backend.nodes.grain_visualization import GrainVisualization + + node = GrainVisualization() + field, mask = _make_test_inputs() + result, labeled = node.process(field, mask, style="inscribed_disc", fill=False) + unique_ids = set(np.unique(labeled.data)) - {0.0} + assert len(unique_ids) == 2, ( + f"Expected 2 unique nonzero grain labels, got {len(unique_ids)}: {unique_ids}" + ) + + +def test_disc_style(): + from backend.nodes.grain_visualization import GrainVisualization + + node = GrainVisualization() + field, mask = _make_test_inputs() + + result_outline, _ = node.process(field, mask, style="inscribed_disc", fill=False) + assert np.any(result_outline > 0), "inscribed_disc outline produced an empty mask" + + result_filled, _ = node.process(field, mask, style="inscribed_disc", fill=True) + assert np.any(result_filled > 0), "inscribed_disc filled produced an empty mask" + + +def test_bounding_box_style(): + from backend.nodes.grain_visualization import GrainVisualization + + node = GrainVisualization() + field, mask = _make_test_inputs() + + result_outline, _ = node.process(field, mask, style="bounding_box", fill=False) + assert np.any(result_outline > 0), "bounding_box outline produced an empty mask" + + result_filled, _ = node.process(field, mask, style="bounding_box", fill=True) + assert np.any(result_filled > 0), "bounding_box filled produced an empty mask" diff --git a/tests/node_tests/logistic_classification.py b/tests/node_tests/logistic_classification.py new file mode 100644 index 0000000..73f55c2 --- /dev/null +++ b/tests/node_tests/logistic_classification.py @@ -0,0 +1,111 @@ +import numpy as np +from backend.execution_context import active_node, execution_callbacks +from tests.node_tests._shared import make_field + + +def test_output_shapes(): + from backend.nodes.logistic_classification import LogisticClassification + node = LogisticClassification() + + data = np.random.default_rng(0).standard_normal((64, 64)) + field = make_field(data=data) + + previews = [] + with execution_callbacks(preview=lambda nid, uri: previews.append(uri)), active_node("test"): + mask, prob = node.process( + field, + use_gaussians=True, + n_gaussians=4, + use_sobel=True, + use_laplacian=True, + regularization=1.0, + max_iter=500, + seed=42, + ) + + assert mask.shape == field.data.shape + assert prob.data.shape == field.data.shape + + +def test_mask_binary(): + from backend.nodes.logistic_classification import LogisticClassification + node = LogisticClassification() + + data = np.zeros((32, 32)) + data[:, 16:] = 1.0 + field = make_field(data=data) + + with execution_callbacks(preview=lambda nid, uri: None), active_node("test"): + mask, _ = node.process( + field, + use_gaussians=True, + n_gaussians=2, + use_sobel=True, + use_laplacian=True, + regularization=1.0, + max_iter=500, + seed=42, + ) + + unique = set(np.unique(mask)) + assert unique <= {0, 255}, f"Mask contains non-binary values: {unique}" + + +def test_probability_range(): + from backend.nodes.logistic_classification import LogisticClassification + node = LogisticClassification() + + data = np.random.default_rng(7).standard_normal((48, 48)) + field = make_field(data=data) + + with execution_callbacks(preview=lambda nid, uri: None), active_node("test"): + _, prob = node.process( + field, + use_gaussians=True, + n_gaussians=3, + use_sobel=True, + use_laplacian=True, + regularization=1.0, + max_iter=500, + seed=42, + ) + + assert prob.data.min() >= 0.0, f"Probability min {prob.data.min()} < 0" + assert prob.data.max() <= 1.0, f"Probability max {prob.data.max()} > 1" + + +def test_with_training(): + from backend.nodes.logistic_classification import LogisticClassification + node = LogisticClassification() + + # Create a field with two distinct regions + data = np.zeros((64, 64)) + data[:, 32:] = 2.0 + data += np.random.default_rng(1).standard_normal((64, 64)) * 0.1 + field = make_field(data=data) + + # Create a training mask marking the right half as positive + training_mask = np.zeros((64, 64), dtype=np.uint8) + training_mask[:, 32:] = 255 + + with execution_callbacks(preview=lambda nid, uri: None), active_node("test"): + mask, prob = node.process( + field, + use_gaussians=True, + n_gaussians=3, + use_sobel=True, + use_laplacian=True, + regularization=1.0, + max_iter=500, + seed=42, + training_mask=training_mask, + ) + + assert mask.dtype == np.uint8 + assert mask.shape == field.data.shape + # The classifier should learn that the right half is positive + right_positive = np.count_nonzero(mask[:, 32:] == 255) + left_positive = np.count_nonzero(mask[:, :32] == 255) + assert right_positive > left_positive, ( + f"Expected more positives on right ({right_positive}) than left ({left_positive})" + ) diff --git a/tests/node_tests/mark_disconnected.py b/tests/node_tests/mark_disconnected.py new file mode 100644 index 0000000..8eb1635 --- /dev/null +++ b/tests/node_tests/mark_disconnected.py @@ -0,0 +1,43 @@ +import numpy as np +import pytest +from tests.node_tests._shared import make_field + + +def test_output_shape(): + from backend.nodes.mark_disconnected import MarkDisconnected + + node = MarkDisconnected() + field = make_field(shape=(64, 64)) + mask, = node.process(field, defect_type="both", radius=5, threshold=0.1) + assert mask.shape == (64, 64) + + +def test_flat_surface_no_defects(): + from backend.nodes.mark_disconnected import MarkDisconnected + + node = MarkDisconnected() + data = np.ones((64, 64)) * 5.0 + field = make_field(data=data) + mask, = node.process(field, defect_type="both", radius=5, threshold=0.1) + assert np.count_nonzero(mask) == 0 + + +def test_spike_detected(): + from backend.nodes.mark_disconnected import MarkDisconnected + + node = MarkDisconnected() + data = np.ones((64, 64), dtype=np.float64) + mean_val = data.mean() + data[32, 32] = mean_val * 100 # large spike + field = make_field(data=data) + mask, = node.process(field, defect_type="positive", radius=3, threshold=0.05) + assert mask[32, 32] == 255 + + +def test_output_is_uint8(): + from backend.nodes.mark_disconnected import MarkDisconnected + + node = MarkDisconnected() + field = make_field(shape=(32, 32)) + mask, = node.process(field, defect_type="negative", radius=5, threshold=0.1) + assert mask.dtype == np.uint8 diff --git a/tests/node_tests/mask_noisify.py b/tests/node_tests/mask_noisify.py new file mode 100644 index 0000000..dacb63c --- /dev/null +++ b/tests/node_tests/mask_noisify.py @@ -0,0 +1,50 @@ +import numpy as np + + +def _make_test_mask(): + mask = np.zeros((64, 64), dtype=np.uint8) + mask[20:40, 20:40] = 255 + return mask + + +def test_output_shape(): + from backend.nodes.mask_noisify import MaskNoisify + node = MaskNoisify() + mask = _make_test_mask() + + result, = node.process(mask, density=0.1, direction="both", + boundaries_only=True, seed=42) + assert result.shape == mask.shape + assert result.dtype == np.uint8 + + +def test_zero_density_unchanged(): + from backend.nodes.mask_noisify import MaskNoisify + node = MaskNoisify() + mask = _make_test_mask() + + result, = node.process(mask, density=0.0, direction="both", + boundaries_only=True, seed=42) + assert np.array_equal(result, mask) + + +def test_density_modifies_mask(): + from backend.nodes.mask_noisify import MaskNoisify + node = MaskNoisify() + mask = _make_test_mask() + + result, = node.process(mask, density=0.5, direction="both", + boundaries_only=True, seed=42) + assert not np.array_equal(result, mask) + + +def test_seed_reproducibility(): + from backend.nodes.mask_noisify import MaskNoisify + node = MaskNoisify() + mask = _make_test_mask() + + result_a, = node.process(mask, density=0.3, direction="both", + boundaries_only=True, seed=123) + result_b, = node.process(mask, density=0.3, direction="both", + boundaries_only=True, seed=123) + assert np.array_equal(result_a, result_b) diff --git a/tests/node_tests/mask_shift.py b/tests/node_tests/mask_shift.py new file mode 100644 index 0000000..fe21996 --- /dev/null +++ b/tests/node_tests/mask_shift.py @@ -0,0 +1,74 @@ +import numpy as np +import pytest + + +def _make_mask(): + """Create a simple test mask: 10x10 block of 255 in a 64x64 field.""" + mask = np.zeros((64, 64), dtype=np.uint8) + mask[10:20, 10:20] = 255 + return mask + + +def test_output_shape(): + from backend.nodes.mask_shift import MaskShift + node = MaskShift() + mask = _make_mask() + + result, = node.process(mask, shift_x=5, shift_y=3, border_mode="zero") + assert result.shape == mask.shape + assert result.dtype == np.uint8 + + result_wrap, = node.process(mask, shift_x=-10, shift_y=7, border_mode="wrap") + assert result_wrap.shape == mask.shape + + result_mirror, = node.process(mask, shift_x=2, shift_y=-4, border_mode="mirror") + assert result_mirror.shape == mask.shape + + +def test_zero_shift_unchanged(): + from backend.nodes.mask_shift import MaskShift + node = MaskShift() + mask = _make_mask() + + result_zero, = node.process(mask, shift_x=0, shift_y=0, border_mode="zero") + assert np.array_equal(result_zero, mask) + + result_wrap, = node.process(mask, shift_x=0, shift_y=0, border_mode="wrap") + assert np.array_equal(result_wrap, mask) + + result_mirror, = node.process(mask, shift_x=0, shift_y=0, border_mode="mirror") + assert np.array_equal(result_mirror, mask) + + +def test_wrap_mode(): + from backend.nodes.mask_shift import MaskShift + node = MaskShift() + mask = _make_mask() + + # Shift block right by 60 pixels — the block at cols 10:20 should wrap + # and appear at cols 70%64=6 to 80%64=16, spanning the boundary. + result, = node.process(mask, shift_x=60, shift_y=0, border_mode="wrap") + assert result.dtype == np.uint8 + # The total number of masked pixels should be preserved in wrap mode + assert np.count_nonzero(result) == np.count_nonzero(mask) + # Original location should not all still be set + # (shift is large enough to move block away from original position) + assert not np.array_equal(result, mask) + + +def test_zero_mode_fills(): + from backend.nodes.mask_shift import MaskShift + node = MaskShift() + mask = _make_mask() + + # Shift right by 5 — left 5 columns should be zeroed + result, = node.process(mask, shift_x=5, shift_y=0, border_mode="zero") + assert np.all(result[:, :5] == 0) + # Block should now be at cols 15:25, rows 10:20 + assert np.all(result[10:20, 15:25] == 255) + + # Shift down by 5 — top 5 rows should be zeroed + result2, = node.process(mask, shift_x=0, shift_y=5, border_mode="zero") + assert np.all(result2[:5, :] == 0) + # Block should now be at rows 15:25, cols 10:20 + assert np.all(result2[15:25, 10:20] == 255) diff --git a/tests/node_tests/neural_classification.py b/tests/node_tests/neural_classification.py new file mode 100644 index 0000000..9734f05 --- /dev/null +++ b/tests/node_tests/neural_classification.py @@ -0,0 +1,72 @@ +import numpy as np +from tests.node_tests._shared import make_field + + +def test_output_shapes(): + from backend.nodes.neural_classification import NeuralClassification + + node = NeuralClassification() + data = np.random.default_rng(0).standard_normal((32, 32)) + field = make_field(data=data) + + mask, prob_field = node.process(field, n_gaussians=3, n_hidden=8, + train_steps=20, seed=7) + assert mask.shape == (32, 32) + assert prob_field.data.shape == (32, 32) + + +def test_mask_is_binary(): + from backend.nodes.neural_classification import NeuralClassification + + node = NeuralClassification() + data = np.random.default_rng(1).standard_normal((24, 24)) + field = make_field(data=data) + + mask, _ = node.process(field, n_gaussians=2, n_hidden=8, + train_steps=10, seed=0) + unique = set(np.unique(mask).tolist()) + assert unique <= {0, 255}, f"Unexpected mask values: {unique}" + + +def test_probability_range(): + from backend.nodes.neural_classification import NeuralClassification + + node = NeuralClassification() + data = np.random.default_rng(2).standard_normal((32, 32)) + field = make_field(data=data) + + _, prob_field = node.process(field, n_gaussians=4, n_hidden=16, + train_steps=50, seed=42) + assert prob_field.data.min() >= 0.0 + assert prob_field.data.max() <= 1.0 + + +def test_with_training_mask(): + from backend.nodes.neural_classification import NeuralClassification + + node = NeuralClassification() + + # Create a field with two distinct height regions + data = np.zeros((48, 48), dtype=np.float64) + data[:, 24:] = 5.0 # right half is elevated + field = make_field(data=data) + + # Training mask: left half = 0 (class A), right half = 255 (class B) + training_mask = np.zeros((48, 48), dtype=np.uint8) + training_mask[:, 24:] = 255 + + mask, prob_field = node.process(field, n_gaussians=4, n_hidden=16, + train_steps=200, seed=42, + training_mask=training_mask) + + assert mask.dtype == np.uint8 + assert mask.shape == (48, 48) + assert prob_field.data.shape == (48, 48) + + # The network should learn to classify the two regions correctly. + # Check that most of the right half is class B and left half is class A. + right_classified = np.count_nonzero(mask[:, 24:] == 255) + left_classified = np.count_nonzero(mask[:, :24] == 0) + total_half = 48 * 24 + assert right_classified > total_half * 0.8, "Right half should mostly be class B" + assert left_classified > total_half * 0.8, "Left half should mostly be class A" diff --git a/tests/node_tests/pixel_classification.py b/tests/node_tests/pixel_classification.py new file mode 100644 index 0000000..2cc9f29 --- /dev/null +++ b/tests/node_tests/pixel_classification.py @@ -0,0 +1,49 @@ +import numpy as np +from tests.node_tests._shared import make_field + + +def test_output_shape(): + from backend.nodes.pixel_classification import PixelClassification + + node = PixelClassification() + field = make_field(shape=(64, 64)) + classified, mask = node.process(field, n_classes=3, feature="height", method="quantile") + assert classified.data.shape == field.data.shape + + +def test_correct_number_of_classes(): + from backend.nodes.pixel_classification import PixelClassification + + node = PixelClassification() + field = make_field(shape=(64, 64)) + for n in (2, 4, 5): + classified, _ = node.process(field, n_classes=n, feature="height", method="quantile") + unique = np.unique(classified.data) + assert len(unique) <= n, f"Expected at most {n} classes, got {len(unique)}" + + +def test_equal_range_method(): + from backend.nodes.pixel_classification import PixelClassification + + node = PixelClassification() + # Linear ramp: equal_range should produce evenly distributed labels + ramp = np.linspace(0, 1, 64 * 64).reshape(64, 64) + field = make_field(data=ramp) + classified, _ = node.process(field, n_classes=4, feature="height", method="equal_range") + labels = classified.data.astype(int) + unique = np.unique(labels) + assert len(unique) == 4 + # Each class should have roughly 25% of pixels + counts = [np.sum(labels == u) for u in unique] + for c in counts: + assert abs(c - 64 * 64 / 4) < 64 * 64 * 0.05 # within 5% + + +def test_mask_output(): + from backend.nodes.pixel_classification import PixelClassification + + node = PixelClassification() + field = make_field(shape=(32, 32)) + _, mask = node.process(field, n_classes=3, feature="height", method="otsu") + assert mask.dtype == np.uint8 + assert set(np.unique(mask)).issubset({0, 255}) diff --git a/tests/node_tests/presentation_ops.py b/tests/node_tests/presentation_ops.py new file mode 100644 index 0000000..0f331d5 --- /dev/null +++ b/tests/node_tests/presentation_ops.py @@ -0,0 +1,43 @@ +import numpy as np +from tests.node_tests._shared import make_field + + +def test_output_shape(): + from backend.nodes.presentation_ops import PresentationOps + + node = PresentationOps() + field = make_field(shape=(32, 32)) + result, = node.process(field, "logscale", blend_factor=0.5) + assert result.data.shape == field.data.shape + + +def test_logscale(): + from backend.nodes.presentation_ops import PresentationOps + + node = PresentationOps() + data = np.array([[1.0, 10.0], [100.0, 1000.0]]) + field = make_field(data=data) + result, = node.process(field, "logscale", blend_factor=0.5) + assert np.all(np.isfinite(result.data)) + # logscale should preserve ordering + assert result.data[0, 0] < result.data[0, 1] < result.data[1, 0] < result.data[1, 1] + + +def test_blend_at_zero(): + from backend.nodes.presentation_ops import PresentationOps + + node = PresentationOps() + field = make_field(data=np.array([[1.0, 2.0], [3.0, 4.0]])) + overlay = make_field(data=np.array([[10.0, 20.0], [30.0, 40.0]])) + result, = node.process(field, "blend", blend_factor=0.0, overlay=overlay) + assert np.allclose(result.data, field.data) + + +def test_blend_at_one(): + from backend.nodes.presentation_ops import PresentationOps + + node = PresentationOps() + field = make_field(data=np.array([[1.0, 2.0], [3.0, 4.0]])) + overlay = make_field(data=np.array([[10.0, 20.0], [30.0, 40.0]])) + result, = node.process(field, "blend", blend_factor=1.0, overlay=overlay) + assert np.allclose(result.data, overlay.data) diff --git a/tests/node_tests/psf_estimation.py b/tests/node_tests/psf_estimation.py new file mode 100644 index 0000000..4a10068 --- /dev/null +++ b/tests/node_tests/psf_estimation.py @@ -0,0 +1,56 @@ +import numpy as np +from scipy.ndimage import gaussian_filter +from tests.node_tests._shared import make_field + + +def _make_test_pair(shape=(64, 64), sigma=2.0): + """Return (measured, ideal) where measured = gaussian_filter(ideal).""" + ideal = make_field(shape=shape) + measured_data = gaussian_filter(ideal.data, sigma=sigma) + measured = make_field(data=measured_data) + return measured, ideal + + +def test_output_shape(): + from backend.nodes.psf_estimation import PSFEstimation + + node = PSFEstimation() + measured, ideal = _make_test_pair() + psf_field, _ = node.process(measured, ideal, "wiener", 0.01, 32) + assert psf_field.data.shape == (32, 32) + + +def test_psf_normalized(): + from backend.nodes.psf_estimation import PSFEstimation + + node = PSFEstimation() + measured, ideal = _make_test_pair() + for method in ("wiener", "least_squares", "gaussian_fit"): + psf_field, _ = node.process(measured, ideal, method, 0.01, 32) + assert abs(psf_field.data.sum() - 1.0) < 1e-6, ( + f"{method}: PSF sum = {psf_field.data.sum()}" + ) + + +def test_gaussian_fit_parameters(): + from backend.nodes.psf_estimation import PSFEstimation + + node = PSFEstimation() + measured, ideal = _make_test_pair() + _, parameters = node.process(measured, ideal, "gaussian_fit", 0.01, 32) + names = {row["quantity"] for row in parameters} + assert "sigma_x" in names + assert "sigma_y" in names + assert "amplitude" in names + + +def test_all_methods_finite(): + from backend.nodes.psf_estimation import PSFEstimation + + node = PSFEstimation() + measured, ideal = _make_test_pair() + for method in ("wiener", "least_squares", "gaussian_fit"): + psf_field, _ = node.process(measured, ideal, method, 0.01, 32) + assert np.isfinite(psf_field.data).all(), ( + f"{method}: PSF contains non-finite values" + ) diff --git a/tests/node_tests/super_resolution.py b/tests/node_tests/super_resolution.py new file mode 100644 index 0000000..9bfd603 --- /dev/null +++ b/tests/node_tests/super_resolution.py @@ -0,0 +1,51 @@ +import numpy as np +from tests.node_tests._shared import make_field + + +def test_output_shape_single(): + """Single input with upscale=2 gives 2x output size.""" + from backend.nodes.super_resolution import SuperResolution + + field = make_field(shape=(32, 32)) + node = SuperResolution() + result, = node.process(field, upscale=2) + assert result.data.shape == (64, 64) + + +def test_output_shape_multi(): + """Multiple inputs still give 2x output size.""" + from backend.nodes.super_resolution import SuperResolution + + rng = np.random.default_rng(0) + f1 = make_field(data=rng.standard_normal((32, 32))) + f2 = make_field(data=rng.standard_normal((32, 32))) + f3 = make_field(data=rng.standard_normal((32, 32))) + node = SuperResolution() + result, = node.process(f1, upscale=2, field2=f2, field3=f3) + assert result.data.shape == (64, 64) + + +def test_finite_values(): + """Output values must all be finite.""" + from backend.nodes.super_resolution import SuperResolution + + rng = np.random.default_rng(1) + f1 = make_field(data=rng.standard_normal((32, 32))) + f2 = make_field(data=rng.standard_normal((32, 32))) + node = SuperResolution() + result, = node.process(f1, upscale=2, field2=f2) + assert np.all(np.isfinite(result.data)) + + +def test_upscale_factor(): + """Output dimensions should equal input dimensions times upscale factor.""" + from backend.nodes.super_resolution import SuperResolution + + field = make_field(shape=(32, 32)) + node = SuperResolution() + for factor in (2, 3, 4): + result, = node.process(field, upscale=factor) + expected = (32 * factor, 32 * factor) + assert result.data.shape == expected, ( + f"upscale={factor}: expected {expected}, got {result.data.shape}" + ) diff --git a/tests/node_tests/tip_shape_estimate.py b/tests/node_tests/tip_shape_estimate.py new file mode 100644 index 0000000..b2542ce --- /dev/null +++ b/tests/node_tests/tip_shape_estimate.py @@ -0,0 +1,86 @@ +import numpy as np +import pytest +from tests.node_tests._shared import make_field + + +def run_tip_shape(field, feature_type="edge", feature_radius=100e-9, n_points=100): + from backend.nodes.tip_shape_estimate import TipShapeEstimate + node = TipShapeEstimate() + tip_shape, parameters = node.process( + field=field, + feature_type=feature_type, + feature_radius=feature_radius, + n_points=n_points, + ) + return tip_shape, parameters + + +# ── Output shape and type ──────────────────────────────────────────────────── + +def test_output_shape(): + """tip_shape must be a 2D DataField.""" + from backend.data_types import DataField + field = make_field(shape=(64, 64), xreal=64e-9, yreal=64e-9) + tip_shape, _ = run_tip_shape(field, feature_type="edge", n_points=33) + assert isinstance(tip_shape, DataField) + assert tip_shape.data.ndim == 2 + assert tip_shape.data.shape[0] == tip_shape.data.shape[1] + + +# ── Parameters table ───────────────────────────────────────────────────────── + +def test_parameters_table(): + """parameters must be a RecordTable containing a tip_radius entry.""" + from backend.data_types import RecordTable + field = make_field(shape=(64, 64), xreal=64e-9, yreal=64e-9) + _, parameters = run_tip_shape(field, feature_type="edge", n_points=33) + assert isinstance(parameters, RecordTable) + quantities = {row["quantity"] for row in parameters} + assert "tip_radius" in quantities + + +# ── Tip apex is maximum ────────────────────────────────────────────────────── + +def test_tip_apex_is_maximum(): + """The centre of the tip shape should be the highest point.""" + field = make_field(shape=(64, 64), xreal=64e-9, yreal=64e-9) + tip_shape, _ = run_tip_shape(field, feature_type="edge", n_points=33) + n = tip_shape.data.shape[0] + ci = n // 2 + assert tip_shape.data[ci, ci] == pytest.approx(tip_shape.data.max(), abs=1e-20) + + +# ── Sphere feature produces valid estimate ─────────────────────────────────── + +def test_sphere_feature(): + """Using a synthetic sphere as input produces a valid tip estimate.""" + from backend.data_types import DataField + + # Create a synthetic sphere on a 64x64 grid. + R = 100e-9 # sphere radius + n = 64 + pixel_size = 10e-9 + xreal = n * pixel_size + Y, X = np.mgrid[:n, :n] + r = np.sqrt((X - 32) ** 2 + (Y - 32) ** 2) * pixel_size + data = np.sqrt(np.maximum(R ** 2 - r ** 2, 0.0)) + + field = make_field(data=data, xreal=xreal, yreal=xreal) + + tip_shape, parameters = run_tip_shape( + field, feature_type="sphere", feature_radius=R, n_points=33, + ) + + # The output must be a valid 2D DataField. + assert isinstance(tip_shape, DataField) + assert tip_shape.data.ndim == 2 + + # Apex should be the maximum. + ci = tip_shape.data.shape[0] // 2 + assert tip_shape.data[ci, ci] == pytest.approx(tip_shape.data.max(), abs=1e-20) + + # Parameters should contain tip_radius. + quantities = {row["quantity"]: row["value"] for row in parameters} + assert "tip_radius" in quantities + # The estimated radius must be a positive finite number. + assert quantities["tip_radius"] > 0