Math utils¶

class surfalize.mathutils.Cylinder(point, direction, radius)¶

Bases: object

Representation of a cylinder in 3d space, defined by its rotation axis (a point on the axis and a unit direction vector) and its radius.

The cylinder is the set of points at perpendicular distance radius from the axis line. This class provides the geometry needed for cylinder form removal: computing the distance of points to the axis, fitting a cylinder to a point cloud and intersecting vertical lines (along z) with the cylinder surface.

Parameters:

pointarray-like[float, float, float]: A point lying on the rotation axis.
directionarray-like[float, float, float]: Direction vector of the rotation axis. It is normalized internally and does not need to be a unit vector.
radiusfloat: Radius of the cylinder.

Methods

`distance_to_axis`(points)	Computes the perpendicular distance of one or more points to the rotation axis.
`from_fit`(points, guess[, radius])	Fits a cylinder to a point cloud by minimizing the sum of squared radial deviations from the cylinder surface using scipy.optimize.least_squares.
`intersect_vertical`(x, y, z_reference)	Intersects vertical lines (parallel to the z-axis) at positions (x, y) with the cylinder surface and returns the z-coordinate of the intersection.
`residuals`(points)	Computes the signed radial deviation of points from the cylinder surface, i.e. the distance to the axis minus the radius.

distance_to_axis(points)¶

Computes the perpendicular distance of one or more points to the rotation axis.

Parameters:

pointsndarray: Array of shape (N, 3) containing the (x, y, z) coordinates of N points.

Returns:

ndarray: Array of shape (N,) containing the distance of each point to the axis.

classmethod from_fit(points, guess, radius=None)¶

Fits a cylinder to a point cloud by minimizing the sum of squared radial deviations from the cylinder surface using scipy.optimize.least_squares.

The axis is parametrized minimally with respect to the dominant component of the initial guess direction to avoid gauge freedom: the dominant direction component is fixed to one while the remaining two components are free, and the axis point is constrained to the plane where the dominant coordinate is zero. This leaves four free parameters for the axis (two for orientation, two for position), plus the radius if it is not fixed.

Parameters:

pointsndarray: Array of shape (N, 3) containing the (x, y, z) coordinates of the point cloud.
guessCylinder: Initial guess for the cylinder. A good guess for the axis orientation, position and radius is required for reliable convergence.
radiusfloat | None, default None: If provided, the radius is held fixed at this value and only the axis is optimized. If None, the radius is optimized as well, starting from guess.radius.

Returns:

Cylinder: The fitted cylinder.

intersect_vertical(x, y, z_reference)¶

Intersects vertical lines (parallel to the z-axis) at positions (x, y) with the cylinder surface and returns the z-coordinate of the intersection. For each (x, y) position there are generally two intersections (the upper and lower side of the cylinder); the one closest to z_reference is returned. Positions whose vertical line does not intersect the cylinder yield NaN.

Parameters:

x, yndarray: Coordinates at which to intersect.
z_referencendarray: Reference z-values used to select between the two intersection points (typically the measured heights).

Returns:

ndarray: z-coordinate of the cylinder surface at each (x, y) position.

residuals(points)¶

Computes the signed radial deviation of points from the cylinder surface, i.e. the distance to the axis minus the radius.

Parameters:

pointsndarray: Array of shape (N, 3) containing the (x, y, z) coordinates of N points.

Returns:

ndarray: Array of shape (N,) containing the radial deviation of each point.

class surfalize.mathutils.Sinusoid(amplitude, period, x0, y0)¶

Bases: object

Generic sinusoid representation.

Parameters:

amplitudefloat: amplitude.
periodfloat: period.
x0float: offset in x.
y0float: offset in y.

Methods

`__call__`(x)	Computes the value of a generic sinusoid at position x.
`first_extremum`()	Compute the position of the first extremum (peak or valley) for x >= 0.
`first_peak`()	Compute the position of the first peak for x >= 0.
`from_fit`(xdata, ydata[, p0, infer_p0])	Fit a general sinusoid to x and y data using scipy.optimize.curve_fit.

first_extremum()¶

Compute the position of the first extremum (peak or valley) for x >= 0.

Returns:

xfefloat

first_peak()¶

Compute the position of the first peak for x >= 0.

Returns:

xfpfloat

classmethod from_fit(xdata, ydata, p0=None, infer_p0=False)¶

Fit a general sinusoid to x and y data using scipy.optimize.curve_fit. Optionally, an initial guess for the sinusoid parameters a, p, x0, y0 can be specified using the p0 keyword argument, where a is the amplitude, p is the period, x0 is the lateral and y0 the vertical offset.

This function uses unbounded fitting, which can result in negative amplitude, because bounded fitting using scipy.optimize.curve_fit invokes a different algorithm which seems to perform worse on this specific fitting problem. Instead, a result with negative amplitude is converted into the positive amplitude equivalent by phase shifting.

Parameters:

xdatalist-like: array of x-data
ydatalist-like: array of y-data
p0list-like[float, float, float, float] | None, defaults to None: Optional initial guess for the parameters a, p, x0, y0.
infer_p0bool, defaults to False: If True, automatically infers starting guesses of the parameters. Any values provided to the p0 keyword argument will be overwritten with the automatically inferred values.

Returns:

Sinusoid

surfalize.mathutils.argclosest(x, xdata)¶

Returns the index of the value in an array that is closest to the value x.

Parameters:

xfloat: value to which closest array value index should be computed
xdataarray_like: array of x-values

Returns:

indexint: Index of the value in xdata that is closest to x.

surfalize.mathutils.argmax_all(arr)¶

Returns all indices where the array reaches its maximum value.

Parameters:

arrarray-like: Input array

Returns:

numpy.ndarray: Array of indices where the maximum value occurs

surfalize.mathutils.argmin_all(arr)¶

Returns all indices where the array reaches its minimum value.

Parameters:

arrarray-like: Input array

Returns:

numpy.ndarray: Array of indices where the minimum value occurs

surfalize.mathutils.closest(x, data)¶

Returns the value in an array that is closest to the value x.

Parameters:

xfloat: value to which closest array value index should be computed
xdataarray_like: array of x-values

Returns:

value: Value in xdata that is closest to x.

surfalize.mathutils.get_period_fft_1d(xdata, ydata)¶

Estimates the dominant period from a 1d periodic profile with uniformly spaced points.

Parameters:

xdatalist-like: array of uniformly spaced xdata.
ydatalist-like: array of ydata with the same size as xdata.

Returns:

periodfloat

surfalize.mathutils.interp1d(xdata, ydata, assume_sorted=False)¶

Creates a function that linearly interpolates the given x- and y-data at any value of x.

Mimics scipy.interpolate.interp1d, since the scipy version is no longer supported and might be removed in future versions. Contrary to the scipy implementation, assume_sorted is False by default, since a default True value would lead to more errors. If the xdata is sorted in ascending order, then assume_sorted can be set to True to gain a performance increase.

Parameters:

xdataarray_like: array of x-values
ydataarray_like: array of y-values
assume_sortedbool, default False: If True, they xdata array must be supplied with ascendingly ordered values and sorting is skipped. If False, the array xdata will be sorted in ascending order and the array ydata will be sorted accordingly.

Returns:

function: Linear interpolation function y(x).

surfalize.mathutils.interpolate_line_on_2d_array(array, start, end, order=3, num_points=100)¶

Interpolates a line between two points on a 2d array using spline interpolation.

Parameters:

array2d array-like
starttuple[int, int]: Index of the start point.
endtuple[int, int]: Index of the end point.
orderint: Order of spline interpolation. Defaults to 3.
num_pointsint: Number of points of the interpolated line. Defaults to 100.

Returns:

np.ndarray

surfalize.mathutils.otsu_threshold(values, nbins=256)¶

Computes the optimal threshold that separates the values into two classes using Otsu’s method, which maximizes the between-class variance. For one-dimensional two-class segmentation this yields the global optimum and is a fast, deterministic alternative to a two-cluster k-means.

Parameters:

valuesarray_like: Input values. Non-finite values (NaN, inf) are ignored.
nbinsint, default 256: Number of histogram bins used to evaluate the threshold.

Returns:

thresholdfloat: Value separating the two classes. Values greater than the threshold belong to the upper class.

surfalize.mathutils.resolve_box(box, border, x_max, y_max)¶

Resolves a rectangle specification from either an explicit box or a border inset.

Exactly one of box or border must be given. box is returned unchanged as an (x0, x1, y0, y1) tuple. border specifies the distance from each edge inwards to the rectangle, in the same units and axis order as box: a scalar insets all four edges equally, while an (x0, x1, y0, y1) tuple insets each edge individually. x_max and y_max are the coordinates of the far edges (width_um/height_um in physical units, or size.x - 1/size.y - 1 in pixels).

Parameters:

boxtuple[float, float, float, float] | None: Explicit rectangle as (x0, x1, y0, y1), or None if border is used instead.
borderfloat | tuple[float, float, float, float] | None: Distance from each edge inwards, as a scalar or (x0, x1, y0, y1) tuple, or None if box is used instead.
x_max, y_maxfloat: Coordinates of the far edges along x and y.

Returns:

boxtuple[float, float, float, float]: The resolved (x0, x1, y0, y1) rectangle.