nabu.estimation.motion

Classes

PairType
MotionEstimation — A class for estimating rigid sample motion during acquisition. The motion is estimated by 1. Measuring the translations between radios (either opposite radios for 360-degrees scan ; or "reference radios") 2. Fitting these measured translations with a displacement model.

source enum PairType(*args, **kwds)

Bases : Enum

Attributes

OPPOSITE
RETURN

source class MotionEstimation(projs_stack1, projs_stack2, angles1_rad, angles2_rad, indices1=None, indices2=None, n_projs_tot=None, n_calc_threads=None, reference='begin', shifts_estimator='phase_cross_correlation', shifts_estimator_kwargs=None)

A class for estimating rigid sample motion during acquisition. The motion is estimated by 1. Measuring the translations between radios (either opposite radios for 360-degrees scan ; or "reference radios") 2. Fitting these measured translations with a displacement model.

For (1), there are various available functions to estimate the shift between projections. The workhorse is phase cross-correlation, but this class allows to use other functions. For (2), you can pick several displacement model. The default one is a low-degree polynomial.

Once the displacement model is computed, you have the displacement in sample reference frame, and you can "project" these displacements in the detector reference frame. Having the displacements converted as detector shifts allows you to (a) assess the fit between displacement model, and measured detector shifts (b) get a list movements to correct during reconstruction, vertical and/or horizontal. In nabu pipeline, this is handled by "translation_movements_file" in [reconstruction] section.

Parameters

projs_stack1 : numpy.ndarray — Stack of projections, with shape (n_projs, n_rows, n_cols). It has to be of the same shape as 'projs_stack2'. Projection number 'i' in 'projs_stack1' will be compared with projection number 'i' in 'projs_stack2'
projs_stack2 : numpy.ndarray — Stack of projections, with shape (n_projs, n_rows, n_cols). It has to be of the same shape as 'projs_stack1'. Projection number 'i' in 'projs_stack1' will be compared with projection number 'i' in 'projs_stack2'
angles1_rad : numpy.ndarray or list of float — Angles (in radians) of each projection, i.e, projection number 'i' in 'projs_stack1' was acquired at angle theta=angles1_rad[i]
angles2_rad : numpy.ndarray or list of float — Angles (in radians) of each projection, i.e, projection number 'i' in 'projs_stack2' was acquired at angle theta=angles2_rad[i]
indices1 : numpy.ndarray or list of int, optional — Indices corresponding to each projection in 'projs_stack1'. It is used to calculate the curvilinear coordinate for the fit.
indices2 : numpy.ndarray or list of int, optional — Indices corresponding to each projection in 'projs_stack2'. It is used to calculate the curvilinear coordinate for the fit. It is mostly important if projections in 'projs_stack2' are "return projections" (see Notes below for what this means)
n_calc_threads : int, optional — Number of threads to use for calculating phase cross correlation on pairs of images. Default is to use all available threads.

Notes

"Return projections" is the name of extra projections that might be acquired at the end of the scan. For example, for a 180-degrees scan with rotation angles [0, ..., 179.8], extra projections can be saved at angles [180, 90, 0] (the rotation stage rewinds to its original angular position). These extra projection serve to check whether the sample/stage moved during the scan. In this case, angles2_rad = np.radians([180, 90, 0]).

This class works by fitting the measured displacements (in detector space) with a model. This model uses a "normalized coordinate" built upon projection angles & indices. Each pair of radios (projs_stack1[i], projs_stack2[i]) has angles (angles1_rad[i], angles2_rad[i]) and normalized coordinates (a1[i], a2[i]) where a1 = normalize_coordinates(angles1_rad) and a2 = normalize_coordinates(angles2_rad)

Methods

normalize_coordinates — Get the "curvilinear coordinates" that are used (instead of projection angles in radians or degrees) for fit. These coordinates depend on: - how we normalize (wrt total number of angles, or wrt angle max) - the reference projection (start or end)
get_images
compute_detector_shifts — This function computes the shifts between two images of all pairs.
get_model_matrix — Compute the model matrix for horizontal components (x, y) For degree 1: M = [cos(theta) * (a^+ - a) ; -sin(theta) * (a^+ - a) ; 2 ] This matrix needs three ingredients: - The angles "theta" for which pairs of radios are compared. We take the "angles1_rad" provided at this class instanciation. - The corresponding normalized coordinate "a" : a = normalize_coordinates(angles1_rad) - Normalized coordinate "a_plus" of each second radio in a pair.
apply_fit_horiz — Apply the fitted parameters to get the sample displacement in (x, y)
apply_fit_vertic
estimate_horizontal_motion — Estimation of the horizontal motion component.
estimate_vertical_motion — Estimation of the motion vertical component.
mask_pairs — Mask pairs to ignore them for the fit.
project_sample_motion_onto_detector — Convert vectors of motion (t_x, t_y, t_z), from the sample domain, to vectors of motion (t_u, t_v) in the detector domain
apply_fit_to_get_detector_displacements
get_max_fit_error — Get the maximum error of the fit of displacement (in pixel units), projected in the detector domain ; i.e compare shifts_u (measured via phase cross correlation) to M.dot(coeffs) (model fit)
plot_detector_shifts — Plot the detector shifts, i.e difference of u-movements between theta and theta^+. This can be used to compare the fit against measured shifts between pairs of projections.
plot_movements — Plot the movements in the sample and detector domain, for a given vector of angles. Mind the difference with plot_sample_shifts(): in this case we plot proj(txy(theta)) whereas plot_sample_shifts() plots proj(txy(theta^+)) - proj(txy(theta)), which can be used to compare with measured shifts between projections.

source method MotionEstimation.normalize_coordinates(angles, part=1)

Get the "curvilinear coordinates" that are used (instead of projection angles in radians or degrees) for fit. These coordinates depend on: - how we normalize (wrt total number of angles, or wrt angle max) - the reference projection (start or end)

Parameters

angles : array — Array with projection angles
part : int, optional — Which part of the scan the provided angles belong to. Using "part=1" (resp. part=2) means that these angles correspond to "angles1_rad" (resp. angles2_rad)

Raises

ValueError

source method MotionEstimation.get_images(idx, flip_if_needed=True)

source method MotionEstimation.compute_detector_shifts()

This function computes the shifts between two images of all pairs.

source method MotionEstimation.get_model_matrix(do_cor_estimation=True, degree=1)

Compute the model matrix for horizontal components (x, y) For degree 1: M = [cos(theta) * (a^+ - a) ; -sin(theta) * (a^+ - a) ; 2 ] This matrix needs three ingredients: - The angles "theta" for which pairs of radios are compared. We take the "angles1_rad" provided at this class instanciation. - The corresponding normalized coordinate "a" : a = normalize_coordinates(angles1_rad) - Normalized coordinate "a_plus" of each second radio in a pair.

In other words, the pair of radios (projs_stack1[i], projs_stack2[i]) have angles (angles1_rad[i], angles2_rad[i]) and normalized coordinates (a[i], a_plus[i])

The resulting matrix strongly depends on how the angles are ordered/normalized.

source method MotionEstimation.apply_fit_horiz(angles=None, angles_normalized=None)

Apply the fitted parameters to get the sample displacement in (x, y)

Parameters

angles : array, optional — Angles the fit is applied onto.
angles_normalized : array, optional — Normalized angles the fit is applied onto. If provided, takes precedence over 'angles' (see notes below)

Returns

txy : array — Array of shape (n_provided_angles, 2) where txy[:, 0] is the motion x-component, and txy[:, 1] is the motion y-component

Notes

The fit is assumed to have been done beforehand on a series of detector shifts measurements. Once the fit is done, coefficients are extracted and stored by this class instance. The parameter 'angles' provided to this function is normalized before applying the fit. For degree 2, applying the fit is roughly: t_x = alpha_x * a^2 + beta_x * a where 'a' is the normalized angle coordinate, and (alpha_x, beta_x) the coefficients extracted from the fit.

Raises

RuntimeError
ValueError

source method MotionEstimation.apply_fit_vertic(angles=None, angles_normalized=None)

Raises

RuntimeError
ValueError

source method MotionEstimation.estimate_horizontal_motion(degree=1, cor=None, recalculate_shifts=False)

Estimation of the horizontal motion component.

Parameters

degree : int (default=1). — Degree of the polynomial model of the motion in the horizontal plane, for both x and y components.
cor : float, optional — Center of rotation, relative to the middle of the image. If None (default), it will be estimated along with the horizontal movement components. If provided (scalar value), use this value and estimate only the horizontal movement components.
recalculate_shifts : bool, optional — Whether to re-calculate detector shifts (usually with phase cross-correlation) if already calculated

source method MotionEstimation.estimate_vertical_motion(degree=1, recalculate_shifts=False)

Estimation of the motion vertical component.

source method MotionEstimation.mask_pairs(indices)

Mask pairs to ignore them for the fit.

source method MotionEstimation.project_sample_motion_onto_detector(t_xy, t_z, angles, cor=None)

Convert vectors of motion (t_x, t_y, t_z), from the sample domain, to vectors of motion (t_u, t_v) in the detector domain

Parameters

t_xy : numpy.ndarray — Sample horizontal movements with shape (n_angles, 2). The first (resp. second) column are the x-movements (resp. y-movements)
t_z : numpy.ndarray — Sample vertical movements, with the size n_angles
angles : numpy.ndarray — Rotation angle (in degree) corresponding to each component
cor : float, optional — Center of rotation

source method MotionEstimation.apply_fit_to_get_detector_displacements()

source method MotionEstimation.get_max_fit_error()

Get the maximum error of the fit of displacement (in pixel units), projected in the detector domain ; i.e compare shifts_u (measured via phase cross correlation) to M.dot(coeffs) (model fit)

source method MotionEstimation.plot_detector_shifts()

Plot the detector shifts, i.e difference of u-movements between theta and theta^+. This can be used to compare the fit against measured shifts between pairs of projections.

The sample movements were inferred from pairs of projections

projs1[i] at angles1[i], and projs2[i] at angles2[i] From these pairs of projections, a model of motion is built and we can generate: - The motion in sample domain (tx(theta), ty(theta), tz(theta)) for arbitrary theta - The motion in detector domain (t_u(theta), t_v(theta)) (for arbitrary theta) which is the parallel projection of the above

What was used for the fit is t_u(theta^+) - t_u(theta). This function will plot this "difference" between t_u(theta^+) and t_u(theta).

source method MotionEstimation.plot_movements(angles_rad=None, gt_xy=None, gt_z=None, cor_offset=0)

Plot the movements in the sample and detector domain, for a given vector of angles. Mind the difference with plot_sample_shifts(): in this case we plot proj(txy(theta)) whereas plot_sample_shifts() plots proj(txy(theta^+)) - proj(txy(theta)), which can be used to compare with measured shifts between projections.

Raises

ImportError
RuntimeError
ValueError