Processing

This module provides data processing utilities for SPM experiments, including distortion correction, general operations, and helper functions.

Distortion Correction

Algorithms for correcting image distortions and aligning SPM data.

hystorian.processing.distortion.compute_hessian(jac)

Computes the Hessian matrix from the Jacobian for ECC optimization.

The Hessian matrix quantifies the second-order partial derivatives of the cost function with respect to the transformation parameters. It is constructed by summing the outer products of the Jacobian across all spatial axes, resulting in a symmetric matrix that is used to update the transformation parameters during optimization.

Parameters:: jac (ndarray) – The Jacobian matrix, typically of shape (num_params, H, W) for 2D images or (num_params, D, H, W) for 3D images.
Returns:: The Hessian matrix, of shape (num_params, num_params).
Return type:: ndarray

Notes

The implementation uses np.tensordot to sum the products of the Jacobian over all spatial axes, which generalizes to higher dimensions.

This is equivalent to the following code:

hessian = np.empty((np.shape(jac)[0], np.shape(jac)[0]))
for i in range(np.shape(jac)[0]):
    hessian[i,:] = np.sum(np.multiply(jac[i,:,:], jac), axis=(1,2))
    for j in range(i+1, np.shape(jac)[0]):
        hessian[i,j] = np.sum(np.multiply(jac[i,:,:], jac[j,:,:]))
        hessian[j,i] = hessian[i,j]

hystorian.processing.distortion.compute_jacobian(grad, xy_grid, warp_matrix, motion_type='affine')

Computes the Jacobian matrix of the warped image with respect to the transformation parameters.

The Jacobian describes how small changes in the transformation parameters affect the pixel intensities of the warped image. Mathematically, the Jacobian is defined as the matrix of partial derivatives of the warped image with respect to the parameter vector, i.e., J(x; p) = ∂w(x; p)/∂p (see Evangelidis & Psarakis, 2008, IEEE TPAMI, p. 1859).

This function supports multiple motion models, including translation, affine, euclidean, homography, and their 3D variants. The appropriate Jacobian computation is selected based on the dimensionality of the input and the specified motion type.

Parameters:

grad (ndarray) – Gradient(s) of the warped image.
xy_grid (list of ndarray) – Meshgrid coordinates corresponding to the image dimensions.
warp_matrix (ndarray) – Current transformation matrix.
motion_type (str, optional) – Type of transformation: “translation”, “affine”, “euclidean”, “homography”. Default is “affine”.

Returns:

The computed Jacobian matrix, with shape depending on the motion model and image dimensionality.

Return type:

ndarray

Notes

The returned Jacobian is used in ECC optimization to update transformation parameters.
For 2D images, the Jacobian shape is (num_params, H, W).
For 3D images, the Jacobian shape is (num_params, D, H, W).
The grad parameter should be a sequence (such as a list or array) containing the spatial gradients of the warped image. For 2D images, this is typically [grad_x, grad_y], and for 3D images, [grad_x, grad_y, grad_z]. Each gradient array should have the same shape as the warped image and represent the rate of change of pixel intensities along the respective axis.

hystorian.processing.distortion.custom_warp(im, mat, motion_type='affine', order=1, cval=0.0)

Applies a geometric transformation to an image using a specified transformation matrix.

Parameters:

im (ndarray) – The input image array to be transformed.
mat (ndarray) – The transformation matrix. For affine, this is typically a 2x3 or 3x3 matrix. For homography, this should be a 3x3 matrix.
motion_type (str, optional) – The type of transformation to apply. Supported values are “affine” and “homography”. Default is “affine”.
order (int, optional) – The order of the spline interpolation used for the transformation. Default is 1 (linear interpolation).

Returns:

The transformed (warped) image.

Return type:

ndarray

Notes

This function supports both affine and homography transformations.
For affine transformations, it uses scipy.ndimage.affine_transform, which applies a linear mapping from input coordinates to output coordinates.
For homography transformations, it uses scipy.ndimage.geometric_transform with a custom coordinate mapping function, allowing for more complex projective transformations.
The function assumes the transformation matrix is correctly formatted for the chosen motion type.

hystorian.processing.distortion.find_transform_ecc(ir, iw, warp_matrix=None, motion_type='affine', number_of_iterations=200, termination_eps=-1.0, gauss_filt_size=5.0, order=1)

Estimates the geometric transformation matrix that best aligns a warped image (iw) to a reference image (ir) by maximizing the Enhanced Correlation Coefficient (ECC).

The function iteratively updates the transformation matrix to maximize the normalized correlation between the reference and warped images. It supports several motion models, including translation, affine, euclidean, and homography. Gaussian filtering is applied to both images before alignment to reduce noise and improve stability.

Parameters:

ir (ndarray) – Reference image (target for alignment).
iw (ndarray) – Warped image to be corrected.
warp_matrix (ndarray, optional) – Initial guess for the transformation matrix. If None, an identity matrix is used. Should be 3x3 for 2D images or 4x4 for 3D images. Default is None.
motion_type (str, optional) – Type of transformation: “translation”, “affine”, “euclidean”, or “homography”. Default is “affine”.
number_of_iterations (int, optional) – Maximum number of iterations before the algorithm stops. Default is 200.
termination_eps (float, optional) – Threshold for the absolute difference between the normalized correlation of successive iterations. If the difference is less than this value, the algorithm stops. Default is -1.0 (no early stopping).
gauss_filt_size (float, optional) – Standard deviation for Gaussian kernel used to blur ir and iw before alignment. Default is 5.0.
order (int, optional) – Interpolation order used in the warp function. Default is 1.

Returns:

warp_matrix – The transformation matrix that best aligns iw to ir.

Return type:

ndarray

Raises:

ValueError –

If the correlation coefficient (rho) becomes NaN, indicating a failure in the algorithm. - If the algorithm stops before convergence, indicating that the images may be uncorrelated or non-overlapping.

hystorian.processing.distortion.project_onto_jacobian(jac, mat)

Projects a matrix onto the Jacobian by computing the sum of element-wise products.

This function multiplies the Jacobian array and the input matrix element-wise and then sums over all spatial axes. The result is a vector where each element represents the projection of the input matrix onto one parameter direction of the Jacobian.

Parameters:

jac (ndarray) – The Jacobian matrix, typically of shape (num_params, H, W) for 2D images or (num_params, D, H, W) for 3D images.
mat (ndarray) – The matrix to project, typically of shape (H, W) for 2D or (D, H, W) for 3D.

Returns:

The projection result, a vector of length num_params.

Return type:

ndarray

Notes

In the original code, the Jacobian was stored as a 2D [K*H, W] array, and the code looped through K, splitting the matrix into K submatrices.
Each submatrix and the input matrix were flattened into vectors, and a dot product was applied. This is equivalent to multiplying the two matrices together element-by-element, then summing the result. Here, the Jacobian is stored as a 3D array [K, H, W], so broadcasting is used to avoid explicit looping. The summation is performed over all spatial axes, which generalizes to higher dimensions.

hystorian.processing.distortion.update_warping_matrix(map_matrix, update, motion_type='affine')

Updates the transformation matrix by applying the parameter increments computed during ECC optimization.

This function supports multiple motion models, including translation, affine, euclidean, homography, and their 3D variants. The update is performed by adding the parameter increments to the appropriate elements of the transformation matrix, according to the selected motion type.

Parameters:

map_matrix (ndarray) – The current transformation matrix to be updated. For 2D images, this is typically a 3x3 matrix. For 3D images, this is typically a 4x4 matrix.
update (ndarray or sequence) – The parameter increments to be applied to the transformation matrix. The length and meaning of this vector depend on the motion type.
motion_type (str, optional) – The type of transformation model. Supported values are “translation”, “affine”, “euclidean”, “homography”, and their 3D variants. Default is “affine”.

Returns:

The updated transformation matrix.

Return type:

ndarray

General Correction

Functions for general image correction and preprocessing.

hystorian.processing.correction.line_fit(line, order=1, box=None)

Do a nth order polynomial line flattening

Parameters:

line (1d array-like)
order (integer)

Returns:

result – same shape as data

Return type:

1d array-like

hystorian.processing.correction.line_flatten_image(data, order=1, axis=0, box=None)

Do a line flattening

Parameters:

data (2d array)
order (integer)
axis (integer) – axis perpendicular to lines

Returns:

result – same shape as data

Return type:

array-like

hystorian.processing.correction.plane_flatten_image(data, order=1, box=None)

Do a plane flattening

Parameters:

data (2d array)
order (integer)

Returns:

result – same shape as data

Return type:

array-like

SPM Utilities

Tools and helpers specific to Scanning Probe Microscopy data.

hystorian.processing.spm.binarize_phase(phase: ndarray[Any, dtype[ScalarType]]) → ndarray[Any, dtype[ScalarType]]

Binarize a phase image by finding the midpoint between two largest peaks in the histogram of phases.

Parameters:: phase (npt.NDArray) – Phase image to binarize
Returns:: Binarized phase image, where values below the midpoint are 0 and above are 1
Return type:: npt.NDArray

hystorian.processing.spm.clean_loop(bias, phase, amp, threshold=None)

Used to determine if a SSPFM loop is good or not by calculating the area encompassed by the hysteresis curve and comparing it to a threshold

Parameters:

bias (nd-array) – 3d-array containing the bias applied to each point of the grid
phase (nd-array) – 3d-array containing the phase applied to each point of the grid
amp (nd-array) – 3d-array containing the amplitude applied to each point of the grid
threshold (int or float, optional) – minimal value of the loop area to be considered a good loop (default: None) if set to None will use threshold = np.mean(area_grid_full) - 2 * np.std(area_grid_full)

Returns:

good_bias (list) – list of the bias corresponding to good loops
good_phase (list) – list of the phase corresponding to good loops
good_amp (list) – list of the amplitudes corresponding to good loops
mask (list) – a 2D mask where a 1 correspond to a good loop and 0 to a bad loop, can be used to mask the input data.

hystorian.processing.spm.pfm_params_map(bias: ndarray[Any, dtype[ScalarType]], phase: ndarray[Any, dtype[ScalarType]]) → tuple[ndarray[Any, dtype[ScalarType]]]

PFM_params_map calculates physically relevant hysteresis parameters from bias and phase channels.

Parameters:

bias (npt.NDArray) – array containing the bias acquired during an SSPFM. This should be the “on” bias computed with extract_hist
phase (npt.NDArray) – array containing the phase acquired during an SSPFM. This should be the “off” phase computed with extract_hist

Returns:

coerc_pos (npt.NDArray) – the positive coercive bias
coerc_neg (npt.NDArray) – the negative coercive bias
step_left (npt.NDArray) – The size of the phase jump (left side)
step_right (npt.NDArray) – The size of the phase jump (right side)
imprint (npt.NDArray) – The imprint of the switching loop
phase_shift (npt.NDArray) – The phase shift of the switching loop

hystorian.processing.spm.shift_and_wrap_phase(phase: ndarray[Any, dtype[ScalarType]], phase_step: int = 1) → ndarray[Any, dtype[ScalarType]]

Finds the smallest shift that minimizes the phase jump in a phase series

Parameters:

phase (a 1D array of consecutive phase measurements)
phase_step (the algorithm will try shifts every phase_step. Increase this to speed up execution.)

Return type:

The phase image shifted such that there is the smallest jumps between consecutive phase points.

Operations

General mathematical and image operations for data processing.

hystorian.processing.operation.gauss_area(x, y)

Determine the area created by the polygon formed by x,y using the Gauss’s area formula (also called shoelace formula)

Parameters:

x (array_like) – values along the x axis
y (array_like) – values along the x axis

Returns:

value corresponding at the encompassed area

Return type:

float