diff --git a/GWYDDION_FEATURE_GAP.md b/GWYDDION_FEATURE_GAP.md index e708040..5b266f2 100644 --- a/GWYDDION_FEATURE_GAP.md +++ b/GWYDDION_FEATURE_GAP.md @@ -28,15 +28,15 @@ Reference for future implementation. Grouped by value to typical SPM workflows. | # | Feature | Gwyddion Source | Description | |---|---------|---------------|-------------| | 15 | Correlation / Pattern Matching | crosscor.c, maskcor.c | Find repeated features or align images via cross-correlation. | -| 16 | Slope Distribution | slope_dist.c | Angular histogram of surface slopes. Characterizes surface texture directionality. | -| 17 | Grain Filtering | grain_filter.c | Remove grains by size, height, or border contact. Refine grain masks post-detection. | -| 18 | Field Arithmetic | arithmetic.c | Add/subtract/multiply/divide two DATA_FIELDs. Useful for difference maps, normalization. | +| ~~16~~ | ~~Slope Distribution~~ | ~~slope_dist.c~~ | ~~Angular histogram of surface slopes. Characterizes surface texture directionality.~~ **DONE** | +| ~~17~~ | ~~Grain Filtering~~ | ~~grain_filter.c~~ | ~~Remove grains by size, height, or border contact. Refine grain masks post-detection.~~ **DONE** | +| ~~18~~ | ~~Field Arithmetic~~ | ~~arithmetic.c~~ | ~~Add/subtract/multiply/divide two DATA_FIELDs. Useful for difference maps, normalization.~~ **DONE** | | 19 | Spot Removal | spotremove.c | Interpolate over selected point defects (dust, spikes). | -| 20 | Tip Modeling / Deconvolution | tip_blind.c, tip_model.c | Estimate tip shape from image, deconvolve to recover true surface. | -| 21 | Radial Profile | rprofile tool | Azimuthally averaged profile from a center point. Good for circular features. | +| ~~20~~ | ~~Tip Modeling / Deconvolution~~ | ~~tip_blind.c, tip_model.c~~ | ~~Estimate tip shape from image, deconvolve to recover true surface.~~ **DONE** | +| ~~21~~ | ~~Radial Profile~~ | ~~rprofile tool~~ | ~~Azimuthally averaged profile from a center point. Good for circular features.~~ **DONE** | | 22 | Wavelet Transform | dwt.c, cwt.c | Discrete/continuous wavelet analysis. Multi-scale roughness decomposition. | -| 23 | Scale / Resample | scale.c, resample.c | Resize fields with interpolation. | -| 24 | Gradient | gradient.c | Compute x/y gradient magnitude maps. | +| ~~23~~ | ~~Scale / Resample~~ | ~~scale.c, resample.c~~ | ~~Resize fields with interpolation.~~ **DONE** | +| ~~24~~ | ~~Gradient~~ | ~~gradient.c~~ | ~~Compute x/y gradient magnitude maps.~~ **DONE** | | 25 | Custom Convolution | convolution_filter.c | User-defined kernel convolution. | | 26 | Local Contrast Enhancement | local_contrast.c | Enhance visibility of local features in images. | @@ -99,3 +99,6 @@ For reference, these Gwyddion equivalents are already covered: | Watershed Segmentation | grains | grain_wshed.c | | Grain Analysis | grains | grain_stat.c | | Preview / 3D View / Print Table | display | Presentation, 3D view | +| Tip Model | tip | tip_model.c, tip.c | +| Tip Deconvolution | tip | tip_blind.c, tip.c (gwy_tip_erosion) | +| Blind Tip Estimate | tip | tip_blind.c, morph_lib.c (gwy_tip_estimate_partial/full + gwy_tip_cmap) | diff --git a/backend/node_menu.py b/backend/node_menu.py index d388df7..74940f5 100644 --- a/backend/node_menu.py +++ b/backend/node_menu.py @@ -41,11 +41,14 @@ MENU_LAYOUT: dict[str, list[str]] = { "CropResizeField", "RotateField", "FlipField", + "Resample", + "FieldArithmetic", ], "Filter": [ "GaussianFilter", "MedianFilter", "EdgeDetect", + "Gradient", "FFTFilter", "ScarRemoval", ], @@ -76,6 +79,8 @@ MENU_LAYOUT: dict[str, list[str]] = { "ACF2D", "ACF1D", "PSDF", + "SlopeDistribution", + "RadialProfile", "Statistics", "Stats", ], @@ -92,6 +97,12 @@ MENU_LAYOUT: dict[str, list[str]] = { "GrainDistanceTransform", "WatershedSegmentation", "GrainAnalysis", + "GrainFilter", + ], + "Tip": [ + "TipModel", + "TipDeconvolution", + "BlindTipEstimate", ], } diff --git a/backend/nodes/__init__.py b/backend/nodes/__init__.py index 4660955..f58dc0a 100644 --- a/backend/nodes/__init__.py +++ b/backend/nodes/__init__.py @@ -58,4 +58,14 @@ from backend.nodes import ( watershed_segmentation, grain_analysis, fft_1d, + # Medium-value nodes (Gwyddion feature-gap) + field_arithmetic, + gradient, + resample, + slope_distribution, + grain_filter, + radial_profile, + tip_model, + tip_deconvolution, + tip_blind_estimate, ) diff --git a/backend/nodes/field_arithmetic.py b/backend/nodes/field_arithmetic.py new file mode 100644 index 0000000..5c347d5 --- /dev/null +++ b/backend/nodes/field_arithmetic.py @@ -0,0 +1,61 @@ +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="Field Arithmetic") +class FieldArithmetic: + @classmethod + def INPUT_TYPES(cls): + return { + "required": { + "field_a": ("DATA_FIELD",), + "field_b": ("DATA_FIELD",), + "operation": (["add", "subtract", "multiply", "divide", "min", "max", "hypot"],), + } + } + + OUTPUTS = ( + ('DATA_FIELD', 'result'), + ) + FUNCTION = "process" + + DESCRIPTION = ( + "Apply a point-wise arithmetic operation to two DATA_FIELDs of the same resolution. " + "add/subtract/multiply/divide/min/max perform element-wise operations; " + "hypot computes sqrt(a² + b²) per pixel. " + "Equivalent to gwy_data_field_sum_fields / subtract_fields / multiply_fields / " + "divide_fields / min_of_fields / max_of_fields / hypot_of_fields in arithmetic.c." + ) + + def process(self, field_a: DataField, field_b: DataField, operation: str) -> tuple: + if field_a.data.shape != field_b.data.shape: + raise ValueError( + f"Fields must have the same resolution: " + f"{field_a.data.shape} vs {field_b.data.shape}" + ) + + a = field_a.data + b = field_b.data + + if operation == "add": + result = a + b + elif operation == "subtract": + result = a - b + elif operation == "multiply": + result = a * b + elif operation == "divide": + result = a / b + elif operation == "min": + result = np.minimum(a, b) + elif operation == "max": + result = np.maximum(a, b) + elif operation == "hypot": + result = np.hypot(a, b) + else: + raise ValueError(f"Unknown operation: {operation!r}") + + return (field_a.replace(data=result),) diff --git a/backend/nodes/gradient.py b/backend/nodes/gradient.py new file mode 100644 index 0000000..0638d55 --- /dev/null +++ b/backend/nodes/gradient.py @@ -0,0 +1,63 @@ +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="Gradient") +class Gradient: + @classmethod + def INPUT_TYPES(cls): + return { + "required": { + "field": ("DATA_FIELD",), + "component": (["magnitude", "x", "y", "azimuth"],), + } + } + + OUTPUTS = ( + ('DATA_FIELD', 'gradient'), + ) + FUNCTION = "process" + + DESCRIPTION = ( + "Compute the spatial gradient using a Sobel operator. " + "'x'/'y' give the physical gradient components (z_unit/xy_unit); " + "'magnitude' gives sqrt(gx²+gy²); " + "'azimuth' gives the local slope direction in radians via atan2(gy, gx). " + "Equivalent to gwy_data_field_filter_sobel in Gwyddion (gradient.c)." + ) + + def process(self, field: DataField, component: str) -> tuple: + from scipy.ndimage import sobel + + data = field.data + # Sobel kernel sums to ±8 over 2-pixel span; divide by 8·dx to get z/xy slope. + gx = sobel(data, axis=1) / (8.0 * field.dx) + gy = sobel(data, axis=0) / (8.0 * field.dy) + + if component == "magnitude": + result = np.hypot(gx, gy) + z = str(field.si_unit_z or "").strip() + xy = str(field.si_unit_xy or "").strip() + out_unit_z = f"{z}/{xy}" if z and xy else (z or xy) + elif component == "x": + result = gx + z = str(field.si_unit_z or "").strip() + xy = str(field.si_unit_xy or "").strip() + out_unit_z = f"{z}/{xy}" if z and xy else (z or xy) + elif component == "y": + result = gy + z = str(field.si_unit_z or "").strip() + xy = str(field.si_unit_xy or "").strip() + out_unit_z = f"{z}/{xy}" if z and xy else (z or xy) + elif component == "azimuth": + # Azimuth: local slope direction, radians, matches Gwyddion's filter_azimuth + result = np.arctan2(gy, gx) + out_unit_z = "rad" + else: + raise ValueError(f"Unknown gradient component: {component!r}") + + return (field.replace(data=result, si_unit_z=out_unit_z),) diff --git a/backend/nodes/grain_filter.py b/backend/nodes/grain_filter.py new file mode 100644 index 0000000..3ed5eaf --- /dev/null +++ b/backend/nodes/grain_filter.py @@ -0,0 +1,72 @@ +from __future__ import annotations + +import numpy as np + +from backend.node_registry import register_node + + +@register_node(display_name="Grain Filter") +class GrainFilter: + @classmethod + def INPUT_TYPES(cls): + return { + "required": { + "mask": ("IMAGE",), + "min_area": ("INT", {"default": 10, "min": 0, "max": 1000000, "step": 1}), + "max_area": ("INT", {"default": 0, "min": 0, "max": 1000000, "step": 1}), + "remove_border": ("BOOLEAN", {"default": False}), + } + } + + OUTPUTS = ( + ('IMAGE', 'filtered_mask'), + ) + FUNCTION = "process" + + DESCRIPTION = ( + "Remove grains from a binary mask based on size and border contact. " + "'min_area': discard grains smaller than this many pixels (removes specks). " + "'max_area': discard grains larger than this many pixels (0 = no limit). " + "'remove_border': discard any grain that touches the image edge. " + "Equivalent to Gwyddion's grain_filter module (grain_filter.c)." + ) + + def process( + self, + mask: np.ndarray, + min_area: int, + max_area: int, + remove_border: bool, + ) -> tuple: + from scipy.ndimage import label + + binary = np.asarray(mask) > 127 + labeled, n_grains = label(binary) + + # Build per-grain keep table (index 0 = background, always False) + keep = np.zeros(n_grains + 1, dtype=bool) + + for gid in range(1, n_grains + 1): + grain = labeled == gid + area = int(grain.sum()) + + if area < min_area: + continue + if max_area > 0 and area > max_area: + continue + if remove_border and _touches_border(grain): + continue + + keep[gid] = True + + result = keep[labeled] + return (result.astype(np.uint8) * 255,) + + +def _touches_border(grain: np.ndarray) -> bool: + return ( + grain[0, :].any() + or grain[-1, :].any() + or grain[:, 0].any() + or grain[:, -1].any() + ) diff --git a/backend/nodes/radial_profile.py b/backend/nodes/radial_profile.py new file mode 100644 index 0000000..0b88eb1 --- /dev/null +++ b/backend/nodes/radial_profile.py @@ -0,0 +1,75 @@ +from __future__ import annotations + +import numpy as np + +from backend.node_registry import register_node +from backend.data_types import DataField, LineData + + +@register_node(display_name="Radial Profile") +class RadialProfile: + @classmethod + def INPUT_TYPES(cls): + return { + "required": { + "field": ("DATA_FIELD",), + "cx": ("FLOAT", {"default": 0.5, "min": 0.0, "max": 1.0, "step": 0.01}), + "cy": ("FLOAT", {"default": 0.5, "min": 0.0, "max": 1.0, "step": 0.01}), + "n_bins": ("INT", {"default": 128, "min": 4, "max": 4096, "step": 1}), + } + } + + OUTPUTS = ( + ('LINE', 'profile'), + ) + FUNCTION = "process" + + DESCRIPTION = ( + "Compute the azimuthally averaged radial profile from a centre point. " + "cx/cy give the centre as a fraction of the field width/height (0.5 = centre). " + "Output x-axis is radius in physical xy units. " + "Equivalent to gwy_data_field_angular_average used by Gwyddion's Radial Profile tool (rprofile.c)." + ) + + def process(self, field: DataField, cx: float, cy: float, n_bins: int) -> tuple: + yres, xres = field.data.shape + + # Centre in physical coordinates (matches Gwyddion: xc = cx*xreal + xoff) + xc_phys = cx * field.xreal + field.xoff + yc_phys = cy * field.yreal + field.yoff + + # Pixel-centre physical coordinates + xs = (np.arange(xres) + 0.5) * field.dx + field.xoff + ys = (np.arange(yres) + 0.5) * field.dy + field.yoff + gx, gy = np.meshgrid(xs, ys) + + r = np.hypot(gx - xc_phys, gy - yc_phys).ravel() + values = field.data.ravel() + + # Maximum radius — farthest pixel from centre + r_max = float(r.max()) + if r_max == 0.0: + r_max = max(field.dx, field.dy) + + # Bin by radius — matches Gwyddion's lineres-bin approach + bin_edges = np.linspace(0.0, r_max, n_bins + 1) + idx = np.clip( + np.floor(n_bins * r / r_max).astype(np.intp), 0, n_bins - 1 + ) + + sums = np.zeros(n_bins) + counts = np.zeros(n_bins, dtype=np.intp) + np.add.at(sums, idx, values) + np.add.at(counts, idx, 1) + + with np.errstate(invalid="ignore"): + profile = np.where(counts > 0, sums / counts, np.nan) + + centers = 0.5 * (bin_edges[:-1] + bin_edges[1:]) + + return (LineData( + data=profile, + x_axis=centers, + x_unit=field.si_unit_xy, + y_unit=field.si_unit_z, + ),) diff --git a/backend/nodes/resample.py b/backend/nodes/resample.py new file mode 100644 index 0000000..95ca32f --- /dev/null +++ b/backend/nodes/resample.py @@ -0,0 +1,52 @@ +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="Resample") +class Resample: + @classmethod + def INPUT_TYPES(cls): + return { + "required": { + "field": ("DATA_FIELD",), + "width": ("INT", {"default": 256, "min": 2, "max": 16384, "step": 1}), + "height": ("INT", {"default": 256, "min": 2, "max": 16384, "step": 1}), + "interpolation": (["linear", "cubic", "nearest"], {"default": "linear"}), + } + } + + OUTPUTS = ( + ('DATA_FIELD', 'resampled'), + ) + FUNCTION = "process" + + DESCRIPTION = ( + "Resample a DATA_FIELD to a new pixel resolution while preserving physical dimensions. " + "Physical size (xreal, yreal) is unchanged; pixel size dx/dy scales accordingly. " + "Equivalent to gwy_data_field_new_resampled with GWY_INTERPOLATION_LINEAR / " + "GWY_INTERPOLATION_CUBIC / GWY_INTERPOLATION_ROUND (scale.c)." + ) + + _ORDERS = {"nearest": 0, "linear": 1, "cubic": 3} + + def process(self, field: DataField, width: int, height: int, interpolation: str) -> tuple: + from scipy.ndimage import zoom + + if interpolation not in self._ORDERS: + raise ValueError( + f"Unknown interpolation {interpolation!r}. " + f"Choose from: {list(self._ORDERS)}" + ) + + old_yres, old_xres = field.data.shape + zy = height / old_yres + zx = width / old_xres + + result = zoom(field.data, (zy, zx), order=self._ORDERS[interpolation]) + + # Physical dimensions stay the same; xres/yres are re-derived from data.shape. + return (field.replace(data=result, xreal=field.xreal, yreal=field.yreal),) diff --git a/backend/nodes/slope_distribution.py b/backend/nodes/slope_distribution.py new file mode 100644 index 0000000..1b7a63e --- /dev/null +++ b/backend/nodes/slope_distribution.py @@ -0,0 +1,120 @@ +from __future__ import annotations + +import numpy as np + +from backend.node_registry import register_node +from backend.data_types import DataField, LineData + + +@register_node(display_name="Slope Distribution") +class SlopeDistribution: + @classmethod + def INPUT_TYPES(cls): + return { + "required": { + "field": ("DATA_FIELD",), + "distribution": (["theta", "phi", "gradient"],), + "n_bins": ("INT", {"default": 90, "min": 10, "max": 1000, "step": 1}), + } + } + + OUTPUTS = ( + ('LINE', 'distribution'), + ) + FUNCTION = "process" + + DESCRIPTION = ( + "Compute the angular slope distribution of a DATA_FIELD surface. " + "'theta' is the inclination angle (0–max°), probability density (1/deg); " + "'phi' is the azimuthal slope direction (0–360°), weighted by slope² (z/xy)²; " + "'gradient' is the gradient magnitude distribution, probability density (1/(z/xy)). " + "Equivalent to Gwyddion's slope_dist module (slope_dist.c)." + ) + + def process(self, field: DataField, distribution: str, n_bins: int) -> tuple: + from scipy.ndimage import sobel + + # Physical slopes in z_unit/xy_unit — matches Gwyddion's gwy_data_field_filter_sobel + gx = sobel(field.data, axis=1) / (8.0 * field.dx) + gy = sobel(field.data, axis=0) / (8.0 * field.dy) + + gx = gx.ravel() + gy = gy.ravel() + n = len(gx) + + z = str(field.si_unit_z or "").strip() + xy = str(field.si_unit_xy or "").strip() + slope_unit = f"{z}/{xy}" if z and xy else (z or xy) + + if distribution == "phi": + return self._phi(gx, gy, n_bins, slope_unit) + elif distribution == "theta": + return self._theta(gx, gy, n_bins) + elif distribution == "gradient": + return self._gradient(gx, gy, n_bins, slope_unit) + else: + raise ValueError(f"Unknown distribution type: {distribution!r}. " + f"Choose from: theta, phi, gradient") + + # ------------------------------------------------------------------ + # phi: azimuthal angle distribution, weighted by slope² (slope_dist.c:430-466) + # phi = atan2(gy, -gx), canonicalized to [0, 2π); bin weight = gx²+gy² + # ------------------------------------------------------------------ + def _phi(self, gx, gy, n_bins, slope_unit): + phi = np.arctan2(gy, -gx) + phi = phi % (2.0 * np.pi) # canonicalize to [0, 2π) + weights = gx ** 2 + gy ** 2 + + bin_edges = np.linspace(0.0, 2.0 * np.pi, n_bins + 1) + counts = np.zeros(n_bins) + idx = np.clip(np.floor(n_bins * phi / (2.0 * np.pi)).astype(int), 0, n_bins - 1) + np.add.at(counts, idx, weights) + + centers_deg = np.degrees(0.5 * (bin_edges[:-1] + bin_edges[1:])) + y_unit = f"({slope_unit})^2" if slope_unit else "" + return (LineData(data=counts, x_axis=centers_deg, x_unit="deg", y_unit=y_unit),) + + # ------------------------------------------------------------------ + # theta: inclination angle in degrees, normalized probability density (slope_dist.c:468-513) + # theta = (180/π)·atan(|gradient|); normalized by size/(nc·max) + # ------------------------------------------------------------------ + def _theta(self, gx, gy, n_bins): + theta = np.degrees(np.arctan(np.hypot(gx, gy))) + max_theta = float(theta.max()) if theta.size > 0 else 90.0 + if max_theta == 0.0: + max_theta = 90.0 + + bin_edges = np.linspace(0.0, max_theta, n_bins + 1) + counts = np.zeros(n_bins) + idx = np.clip(np.floor(n_bins * theta / max_theta).astype(int), 0, n_bins - 1) + np.add.at(counts, idx, 1) + + nc = len(theta) + if nc > 0 and max_theta > 0: + counts = counts * (n_bins / (nc * max_theta)) + + centers = 0.5 * (bin_edges[:-1] + bin_edges[1:]) + return (LineData(data=counts, x_axis=centers, x_unit="deg", y_unit="1/deg"),) + + # ------------------------------------------------------------------ + # gradient: magnitude distribution, normalized probability density (slope_dist.c:515-560) + # normalized by size/(nc·max) + # ------------------------------------------------------------------ + def _gradient(self, gx, gy, n_bins, slope_unit): + grad = np.hypot(gx, gy) + max_grad = float(grad.max()) if grad.size > 0 else 1.0 + if max_grad == 0.0: + max_grad = 1.0 + + bin_edges = np.linspace(0.0, max_grad, n_bins + 1) + counts = np.zeros(n_bins) + idx = np.clip(np.floor(n_bins * grad / max_grad).astype(int), 0, n_bins - 1) + np.add.at(counts, idx, 1) + + nc = len(grad) + if nc > 0 and max_grad > 0: + counts = counts * (n_bins / (nc * max_grad)) + + centers = 0.5 * (bin_edges[:-1] + bin_edges[1:]) + y_unit = f"1/({slope_unit})" if slope_unit else "" + return (LineData(data=counts, x_axis=centers, x_unit=slope_unit, y_unit=y_unit),) diff --git a/backend/nodes/tip_blind_estimate.py b/backend/nodes/tip_blind_estimate.py new file mode 100644 index 0000000..39c82ed --- /dev/null +++ b/backend/nodes/tip_blind_estimate.py @@ -0,0 +1,582 @@ +from __future__ import annotations + +""" +Blind tip estimation following Gwyddion's morph_lib.c / tip.c. + +Reference: J. S. Villarrubia, J. Res. Natl. Inst. Stand. Technol. 102 (1997) 425. + +─── High-level algorithm ──────────────────────────────────────────────────── + +In AFM the measured image is NOT the true surface — it is the dilation of the +true surface by the tip shape. Blind estimation works backwards from the +measured image to recover an upper bound on the tip shape, without knowing +the true surface in advance. + +The key insight (Villarrubia 1997): wherever the measured image has a sharp +local maximum, the tip apex must have been in contact with the surface at that +point. The shape of the image near that peak therefore constrains the tip +geometry. Iterating over all such peaks and taking the intersection of all +constraints converges to the sharpest tip consistent with the measured data. + +─── Integer quantisation (matches Gwyddion exactly) ──────────────────────── + +All internal arithmetic uses int32 to avoid floating-point drift in +comparisons. The conversion is: + + step = (surface_max - surface_min) / 10 000 (MORPH_LIB_N = 10 000) + isurf[y,x] = floor( (surf[y,x] - surf_min) / step ) → values in [0, 10000] + itip[y,x] = floor( (tip[y,x] - tip_max) / step ) → values in [-10000, 0] + +Tips are stored max-zeroed: the apex pixel is always 0, all other pixels ≤ 0. +This is the convention used by every function in morph_lib.c. + +─── Tip initialisation ────────────────────────────────────────────────────── + +We start with itip0 = all zeros. In the max-zeroed convention this means +every pixel has the same height as the apex — i.e. the tip is perfectly flat +(infinitely wide). The estimation can only decrease pixel values, so the +flat start is the most conservative (widest) possible initial assumption. +Every iteration narrows the tip to fit the image data more tightly. + +─── Two estimation modes ──────────────────────────────────────────────────── + +partial (_itip_estimate_partial / _gwy_morph_lib_itip_estimate0): + Collects local maxima in the image and iterates refinement over those + peaks only. Fast; works well for images with sharp, isolated features. + +full (_itip_estimate_full / _gwy_morph_lib_itip_estimate): + Each iteration computes the morphological opening of the image (erosion + then re-dilation with the current tip). All pixels where the measured + image exceeds the opening by more than the threshold are processed. + Slower; more robust for images without distinct isolated peaks. + +─── Certainty map ─────────────────────────────────────────────────────────── + +After estimation the certainty map identifies which surface pixels are +reliably reconstructed. A surface pixel (sy, sx) is "certain" if exactly +ONE scanning position (i, j) in the image could have produced it — i.e. only +one tip placement was consistent with the measured height at that point. +Pixels with two or more possible tip placements are ambiguous and are left +unmarked (0). +""" + +import numpy as np +from scipy.ndimage import grey_erosion, grey_dilation + +from backend.node_registry import register_node +from backend.data_types import DataField + +# Number of quantisation levels. Gwyddion defines this as MORPH_LIB_N = 10000. +_MORPH_LIB_N = 10000.0 + +# Sentinel for "no contribution yet" in the dilation accumulator. +_INT_MIN = np.iinfo(np.int32).min + + +# --------------------------------------------------------------------------- +# Quantisation helpers +# --------------------------------------------------------------------------- + +def _quantize_surface(data: np.ndarray, surf_min: float, step: float) -> np.ndarray: + """ + Convert a float surface to an integer array in [0, MORPH_LIB_N]. + + Matches Gwyddion's i_datafield_to_field(surface, maxzero=FALSE, min=surf_min, step). + """ + return np.floor((data - surf_min) / step).astype(np.int32) + + +def _quantize_tip(tip_data: np.ndarray, tip_max: float, step: float) -> np.ndarray: + """ + Convert a float tip to a max-zeroed integer array (values ≤ 0). + + Matches Gwyddion's i_datafield_to_field(tip, maxzero=TRUE, min=tip_min, step). + With tip_min = 0 (our tips are always shifted so min = 0) this simplifies to: + itip[y,x] = floor( (tip[y,x] − tip_max) / step ) + """ + return np.floor((tip_data - tip_max) / step).astype(np.int32) + + +def _dequantize(idata: np.ndarray, offset: float, step: float) -> np.ndarray: + """Reverse of quantisation: float = idata * step + offset.""" + return idata.astype(np.float64) * step + offset + + +# --------------------------------------------------------------------------- +# _useit — is pixel (x, y) a local maximum suitable for tip refinement? +# Matches Gwyddion's static useit() in morph_lib.c +# --------------------------------------------------------------------------- + +def _useit(image: np.ndarray, x: int, y: int, delta: int) -> bool: + """ + Return True if (x, y) should be used as a refinement point. + + A pixel qualifies when: + 1. It is the maximum value within a (2·delta + 1)² neighbourhood. + 2. Fewer than 1/5 of the neighbourhood pixels share that maximum + (to exclude flat plateaux, which give no tip shape information). + + delta = max(max(txres, tyres) // 10, 1) so larger tips use a wider + neighbourhood, matching Gwyddion's choice. + """ + yres, xres = image.shape + xmin, xmax = max(x - delta, 0), min(x + delta, xres - 1) + ymin, ymax = max(y - delta, 0), min(y + delta, yres - 1) + region = image[ymin:ymax + 1, xmin:xmax + 1] + max_val = int(region.max()) + if int(image[y, x]) < max_val: + return False # not the maximum + count = int((region >= max_val).sum()) + return count <= (2 * delta + 1) ** 2 // 5 # not a flat plateau + + +# --------------------------------------------------------------------------- +# _estimate_point_interior — refine tip from a single interior image pixel +# Matches Gwyddion's itip_estimate_point() interior branch in morph_lib.c +# --------------------------------------------------------------------------- + +def _estimate_point_interior( + ixp: int, jxp: int, + image: np.ndarray, + txres: int, tyres: int, + xc: int, yc: int, + tip0: np.ndarray, + thresh: int, +) -> int: + """ + Use the image pixel at (ixp, jxp) to narrow the tip estimate tip0. + + Returns the number of tip pixels that were updated (refined). + + ── Geometric meaning ────────────────────────────────────────────────── + + Imagine placing the tip so its apex (xc, yc) is directly above image + pixel (ixp, jxp). For each tip pixel (jd, id) offset from the apex: + + The tip "touches" the image at reference pixel (jxp+yc-jd, ixp+xc-id). + For the tip to be inside the surface (physically valid), the image + height at that reference pixel must be ≥ image_p − tip0[jd, id]. + + Pixels that pass this check are called "good pixels" — they are the tip + positions that are geometrically consistent with the apex sitting at + (ixp, jxp). + + ── Computing the dilation at each tip pixel ─────────────────────────── + + For each tip output pixel (jx, ix) we compute the maximum height that + the tip could produce at any scanning position consistent with the good + pixels. This is: + + dil[jx, ix] = max over good (jd, id) of: + image[jxp + jx - jd, ixp + ix - id] + tip0[jd, id] - image_p + + Intuitively: if the tip pixel (jd, id) was in contact with the surface + at the reference point, then the image at the laterally shifted position + (jxp+jx-jd, ixp+ix-id) constrains how high tip pixel (jx, ix) can be. + + ── Updating the tip ─────────────────────────────────────────────────── + + If dil[jx, ix] < tip0[jx, ix] - thresh, the new constraint is tighter + than the current estimate, so we update: + tip0[jx, ix] = dil[jx, ix] + thresh + + The +thresh prevents over-sharpening due to noise: we only accept a + refinement if it exceeds the current estimate by more than the threshold. + + ── Vectorisation ────────────────────────────────────────────────────── + + For each good pixel (jd, id) we extract a (tyres x txres) window from + the image, shifted so that image[jxp-jd, ixp-id] aligns with tip[0, 0]. + Taking the element-wise max across all good pixels yields dil efficiently + without the O(tyres² x txres²) double loop of the original C code. + """ + image_p = int(image[jxp, ixp]) + + # Build image_at_ref[jd, id] = image[jxp+yc-jd, ixp+xc-id]. + # We want rows jxp+yc down to jxp+yc-tyres+1 (reversed), which numpy + # gives as np.flip( image[jxp+yc-tyres+1 : jxp+yc+1, ...] ). + image_at_ref = np.flip( + image[jxp + yc - tyres + 1:jxp + yc + 1, + ixp + xc - txres + 1:ixp + xc + 1], + axis=(0, 1), + ) # shape (tyres, txres); integer values + + # Good-pixel mask: image_p - image_at_ref[jd,id] <= tip0[jd,id] + # Rearranged: image_at_ref[jd,id] >= image_p - tip0[jd,id] + # Since tip0 ≤ 0, the right-hand side is ≥ image_p. + # A good pixel is one where the reference height is high enough that + # placing the tip apex at (ixp, jxp) is geometrically consistent. + diff = image_p - image_at_ref # positive where the reference is low + good_mask = diff <= tip0 # True for geometrically valid placements + good_jd, good_id = np.where(good_mask) + + if good_jd.size == 0: + # No valid tip placement found at this image point; nothing to refine. + return 0 + + # Accumulate dil[jx, ix] = max over good (jd, id) of the window contribution. + # Initialise to _INT_MIN so we can detect pixels with no contribution. + dil = np.full((tyres, txres), _INT_MIN, dtype=np.int64) + + for k in range(good_jd.size): + jd = int(good_jd[k]) + id_ = int(good_id[k]) + # For this good pixel (jd, id), the contribution to dil[jx, ix] is: + # image[jxp + jx - jd, ixp + ix - id] + tip0[jd, id] - image_p + # Extracting the (tyres × txres) window starting at (jxp-jd, ixp-id) + # gives all (jx, ix) contributions simultaneously. + row0 = jxp - jd + col0 = ixp - id_ + window = image[row0:row0 + tyres, col0:col0 + txres].astype(np.int64) + contrib = window + int(tip0[jd, id_]) - image_p + np.maximum(dil, contrib, out=dil) + + # Update tip0 where the new constraint is strictly tighter. + # Only pixels with at least one good-pixel contribution are candidates. + update = (dil > _INT_MIN) & (dil < tip0.astype(np.int64) - thresh) + count = int(update.sum()) + if count: + tip0[update] = (dil[update] + thresh).astype(np.int32) + return count + + +# --------------------------------------------------------------------------- +# Partial estimation — iterate over local maxima only +# Matches _gwy_morph_lib_itip_estimate0() in morph_lib.c +# --------------------------------------------------------------------------- + +def _itip_estimate_partial( + image: np.ndarray, + txres: int, tyres: int, + xc: int, yc: int, + tip0: np.ndarray, + thresh: int, + use_edges: bool, +) -> int: + """ + Refine tip0 using the partial (local-maxima) method. + + Steps: + 1. Scan the image for pixels that are local maxima in a ±delta + neighbourhood and are not on a flat plateau (useit criterion). + These are the image features most likely to have been produced by + the tip apex — they constrain the tip shape most tightly. + 2. Iterate _estimate_point_interior over every qualifying pixel until + either no pixel produces a refinement, or fewer than 20 pixels do + (convergence criterion matching Gwyddion's maxcount = 20). + + Returns the total number of tip-pixel updates across all iterations. + """ + yres, xres = image.shape + + # delta: neighbourhood radius for the local-maximum test. + # Larger tips look over a wider window to avoid spurious maxima. + delta = max(max(txres, tyres) // 10, 1) + maxcount = 20 # convergence: stop when ≤ maxcount pixels updated per pass + + # Step 1: collect all valid local maxima in the interior. + # Interior constraint: the full tip must fit inside the image at every + # candidate point, so we stay ≥ (tip_half_size) from each edge. + maxima_x, maxima_y = [], [] + for jxp in range(tyres - 1 - yc, yres - yc): + for ixp in range(txres - 1 - xc, xres - xc): + if _useit(image, ixp, jxp, delta): + maxima_x.append(ixp) + maxima_y.append(jxp) + + if not maxima_x: + return 0 # no usable features found; tip stays flat + + # Step 2: iterate until convergence. + sumcount = 0 + while True: + count = 0 + for ixp, jxp in zip(maxima_x, maxima_y): + # Restrict to strictly interior points (tip fits fully inside image). + if (jxp >= tyres - 1 and jxp <= yres - tyres + and ixp >= txres - 1 and ixp <= xres - txres): + count += _estimate_point_interior( + ixp, jxp, image, txres, tyres, xc, yc, tip0, thresh) + sumcount += count + if count == 0 or count <= maxcount: + break # converged + + return sumcount + + +# --------------------------------------------------------------------------- +# Full estimation — iterate over all points above the morphological opening +# Matches _gwy_morph_lib_itip_estimate() in morph_lib.c +# --------------------------------------------------------------------------- + +def _itip_estimate_full( + image: np.ndarray, + txres: int, tyres: int, + xc: int, yc: int, + tip0: np.ndarray, + thresh: int, + use_edges: bool, +) -> int: + """ + Refine tip0 using the full estimation method. + + Each iteration: + 1. Computes the morphological opening of the image with the current tip + estimate: opening = dilate(erode(image, tip0), tip0). + The opening removes features narrower than the tip; any pixel where + image > opening indicates that the tip apex was directly above a + surface feature that is sharper than the current tip estimate. + 2. Calls _estimate_point_interior for every such pixel. + 3. Repeats until no pixel produces a refinement (count == 0). + + This is slower than the partial method but does not require isolated + sharp peaks — it converges from any image with sufficient surface relief. + """ + yres, xres = image.shape + sumcount = 0 + + while True: + # Morphological opening with the current integer tip estimate. + # scipy grey_erosion/dilation treat the structure as centred at + # its middle pixel, matching Gwyddion's convention of apex at (yc, xc). + tip_float = tip0.astype(np.float64) + eroded = grey_erosion(image.astype(np.float64), structure=tip_float) + opened = grey_dilation(eroded, structure=tip_float) + opened_int = opened.astype(np.int64) + + count = 0 + for jxp in range(tyres - 1, yres - tyres + 1): + for ixp in range(txres - 1, xres - txres + 1): + # Only process pixels where the image exceeds its opening by + # more than the noise threshold. These are genuine sharp + # features that cannot be explained by the current tip shape. + if int(image[jxp, ixp]) - int(opened_int[jxp, ixp]) > thresh: + count += _estimate_point_interior( + ixp, jxp, image, txres, tyres, xc, yc, tip0, thresh) + sumcount += count + if count == 0: + break # no more refinements possible: converged + + return sumcount + + +# --------------------------------------------------------------------------- +# _certainty_map_fast — vectorised certainty map +# Matches _gwy_morph_lib_icmap() in morph_lib.c and the floating-point +# certainty_map() helper in tip.c (disabled #if 0 branch, which is cleaner). +# --------------------------------------------------------------------------- + +def _certainty_map_fast( + image: np.ndarray, + tip: np.ndarray, + rsurf: np.ndarray, + xc: int, yc: int, + cmap_thresh: float, +) -> np.ndarray: + """ + Compute the certainty map for the reconstructed surface. + + A surface pixel (sy, sx) is marked certain (= 1) if there is exactly ONE + scanning position (image pixel i, j) at which the tip was unambiguously + in single-point contact with the surface at (sy, sx). + + ── Contact condition ──────────────────────────────────────────────────── + + For tip pixel (ti, tj) and its reflected counterpart tip_flip[ti, tj]: + + image[i, j] − tip_flip[ti, tj] − rsurf[sy, sx] ≈ 0 + + means that placing the tip so that tip pixel (ti, tj) touches the surface + at (sy, sx) exactly accounts for the measured height at image pixel (i, j). + The reflected tip is used because dilation/erosion use the tip flipped + 180°. Coordinates: + sy = ti + i − ryc, sx = tj + j − rxc + where ryc = tyres−1−yc and rxc = txres−1−xc are the reflected-apex + offsets. + + ── Why "exactly one" contact? ─────────────────────────────────────────── + + If two or more tip placements both satisfy the contact condition, the + reconstruction is ambiguous — either placement could have produced the + observed image value. Only when a single tip placement is consistent is + the surface height at (sy, sx) uniquely determined. + + ── Vectorisation strategy ─────────────────────────────────────────────── + + Instead of looping over image pixels (O(yres × xres)) and for each one + checking all tip pixels (O(tyres × txres)), we loop over tip pixels + (typically much smaller than the image) and for each tip pixel accumulate + a contact-count array over the surface. This turns the O(N² × T²) scalar + loop into an O(T² × N²/T²) = O(N²) numpy operation per tip pixel. + """ + yres, xres = image.shape + tyres, txres = tip.shape + # Reflected-apex offsets: ryc = tyres-1-yc, rxc = txres-1-xc. + # For a centred tip (yc = tyres//2), ryc = tyres//2 as well. + ryc = tyres - 1 - yc + rxc = txres - 1 - xc + + # The certainty map algorithm uses the tip rotated 180° (reflected both + # axes), because the erosion that produced rsurf used the flipped tip. + tip_flip = np.flipud(np.fliplr(tip)).astype(np.float64) + + # count_arr[sy, sx] = number of tip placements that produce a contact at (sy, sx). + count_arr = np.zeros((yres, xres), dtype=np.int32) + + for ti in range(tyres): + for tj in range(txres): + tf_val = tip_flip[ti, tj] + + # Determine which surface pixels (sy, sx) are reachable from valid + # interior image pixels (i, j) when tip pixel (ti, tj) is in contact: + # sy = ti + i - ryc → sy in [ti, yres - tyres + ti] (clipped to image) + # sx = tj + j - rxc → sx in [tj, xres - txres + tj] + sy_lo = max(ti, 0) + sy_hi = min(yres - tyres + ti, yres - 1) + sx_lo = max(tj, 0) + sx_hi = min(xres - txres + tj, xres - 1) + if sy_lo > sy_hi or sx_lo > sx_hi: + continue # this tip pixel is entirely out of range + + sy_range = np.arange(sy_lo, sy_hi + 1) + sx_range = np.arange(sx_lo, sx_hi + 1) + + # Corresponding image pixels for this tip placement: + i_range = sy_range + ryc - ti + j_range = sx_range + rxc - tj + + # Check the contact condition for every (sy, sx) / (i, j) pair: + # |image[i, j] - tip_flip[ti, tj] - rsurf[sy, sx]| <= cmap_thresh + img_block = image[np.ix_(i_range, j_range)] + rs_block = rsurf[np.ix_(sy_range, sx_range)] + contact = np.abs(img_block - tf_val - rs_block) <= cmap_thresh + + # Accumulate into the contact-count array. + count_arr[np.ix_(sy_range, sx_range)] += contact.astype(np.int32) + + # Surface pixels with exactly one contact are certain; all others are not. + return (count_arr == 1).astype(np.float64) + + +# --------------------------------------------------------------------------- +# Node +# --------------------------------------------------------------------------- + +@register_node(display_name="Blind Tip Estimate") +class BlindTipEstimate: + @classmethod + def INPUT_TYPES(cls): + return { + "required": { + "field": ("DATA_FIELD",), + "n_pixels": ("INT", {"default": 33, "min": 3, "max": 129, "step": 2}), + "threshold": ("FLOAT", { + "default": 0.0, "min": 0.0, "max": 1e-6, "step": 1e-10, + }), + "method": (["partial", "full"], {"default": "partial"}), + "use_edges": ("BOOLEAN", {"default": False}), + } + } + + OUTPUTS = ( + ('DATA_FIELD', 'tip'), + ('DATA_FIELD', 'certainty'), + ) + FUNCTION = "process" + + DESCRIPTION = ( + "Blind tip estimation from a measured SPM image (Villarrubia algorithm). " + "Iteratively refines the tip shape using only the sharpest image features. " + "partial: uses local maxima only (faster, needs sharp isolated features). " + "full: uses all points above morphological opening (slower, more robust). " + "threshold: noise floor in metres — start at 0 and increase if tip is too sharp. " + "Output tip has apex=max, edges=0 (same convention as TipModel). " + "Certainty map marks surface pixels where the tip was in unambiguous single contact. " + "Equivalent to gwy_tip_estimate_partial / gwy_tip_estimate_full + gwy_tip_cmap (tip.c)." + ) + + def process( + self, + field: DataField, + n_pixels: int, + threshold: float, + method: str, + use_edges: bool, + ) -> tuple: + # ── Tip grid setup ──────────────────────────────────────────────── + # Ensure odd size so the apex pixel is exactly centred. + n = n_pixels if n_pixels % 2 == 1 else n_pixels + 1 + txres = tyres = n + xc = yc = n // 2 # apex is at the centre pixel + + # ── Integer quantisation ────────────────────────────────────────── + # All estimation arithmetic is performed in integer units to match + # Gwyddion's morph_lib.c exactly. step is chosen so the full surface + # height range maps to MORPH_LIB_N = 10 000 integer levels. + surf = field.data + surf_min = float(surf.min()) + surf_max = float(surf.max()) + surf_range = surf_max - surf_min + if surf_range == 0.0: + surf_range = 1.0 # guard against featureless (flat) images + + step = surf_range / _MORPH_LIB_N + isurf = _quantize_surface(surf, surf_min, step) # int32 in [0, 10000] + + # ── Initial tip: flat (all zeros) ──────────────────────────────── + # In the max-zeroed convention, all zeros means every tip pixel is at + # the same height as the apex — a perfectly flat (infinitely wide) tip. + # This is the most conservative starting point; the algorithm can only + # decrease tip values (narrow the tip), never increase them. + itip0 = np.zeros((tyres, txres), dtype=np.int32) + + # ── Noise threshold in integer units ───────────────────────────── + # A refinement is only accepted if the new constraint is more than + # thresh_int quantisation levels below the current estimate. + # Gwyddion recommends starting at 0 and increasing if the result is + # too sharp (noise is interpreted as sharp features). + thresh_int = max(int(threshold / step), 0) + + # ── Run estimation ──────────────────────────────────────────────── + if method == "partial": + _itip_estimate_partial(isurf, txres, tyres, xc, yc, itip0, thresh_int, use_edges) + else: + _itip_estimate_full(isurf, txres, tyres, xc, yc, itip0, thresh_int, use_edges) + + # ── Dequantise tip back to physical units ───────────────────────── + # Matches Gwyddion's i_field_to_datafield then gwy_data_field_add(-min): + # tip_data[y,x] = itip0[y,x] * step (offset = tipmin = 0 for flat start) + # then shift so minimum = 0 (apex = maximum after shift). + tip_data = itip0.astype(np.float64) * step + tip_data -= tip_data.min() # shift: minimum → 0, apex → maximum + + pixel_size = (field.dx + field.dy) * 0.5 + xreal = n * pixel_size + + tip_field = DataField( + data=tip_data, + xreal=xreal, + yreal=xreal, + si_unit_xy=field.si_unit_xy, + si_unit_z=field.si_unit_z, + ) + + # ── Certainty map ───────────────────────────────────────────────── + # Reconstruct the surface by eroding the measured image with the + # estimated tip (this is gwy_tip_erosion / TipDeconvolution). + # tip_max_zeroed converts back to the erosion convention (max = 0). + tip_max_zeroed = tip_data - tip_data.max() # values ≤ 0, apex = 0 + rsurf = grey_erosion(surf, structure=tip_max_zeroed) + + # cmap_thresh: tolerance for the "exact contact" comparison. + # Gwyddion uses 50 × step (the icmap integer branch uses equality, + # which corresponds to ±0.5 step; 50 steps gives a looser float match). + cmap_thresh = 50.0 * step + cmap_data = _certainty_map_fast(surf, tip_data, rsurf, xc, yc, cmap_thresh) + + cmap_field = field.replace( + data=cmap_data, + si_unit_z="", # certainty is dimensionless + ) + + return (tip_field, cmap_field) diff --git a/backend/nodes/tip_deconvolution.py b/backend/nodes/tip_deconvolution.py new file mode 100644 index 0000000..2a28b34 --- /dev/null +++ b/backend/nodes/tip_deconvolution.py @@ -0,0 +1,44 @@ +from __future__ import annotations + +import numpy as np +from scipy.ndimage import grey_erosion + +from backend.node_registry import register_node +from backend.data_types import DataField + + +@register_node(display_name="Tip Deconvolution") +class TipDeconvolution: + @classmethod + def INPUT_TYPES(cls): + return { + "required": { + "field": ("DATA_FIELD",), + "tip": ("DATA_FIELD",), + } + } + + OUTPUTS = ( + ('DATA_FIELD', 'surface'), + ) + FUNCTION = "process" + + DESCRIPTION = ( + "Reconstruct the true surface from a tip-broadened measured image. " + "Uses morphological grey erosion (Villarrubia algorithm): " + " mytip = flip(tip) − max(flip(tip)) [max shifted to 0] " + " surface[y,x] = min_{dy,dx}[image[y+dy, x+dx] − mytip[dy,dx]] " + "Connect the tip output from a TipModel node. " + "The tip pixel size must match the image pixel size. " + "Equivalent to gwy_tip_erosion (tip.c)." + ) + + def process(self, field: DataField, tip: DataField) -> tuple: + # Gwyddion gwy_tip_erosion: + # mytip = flip(tip) − max(flip(tip)) (values ≤ 0, apex = 0) + # result[y,x] = min_{ty,tx}[surface[y+ty, x+tx] − mytip[ty,tx]] + tip_flipped = np.flipud(np.fliplr(tip.data)) + mytip = tip_flipped - tip_flipped.max() # shift so max = 0 + + result = grey_erosion(field.data, structure=mytip) + return (field.replace(data=result),) diff --git a/backend/nodes/tip_model.py b/backend/nodes/tip_model.py new file mode 100644 index 0000000..57b1444 --- /dev/null +++ b/backend/nodes/tip_model.py @@ -0,0 +1,109 @@ +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="Tip Model") +class TipModel: + @classmethod + def INPUT_TYPES(cls): + return { + "required": { + "field": ("DATA_FIELD",), + "shape": (["parabola", "cone", "sphere"], {"default": "parabola"}), + "radius": ("FLOAT", { + "default": 10e-9, "min": 1e-10, "max": 1e-6, "step": 1e-10, + }), + "half_angle": ("FLOAT", { + "default": 20.0, "min": 1.0, "max": 89.0, "step": 0.5, + }), + "n_pixels": ("INT", {"default": 65, "min": 3, "max": 511, "step": 2}), + } + } + + OUTPUTS = ( + ('DATA_FIELD', 'tip'), + ) + FUNCTION = "process" + + DESCRIPTION = ( + "Generate a synthetic AFM tip model DATA_FIELD. " + "The input field sets the pixel size for the tip. " + "The apex (centre pixel) is the maximum value; edges are shifted to zero. " + "Shapes: parabola — paraboloid with apex radius R; " + "cone — sphere-capped cone (radius R, half_angle from tip axis in degrees); " + "sphere — ball-on-stick (sphere cap only). " + "Equivalent to gwy_tip_model_preset_create (tip.c / tip_model.c)." + ) + + def process( + self, + field: DataField, + shape: str, + radius: float, + half_angle: float, + n_pixels: int, + ) -> tuple: + pixel_size = (field.dx + field.dy) * 0.5 + + # Ensure odd size so the centre pixel is well-defined + n = n_pixels if n_pixels % 2 == 1 else n_pixels + 1 + ci = n // 2 + + # Physical offsets from the centre pixel + offsets = (np.arange(n) - ci) * pixel_size + gx, gy = np.meshgrid(offsets, offsets) + r2 = gx**2 + gy**2 + r = np.sqrt(r2) + + if shape == "parabola": + # Gwyddion parabola(): a = 1/(2R), z0 = 2a·x_half² + # z[y,x] = z0 − a·r² → min at corners is exactly 0 + a = 0.5 / radius + x_half = ci * pixel_size # half-width in physical units + z0 = 2.0 * a * x_half**2 + data = z0 - a * r2 + data -= data.min() # shift so min = 0 + + elif shape == "cone": + # Gwyddion cone(): + # angle = half-angle from horizontal = π/2 − half_angle_from_axis + # z0 = R/sin(angle) + # r_cross² = (R·cos(angle))² + # inner (r² < r_cross²): z = sqrt(R² − r²) ← spherical cap + # outer: z = z0 − r/tan(angle) + angle = np.radians(90.0 - half_angle) # slope from horizontal + z0 = radius / np.sin(angle) + br2 = (radius * np.cos(angle))**2 + ta = 1.0 / np.tan(angle) + inner = np.sqrt(np.maximum(0.0, radius**2 - r2)) + outer = z0 - ta * r + data = np.where(r2 < br2, inner, outer) + data -= data.min() + + elif shape == "sphere": + # Gwyddion sphere(): cone with near-vertical sides (angle ≈ 0 from horizontal) + angle = 1e-5 # radians — essentially vertical cone + z0 = radius / np.sin(angle) + br2 = (radius * np.cos(angle))**2 + ta = 1.0 / np.tan(angle) + inner = np.sqrt(np.maximum(0.0, radius**2 - r2)) + outer = z0 - ta * r + data = np.where(r2 < br2, inner, outer) + data -= data.min() + + else: + raise ValueError(f"Unknown tip shape {shape!r}. Choose: parabola, cone, sphere.") + + xreal = n * pixel_size + tip_field = DataField( + data=data, + xreal=xreal, + yreal=xreal, + si_unit_xy=field.si_unit_xy, + si_unit_z=field.si_unit_z, + ) + return (tip_field,) diff --git a/frontend/src/CustomNode.jsx b/frontend/src/CustomNode.jsx index bed5e7b..6f3c5cb 100644 --- a/frontend/src/CustomNode.jsx +++ b/frontend/src/CustomNode.jsx @@ -254,6 +254,51 @@ class PreviewBoundary extends React.Component { } } +// ── SI prefix helpers ───────────────────────────────────────────────── + +const _SI_PREFIXES = [ + { prefix: 'T', factor: 1e12 }, + { prefix: 'G', factor: 1e9 }, + { prefix: 'M', factor: 1e6 }, + { prefix: 'k', factor: 1e3 }, + { prefix: '', factor: 1 }, + { prefix: 'm', factor: 1e-3 }, + { prefix: 'μ', factor: 1e-6 }, + { prefix: 'n', factor: 1e-9 }, + { prefix: 'p', factor: 1e-12 }, + { prefix: 'f', factor: 1e-15 }, +]; + +// Map of suffix characters → multiplier (accept both 'u' and 'μ' for micro) +const _SI_PARSE_MAP = { T:1e12, G:1e9, M:1e6, k:1e3, m:1e-3, u:1e-6, μ:1e-6, n:1e-9, p:1e-12, f:1e-15 }; + +function formatSI(v, prec) { + if (!Number.isFinite(v)) return String(v); + if (v === 0) return prec != null ? `0.${'0'.repeat(prec)}` : '0'; + const abs = Math.abs(v); + // Pick the largest SI prefix whose factor is ≤ |v| (gives value in [1, 1000)) + let chosen = _SI_PREFIXES[_SI_PREFIXES.length - 1]; + for (const p of _SI_PREFIXES) { + if (abs >= p.factor * (1 - 1e-10)) { chosen = p; break; } + } + const scaled = v / chosen.factor; + return (prec != null ? scaled.toFixed(prec) : String(scaled)) + chosen.prefix; +} + +// Parse a string that may carry an SI suffix (e.g. "20n", "1.5μ", "500p") +// Falls back to standard parseFloat for plain numbers and scientific notation. +function parseSI(text) { + const t = (text || '').trim(); + if (!t) return NaN; + const lastChar = t.slice(-1); + const factor = _SI_PARSE_MAP[lastChar]; + if (factor != null) { + const num = parseFloat(t.slice(0, -1)); + if (!isNaN(num)) return num * factor; + } + return parseFloat(t); +} + // ── Draggable number input ──────────────────────────────────────────── function DraggableNumber({ value, step, min, max, precision, onChange }) { @@ -262,7 +307,7 @@ function DraggableNumber({ value, step, min, max, precision, onChange }) { const dragState = useRef(null); const elRef = useRef(null); - const display = precision != null ? Number(value).toFixed(precision) : String(value); + const display = precision != null ? formatSI(Number(value), precision) : String(value); const clamp = useCallback((v) => { if (min != null && v < min) v = min; @@ -282,9 +327,7 @@ function DraggableNumber({ value, step, min, max, precision, onChange }) { const dx = e.clientX - dragState.current.startX; const delta = dx * (step || 0.01); const raw = dragState.current.startVal + delta; - const rounded = precision != null - ? parseFloat(raw.toFixed(precision)) - : Math.round(raw); + const rounded = precision != null ? raw : Math.round(raw); onChange(clamp(rounded)); }, [step, precision, clamp, onChange]); @@ -308,16 +351,14 @@ function DraggableNumber({ value, step, min, max, precision, onChange }) { const delta = (e.deltaY < 0 ? 1 : -1) * baseStep * multiplier; const startVal = Number(value); const raw = (Number.isFinite(startVal) ? startVal : 0) + delta; - const rounded = precision != null - ? parseFloat(raw.toFixed(precision)) - : Math.round(raw); + const rounded = precision != null ? raw : Math.round(raw); onChange(clamp(rounded)); }, [editing, step, value, precision, onChange, clamp]); const commitEdit = useCallback(() => { setEditing(false); - const parsed = parseFloat(editText); - if (!isNaN(parsed)) onChange(clamp(precision != null ? parseFloat(parsed.toFixed(precision)) : Math.round(parsed))); + const parsed = parseSI(editText); + if (!isNaN(parsed)) onChange(clamp(precision != null ? parsed : Math.round(parsed))); }, [editText, precision, clamp, onChange]); if (editing) { diff --git a/tests/node_tests/acf_1d.py b/tests/node_tests/acf_1d.py index 63bc8f6..e53bd11 100644 --- a/tests/node_tests/acf_1d.py +++ b/tests/node_tests/acf_1d.py @@ -24,9 +24,12 @@ def test_acf_1d(): assert isinstance(acf, LineData) assert isinstance(measurement, RecordTable) - # ACF should be symmetric about zero lag - center = len(acf) // 2 - assert np.allclose(acf.data, acf.data[::-1], atol=1e-10) + # Only positive lags are output + assert acf.x_axis is not None + assert np.all(acf.x_axis > 0) + + # ACF should be monotonically decreasing at the very start (falls from the zero-lag peak) + assert acf.data[0] > 0 # Peak period should be close to the input period in metres expected_period_m = period * 1e-9 @@ -35,13 +38,6 @@ def test_acf_1d(): assert abs(measurement[0]["value"] - expected_period_m) / expected_period_m < 0.1 assert measurement[0]["unit"] == "m" - # x_axis should be centred on zero - assert acf.x_axis is not None - assert acf.x_axis[center] == 0.0 or abs(acf.x_axis[center]) < 1e-15 - - # ACF at zero lag should equal variance (signal is mean-subtracted) - assert acf.data[center] > 0 - def test_acf_1d_no_peak(): from backend.nodes.acf_1d import ACF1D @@ -68,5 +64,5 @@ def test_acf_1d_level_none(): profile = LineData(data=data, x_axis=np.arange(32, dtype=np.float64)) acf, _ = node.process(profile, level="none") - # ACF of a constant is a constant - assert acf.data[len(acf) // 2] > 0 + # ACF of a constant is a constant — first positive-lag value should be positive + assert acf.data[0] > 0 diff --git a/tests/node_tests/field_arithmetic.py b/tests/node_tests/field_arithmetic.py new file mode 100644 index 0000000..05e660b --- /dev/null +++ b/tests/node_tests/field_arithmetic.py @@ -0,0 +1,75 @@ +import numpy as np +import pytest +from tests.node_tests._shared import make_field + + +def test_field_arithmetic_basic(): + from backend.nodes.field_arithmetic import FieldArithmetic + + node = FieldArithmetic() + a = make_field(data=np.array([[1.0, 2.0], [3.0, 4.0]])) + b = make_field(data=np.array([[1.0, 1.0], [1.0, 1.0]])) + + result, = node.process(a, b, "add") + assert np.allclose(result.data, [[2.0, 3.0], [4.0, 5.0]]) + + result, = node.process(a, b, "subtract") + assert np.allclose(result.data, [[0.0, 1.0], [2.0, 3.0]]) + + result, = node.process(a, b, "multiply") + assert np.allclose(result.data, [[1.0, 2.0], [3.0, 4.0]]) + + result, = node.process(a, b, "divide") + assert np.allclose(result.data, [[1.0, 2.0], [3.0, 4.0]]) + + result, = node.process(a, b, "min") + assert np.allclose(result.data, [[1.0, 1.0], [1.0, 1.0]]) + + result, = node.process(a, b, "max") + assert np.allclose(result.data, [[1.0, 2.0], [3.0, 4.0]]) + + +def test_field_arithmetic_hypot(): + from backend.nodes.field_arithmetic import FieldArithmetic + + node = FieldArithmetic() + a = make_field(data=np.array([[3.0, 0.0], [0.0, 5.0]])) + b = make_field(data=np.array([[4.0, 5.0], [3.0, 12.0]])) + + result, = node.process(a, b, "hypot") + assert np.allclose(result.data, [[5.0, 5.0], [3.0, 13.0]]) + + +def test_field_arithmetic_metadata_inherited(): + from backend.nodes.field_arithmetic import FieldArithmetic + + node = FieldArithmetic() + a = make_field(data=np.ones((4, 4))) + b = make_field(data=np.ones((4, 4))) + + result, = node.process(a, b, "add") + assert result.xreal == a.xreal + assert result.si_unit_xy == a.si_unit_xy + assert result.si_unit_z == a.si_unit_z + + +def test_field_arithmetic_shape_mismatch(): + from backend.nodes.field_arithmetic import FieldArithmetic + + node = FieldArithmetic() + a = make_field(data=np.ones((4, 4))) + b = make_field(data=np.ones((3, 3))) + + with pytest.raises(ValueError, match="resolution"): + node.process(a, b, "add") + + +def test_field_arithmetic_unknown_operation(): + from backend.nodes.field_arithmetic import FieldArithmetic + + node = FieldArithmetic() + a = make_field(data=np.ones((4, 4))) + b = make_field(data=np.ones((4, 4))) + + with pytest.raises(ValueError): + node.process(a, b, "power") diff --git a/tests/node_tests/gradient.py b/tests/node_tests/gradient.py new file mode 100644 index 0000000..e930ca7 --- /dev/null +++ b/tests/node_tests/gradient.py @@ -0,0 +1,91 @@ +import numpy as np +import pytest +from tests.node_tests._shared import make_field + + +def test_gradient_flat_field(): + from backend.nodes.gradient import Gradient + + node = Gradient() + field = make_field(data=np.zeros((16, 16))) + + for component in ("magnitude", "x", "y"): + result, = node.process(field, component) + assert np.allclose(result.data, 0.0) + assert result.data.shape == field.data.shape + + +def test_gradient_x_ramp(): + """Linear ramp in x should give uniform x-gradient and zero y-gradient.""" + from backend.nodes.gradient import Gradient + + node = Gradient() + data = np.tile(np.linspace(0.0, 1.0, 32), (32, 1)) + field = make_field(data=data) + + gx, = node.process(field, "x") + gy, = node.process(field, "y") + + # Interior pixels should have nearly identical x-gradient values + interior_gx = gx.data[1:-1, 1:-1] + assert np.all(interior_gx > 0) + assert np.std(interior_gx) < np.mean(interior_gx) * 0.01 + + # y-gradient should be zero everywhere in the interior + assert np.allclose(gy.data[1:-1, 1:-1], 0.0, atol=1e-10) + + +def test_gradient_magnitude_non_negative(): + from backend.nodes.gradient import Gradient + + node = Gradient() + field = make_field() + + result, = node.process(field, "magnitude") + assert np.all(result.data >= 0.0) + + +def test_gradient_azimuth_range(): + from backend.nodes.gradient import Gradient + + node = Gradient() + field = make_field() + + result, = node.process(field, "azimuth") + assert np.all(result.data >= -np.pi - 1e-10) + assert np.all(result.data <= np.pi + 1e-10) + assert result.si_unit_z == "rad" + + +def test_gradient_units(): + from backend.nodes.gradient import Gradient + + node = Gradient() + field = make_field() # si_unit_z="m", si_unit_xy="m" + + result, = node.process(field, "magnitude") + assert result.si_unit_z == "m/m" + + result, = node.process(field, "x") + assert result.si_unit_z == "m/m" + + +def test_gradient_shape_preserved(): + from backend.nodes.gradient import Gradient + + node = Gradient() + field = make_field(shape=(48, 64)) + + for component in ("magnitude", "x", "y", "azimuth"): + result, = node.process(field, component) + assert result.data.shape == (48, 64) + + +def test_gradient_unknown_component(): + from backend.nodes.gradient import Gradient + + node = Gradient() + field = make_field() + + with pytest.raises(ValueError): + node.process(field, "diagonal") diff --git a/tests/node_tests/grain_filter.py b/tests/node_tests/grain_filter.py new file mode 100644 index 0000000..699fd9c --- /dev/null +++ b/tests/node_tests/grain_filter.py @@ -0,0 +1,91 @@ +import numpy as np +import pytest + + +def _make_mask(*rects): + """Create a 30x30 uint8 mask with 255 in the given (row_start, row_end, col_start, col_end) rects.""" + m = np.zeros((30, 30), dtype=np.uint8) + for r0, r1, c0, c1 in rects: + m[r0:r1, c0:c1] = 255 + return m + + +def test_grain_filter_min_area(): + from backend.nodes.grain_filter import GrainFilter + + node = GrainFilter() + # Small grain: 2×2 = 4 px; large grain: 6×6 = 36 px + mask = _make_mask((1, 3, 1, 3), (15, 21, 15, 21)) + + # Remove grains smaller than 10 px → only large grain survives + result, = node.process(mask, min_area=10, max_area=0, remove_border=False) + assert result[1:3, 1:3].sum() == 0 # small grain removed + assert result[15:21, 15:21].sum() > 0 # large grain kept + + +def test_grain_filter_max_area(): + from backend.nodes.grain_filter import GrainFilter + + node = GrainFilter() + mask = _make_mask((1, 3, 1, 3), (15, 21, 15, 21)) + + # Remove grains larger than 10 px → only small grain survives + result, = node.process(mask, min_area=1, max_area=10, remove_border=False) + assert result[1:3, 1:3].sum() > 0 # small grain kept + assert result[15:21, 15:21].sum() == 0 # large grain removed + + +def test_grain_filter_remove_border(): + from backend.nodes.grain_filter import GrainFilter + + node = GrainFilter() + # Border grain (touches top-left corner) and interior grain + mask = _make_mask((0, 3, 0, 3), (10, 15, 10, 15)) + + result, = node.process(mask, min_area=1, max_area=0, remove_border=True) + assert result[0:3, 0:3].sum() == 0 # border grain removed + assert result[10:15, 10:15].sum() > 0 # interior grain kept + + +def test_grain_filter_no_grains(): + from backend.nodes.grain_filter import GrainFilter + + node = GrainFilter() + mask = np.zeros((20, 20), dtype=np.uint8) + + result, = node.process(mask, min_area=1, max_area=0, remove_border=False) + assert result.sum() == 0 + + +def test_grain_filter_all_grains_kept(): + from backend.nodes.grain_filter import GrainFilter + + node = GrainFilter() + mask = _make_mask((5, 9, 5, 9), (15, 19, 15, 19)) + + # Permissive thresholds: no grains should be removed + result, = node.process(mask, min_area=1, max_area=0, remove_border=False) + assert result[5:9, 5:9].sum() > 0 + assert result[15:19, 15:19].sum() > 0 + + +def test_grain_filter_max_area_zero_means_no_limit(): + from backend.nodes.grain_filter import GrainFilter + + node = GrainFilter() + # Large grain 10×10 = 100 px; max_area=0 means no upper limit + mask = _make_mask((5, 15, 5, 15)) + + result, = node.process(mask, min_area=1, max_area=0, remove_border=False) + assert result[5:15, 5:15].sum() > 0 + + +def test_grain_filter_output_dtype(): + from backend.nodes.grain_filter import GrainFilter + + node = GrainFilter() + mask = _make_mask((5, 10, 5, 10)) + + result, = node.process(mask, min_area=1, max_area=0, remove_border=False) + assert result.dtype == np.uint8 + assert set(result.flat).issubset({0, 255}) diff --git a/tests/node_tests/radial_profile.py b/tests/node_tests/radial_profile.py new file mode 100644 index 0000000..0ab151f --- /dev/null +++ b/tests/node_tests/radial_profile.py @@ -0,0 +1,85 @@ +import numpy as np +import pytest +from tests.node_tests._shared import make_field +from backend.data_types import LineData + + +def test_radial_profile_constant_field(): + """Constant field: profile should equal the constant at every finite bin.""" + from backend.nodes.radial_profile import RadialProfile + + node = RadialProfile() + field = make_field(data=np.full((64, 64), 2.5)) + + result, = node.process(field, cx=0.5, cy=0.5, n_bins=32) + assert isinstance(result, LineData) + assert len(result.data) == 32 + + finite = result.data[np.isfinite(result.data)] + assert finite.size > 0 + assert np.allclose(finite, 2.5, atol=1e-10) + + +def test_radial_profile_units(): + from backend.nodes.radial_profile import RadialProfile + + node = RadialProfile() + field = make_field() + + result, = node.process(field, cx=0.5, cy=0.5, n_bins=32) + assert result.x_unit == field.si_unit_xy + assert result.y_unit == field.si_unit_z + + +def test_radial_profile_x_axis_monotone(): + """Radius axis must be strictly increasing and start near zero.""" + from backend.nodes.radial_profile import RadialProfile + + node = RadialProfile() + field = make_field() + + result, = node.process(field, cx=0.5, cy=0.5, n_bins=64) + assert result.x_axis[0] >= 0.0 + assert np.all(np.diff(result.x_axis) > 0) + + +def test_radial_profile_off_centre(): + """Off-centre origin produces a valid profile with the same number of bins.""" + from backend.nodes.radial_profile import RadialProfile + + node = RadialProfile() + field = make_field(data=np.ones((64, 64))) + + result, = node.process(field, cx=0.0, cy=0.0, n_bins=32) + assert len(result.data) == 32 + finite = result.data[np.isfinite(result.data)] + assert np.allclose(finite, 1.0, atol=1e-10) + + +def test_radial_profile_radial_symmetry(): + """A radially symmetric field should give a smooth, non-constant profile.""" + from backend.nodes.radial_profile import RadialProfile + + node = RadialProfile() + yres, xres = 64, 64 + ys, xs = np.mgrid[0:yres, 0:xres] + r = np.hypot(xs - xres / 2.0, ys - yres / 2.0) + data = np.cos(r * np.pi / (xres / 2.0)) + field = make_field(data=data) + + result, = node.process(field, cx=0.5, cy=0.5, n_bins=32) + finite = result.data[np.isfinite(result.data)] + # The profile should vary (not constant) + assert np.std(finite) > 0.01 + + +def test_radial_profile_n_bins(): + from backend.nodes.radial_profile import RadialProfile + + node = RadialProfile() + field = make_field() + + for n in (16, 64, 256): + result, = node.process(field, cx=0.5, cy=0.5, n_bins=n) + assert len(result.data) == n + assert len(result.x_axis) == n diff --git a/tests/node_tests/resample.py b/tests/node_tests/resample.py new file mode 100644 index 0000000..5bfdfaa --- /dev/null +++ b/tests/node_tests/resample.py @@ -0,0 +1,83 @@ +import numpy as np +import pytest +from tests.node_tests._shared import make_field + + +def test_resample_upsample(): + from backend.nodes.resample import Resample + + node = Resample() + field = make_field(shape=(32, 32)) + result, = node.process(field, width=64, height=64, interpolation="linear") + + assert result.data.shape == (64, 64) + assert np.isclose(result.xreal, field.xreal) + assert np.isclose(result.yreal, field.yreal) + # dx should halve since physical size is same but twice the pixels + assert np.isclose(result.dx, field.dx / 2, rtol=1e-9) + + +def test_resample_downsample(): + from backend.nodes.resample import Resample + + node = Resample() + field = make_field(shape=(64, 64)) + result, = node.process(field, width=32, height=32, interpolation="linear") + + assert result.data.shape == (32, 32) + assert np.isclose(result.xreal, field.xreal) + + +def test_resample_non_square(): + from backend.nodes.resample import Resample + + node = Resample() + field = make_field(shape=(32, 64)) + result, = node.process(field, width=128, height=64, interpolation="nearest") + + assert result.data.shape == (64, 128) + + +def test_resample_interpolation_modes(): + from backend.nodes.resample import Resample + + node = Resample() + field = make_field(shape=(32, 32)) + + for interp in ("linear", "cubic", "nearest"): + result, = node.process(field, width=64, height=64, interpolation=interp) + assert result.data.shape == (64, 64) + + +def test_resample_constant_field(): + """Resampling a constant field should remain constant regardless of interpolation.""" + from backend.nodes.resample import Resample + + node = Resample() + field = make_field(data=np.full((16, 16), 3.14)) + + for interp in ("linear", "cubic", "nearest"): + result, = node.process(field, width=32, height=32, interpolation=interp) + assert np.allclose(result.data, 3.14, atol=1e-6) + + +def test_resample_metadata_preserved(): + from backend.nodes.resample import Resample + + node = Resample() + field = make_field() + result, = node.process(field, width=64, height=64, interpolation="linear") + + assert result.si_unit_xy == field.si_unit_xy + assert result.si_unit_z == field.si_unit_z + assert np.isclose(result.xreal, field.xreal) + + +def test_resample_unknown_interpolation(): + from backend.nodes.resample import Resample + + node = Resample() + field = make_field() + + with pytest.raises(ValueError, match="interpolation"): + node.process(field, width=64, height=64, interpolation="lanczos") diff --git a/tests/node_tests/slope_distribution.py b/tests/node_tests/slope_distribution.py new file mode 100644 index 0000000..82b922e --- /dev/null +++ b/tests/node_tests/slope_distribution.py @@ -0,0 +1,98 @@ +import numpy as np +import pytest +from tests.node_tests._shared import make_field +from backend.data_types import LineData + + +def test_slope_distribution_output_shape(): + from backend.nodes.slope_distribution import SlopeDistribution + + node = SlopeDistribution() + field = make_field() + + for dist in ("theta", "phi", "gradient"): + result, = node.process(field, dist, 90) + assert isinstance(result, LineData) + assert len(result.data) == 90 + assert len(result.x_axis) == 90 + + +def test_slope_distribution_theta_non_negative_and_normalized(): + from backend.nodes.slope_distribution import SlopeDistribution + + node = SlopeDistribution() + field = make_field() + + result, = node.process(field, "theta", 90) + assert np.all(result.data >= 0.0) + assert result.x_unit == "deg" + # Probability density: integral ≈ 1 (bin_width × sum ≈ 1) + bin_width = result.x_axis[1] - result.x_axis[0] + integral = float(np.sum(result.data) * bin_width) + assert abs(integral - 1.0) < 0.05 + + +def test_slope_distribution_gradient_normalized(): + from backend.nodes.slope_distribution import SlopeDistribution + + node = SlopeDistribution() + field = make_field() + + result, = node.process(field, "gradient", 100) + assert np.all(result.data >= 0.0) + # Probability density: integral ≈ 1 + bin_width = result.x_axis[1] - result.x_axis[0] + integral = float(np.sum(result.data) * bin_width) + assert abs(integral - 1.0) < 0.05 + + +def test_slope_distribution_phi_units(): + from backend.nodes.slope_distribution import SlopeDistribution + + node = SlopeDistribution() + field = make_field() + + result, = node.process(field, "phi", 180) + assert result.x_unit == "deg" + assert np.all(result.data >= 0.0) + # x-axis must span [0, 360) + assert result.x_axis[0] >= 0.0 + assert result.x_axis[-1] < 360.0 + + +def test_slope_distribution_flat_field_theta(): + """Flat field: all gradients are zero → theta=0 everywhere, first bin gets all weight.""" + from backend.nodes.slope_distribution import SlopeDistribution + + node = SlopeDistribution() + field = make_field(data=np.zeros((32, 32))) + + result, = node.process(field, "theta", 90) + # With max_theta=0 we use fallback 90; all pixels land in bin 0 + assert result.data[0] == pytest.approx(result.data.sum(), abs=1e-10) + + +def test_slope_distribution_x_ramp_phi(): + """X-only ramp: slope direction should be concentrated near phi=180° (ascending left, or 0°/360°).""" + from backend.nodes.slope_distribution import SlopeDistribution + + node = SlopeDistribution() + data = np.tile(np.linspace(0.0, 1.0, 64), (64, 1)) + field = make_field(data=data) + + result, = node.process(field, "phi", 360) + # Peak should be within first or last few bins (phi ≈ 0 or 2π, i.e., ascending in +x direction) + # atan2(0, -gx) with gx>0 gives atan2(0,-gx) = π, so peak near 180° + peak_idx = int(np.argmax(result.data)) + peak_deg = float(result.x_axis[peak_idx]) + assert 150.0 < peak_deg < 210.0, f"Unexpected peak at {peak_deg}°" + + +def test_slope_distribution_unknown_distribution(): + from backend.nodes.slope_distribution import SlopeDistribution + + node = SlopeDistribution() + field = make_field() + + with pytest.raises(ValueError): + node.process(field, "azimuth", 90) diff --git a/tests/node_tests/tip_blind_estimate.py b/tests/node_tests/tip_blind_estimate.py new file mode 100644 index 0000000..21804ea --- /dev/null +++ b/tests/node_tests/tip_blind_estimate.py @@ -0,0 +1,168 @@ +import numpy as np +import pytest +from tests.node_tests._shared import make_field +from backend.data_types import DataField + + +def run_blind(field, n_pixels=17, threshold=0.0, method="partial", use_edges=False): + from backend.nodes.tip_blind_estimate import BlindTipEstimate + node = BlindTipEstimate() + tip, certainty = node.process( + field=field, n_pixels=n_pixels, threshold=threshold, + method=method, use_edges=use_edges, + ) + return tip, certainty + + +# ── Output types and dimensions ────────────────────────────────────────────── + +def test_outputs_are_data_fields(): + field = make_field(shape=(32, 32), xreal=32e-9, yreal=32e-9) + tip, certainty = run_blind(field, n_pixels=9) + assert isinstance(tip, DataField) + assert isinstance(certainty, DataField) + + +def test_tip_output_shape(): + """Tip must be (n_pixels, n_pixels) — odd, bumped if even.""" + field = make_field(shape=(32, 32), xreal=32e-9, yreal=32e-9) + for n in (7, 11, 17): + tip, _ = run_blind(field, n_pixels=n) + assert tip.data.shape == (n, n), f"Expected ({n},{n}), got {tip.data.shape}" + + +def test_tip_n_pixels_even_bumped(): + field = make_field(shape=(32, 32), xreal=32e-9, yreal=32e-9) + tip, _ = run_blind(field, n_pixels=16) + assert tip.data.shape[0] == 17 + + +def test_certainty_output_matches_field_shape(): + field = make_field(shape=(48, 64)) + _, certainty = run_blind(field, n_pixels=9) + assert certainty.data.shape == field.data.shape + + +def test_certainty_is_binary(): + """Certainty map values must all be 0.0 or 1.0.""" + field = make_field(shape=(32, 32), xreal=32e-9, yreal=32e-9) + _, certainty = run_blind(field, n_pixels=9) + vals = np.unique(certainty.data) + for v in vals: + assert v in (0.0, 1.0), f"Non-binary certainty value: {v}" + + +# ── Tip conventions ─────────────────────────────────────────────────────────── + +def test_tip_min_is_zero(): + field = make_field(shape=(32, 32), xreal=32e-9, yreal=32e-9) + tip, _ = run_blind(field, n_pixels=9) + assert tip.data.min() >= -1e-15 + + +def test_tip_max_at_centre(): + """Apex (centre pixel) must be the maximum of the estimated tip.""" + field = make_field(shape=(32, 32), xreal=32e-9, yreal=32e-9) + tip, _ = run_blind(field, n_pixels=9) + ci = tip.data.shape[0] // 2 + assert tip.data[ci, ci] == pytest.approx(tip.data.max(), abs=1e-20) + + +def test_tip_units_inherited(): + field = make_field(xreal=1e-6, yreal=1e-6) + field.si_unit_xy = "nm" + field.si_unit_z = "V" + tip, _ = run_blind(field, n_pixels=9) + assert tip.si_unit_xy == "nm" + assert tip.si_unit_z == "V" + + +# ── Flat field gives flat (zero) tip ───────────────────────────────────────── + +def test_flat_field_gives_zero_tip(): + """A flat image has no features → blind estimation cannot refine the tip → stays flat.""" + flat = make_field(data=np.full((32, 32), 3.14)) + tip, _ = run_blind(flat, n_pixels=7) + # Tip should be all zeros (no refinement possible) + assert np.allclose(tip.data, 0.0, atol=1e-12) + + +# ── Single spike → estimated tip matches expected shape ────────────────────── + +def test_spike_gives_sharp_tip(): + """ + A single sharp spike dilated by a known paraboloid tip gives a broadened image. + Blind estimation on the broadened image should recover a tip that is ≤ the true tip + everywhere (blind estimation is an upper bound on the true tip shape). + """ + from scipy.ndimage import grey_dilation + from backend.nodes.tip_model import TipModel + + pixel_size = 1e-9 + n = 64 + field_ref = make_field(shape=(n, n), xreal=n * pixel_size, yreal=n * pixel_size) + + # Create true parabolic tip (33×33, radius=20nm) + true_tip_node = TipModel() + true_tip, = true_tip_node.process( + field=field_ref, shape="parabola", radius=20e-9, + half_angle=20.0, n_pixels=33, + ) + + # Spike surface + surface = np.zeros((n, n)) + surface[n // 2, n // 2] = 1e-9 # 1 nm spike + + # Dilated (measured) image + dil_struct = true_tip.data - true_tip.data.max() + measured_data = grey_dilation(surface, structure=dil_struct) + measured = make_field(data=measured_data, xreal=n * pixel_size, yreal=n * pixel_size) + + # Blind estimation + est_tip, certainty = run_blind(measured, n_pixels=33, threshold=0.0, method="partial") + + # The estimated tip must be ≤ the true tip everywhere (blind est. is upper bound) + # Both are min=0, max=apex. Compare normalised shapes. + true_norm = true_tip.data / true_tip.data.max() if true_tip.data.max() > 0 else true_tip.data + est_norm = est_tip.data / est_tip.data.max() if est_tip.data.max() > 0 else est_tip.data + # After normalisation: est ≤ true everywhere (blind is conservative) + assert np.all(est_norm <= true_norm + 1e-6), \ + "Blind estimate exceeded true tip (not a valid upper bound)" + + +# ── Partial vs full ─────────────────────────────────────────────────────────── + +def test_full_method_runs(): + """Full estimation should run without error on a small image.""" + field = make_field(shape=(24, 24), xreal=24e-9, yreal=24e-9) + tip, certainty = run_blind(field, n_pixels=7, method="full") + assert isinstance(tip, DataField) + assert isinstance(certainty, DataField) + + +# ── Certainty increases with sharp features ─────────────────────────────────── + +def test_certainty_nonzero_for_sharp_image(): + """An image with distinct features should produce some certain pixels.""" + from scipy.ndimage import grey_dilation + from backend.nodes.tip_model import TipModel + + pixel_size = 1e-9 + n = 64 + field_ref = make_field(shape=(n, n), xreal=n * pixel_size, yreal=n * pixel_size) + true_tip_node = TipModel() + true_tip, = true_tip_node.process( + field=field_ref, shape="parabola", radius=20e-9, + half_angle=20.0, n_pixels=17, + ) + + # Create a sharp pyramid surface + cx, cy = n // 2, n // 2 + ys, xs = np.mgrid[0:n, 0:n] + surface = np.maximum(0, 5e-9 - 0.3e-9 * np.maximum(np.abs(xs - cx), np.abs(ys - cy))) + dil_struct = true_tip.data - true_tip.data.max() + measured_data = grey_dilation(surface, structure=dil_struct) + measured = make_field(data=measured_data, xreal=n * pixel_size, yreal=n * pixel_size) + + _, certainty = run_blind(measured, n_pixels=17, method="partial") + assert certainty.data.sum() > 0, "No certain pixels found for a sharp image" diff --git a/tests/node_tests/tip_deconvolution.py b/tests/node_tests/tip_deconvolution.py new file mode 100644 index 0000000..89d6093 --- /dev/null +++ b/tests/node_tests/tip_deconvolution.py @@ -0,0 +1,120 @@ +import numpy as np +import pytest +from tests.node_tests._shared import make_field +from backend.data_types import DataField + + +def run_deconv(field, tip): + from backend.nodes.tip_deconvolution import TipDeconvolution + node = TipDeconvolution() + result, = node.process(field=field, tip=tip) + return result + + +def make_tip(shape="parabola", radius=10e-9, half_angle=20.0, n_pixels=33): + from backend.nodes.tip_model import TipModel + field = make_field(shape=(64, 64), xreal=64e-9, yreal=64e-9) + node = TipModel() + result, = node.process(field=field, shape=shape, radius=radius, + half_angle=half_angle, n_pixels=n_pixels) + return result + + +# ── Output type and shape ──────────────────────────────────────────────────── + +def test_deconv_output_is_data_field(): + field = make_field() + tip = DataField(data=np.zeros((1, 1)), xreal=1e-9, yreal=1e-9) + result = run_deconv(field, tip) + assert isinstance(result, DataField) + + +def test_deconv_output_shape_matches_input(): + field = make_field(shape=(64, 64)) + tip = make_tip(n_pixels=11) + result = run_deconv(field, tip) + assert result.data.shape == field.data.shape + + +def test_deconv_preserves_physical_dimensions(): + field = make_field(xreal=2e-6, yreal=3e-6) + tip = DataField(data=np.zeros((1, 1)), xreal=1e-9, yreal=1e-9) + result = run_deconv(field, tip) + assert result.xreal == field.xreal + assert result.yreal == field.yreal + + +def test_deconv_preserves_units(): + field = make_field() + field.si_unit_xy = "nm" + field.si_unit_z = "V" + tip = DataField(data=np.zeros((1, 1)), xreal=1e-9, yreal=1e-9) + result = run_deconv(field, tip) + assert result.si_unit_xy == "nm" + assert result.si_unit_z == "V" + + +# ── Identity with a point tip ──────────────────────────────────────────────── + +def test_deconv_point_tip_is_identity(): + """A 1×1 tip with value 0 is the identity: erosion with structure [[0]] = original.""" + field = make_field(data=np.random.default_rng(0).standard_normal((32, 32)) + 10) + tip = DataField(data=np.zeros((1, 1)), xreal=1e-9, yreal=1e-9) + result = run_deconv(field, tip) + assert np.allclose(result.data, field.data, atol=1e-12) + + +# ── Flat field invariance ──────────────────────────────────────────────────── + +@pytest.mark.parametrize("shape", ["parabola", "cone", "sphere"]) +def test_deconv_flat_field_stays_flat(shape): + """Deconvolution of a constant-valued field must remain constant.""" + flat_data = np.full((64, 64), 5.0) + field = make_field(data=flat_data) + tip = make_tip(shape=shape, n_pixels=15) + result = run_deconv(field, tip) + interior = result.data[8:-8, 8:-8] # avoid border effects + assert np.allclose(interior, 5.0, atol=1e-10) + + +# ── Deconvolution sharpens features ───────────────────────────────────────── + +def test_deconv_sharpens_broadened_image(): + """ + Forward tip dilation broadens a spike; deconvolution (erosion) should remove + the broadening. Result must be ≤ input everywhere (erosion is a lower bound). + """ + from scipy.ndimage import grey_dilation + + # Build a field with a single spike + data = np.zeros((64, 64)) + data[32, 32] = 1.0 + field = make_field(data=data) + + # Create a small parabolic tip + tip = make_tip(shape="parabola", radius=50e-9, n_pixels=11) + tip_data = tip.data + + # Simulate measured image via tip dilation (Gwyddion gwy_tip_dilation): + # dilation_tip = tip - max(tip) (max shifted to 0, values ≤ 0) + # measured[y,x] = max_{ty,tx}[surface[y−ty, x−tx] + dilation_tip[ty,tx]] + dilation_struct = tip_data - tip_data.max() + measured_data = grey_dilation(data, structure=dilation_struct) + measured = make_field(data=measured_data) + + # Deconvolve + result = run_deconv(measured, tip) + + # Erosion is a lower bound on the input + assert np.all(result.data <= measured_data + 1e-10) + + # The spike should be recovered: result at spike position ≥ most surroundings + assert result.data[32, 32] > result.data[32, 36] + + +def test_deconv_erosion_never_exceeds_input(): + """Grey erosion is always ≤ the input (fundamental morphological property).""" + field = make_field(data=np.abs(np.random.default_rng(7).standard_normal((32, 32)))) + tip = make_tip(shape="parabola", n_pixels=7) + result = run_deconv(field, tip) + assert np.all(result.data <= field.data + 1e-10) diff --git a/tests/node_tests/tip_model.py b/tests/node_tests/tip_model.py new file mode 100644 index 0000000..af9c675 --- /dev/null +++ b/tests/node_tests/tip_model.py @@ -0,0 +1,153 @@ +import numpy as np +import pytest +from tests.node_tests._shared import make_field + + +def make_tip(shape="parabola", radius=10e-9, half_angle=20.0, n_pixels=65, field=None): + from backend.nodes.tip_model import TipModel + if field is None: + # 1 nm/pixel reference field + field = make_field(shape=(128, 128), xreal=128e-9, yreal=128e-9) + node = TipModel() + result, = node.process(field=field, shape=shape, radius=radius, + half_angle=half_angle, n_pixels=n_pixels) + return result + + +# ── Shape and dimensionality ──────────────────────────────────────────────── + +def test_tip_output_is_data_field(): + from backend.data_types import DataField + tip = make_tip() + assert isinstance(tip, DataField) + + +def test_tip_shape_matches_n_pixels(): + """Output grid must be (n_pixels, n_pixels).""" + for n in (5, 9, 33, 65): + tip = make_tip(n_pixels=n) + assert tip.data.shape == (n, n) + + +def test_tip_n_pixels_even_forced_odd(): + """Even n_pixels should be silently bumped to n+1.""" + tip = make_tip(n_pixels=64) + assert tip.data.shape[0] % 2 == 1 + assert tip.data.shape[0] == 65 + + +def test_tip_xreal_equals_n_times_pixel_size(): + """xreal must equal n_pixels * pixel_size of the reference field.""" + pixel_size = 1e-9 + field = make_field(shape=(128, 128), xreal=128 * pixel_size, yreal=128 * pixel_size) + tip = make_tip(n_pixels=65, field=field) + assert abs(tip.xreal - 65 * pixel_size) < 1e-25 + + +def test_tip_units_inherited_from_field(): + field = make_field() + field.si_unit_xy = "nm" + field.si_unit_z = "V" + tip = make_tip(field=field) + assert tip.si_unit_xy == "nm" + assert tip.si_unit_z == "V" + + +# ── Apex and minimum conventions ──────────────────────────────────────────── + +@pytest.mark.parametrize("shape", ["parabola", "cone", "sphere"]) +def test_tip_min_is_zero(shape): + """All shapes must shift the data so the minimum is 0.""" + tip = make_tip(shape=shape) + assert tip.data.min() >= -1e-15 + + +@pytest.mark.parametrize("shape", ["parabola", "cone", "sphere"]) +def test_tip_apex_at_centre(shape): + """Centre pixel must equal the maximum for all shapes.""" + tip = make_tip(shape=shape) + n = tip.data.shape[0] + ci = n // 2 + assert tip.data[ci, ci] == pytest.approx(tip.data.max(), abs=1e-20) + + +@pytest.mark.parametrize("shape", ["parabola", "cone", "sphere"]) +def test_tip_radial_symmetry(shape): + """All shapes must be left-right and up-down symmetric.""" + tip = make_tip(shape=shape) + data = tip.data + assert np.allclose(data, np.flipud(data), atol=1e-12) + assert np.allclose(data, np.fliplr(data), atol=1e-12) + + +# ── Parabola formula verification ──────────────────────────────────────────── + +def test_parabola_apex_formula(): + """Parabola apex height = (half_width)² / R (Gwyddion: z0 = 2a·x_half²).""" + radius = 20e-9 + n = 65 + pixel_size = 1e-9 + field = make_field(shape=(256, 256), xreal=256 * pixel_size, yreal=256 * pixel_size) + tip = make_tip(shape="parabola", radius=radius, n_pixels=n, field=field) + + x_half = (n // 2) * pixel_size # = 32 nm + expected_apex = x_half**2 / radius # = 1024e-18 / 20e-9 = 51.2 nm + assert abs(tip.data[n // 2, n // 2] - expected_apex) < 1e-12 * expected_apex + + +def test_parabola_decreases_monotonically_from_centre(): + """Parabola row profile must be monotonically decreasing from centre outward.""" + tip = make_tip(shape="parabola", n_pixels=33) + ci = 16 + row = tip.data[ci, :] + assert np.all(np.diff(row[:ci + 1]) >= 0) # left half increases toward centre + assert np.all(np.diff(row[ci:]) <= 0) # right half decreases from centre + + +def test_parabola_corner_value_is_zero(): + """Gwyddion constructs z0 so that z at (±half, ±half) = 0 exactly.""" + n = 33 + tip = make_tip(shape="parabola", n_pixels=n) + corner = tip.data[0, 0] + # The corner is by construction the minimum; after the min-shift it should be ~0 + assert abs(corner) < 1e-15 * tip.data[n // 2, n // 2] + + +# ── Cone shape ────────────────────────────────────────────────────────────── + +def test_cone_apex_greater_than_edges(): + tip = make_tip(shape="cone", half_angle=20.0, n_pixels=33) + ci = 16 + assert tip.data[ci, ci] > tip.data[0, 0] + assert tip.data[0, 0] >= -1e-15 # edges should be ~0 + + +def test_cone_decreases_away_from_centre(): + """Cone must be monotonically decreasing along any radial line from centre.""" + tip = make_tip(shape="cone", half_angle=20.0, n_pixels=33) + ci = 16 + row = tip.data[ci, :] + assert np.all(np.diff(row[ci:]) <= 0) + + +def test_cone_different_half_angles_differ(): + """Tips with different half-angles must produce different shapes.""" + tip_wide = make_tip(shape="cone", half_angle=40.0, n_pixels=33) + tip_narrow = make_tip(shape="cone", half_angle=10.0, n_pixels=33) + assert not np.allclose(tip_wide.data, tip_narrow.data) + + +# ── Sphere shape ───────────────────────────────────────────────────────────── + +def test_sphere_curvature_near_apex(): + """Near the apex the sphere cap gives z ≈ sqrt(R²-r²), so + z[ci,ci] - z[ci, ci+1] ≈ pixel_size² / (2R) for r << R.""" + radius = 100e-9 + pixel_size = 1e-9 + field = make_field(shape=(256, 256), xreal=256 * pixel_size, yreal=256 * pixel_size) + tip = make_tip(shape="sphere", radius=radius, n_pixels=33, field=field) + ci = 16 + delta_z = tip.data[ci, ci] - tip.data[ci, ci + 1] + expected = pixel_size**2 / (2 * radius) + # relative tolerance 0.1% (small-angle approximation holds well for r << R) + assert abs(delta_z - expected) / expected < 1e-3