matei/tono
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

tip modelling and deconvolution

Browse Source
This commit is contained in:
matei jordache
2026-03-29 21:49:17 -07:00
parent 24b2c55f2a
commit 1df4df2811
23 changed files with 2231 additions and 28 deletions
Download Patch File Download Diff File Expand all files Collapse all files

11
backend/node_menu.py
View File

@@ -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",
],
}

10
backend/nodes/__init__.py
View File

@@ -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,
)

61
backend/nodes/field_arithmetic.py Normal file
View File

@@ -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),)

63
backend/nodes/gradient.py Normal file
View File

@@ -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),)

72
backend/nodes/grain_filter.py Normal file
View File

@@ -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()
)

75
backend/nodes/radial_profile.py Normal file
View File

@@ -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,
),)

52
backend/nodes/resample.py Normal file
View File

@@ -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),)

120
backend/nodes/slope_distribution.py Normal file
View File

@@ -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),)

582
backend/nodes/tip_blind_estimate.py Normal file
View File

@@ -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)

44
backend/nodes/tip_deconvolution.py Normal file
View File

@@ -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),)

109
backend/nodes/tip_model.py Normal file
View File

@@ -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,)